Random-mass disorder in the critical Gross-Neveu-Yukawa models

An important yet largely unsolved problem in the statistical mechanics of disordered quantum systems is to understand how quenched disorder affects quantum phase transitions in systems of itinerant fermions. In the clean limit, continuous quantum phase transitions of the symmetry-breaking type in Dirac materials such as graphene and the surfaces of topological insulators are described by relativistic (2+1)-dimensional quantum field theories of the Gross-Neveu-Yukawa (GNY) type. We study the universal critical properties of the chiral Ising, XY, and Heisenberg GNY models perturbed by quenched random-mass disorder, both uncorrelated or with long-range power-law correlations. Using the replica method combined with a controlled triple epsilon expansion below four dimensions, we find a variety of new finite-randomness critical and multicritical points with nonzero Yukawa coupling between low-energy Dirac fields and bosonic order parameter fluctuations, and compute their universal critical exponents. Analyzing bifurcations of the renormalization-group flow, we find instances of the fixed-point annihilation scenario---continuously tuned by the power-law exponent of long-range disorder correlations and associated with an exponentially large crossover length---as well as the transcritical bifurcation and the supercritical Hopf bifurcation. The latter is accompanied by the birth of a stable limit cycle on the critical hypersurface, which represents the first instance of fermionic quantum criticality with emergent discrete scale invariance.


I. INTRODUCTION
Understanding the effect of quenched disorder on continuous quantum phase transitions is a question of enduring interest [1], motivated by the ubiquitous presence of imperfections in the condensed matter systems that exhibit such transitions. In the clean limit, the divergence of the correlation length at criticality produces universal critical phenomena that are controlled by renormalization-group (RG) fixed points of a translationally invariant continuum quantum field theory. A given disorder configuration manifestly breaks translation symmetry even on long length scales, and thus produces behavior very different from that of a translationally invariant system. However, that symmetry is restored in physical properties upon averaging over such configurations. The cases of main interest, then, are those in which disorder qualitatively affects the long-distance physics even after disorder averaging. For quantum critical points (QCPs) described at long distances by a (2+1)D strongly interacting conformal field theory in the clean limit, such as many QCPs of interest in condensed matter physics [2], determining the fate of the system in the infrared after disorder averaging is a problem fraught with technical difficulties. For example, thermodynamic properties are in principle determined by first computing the partition function of a strongly coupled quantum field theory with spatially random couplings, then averaging its logarithm over some chosen probability distribution. An approach better suited to determining the long-distance behavior of the system, our only concern here, is to investigate the RG flow of disorder-averaged observables. It was recently shown [3,4] that this is equivalent to studying the RG flow of an effective theory with disorder-induced translationally-invariant interactions, derived using the standard replica trick [5], despite the formally nonlocal nature of such theories and oft-invoked concerns about the validity of analytically continuing the number of replicas to zero.
A situation of particular interest is one in which disorder produces RG flows on the critical hypersurface that connect the clean fixed point (CFP) describing the transition in the absence of disorder to fixed points characterized by a nonzero value of the effective disorder coupling(s). Such disordered fixed points (DFPs) exhibit scaling behavior but, by contrast with CFPs, no selfaveraging in the thermodynamic limit [6]. We will be exclusively concerned with random-mass disorder, also known as random-T c disorder in the context of classical (thermal) phase transitions, every configuration of which preserves those symmetries of the system that are broken spontaneously at the transition. (In d = 2 spatial dimensions, the focus of this paper, random-field disorder-which violates those symmetries-precludes long-range order, and thus the possibility of a sharp transition [7][8][9][10].) The standard scenario is one in which random-mass disorder is a relevant perturbation at the CFP [ Fig. 1(a)], and drives a direct RG flow to a DFP. For shortrange correlated disorder, this occurs when the correlation length exponent ν CFP of the CFP obeys the Harris inequality ν CFP < 2/d [11]. Examples include the superfluid-Mott glass transition of bosons with particle-hole symmetry in d = 2 [12] and d = 3 [13], described by the O(2) vector model with random-mass disorder in (2+1)D (ν CFP ≈ 0.67 < 1) and (3+1)D (ν CFP = 1/2 < 2/3), respectively. The true correlation length exponent ν in the presence of disorder, i.e., its value at the DFP, obeys the Chayes inequality ν 2/d [14]; the dynamic critical exponent z changes from its Lorentz-invariant value z = 1 at the conformally invariant O(2) Wilson-Fisher fixed point to some noninteger but equally universal (and finite) value z > 1 at the DFP (see Table I). Finiterandomness DFPs in the O(n) vector model are in principle accessible via perturbative RG analyses of the disorder-averaged effective field theory combined with (double) epsilon [15][16][17][18][19] or 1/n expansions [20]. Infinite-randomness DFPs (for which z = ∞) are also possible, such as those describing the random-bond transverse-field Ising model at criticality in (1+1)D [21] (ν CFP = 1 < 2) and (2+1)D [22] (ν CFP ≈ 0.63 < 1). These, however, are not adequately captured by perturbative RG analyses of a disorder-averaged continuum field theory, given the runaway flow to infinite disorder [ Fig. 1(b)]. Rather, they can be quantitatively studied using strong-disorder real-space RG methods [21][22][23][24] which, in spatial dimensions d 2 at least, must be implemented numerically in microscopic lattice models.
If disorder is Harris-irrelevant, ν CFP > 2/d, the standard lore is that disorder has no effect on the phase transition at long distances [ Fig. 1(c)]. However, as observed in Ref. [25], an irrelevant perturbation with a finite coefficient can have nontrivial consequences on the RG flow finitely away from the CFP, just as formally irrelevant interactions at a stable noninteracting fixed point can eventually trigger a phase transition and produce a critical fixed point. The simplest possible RG flow leading to a DFP in the case of Harris-irrelevant disorder is illustrated in Fig. 1(d), and was recently found in a double epsilon-expansion study of the random-mass chiral XY Gross-Neveu-Yukawa (GNY) model [26], a fermionic analog of the O(2) vector model that, absent disorder, describes the quantum phase transition between a Dirac semimetal and a gapped superconductor (Ref. [27,28], and also see Sec. II). Below a separatrix line controlled by a disordered saddle-type fixed point (DFP 1 ), the transition is in the same universality class as the clean sytem, while above that separatrix line, the transition is governed by a disordered critical point (DFP 2 ). A similar RG flow is found in the classical 2D Ising model with binary (±J) random-bond disorder [29]. For a weak concentration of antiferromagnetic bonds randomly distributed amidst ferromagnetic bonds, the paramagnetic-ferromagnetic critical behavior is controlled by the clean 2D Ising fixed point, consistent with the fact that random-mass disorder is (marginally) Harris irrelevant at that fixed point [30,31]. For sufficiently strong disorder, however, the clean critical behavior gives way to critical behavior controlled by a zero-temperature disordered fixed point (spin-glass critical point) via an intervening disordered multicritical point, the Nishimori point [32,33].
Coming back to the RG flow of the random-mass chiral XY GNY model [26], depending on the number of fermion flavors (see Sec. II) the disordered critical point (DFP 2 ) is found to be either a standard sink-type fixed point, as illustrated in Fig. 1(d), or a fixed point of stable-focus type. In the latter case, RG trajectories asymptotically spiral towards the fixed point, implying oscillatory corrections to scaling. Stable-focus fixed points have been found before in replica RG studies of both classical [34] and quantum [15][16][17][18][19]35] disordered systems, and are sometimes considered an artefact of perturbative replica-based RG. However, such flows cannot be ruled out as a matter of principle, since DFPs are in general non-unitary, and real, non-unitary, scale-invariant quantum field theories can have pairs of scaling fields with complex-conjugate dimensions [4,36]. Furthermore, oscillations in scaling laws, characteristic of spiraling or cyclic RG flows, have also been found in numerical studies of disordered holographic models [37] which rely neither on the replica trick nor on perturbation theory (in either the interaction or disorder strengths). A recent Monte Carlo study of classically frustrated 3D Heisenberg antiferromagnets also supports the existence of a stable-focus critical point (in this case, in a clean system) [38].
In the present work, we extend the study of Ref. [26] in two directions. First, in Ref. [26], only short-range correlated (or equivalently at long distances, uncorrelated) random-mass disorder was considered. Here we additionally consider random-mass disorder with correlations between two spatial points x, x that decay asymptotically as a power law, ∼ |x−x | −α , with α < d. (For α > d, the correlations are short range, as the disorder correlation function in momentum space remains finite in the long-wavelength limit.) A clean critical point with correlation length exponent ν CFP is perturbatively stable against such long-range correlated disorder if ν CFP > 2/ min(d, α) [39]; this type of disorder thus generally has a stronger effect at phase transitions than uncorrelated disorder. Second, Ref. [26] only studied the chiral XY GNY model. Here, we perform a comprehensive study of the effect of random-mass disorder in the three standard families of critical GNY models: the chiral Ising, XY, and Heisenberg models [40,41], fermionic analogs of the Ising, XY, and Heisenberg Wilson-Fisher universality classes, respectively. As we briefly review in Sec. II, these chiral GNY models describe a variety of QCPs in condensed matter systems [42].
Our main results are summarized as follows. For the chiral Ising GNY model, we find new disordered multicritical points, and for the chiral XY and Heisenberg GNY models, new disordered critical and multicritical points. As in Ref. [26], some of the disordered QCPs found exhibit usual sink-type RG flows, while others are of stable-focus type. We also explore how the structure of the RG flow on the critical hypersurface evolves upon tuning RG-invariant system parameters, here the number N of fermion flavors and the exponent α describing disorder correlations. We are particularly interested in bifurcations of these RG flows [43], where the number or stability properties of fixed points suddenly change as a function of N and α, called control parameters in bifurcation theory. We find and analyze instances of the saddle-node bifurcation, also known as the fixed-point annihilation scenario [44], at which a repulsive fixed point and an attractive fixed point coalesce and disappear into the complex plane. This type of bifurcation appears or has been argued to appear in RG flows in a variety of problems of current interest in both high-energy physics [36,[44][45][46][47][48][49][50] and condensed matter physics/statistical mechanics [51][52][53][54][55][56][57][58][59]. The characteristic phenomenology associated with it includes Berezinskii-Kosterlitz-Thouless/Miransky scaling, walking/pseudo-critical behavior, and weakly first-order transitions. In our particular problem, it manifests itself in the existence of an anomalously (i.e., exponentially) large length scale L * that governs the crossover between two distinct universality classes of critical behavior. In much previous work, the saddle-node bifurcation is tuned by a parameter such as space(time) dimensionality d or the integer number N of components of a fermionic or bosonic field, and thus cannot be approached continuously in practice. Here, for fixed d and N the bifurcation can be approached by continuously tuning the exponent α for disorder correlations.
Besides the saddle-node bifurcation, we also discover instances of more exotic bifurcations [43]: the transcritical bifurcation, at which two fixed points exchange their stability properties without annihilating, and the supercritical Hopf (or Poincaré-Andronov-Hopf) bifurcation [60]. The latter is a bifurcation at which a stable-focus QCP loses its stability by giving birth to a stable limit cycle, which then controls the asymptotic critical behavior. A possibility first considered by Wilson [61], stable RG limit cycles lead to log-periodic scaling behavior [62], i.e., discrete scale invariance (as opposed to log-periodic behavior of corrections to scaling at stable-focus points). Hopf bifurcations in RG flows were found in classical disordered O(n) models [39,63], but only the subcritical Hopf bifurcation [60] was found, where an unstable-focus fixed point becomes stable and gives birth to an unstable limit cycle. As a result, the models studied in Refs. [39,63] did not exhibit log-periodic critical scaling behavior in the long-distance limit.
The rest of the paper is structured as follows. In Sec. II, we briefly describe the chiral GNY models with long-range correlated random-mass quenched disorder. In Sec. III, we describe the perturbative RG scheme used to derive beta functions on the critical hypersurface. By contrast with Ref. [26], where the double epsilon [16][17][18][19] expansion was sufficient to tame RG flows in the presence of uncorrelated disorder, here we use a controlled triple epsilon expansion [64] at one-loop order that allows us to tame the flow of both interaction and correlated disorder strengths. In Sec. IV, we investigate the fixed points of the RG beta functions derived in Sec. III, focusing on DFPs and analyzing their linear stability. We compute critical exponents and anomalous dimensions at all DFPs. In Sec. V, we discuss qualitative features of the RG flow, including various bifurcations that occur under changes of the control parameters N and α, and their consequences for critical properties. We conclude in Sec. VI with a summary of our main results and a few directions for further research. Three appendices (App. A-C) contain the details of some calculations.

II. THE RANDOM-MASS GNY MODELS
Our starting point is the family of chiral O(n) GNY models in 2+1 dimensions at zero temperature, described by the Euclidean action: where x denotes spatial coordinates, and τ is imaginary time. The model consists of a real ncomponent scalar field φ = (φ 1 , . . . , φ n ), the order parameter, governed by the Lagrangian: where It is coupled to a Dirac fermion field ψ, described by the Lagrangian: The scalar mass squared r in Eq. (2) tunes the model through criticality: r < 0 gives a phase with spontaneously broken O(n) symmetry, r > 0 is the symmetric phase, and r = 0 is the critical point. The parameter λ 2 describes self-interactions of the order parameter. We define the Dirac adjoint in Eq. (3) as ψ = −iψ † γ 0 . We denote γ = (γ 1 , γ 2 ), and γ µ , µ = 0, 1, 2 are Hermitian Dirac matrices obeying the SO(3) Clifford algebra {γ µ , γ ν } = 2δ µν . In the ordinary GNY model, Lorentz invariance (exact or emergent at criticality [65]) demands that the fermion c f and boson c b velocities be equal, but in the presence of quenched disorder, to be introduced below, the ratio c = c f /c b will flow under RG transformations. We perform perturbative calculations near four dimensions at one-loop order in the context of a particular epsilon-expansion scheme to be explained below, but we are ultimately interested in (2+1)D physics. As is customary for these types of problems (see, e.g., Ref. [41]), we adopt a naive dimensional-regularization prescription according to which all Dirac matrices anticommute [66] and spinor traces over products of an odd number of Dirac matrices vanish [67]. In addition to a spinor index, the field ψ carries a flavor index. With the dimensional-regularization prescription just mentioned, perturbative results only depend on the total number of (complex) fermionic degrees of freedom, i.e., the dimension of the chosen representation of the Dirac algebra, times the number of flavors. We will present our results in terms of the number N of flavors of two-component Dirac fermions (i.e., the number of linear band crossing points at the Fermi level in a condensed matter system), but they can alternatively be interpreted as pertaining to N f = N/2 flavors of four-component Dirac fermions when N is even.
We consider the cases n = 1, 2, 3, corresponding to the chiral Ising, XY, and Heisenberg GNY models, respectively [40,41]. The form of the Yukawa coupling L ψφ in Eq. (1) differs in each case. In the chiral Ising GNY model [68], a single real scalar φ couples to the fermion mass iψψ, with coupling strength h. The Yukawa coupling in the chiral XY GNY model can be formulated in different but equivalent ways, depending on the choice of spinor representation. In the fourcomponent representation, the Yukawa coupling can be written as a coupling to both the ordinary mass iψψ and an axial mass ψγ 5 ψ, and is equivalent to the Nambu-Jona-Lasinio model [69]. Here, one utilizes a four-dimensional representation γ µ , µ = 0, 1, 2, 3 of the SO(4) Clifford algebra, and γ 5 = γ 0 γ 1 γ 2 γ 3 . In a different spinor representation [70], the model can be written as a coupling to a Majorana mass, where the O(2) order parameter φ = (φ 1 , φ 2 ) is expressed as a complex scalar field φ = φ 1 + iφ 2 . Finally, the Yukawa coupling in the chiral Heisenberg GNY model is: where σ = (σ 1 , σ 2 , σ 3 ) forms a spin-1/2 representation of the SU (2) algebra.
For different values of N , the O(n) GNY models introduced above describe a variety of quantum phase transitions in (2+1)D condensed matter systems [42]. For N = 4 (spinful fermions) and N = 2 (spinless fermions), the chiral Ising GNY model (n = 1) describes a transition from a Dirac semimetal to an insulator with charge-density-wave order on the honeycomb lattice [71]. For N = 1, the model describes a ferromagnetic transition on the surface of a 3D topological insulator [72]. For N = 1/2, which can be interpreted as a model containing a single flavor of two-component Majorana fermions, the model describes the time-reversal symmetry-breaking transition on the surface of a 3D topological superconductor [73], which exhibits an emergent N = 1 supersymmetry [73][74][75]. Turning to the chiral XY GNY model (n = 2), the cases N = 4 and N = 2 describe a quantum phase transition from a Dirac semimetal (spinful or spinless, respectively) to an insulator with Kekulé valence-bond-solid (VBS) order on the honeycomb lattice [27,76], or to an insulator with columnar VBS order on the π-flux square lattice [77]. The spontaneously broken symmetries in those examples are discrete Z 3 and Z 4 point group symmetries, respectively, but those anisotropies are irrelevant perturbations at the O(2)-symmetric GNY fixed point, at least in the large-N limit [78,79]. However, in those VBS realizations of chiral XY GNY criticality, spatial randomness necessarily couples linearly to the VBS order parameter: it thus acts as random-field disorder, which destroys the d = 2 critical point [80]. Alternatively, the chiral XY GNY model also describes a semimetal-superconductor transition in a system with N two-component Dirac fermions (N = 4 for spinful fermions on the honeycomb lattice [27]), in which case the U (1) ∼ = SO(2) symmetry is exact and random-field disorder is forbidden by conservation of particle number. For N = 1, the model describes a superconducting transition on the surface of a 3D topological insulator, and exhibits an emergent N = 2 supersymmetry [27,28,73,75,81,82]. Finally, for N = 4 the chiral Heisenberg GNY model (n = 3) describes the transition from a Dirac semimetal to an insulator with antiferromagnetic spin-density-wave order on the honeycomb lattice [71].
We model quenched random-mass disorder by randomness in the scalar mass squared, r(x) = r 0 + δr(x), where δr(x) is a Gaussian random variable of zero mean and correlation function [39]: where · · · denotes disorder averaging. The uniform part r 0 is the tuning parameter for the transition, and ∆ and v are the short-range and long-range correlated disorder strengths, respectively. Even when considering initial conditions for the RG with only long-range correlated disorder, ∆ = 0, short-range correlated disorder is generated perturbatively already at one-loop order, see Eq. (42), and should be kept in the space of couplings. By contrast, long-range correlated disorder cannot be generated perturbatively from short-range correlated disorder, see Eq. (43). We use the replica trick to average over disorder [2], which induces an effective two-body interaction, where a, b = 1, . . . , m are replica indices, and the replica limit m → 0 is to be taken at the end of the calculation. As for the superfluid-Mott glass transition [83], randomness in the scalar mass squared preserves the exact particle-hole symmetry of the clean GNY action (1).

III. RG IN THE TRIPLE EPSILON EXPANSION
We first briefly recapitulate the idea of the double epsilon expansion for QCPs perturbed by quenched short-range correlated disorder, first focusing on the purely bosonic random-mass O(n) vector model [16][17][18][19]. In d = 4 − spatial and τ imaginary time dimensions, the order parameter field φ has engineering dimension ∆ φ = (2 − + τ )/2. The couplings λ 2 and ∆ thus have mass dimension − τ and , respectively, and a controlled perturbative RG analysis can be performed by treating and τ as small parameters. For n > 1, a stable DFP with λ 2 * ∼ O( , τ ), ∆ * ∼ O( , τ ) on the critical hypersurface r = 0 is found at one-loop order, with critical exponents [84]: For n = 2, and extrapolating τ to 1 and to 2 or 1, relevant to the boson superfluid-Mott glass transition in (2+1)D and (3+1)D, respectively, one obtains exponents in reasonable agreement with those found in numerical Monte Carlo (MC) simulations (Table I). Ref. [26] observed that the double epsilon expansion can also be applied to short-range correlated random-mass GNY models: the fermion field ψ has engineering dimension ∆ ψ = (3 − + τ )/2, thus the Yukawa coupling h has mass dimension ( − τ )/2 and can also be treated perturbatively. To the difference of the bosonic model, however, one must enlarge the space of running couplings to include the relative velocity c = c f /c b , and ensure that the beta function for this parameter also vanishes at the DFP. (For a disordered system with a single field, the flow of the velocity, e.g., c b for the bosonic O(n) model, can be absorbed in the definition of z, provided the disorder strength flows to a fixed-point value ∆ * [3,4].) In the presence of long-range correlated disorder, we see from Eq. (9) that the coupling constant v has mass dimension 4 − α at the Gaussian fixed point. While for generic α < d < 4 this coupling is strongly relevant, if we set α = 4 − δ and treat δ as a small parameter long-range correlated disorder is only slightly relevant and can be treated perturbatively [39]. This forms the basis of a triple expansion in , τ , δ [64], which thus far has only been applied to bosonic systems. Below we employ this triple epsilon expansion to study the GNY models with both short-range and long-range correlated random-mass disorder.
In the presence of three epsilon-like parameters, the nature of the RG fixed points and their stability depend on two ratios, e.g., / τ and δ/ τ . We restrict our consideration to / τ = 2, which in the limit τ → 1 corresponds to (2+1)D systems. Regarding the δ/ τ ratio, we consider the range 0 < δ/ τ < 4. For δ < 0, long-range correlated disorder is irrelevant at the Gaussian fixed point, and for δ/ τ > 4, the long-range disorder correlations (8) with α = 4 − δ would have the unphysical feature of increasing rather than decaying with distance in the limit τ → 1.

A. Bare vs renormalized actions
We now outline the basic steps of the RG procedure using as example the chiral XY GNY model studied in Ref. [26], but with long-range correlated disorder (9). For the chiral Ising and Heisenberg GNY models, the number of components of the order parameter and the form of the Yukawa coupling change [see Eqs. (4-7)], but the relations (18) between bare and renormalized couplings, and the formal expressions (19)(20)(21)(22)(23)(24) for the beta functions in terms of the anomalous dimensions (14), remain the same.
As in Refs. [85,86], we rescale the time coordinate as well as the fermion and boson fields, and redefine the couplings in the action (1)(2)(3)(4)(5)(6)(7)(8)(9), to eliminate the velocities c f and c b in favor of the dimensionless ratio c 2 = (c f /c b ) 2 , which then appears in front of the time derivative term for the boson field. The replicated bare action for the random-mass chiral XY GNY model is then: where a, b = 1, . . . , m are replica indices, and the corresponding renormalized action is: where µ is a renormalization scale. Due to the anisotropy between space and time, we set x B = x and τ B = ητ , and matching the bare and renormalized kinetic terms for the fermion we find that η = Z 2 /Z 1 . Defining the anomalous dimensions: we find that the dynamic critical exponent z = µ(d ln τ /dµ) [87] is given by: The fermion and boson fields are multiplicatively renormalized, and the fermion and boson anomalous dimensions, η ψ = µ(d ln Z ψ /dµ) and η φ = µ(d ln Z φ /dµ), are given by: Comparing Eqs. (12) and (13), we obtain relations between the bare and (dimensionless) renormalized couplings, Using the fact that the bare couplings do not depend on the renormalization scale µ, we find the RG beta functions β g ≡ µ(dg/dµ), g ∈ {c 2 , λ 2 , h 2 , ∆, v}, to be: From Eq. (24), we find the inverse correlation length exponent [88],

B. Renormalization constants
We calculate the renormalization constants Z i , i = 1, . . . , 8, r at one-loop order in the modified minimal subtraction (MS) scheme with dimensional regularization in 4 − space and τ time dimensions. The relevant Feynman rules and diagrams are shown schematically in Figs. 2 and 3, respectively. The fermion and boson propagators are given by: where I, J = 1, . . . , N and i, j = 1, . . . , n are fermion flavor and O(n) indices, respectively, and For the chiral XY GNY model (n = 2), the diagrams in the clean limit or containing only shortrange correlated disorder vertices were already computed in Ref. [26]; these results are also easily adapted to n = 1 and n = 3. The new diagrams containing long-range correlated disorder vertices are evaluated explicitly in Appendix A for n = 1, 2, 3. Unlike the standard epsilon expansion in 4 − dimensions, in the triple epsilon expansion one-loop diagrams contain simple poles not only in , but also in − τ , δ, and 2δ − . We obtain the following renormalization constants: We have rescaled the couplings according to g/(4π) 2 → g, g ∈ {λ 2 , h 2 , ∆, v, r}, and, as in Ref. [26], we define the dimensionless function, plotted in Fig. 4. At one-loop order there is no renormalization of the Yukawa vertex for the chiral XY GNY model, i.e., the diagram in Fig. 3(f) vanishes for n = 2 [see Eq. (33)], which is easily seen from the form (6) of the Yukawa coupling. We also see from the last term in Eq. (34) that short-range correlated disorder is generated at one-loop order from long-range correlated disorder, via the diagram in Fig. 3(m). By contrast, long-range correlated disorder cannot be generated perturbatively from short-range correlated disorder.
Since h 2 * will be O( , τ ) at one-loop order, as can already be seen from Eq. (41), for a consistent treatment we have to discard the τ (z − 1) terms in the fermion and boson anomalous dimensions.

IV. FIXED POINTS AND CRITICAL EXPONENTS
In Sec. IV A, we discuss the fixed points of the flow equations (39)(40)(41)(42)(43). Depending on their stability, which is analyzed in Sec. IV B, these are bona fide critical points (no relevant direction) or multicritical points (one or more relevant directions). Here, the number of relevant directions refers to the number of such directions on the critical hypersurface, since the tuning parameter r for the transition (see Sec. II) is a relevant direction at all fixed points. As mentioned in Sec. III, we fix = 2 τ , with the extrapolation τ → 1 corresponding to 2+1 dimensions. Throughout the paper, we evaluate quantities such as fixed-point couplings, RG eigenvalues, and critical exponents as a function of the control parameters N 1 and δ = 4 − α ∈ [0, 4], where the latter parameter is to be understood as the ratio δ/ τ evaluated at τ = 1.
We next turn to DFPs, for which ∆ * and/or v * are nonzero. To be physical, all fixed points must obey the following conditions [39]: At fermionic DFPs with h 2 * > 0, the condition β c 2 = 0 together with Eq. (49) further implies that c 2 * > 1. From Eq. (39), we find that at a fermionic fixed point, Equation (49) implies that the right-hand side of this equation is positive. From Fig. 4 and Eq. (37), Thus for the left-hand side of Eq. (50) to be positive also we must have c 2 * > 1. (At a clean fermionic fixed point, the left-hand side must vanish, which can only happen for c 2 * = 1.)

Fixed points with short-range correlated disorder
We first focus on DFPs with ∆ * = 0 and v * = 0, which we term short-range disordered fixed points (SDFPs). From Eq. (41) we find that h 2 . When the fixed-point value of the Yukawa coupling is zero, we reproduce the results of Refs. [17][18][19] for the purely bosonic O(n) vector model with random-mass disorder. For n = 1, there is an accidental degeneracy in the system of equations β λ 2 = 0, β ∆ = 0. The degeneracy is lifted at two-loop order, giving rise to a DFP with λ 2 * , ∆ * ∼ O( √ τ ), for a finite ratio / τ [17]. Our focus, however, is on fermionic DFPs with nonzero h 2 * . We find two fermionic SDFPs for n = 2, 3: where . The chiral XY case (n = 2) was discussed in our earlier work [26]: the fixed-point couplings λ 2 * , h 2 * , and ∆ * are nonnegative, and thus physical, for all N 1. At N = 1, SDFP2 merges with the clean GNY fixed point (48), while SDFP1 runs off to infinity as it is impossible to satisfy β c 2 = 0. (Note that for n = 2, SDFP1,2 here correspond to DFP1,2 in Ref. [26] for N < 4 and to DFP2,1 for N > 4.) In the chiral Heisenberg case (n = 3), the discriminant D S 0 for N N D ≈ 27.856, and the SDFPs (51) are physical only for N > N D .
In the chiral Ising case (n = 1), as previously mentioned the RG equations for λ 2 and ∆ become degenerate for zero Yukawa coupling, and we find only one solution at order O( , τ ) for h 2 * = 0: This SDFP is physical for N 6, and merges with the clean GNY fixed point at N = 6. There is in principle the possibility of an additional SDFP at two-loop order with λ 2 * , ∆ * ∼ O( √ τ ), as in the bosonic case, and h 2 * ∼ O( τ ). We show in Appendix B that this cannot happen, because it is impossible to satisfy the equation β c 2 = 0. We also note that this excludes the possibility of a physical SDFP for the N = 1/2 chiral Ising GNY model, which in the clean limit flows to a conformal field theory with emergent supersymmetry [74,75], the N = 1 Wess-Zumino model.  Fig. 8).

Fixed points with long-range correlated disorder
We now turn to DFPs with v * = 0, which we dub long-range disordered fixed points (LDFPs). For vanishing h 2 * , the purely bosonic random-mass O(n) vector model for n > 1 was studied in the triple epsilon expansion in Ref. [64], where LDFPs were found. For n = 1, long-range correlated disorder lifts the previously mentioned degeneracy in the system of fixed-point equations. For nonzero h 2 * = τ /(N + 4 − n), we find two fermionic LDFPs in all three GNY universality classes, n = 1, 2, 3: where The discriminant D L is nonnegative, and thus the fixed-point couplings real, for either: or: In addition to being real, the fixed-point couplings (53-55) must obey the conditions (49). By contrast with the SDFPs (51-52), which are physical above a certain critical value of N that is independent of δ, the LDFPs are physical only in complicated regions of the N -δ plane that possess several disconnected components and/or curved boundaries. Since the fixed-point couplings (53)(54)(55) do not depend explicitly on c 2 * , we first assume a physical solution for c 2 * exists, and discuss how the remaining conditions delimit those nontrivial regions.
• ∆ * + v * 0: For LDFP1, i.e., Eqs. (53)(54)(55) with + √ D L , the condition is satisfied for different regions of the N -δ plane depending on n: For LDFP2, i.e., Eqs. (53)(54)(55) with − √ D L , we have: Here, and N 2 is the value of N , which depends on n, at which δ 1 = δ D . For N < N < N 2 , δ D < 0, in which case [0, δ D ] denotes the empty set. We use the same notational convention whenever the left limit of the interval is greater than the right one.
• v * 0: For LDFP1, we have the following constraints depending on the value of n: For LDFP2, we have: We further define and N 3 is the n-dependent value of N at which δ D = δ 4 .
For a given GNY symmetry class n, the intersection of all those conditions defines regions in the N -δ plane in which the various fixed points discussed are physical, and over which fixed-point properties are plotted throughout the paper. We now return to the question of whether a physical solution c 2 * to the nonlinear equation β c 2 = 0 exists for the LDFPs (53-55). We solve this equation numerically. For n = 1 and n = 3, we find a unique solution everywhere in the physical regions of the N -δ plane. For n = 2, we likewise find a unique physical solution in the physical regions, but for LDFP1 computations become increasingly difficult upon approach to the point N = 1, δ = 4, where c 2 * grows rapidly. Since exactly at this point LDFP1 coincides with SDFP2, and SDFP2 does not admit a solution to β c 2 = 0 for N = 1 [26], we conjecture that c 2 * gradually runs off to infinity as the point N = 1, δ = 4 is approached. Summarizing, we thus find that for all three GNY symmetry classes, a unique solution c 2 * > 1 exists for the LDFPs (53-55) everywhere inside the physical regions (49) of the N -δ plane. As mentioned previously, h 2 * and c 2 * together determine the dynamic critical exponent z at those fixed points (Sec. IV C, Figs. 9-11).

B. Linear stability analysis
We now investigate the stability properties of the physical fixed points. All bosonic fixed points (i.e., with h 2 * = 0) are unstable with respect to the h 2 direction. Additionally, for all models, the Gaussian fixed points are unstable with respect to all other directions, and the Wilson-Fisher fixed points are unstable with respect to both short-range and long-range correlated disorder. The stability properties of the bosonic DFPs in the absence of Yukawa coupling have been discussed previously in Refs. [17][18][19]64].
At all fermionic fixed points (i.e., with h 2 * = 0), the h 2 direction is irrelevant. Additionally, we find that ∂β c 2 /∂c 2 is positive at all such fixed points. Since β c 2 is the only beta function in which c 2 appears, this means c 2 is also an irrelevant direction. We can thus exclude h 2 and c 2 from RG flow considerations and investigate stability within the three-dimensional subspace with fixed h 2 * and c 2 * of the full five-dimensional space of couplings. We compute the eigenvalues y of the stability matrix M gg ≡ −∂β g /∂g , g, g ∈ {λ 2 , ∆, v}, defined such that y > 0 (y < 0) corresponds to a relevant (irrelevant) direction.

Stability of the clean fixed point
We first focus on the clean GNY fixed point (48), which for the rest of the paper we refer to as the CFP. The RG eigenvalues at the CFP are: and are associated with eigenvectors with nonzero projections along the λ 2 , ∆, and v directions, respectively. The eigenvalue y 1 is negative and thus irrelevant for all n and N . For the flow of short-range correlated disorder (y 2 ), we discuss the three GNY symmetry classes in turn.
Finally, long-range correlated disorder (y 3 ) is irrelevant for δ less than δ 1 , which is defined in Eq. (62). At generic points along the curve δ = δ 1 in the N -δ plane, one of the LDFPs merges with the CFP, and long-range correlated disorder crosses marginality. At the special point N = N 2 along this curve, the two LDFPs (53-55) coincide with one another (and with the CFP).

Stability of short-range disordered fixed points
We now consider the SDFPs of Sec. IV A 1. We begin with the unique SDFP (52) in the chiral Ising class (n = 1), which is physical only for N 6. Long-range correlated disorder is irrelevant at this fixed point provided that δ is less than δ 5 , which is defined in Eq. (71). Along the curve  Similarly to the chiral Ising case, long-range correlated disorder is irrelevant at SDFP1 (SDFP2) provided that δ < δ 3 (δ < δ 4 ), with δ 3 , δ 4 defined in Eqs. (69)(70). The curves δ = δ 3 and δ = δ 4 correspond to the merger of the corresponding SDFP with one of the LDFPs. When δ 3 = δ 4 , the Besides long-range correlated disorder, the other two directions are irrelevant at SDFP1, thus it is a genuine critical point for δ < δ 3 . By contrast, one of those two directions is relevant at SDFP2, thus the latter is a multicritical point. For the chiral XY and Heisenberg models, and for sufficiently large N , the two irrelevant eigenvalues at SDFP1 with eigenvectors in the λ 2 -∆ plane form a complex conjugate pair. SDFP1 is then a fixed point of focus type, with spiraling flows near the fixed point. In the XY case, this happens for N > 32 5 = 6.4, while for the Heisenberg case, this happens for N > 28.087. Critical properties in this case are subject to oscillatory corrections to scaling [15,26].

Stability of long-range disordered fixed points
We finally turn to the stability of the LDFPs of Sec. IV A 2. The eigenvalues of the stability matrix depend on N and δ in a complicated way, and we compute them numerically. In Figs. 5-7, we characterize the stability of the two LDFPs in terms of their number of relevant/irrelevant eigenvalues, for each GNY symmetry class. Eigenvalues are real unless otherwise specified; since the stability matrix is real, complex eigenvalues necessarily appear in complex-conjugate pairs, and imply focus-type behavior as discussed above. For all three GNY symmetry classes, the two LDFPs merge along the curve δ = δ D in the N -δ plane, where the discriminant D L vanishes. In the Ising case (Fig. 5), both LDFPs have at least one relevant eigenvalue on the critical hypersurface and are thus multicritical points (for N = 1/2, only LDFP1 is physical, for δ 1 ≈ 1.143 < δ < 2). In the XY and Heisenberg cases (Figs. 6-7), LDFP1 exists in regions (V and VI) in the N -δ plane with no relevant eigenvalues, and is thus a bona fide critical point in those regions. LDFP2 is always multicritical.

C. Critical exponents
Universal critical exponents at the newly found fermionic DFPs can be computed from Eqs. (44)(45)(46)(47) using the fixed-point couplings found in Sec. IV A 1 and Sec. IV A 2. At the present one-loop order, the fermion η ψ and boson η φ anomalous dimensions depend only on h 2 * , which is the same at all fermionic fixed points. Thus their values at the DFPs are the same as those for the clean chiral GNY universality classes [41]: η ψ = n τ /[2(N + 4 − n)] and η φ = N τ /(N + 4 − n). At higher loop order the anomalous dimensions are expected to differ at the different fermionic fixed points.
Using Eq. (45), the dynamic critical exponent z at the fermionic DFPs is given by and thus depends on the fixed-point velocity parameter c 2 * . The latter is a universal function of N and δ for a given DFP but must be computed numerically; we plot the resulting value of z extrapolated to 2+1 dimensions ( τ → 1) in Fig. 8 for the SDFPs and in Figs. 9-11 for the LDFPs. Since c 2 * > 1, and thus f (c 2 * ) > 1 2 , at all fermionic DFPs (see Sec. IV A), such DFPs necessarily have z > 1. This is in agreement with the general expectation that weak disorder increases z [89]; Refs. [3,4] also derive the leading-order result z − 1 ∝ ∆ * > 0 at SDFPs obtained by perturbing a conformally invariant QCP with weak short-range correlated disorder. Here we find z > 1 at LDFPs as well. The inverse correlation length exponent ν −1 , determined from Eq. (44), is the RG eigenvalue associated with the relevant direction r which tunes across the symmetry-breaking transition. For a bona fide critical point, ν controls the divergence of the correlation length ξ at the transition r = 0 via ξ ∼ r −ν . For multicritical points with additional relevant directions g 1 , g 2 , . . . on the critical hypersurface with real, positive eigenvalues y 1 , y 2 , . . ., the correlation length behaves near the transition as ξ(r, g 1 , g 2 , . . .) = r −ν ξ(g 1 /r νy 1 , g 2 /r νy 2 , . . .), where ξ(x 1 , x 2 , . . .) is a universal scaling function [90]. Complex-conjugate eigenvalues produce a scaling function with oscillatory behavior. At all LDFPs in all three GNY symmetry classes, we find ν −1 = 2 − 1 2 δ, which alternatively can be written as ν = 2/α, with α = 4 − δ the exponent controlling long-range disorder correlations in Eq. (8). This superuniversal behavior was also found at long-range correlated bosonic DFPs and explained by Weinrib and Halperin [39]. Consider a LDFP with correlation length exponent ν(α) in a system with disorder of the type (8). If one further perturbs this fixed point with disorder correlated according to |x − x | −α + such that α + > α, the original asymptotic critical behavior should remain the same, as we expect it is controlled by the longest-range part of the disorder. Conversely, if the perturbation is of the form |x − x | −α − with α − < α, this falls off more slowly than the original disorder, and the original critical behavior should be unstable. Assuming α, α ± < d and applying the modified Harris criterion for long-range correlated disorder, we find ν(α) > 2/α + and ν(α) < 2/α − , for all α − < α < α + . Choosing α ± = α ± ε and taking the limit ε → 0 + , we obtain ν(α) = 2/α. The exponent ν for the SDFPs can likewise be calculated directly from Eq. (44), and we obtain ν −1 = 2 − 1 2 δ 5 for the chiral Ising SDFP, with δ 5 defined in Eq. (71). In light of the result above for ν −1 at LDFPs, this is consistent with the fact that the n = 1 SDFP coalesces with one of the LDFPs at δ = δ 5 . Similarly, for both the chiral XY and Heisenberg models we find that SDFP1 has ν −1 = 2 − 1 2 δ 3 and SDFP2 has ν −1 = 2 − 1 2 δ 4 , with δ 3,4 defined in Eqs. (69)(70). As previously mentioned, the curves δ = δ 3 (δ = δ 4 ) correspond to the merger of SDFP1 (SDFP2) with a LDFP. We plot ν −1 at SDFPs for all three GNY models in Fig. 12, including ν −1 at the clean GNY critical point for comparison. and h 2 assuming their fixed-point values. Since in the chiral Ising case all physical fixed points are multicritical, and for the sake of simplicity, we restrict our attention to the chiral XY and Heisenberg symmetry classes, which exhibit the most interesting phenomena.

A. Transcritical and saddle-node bifurcations
We have already mentioned a number of instances in which two fixed points collide as N or δ are varied. We observe two distinct kinds of bifurcations associated with a collision of two fixed points: the transcritical bifurcation and the saddle-node bifurcation.
The transcritical bifurcation [ Fig. 13(a)] is a bifurcation at which a stable fixed point and an unstable fixed point pass through each other, exchanging their stability properties, but without annihilating [43]. An example of this bifurcation is the merging of the two chiral XY SDFPs (51) as N is varied through N = 4. (There is "exchange" of fixed points provided we track individual fixed points on smooth trajectories, as opposed to their arbitrary definition as SDFP1 and SDFP2 in Eq. (51).) Unlike the saddle-node bifurcation discussed below, the two fixed points remain real before and after the bifurcation. At the transcritical bifurcation, the beta function (and associated RG flow) is not only marginal, but its derivative with respect to the control parameter, here N , must vanish as well. Other examples of this bifurcation include the collision of SDFPs with the CFP (at N = 1 for the chiral XY SDFP2), of LDFPs with the CFP (along the curve δ = δ 1 in the N -δ plane), or of SDFPs with LDFPs (curves δ = δ 3 and δ = δ 4 ). At these latter bifurcations, one of the DFPs becomes unphysical, by either ∆ * , v * , or ∆ * + v * going through zero and becoming negative. However, since the other fixed point remains physical and thus real, this unphysical fixed point necessarily remains real also (for another RG example of this scenario, see Ref. [91]). Thus the bifurcation is distinct from the saddle-node bifurcation, which we now discuss.
The saddle-node bifurcation [ Fig. 13(b)] is a bifurcation at which a stable fixed point and an unstable fixed point merge, leading to marginal behavior as above, but subsequently disappear into the complex plane. This typically happens for a pair of fixed points with critical couplings g * ± ∝ A ± √ D, such that the discriminant D continuously goes through zero at the bifurcation and then becomes negative. Both pairs SDFP1,2 and LDFP1,2 are of this type. The saddle-node bifurcation is accompanied by the characteristic phenomenology of walking RG or quasi-critical behavior [44]; we now explain how this manifests itself in the current problem. Focusing on the example above of the annihilation of LDFPs in the chiral XY and Heisenberg GNY models, we first consider a situation where δ is slightly above δ D . Small regions in the N -δ plane exist such that both LDFPs are physical, with LDFP1 a stable sink-type fixed point (region V) and LDFP2 a multicritical point with one relevant direction (region I). LDFP2 is only physical provided δ < δ 1 [see Eq. (68)], which implies that the CFP is stable (Sec. IV B 1). For this type of region, numerical studies of the RG flow show that RG trajectories with initial conditions near LDFP2 end up at either LDFP1 or the CFP. We thus consider a curvilinear coordinate system such that one of these coordinates, g, passes through all three fixed points [ Fig. 14(a)]. In this section only, we define the infrared (Wilsonian) beta function β(g) ≡ dg/d , where grows towards the infrared. Denoting by g * the common fixed-point coupling of LDFP1 and LDFP2 at the bifurcation δ = δ D , we assume that for δ near δ D and g near g * , β(g) can be well approximated by a quadratic function, β(g) ≈ A(δ) + B(δ)(g − g * ) + C(δ)(g − g * ) 2 . Since β(g * ) = ∂β(g * )/∂g = 0 and ∂ 2 β(g * )/∂g 2 < 0 at δ = δ D , we have A(δ D ) = B(δ D ) = 0 and C(δ D ) ≡ −κ < 0. For δ = δ D + ε with ε small, β(g) should have two real zeros that approach g * as ε → 0 + . Expanding A(δ), B(δ), and C(δ) in powers of ε, we find at leading order a pair of zeros of the form g * ± bε/κ with b ≡ A (δ D ), which are real provided that b > 0, and form a complex-conjugate pair when ε < 0 (δ < δ D ). The beta function thus approximately assumes the form β(g) ≈ b(δ − δ D ) − κ(g − g * ) 2 , illustrated in Fig. 14(b), and considered in Ref. [44].
We now take δ = δ D − ε with ε > 0 small, and consider an RG trajectory with initial coupling g UV > g * and "flow velocity" β(g UV ), which is generically not small. As g approaches g * from above, the flow velocity decreases considerably (i.e., the running coupling "walks"), since β(g * ) ≈ −bε is small. This walking behavior persists until g * − g becomes on the order of bε/κ, after which the coupling starts "running" again. This determines a characteristic RG time ∆ insensitive to the initial condition g UV of the flow. Approximating β(g) ≈ β(g * ) ≈ −bε as constant during the walk, we have β(g * ) ≈ ∆g/∆ ∼ bε/κ/∆ , and thus ∆ ∼ 1/ √ κbε. Alternatively, we may integrate the equation dg/d = β(g) from g UV at UV to g IR < g * at IR . Under the condition |g UV,IR − g * | bε/κ, the result of this integration is insensitive to the precise values of g UV and g IR , and we obtain ∆ ≡ IR − UV = π/ √ κbε. In turn, this RG time determines a characteristic infrared length scale L * = L IR = L UV e ∆ , where we can take L UV ∼ a to be on the order of a microscopic lattice constant a. We obtain: as δ approaches δ D from below. The exponential inverse-square-root divergence, reminiscent of the divergence of the correlation length at the Kosterlitz-Thouless transition [92], is characteristic of the saddle-node bifurcation [44]. The existence of this exponentially large length scale L * a allows for a crossover between two distinct physical regimes [ Fig. 14(c)]. On intermediate length scales a L L * , RG trajectories dwell for an extended period of RG time near g = g * , and we have quasi-critical behavior controlled by a complex pair of LDFPs with real part near g * . This quasi-critical regime is characterized by approximate power-law scaling and drifting (i.e., scaledependent) exponents [55]. On the largest length scales a L * L, the transition is controlled by the true infrared fixed point, the CFP, with genuine scale invariance.

B. Supercritical Hopf bifurcation and limit-cycle fermionic quantum criticality
The third type of bifurcation we observe is the supercritical Hopf bifurcation [ Fig. 13(c)]. This bifurcation occurs as one passes from region VI (blue region) to region IV (purple region) in both the chiral XY [ Fig. 6(b)] and Heisenberg [ Fig. 7(a)] models. For instance, one can consider keeping N fixed and tuning δ (black arrow in those figures). In region VI (δ < δ c,1 ), LDFP1 is a stable-focus fixed point with two complex-conjugate irrelevant eigenvalues, i.e., complex-conjugate eigenvalues with a negative real part [solid red line on left part of Fig. 13(c)]. At the bifurcation (δ = δ c,1 ), the real part of those eigenvalues goes through zero and becomes positive for δ > δ c,1 . LDFP1 thus loses its stability and becomes an unstable-focus fixed point [dashed blue line on the right part of Fig. 13(c)]. At the same type, a stable limit cycle is born [solid red line on the right part of Fig. 13(c)], towards which the spiraling RG trajectories coming out of LDFP1 asymptote, and which controls the critical behavior up to a second threshold value δ c,2 to be discussed shortly. (Trajectories outside the limit cycle also spiral and asymptote to it.) To our knowledge, this is the first instance in the context of quantum phase transitions where the supercritical Hopf bifurcation [60] appears. After Ref. [37], which studied a holographic model of a critical scalar field perturbed by disorder, our result is the second example of quantum phase transition governed by a stable limit cycle; to our knowledge, it is the first example for fermionic systems. The subcritical Hopf bifurcation [60], where an unstable-focus fixed point becomes stable by giving birth to an unstable limit cycle, has been reported previously in RG studies of classical disordered systems [39,63]. The general phenomenology of critical behavior controlled by a stable limit cycle was explored in Ref. [62]. For a stable-focus critical point, spiraling trajectories manifest themselves as oscillatory corrections to scaling [15,26]. By contrast, for a transition governed by a stable limit cycle, thermodynamic quantities exhibit log-periodic scaling behavior at leading order, i.e., discrete scale invariance. For instance, we show in Appendix C that the order parameter susceptibility χ obeys the approximate scaling form: where F is a periodic function. Here ν LC and γ LC = (2 − η φ )ν LC are effective correlation-length and susceptibility exponents for the limit cycle, r is the tuning parameter for the transition, and r 0 is a nonuniversal constant. As δ is further increased past δ c,1 , the limit cycle eventually disappears at a second critical value δ c,2 , but in different ways for the chiral XY and Heisenberg GNY models. In the Heisenberg case, the Hopf bifurcation of Fig. 13(c) occurs again but in reverse: the limit cycle shrinks to a point, which becomes the stable-focus LDFP1 of region VI. In the XY case, our numerical studies suggest that at least for some values of N , the limit cycle is destroyed at δ = δ c,2 (still within region IV) by colliding with the CFP and SDFP2, which are both saddle points in this regime [see Fig. 15(c)]. This is a possible example of heteroclinic bifurcation [93], whose detailed study we reserve for future work.

C. Schematic phase diagrams
From the knowledge of the stability properties of the various fixed points and limit cycles, and numerical investigation of the RG flow connecting those different critical manifolds, schematic In the Heisenberg case, the transition is controlled by a stable limit cycle (LC) for generic bare values of the short-range correlated (∆) and long-range correlated (v) disorder strengths. In the XY case, the transition is controlled by the limit cycle for weak short-range disorder and by a disordered fixed point (SDFP1) for strong short-range disorder.
phase diagrams can be constructed analogously to those in Ref. [26]. For given values of N and δ, we focus on the critical hypersurface (r = 0) and ask how the universality class of the transition depends on the bare couplings in the Lagrangian, which determine the initial conditions for the infrared RG flow. We consider a scenario in which the interaction parameters h and λ are fixed, and vary the two types of disorder, ∆ and v. Since the number of possibilities is very large, given the complexity of the stability/physicality regions, we focus on the two most interesting regions: those which contain the instances of limit-cycle quantum criticality discussed in the previous section.
We first focus on region IV in the chiral Heisenberg GNY model [see Fig. 7(a)]. For generic points in this region (e.g., for δ c,1 < δ < δ c,2 ), one has δ > δ D and δ > δ 1 . Furthermore, we assume N < N D ≈ 27.856. From Sec. IV B 1, we conclude that the CFP has two irrelevant directions in the λ 2 -∆ plane, but that long-range correlated disorder v is relevant, since δ > δ 1 . SDFP1,2 are both unphysical, since N < N D , and LDFP2 is unphysical as well. As seen in the previous section, LDFP1 is of unstable-focus type, with spiraling flow towards a stable limit cycle. The resulting RG flow is illustrated schematically in Fig. 15(a). Consequently, at least for sufficiently small bare values of the disorder, the transition is controlled by limit-cycle quantum criticality for generic disorder [ Fig. 15(b)]. If long-range correlated disorder is turned off completely, the transition reverts back to the clean chiral Heisenberg GNY universality class.
We now turn to region IV in the chiral XY GNY model [see Fig. 6(b)], assuming δ c,1 < δ < δ c,2 . As in the previous case, we generically have δ > δ D , δ > δ 1 , and also δ < δ 4 . As in the Heisenberg case, the CFP has two irrelevant directions in the λ 2 -∆ plane, but v is relevant. There are now nontrivial SDFPs, whose stability was discussed in Sec. IV B 2. For SDFP1, λ 2 and ∆ are both irrelevant, and v is irrelevant as well, since δ < δ 4 < δ 3 . For SDFP2, v is irrelevant since δ < δ 4 , but there is one relevant direction with nonzero ∆ projection. LDFP2 is unphysical, and LDFP1 is an unstable focus with flow towards a stable limit cycle. The resulting RG flow is schematized in Fig. 15(c), and the corresponding phase diagram in Fig. 15(d). For weak ∆, the transition is governed by the limit cycle, but for sufficiently strong ∆, the transition is controlled by a disordered fixed point, SDFP1. CFP and SDFP2 appear as multicritical points.

VI. CONCLUSION
In summary, we have performed a comprehensive study of the three classes of chiral GNY models most relevant for symmetry-breaking quantum phase transitions in (2+1)D gapless Dirac matter-the chiral Ising, XY, and Heisenberg GNY models-in the presence of quenched shortrange and long-range correlated random-mass disorder. Using a controlled triple epsilon expansion below the upper critical dimension for these models, we have found several disordered infrared fixed points characterized by finite short-range and/or long-range correlated randomness, and for which we computed critical exponents. The Boyanovsky-Cardy and quantum Weinrib-Halperin fixed points, while present, are destabilized by the Yukawa interaction in favor of new disordered fermionic QCPs, at which the strength of this interaction remains nonzero in the infrared. Besides local stability, using numerical and analytical approaches we analyzed bifurcations of the RG flow. We found instances of the familiar fixed-point annihilation scenario, which can here be tuned by a genuinely continuous variable-the exponent controlling the algebraic decay of disorder correlations-and with which is associated a parametrically large crossover length scale L * that separates a disordered quasi-critical regime (L L * ) from a clean regime in the deep infrared (L L * ). We also uncovered instances of the transcritical bifurcation, at which fixed points exchange their stability, and the more exotic supercritical Hopf bifurcation. The latter was accompanied by the emergence of a stable limit cycle on the critical hypersurface, thus producing the first instance of fermionic quantum criticality with discrete scale invariance.
Several avenues present themselves for future research. The relative paucity of disordered fixed points found in the chiral Ising class as compared to its continuous-symmetry counterparts, and in fact, the complete absence of bona fide critical points in this class, is in agreement with the conjecture by Motrunich et al. [22] that all discrete symmetry-breaking transitions in (2+1)D dis-ordered systems should fall in the infinite-randomness universality class. Since infinite-randomness fixed points are not accessible to perturbative RG methods, nonperturbative numerical studies of Ising transitions of interacting Dirac fermions with quenched randomness are desirable, e.g., using quantum Monte Carlo methods [94] or, possibly, incorporating fermions into (2+1)-dimensional adaptations of the strong-disorder RG method [95]. In the presence of gapless Dirac fermions strongly coupled to bosonic order parameter fluctuations, rare-region effects [96,97]-which dominate the low-energy physics at infinite-randomness fixed points-may however lead to a different strong-disorder phenomenology than that found in local bosonic models [98].
Besides the pure GNY universality classes, relevant to symmetry-breaking transitions in systems of itinerant Dirac electrons, our method of analysis may also provide a point of entry to study the effect of quenched disorder on more exotic transitions, such as those involving fractionalized phases. The algebraic or Dirac spin liquid [99][100][101][102], a quantum-disordered paramagnet with fractionalized spinon excitations, is described at low energies by (2+1)D quantum electrodynamics (QED 3 ) with N = 4 flavors of two-component gapless Dirac fermions. The effect of quenched disorder on QED 3 itself was studied recently [87,[103][104][105][106]; using the methods presented here, one could additionally study the effect of quenched disorder on quantum phase transitions out of the algebraic spin liquid [107]. Transitions towards conventional phases such as VBS states [79,108,109] or antiferromagnets [110][111][112], or transitions towards gapped chiral [113][114][115] or Z 2 spin liquids [91], are described by GNY theories in all three (Ising, XY, Heisenberg) symmetry classes, augmented by a coupling to fluctuating U (1) gauge fields. The effect of random-mass disorder on the critical fixed points of such QED 3 -GNY theories is an interesting topic for future research. correlated disorder.

Boson two-point function
Four diagrams contribute: Fig. 3(a-d). Diagrams (a-b) appear in the pure GNY models, and have been well studied [40,41]. Diagram (c) appears in the purely bosonic random-mass O(n) vector model [17][18][19] and contributes to δZ 3 and δZ r . Diagram (d) also contributes to δZ 3 and δZ r , and we compute it here. Its contribution to the divergent part of the effective action is: where d D k = d τ k 0 d d k, and we have discarded a term that vanishes in the replica m → 0 limit. Since we anticipate a renormalization of both the time-derivative term [3,4] and the scalar mass term, the latter being necessary to compute the correlation length exponent, we must keep the "mass squared" c 2 k 2 0 + r in the denominator. Such massive Feynman integrals can be evaluated using the Mellin-Barnes representation of hypergeometric functions [116,117]. We have: where 2 F 1 (a, b; c; z) is the Gauss hypergeometric function, and S d = 2/[(4π) d/2 Γ(d/2)]. Taking the limit δ, → 0, the hypergeometric function evaluates to a constant: 2 F 1 (0, −1; 2; z) = 1. The only divergent factor in this limit is Γ(−1 + δ 2 ) → −2/δ, and we obtain: After rescaling the couplings by (4π) 2 , we thus obtain:

Boson self-interaction
Diagrams (g) and (h) are the same as in the pure GNY models, and diagram (i) only involves short-range correlated disorder. Diagram (j) contributes to the boson self-interaction vertex: .
Since we are looking for the correction to a local four-point vertex, we can set the external momenta k, k , k to zero in the integral over the loop momentum p. Using standard Euclidean integrals, we then have, in the limit , δ → 0, Rescaling v by (4π) 2 , we obtain:

Short-range correlated disorder strength
Diagram (k) contributes to both the boson self-interaction vertex and the short-range correlated disorder vertex, and was computed before. Diagrams (l) and (m), which involve long-range correlated disorder, both contribute to the renormalization of the short-range disorder strength.
Diagrams of the type (l) give two distinct contributions, of the form: where , (A.10) . (A.11) As in the previous section, we can set k = k = 0, k = 0 in those loop integrals, which then simply reduce to Eq. (A.7). With v rescaled by (4π) 2 as before, we then obtain: Diagram (m) illustrates that long-range correlated disorder perturbatively generates short-range correlated disorder. We obtain: , (A. 13) an expression analogous to Eqs. (A.9-A.10), but with an additional factor |k − k + p| −δ in the loop integral. Again, the loop integral can be evaluated in the limit of vanishing external momenta. Using Eq. (A.6), we obtain: , (A.14) in the limit , δ → 0, and the corresponding renormalization constant is: 4. Long-range correlated disorder strength Diagrams on Fig. 3(j,l,m) also contribute to the long-range disorder coupling renormalization. Diagram (j) gives: .