Cascade of phase transitions in a planar Dirac material

We investigate a model of interacting Dirac fermions in $2+1$ dimensions with $M$ flavors and $N$ colors having the $\mathrm{U}(M)\times \mathrm{SU}(N)$ symmetry. In the large-$N$ limit, we find that the $\mathrm{U}(M)$ symmetry is spontaneously broken in a variety of ways. In the vacuum, when the parity-breaking flavor-singlet mass is varied, the ground state undergoes a sequence of $M$ first-order phase transitions, experiencing $M+1$ phases characterized by symmetry breaking $\mathrm{U}(M)\to \mathrm{U}(M-k)\times \mathrm{U}(k)$ with $k\in\{0,1,2,\cdots,M\}$, bearing a close resemblance to the vacuum structure of three-dimensional QCD. At finite temperature and chemical potential, a rich phase diagram with first and second-order phase transitions and tricritical points is observed. Also exotic phases with spontaneous symmetry breaking of the form as $\mathrm{U}(3)\to \mathrm{U}(1)^3$, $\mathrm{U}(4)\to \mathrm{U}(2)\times \mathrm{U}(1)^2$, and $\mathrm{U}(5)\to \mathrm{U}(2)^2\times \mathrm{U}(1)$ exist. For a large flavor-singlet mass, the increase of the chemical potential $\mu$ brings about $M$ consecutive first-order transitions that separate the low-$\mu$ phase diagram with vanishing fermion density from the high-$\mu$ region with a high fermion density.

Recently QCD 3 has experienced a flurry of revived attention [44][45][46][47][48][49][50][51][52]. In [49] the present authors have proposed a new random matrix theory (RMT) which, when random matrix elements are integrated out, reduces to a four-fermion model that spontaneously breaks symmetries in exactly the same way as does QCD 3 with a Chern-Simons term [44], thus extending the previous work [53]. Although RMT is a zero-dimensional theory with no gauge interactions, it provides exact descriptions of the low-lying Dirac spectrum owing to the universality of the microscopic domain [54][55][56][57][58].
In this work, we study thermodynamics and symmetry breaking of an unconventional interacting model of Dirac fermions in 2 + 1 dimensions at finite temperature and chemical potential in the large-N limit, where N denotes the number of "colors." Each fermion comes in M different flavors. This model can be viewed as a generalization of the RMT proposed in [49]. The model has three key ingredients: a repulsive interaction, an attractive interaction, and a flavor-symmetric parity-breaking mass term. Their interplay leads to a surprisingly rich phase diagram. At zero temperature and zero density, the model exhibits a spontaneous symmetry breaking patterns U(M ) → U(M − k) × U(k) with various k and experiences a sequence of first-order phase transitions, bearing a close resemblance to three-dimensional QCD [44,51]. The model reduces to a sigma model on a complex Grassmannian at low energy. At nonzero temperature or chemical potential, there appear even more exotic phases where the symmetry is broken as U(3) → U(1) × U(1) × U(1), U(4) → U(2) × U(1) × U(1), and U(5) → U(2) × U(2) × U(1), to name but a few. All these patterns show up in a single model with a few adjustable parameters.
The present work is structured as follows. In section 2, the model is defined and the thermodynamic potential is derived. In section 3, the ground state at zero temperature and density is analyzed. In section 4 the effect of nonzero temperature is considered. In section 5, a nonzero chemical potential is introduced, and the fermion number density is calculated. In section 6, phases at nonzero temperature and density are studied. It is shown that the phase structure changes dramatically, depending on the interaction strength and the flavor-singlet mass. We conclude in section 7, and technical details are worked out in several appendices. Throughout this article we will work in the natural units where = c = k B = 1 and with Einstein's summation convention where we sum over repeated indices.

Planar four-fermion model
We consider a system of two-component Dirac fermions ψ i sα in 2 + 1 dimensions. Here α = 1, 2 are spinor indices, i = 1, · · · , N are color indices and s = 1, · · · , M are flavor indices. The Lagrangian in the Euclidean spacetime is given by transformations of ψ. 1 The mass term κψψ breaks parity symmetry, and µ is the baryon chemical potential. The four-fermion interactions of the form (2.1) arise in the random matrix model proposed in [49] which also gives the sign of the interaction terms. We underline that these signs are essential for the results of the present work.

(2.4)
Next, we introduce a shifted field Φ ≡ Φ + ig 1 g 2 φ1 M + κ Assuming that the condition g 2 2 > M g 2 1 is fulfilled, the integral over the φ field can be carried out and leads to the result 1 A three-dimensional four-fermion model having this symmetry was investigated in [59,60]. We thank K. G. Klimenko for bringing these references to our attention.
with m ≡ g 2 κ 2g 2 1 . (2.7) After substituting this into (2.4) and the shifted field Φ ≡ Φ + ig 1 g 2 φ1 M + κ 2g 2 1 M we arrive at (2.6). Alternatively, the Gaussian integral over φ in equation (2.4) can also be evaluated from the saddle point equation in φ with the saddle point φ = (ig 1 /g 2 ) Tr Φ. In the large-N limit the partition function is dominated by saddle points of the effective potential where L is the linear extent of the plane. Assuming a constant field Φ (τ, This is the main result of this section. The momentum integral for the zero-temperature part is UV divergent and we regularize it by a cutoff Λ. What happens if we switch off the Gross-Neveu-type interaction by letting g 1 → 0? In this limit the potential becomes where a divergent constant independent of E has been dropped. For any κ = 0 the origin of the E k is unstable due to the presence of a linear term and E k develops a nonzero vacuum expectation value (VEV) E ∝ 1 M at all temperatures. There is no spontaneous breaking of U(M ) symmetry. As will be shown in the following sections, the situation is dramatically different for g 1 = 0; we will see a rich pattern of symmetry breaking taking place. According to the Coleman-Mermin-Wagner-Hohenberg theorem, in the absence of long range interactions, continuous symmetries cannot be broken spontaneously in two-dimensions which includes 2+1 dimensions at nonzero temperature. However, fluctuations that destroy the condensate, are suppressed for N → ∞ and spontaneously symmetry breaking is possible also at nonzero temperature. Below we always analyze the large N limit, but in appendix C we argue that even at finite N the same phase transitions still may be observed in particular if they are of first order.

Numerical results
We begin our discussion with the vacuum, T = µ = 0, to see what the underlying phases are. The momentum integral can be done analytically and yields the effective potential It is convenient to introduce dimensionless variables which will be used primarily for the numerical results in this paper. For the analytical results we stick to a different notation, see next subsection, to simplify mathematical manipulations.
The λ-dependence of the minimum of V eff (E) for g 1 = 1 and g 2 = 3 at T = µ = 0. Each e k is represented by a different color. The plots are vertically displaced slightly so that they do not overlap exactly.
The dimensionless potential reads When g 2 2 > π, so that the potential v(e) is confining, v(e) has two minima. For large negative λ, all the e k are degenerate and negative. With increasing λ, one of the eigenvalues jumps to a positive value. After increasing λ further, another eigenvalue jumps. This continues until all eigenvalues become degenerate and positive.
To understand this phenomenon quantitatively, we performed numerical minimization of the potential. In figure 1 we display the λ dependence of the e k at the minimum of the potential V eff (E). The e k jump M times, marking M sequential first-order phase transitions. Hence, there are in total M + 1 different vacuum states. The ground state energy density shown in figure 2 exhibits M sharp kinks associated with the jump of the eigenvalues. In   Figure 2. The ground state energy density for g 1 = 1 and g 2 = 3 as a function of λ.
The blobs denote first-order phase transitions. Each phase is labeled with its unbroken symmetry group. The phase structure shown here generalizes to higher M in an obvious manner.
The above numerical findings will be rigorously derived in the next subsection.  [63][64][65]. However, the sign of the Chern-Simons term can nullify this phase so that subleading saddle points become dominant resulting in symmetry breaking patterns U(M ) → U(j)×U(M −j) and a cascade of phase transitions [44]. Recently it was proposed that QCD 3 with a large number of colors would undergo a sequence of first-order transitions when the flavor-singlet mass is varied, in exactly the same fashion as figure 3 [51]. Although we may have local minima leading to the symmetry breaking pattern U(M ) → U(j) × U(M − j − 1) × U(1), saddle points with even less symmetry are unnatural and unlikely to be global minima of the free energy. It may require fine tuning of the parameters if they exist. This is indeed the case for the model analyzed in the present work and makes our model a fascinating theoretical laboratory of ideas and methods for QCD 3 .
Then, the potential takes the form with three parameters γ 1 , γ 2 and λ determining the phases. They essentially correspond to the relative strength of the two quartic interactions, the inverse strength of the interaction proportional to g 2 and the mass term, respectively, cf. In Appendix A, we study the solutions and some properties of the corresponding phase diagrams for general g(e) and illustrate it with a simple but non-trivial toy model. To solve ∂ e k V eff (e) = 0 for a fixed s, we note that g(e) is a strictly monotonously increasing function if γ 2 ≥ 3/2 (Λ < π/g 2 2 ), see the red dashed curve in figure 4, as can be seen from its derivative (1 + 2e 2 ) 2 + e 2 (4|e| −1 + 1 + 4|e|) . Hence, in the regime γ 2 ≥ 3/2 we have only a single real solution for e k = e (0) with unbroken flavor symmetry. Explicit expressions for e (0) are derived in Appendix B. Let us emphasize that this part of the phase diagram will be avoided when the cutoff Λ is chosen large enough.
When γ 2 < 3/2, there may be two minima e − (s) < 0 < e + (s) and a maximum e 0 (s) of the confining potential v(e) for each e k when s is fixed because the derivative g(e) = v (e) has the form of the blue curve sketched in For M > 2, γ 1 > 0 and γ 2 < 3/2, we can exploit insights from the previous work [49] and the discussion in Appendix A. As is the case in the random matrix theory [49] on which the present model is based, there are always the phases corresponding to the symmetry breaking pattern U(M ) → U(j) × U(M − j) with j = 0, . . . , M (see Appendix A). The integer j is a monotonously increasing function of λ when the solution is given by e = (e a 1 M −j , e b 1 j ) with e a < e b . The phase transition between the phase j and j + 1 with j = 1, . . . , M − 2 evidently has to be of first order because for a second order phase transition to happen we need e a = e b at the transition point.
As we have seen in Appendix A.3.1 also the phase transitions between j = 0 and j = 1 as well as j = M − 1 and j = M are of first order. This time the Hessian is positive semi-definite so that a second order phase transition might be possible, but there is always a direction, namely the eigenvector of the Hessian with the eigenvalue 0, whose leading term in the Taylor expansion becomes negative (see Appendix A.3.1).
Although general arguments indicate [66,67] that the flavor breaking ground state should have maximum symmetry, more exotic symmetry breaking patterns such as  The effective potential (2.13) at T > 0 and µ = 0 can be evaluated analytically. Absorbing where Li s (z) = ∞ k=1 z k k s is the polylogarithm function. 3 In the notation of (3.6), we find i) The vanishing of the slope at e = 0 (blue line in figure 6), Because the asymptotic behavior of g(e) ≈ 2γ 2 e, if the slope at 0 is negative g(e) cannot be a monotonic function, and the equation g(e) = −2s can have three solutions for T < (3 − 2γ 2 )/(6 log 2).
ii) The curve in the (γ 2 , T ) plane (red curve in figure 6) with g(e) = g (e) = 0 (4.5) separates region II and region IV. At those points a minimum of g(e) touches the e-axes, and because g(e) is an odd function, the equation g(e) = −2s can have 5 possible real solutions in the region IV.
iii) The curve in the (g 2 , T ) plane when a minimum and a maximum of g(e) coincide (green curve in figure 6) is given by It indicates definitely a phase transition that splits the region above curve i) and ii) into regions I and II. In region I the function g(e) increases monotonously, while in region II the equation g(e) = −2s can have at most three solutions despite g(e) has two local minima and two local maxima.
At the tricritical point, the potential which had three minima in the region V, joins the potential in regions I and II, with one and two minima, respectively. Because the potential is even, this has to happen at e = 0. The condition for the tricritical point is thus which is solved by A second special point in the (γ 2 , T ) plane is the point on the curve g (e) = g (e) = 0 where This point is at γ cr 2 = 0.3278 with T cr = (3 − 2γ 2 )/(6 log 2). For γ 2 < γ cr 2 the system always experiences a cascade of phase transitions when varying λ.

High temperature regime
At sufficiently high temperatures and fixed γ 2 > 0, see region II in figure 6, the curve g(e) shows a "wiggle" (a local maximum followed by a minimum) for large |e|. Taking into account the arguments of [49] and the discussion in Appendix A.3, we expect a cascade of phase transitions for sufficiently large | λ|. This is indeed observed numerically, see the plots in figures 9 as well as 10 for M = 2, 3, 4. It shows as a strip which obeys approximately a linear relation between T and λ. The cascade of symmetry breaking patterns are those of In Appendices A.3 and A.3.1 we have argued that all phase transitions for a locally double well shaped potential have to be of first order for M ≥ 3. For M = 2, a second order phase transition is possible, but our numerics confirm that at high temperature all transitions are first order. As we will see in the next section, a second order phase transition does occur for M = 2 at lower temperatures.

Low temperature regime
At low temperature T < (3 − 2γ 2 )/(6 log 2) and γ 2 < 3 2 (region III in figure 6) we find a g(e) in the shape of a wiggle, this time about the origin. When increasing the temperature we encounter three scenarios depending on the value of γ 2 which will be discussed in the next three subsections.

Low temperature regime with
√ π < | g 2 | < g tri 2 (γ tri 2 < γ 2 < 3 2 ) For γ tri 2 < γ 2 < 3 2 , the function g(e) becomes a strict monotonously increasing function for T > (3 − 2γ 2 )/(6 log 2), as we already have seen at T = 0 for γ 2 > 3 2 (or | g 2 | <  figure 7). This is the starting point of a strip of two close first order phase transitions in the ( λ, T ) plane and begins at λ = 6.2. The first order transitions are from a phase with unbroken flavor symmetry to a phase with U(1) × U(1) breaking and back to the unbroken phase. Adopting the order parameter M k=2 (−1) k (e k−1 − e k ) (assuming e 1 ≥ e 2 ≥ · · · e M ) the phase diagram for M = 2 in the ( g 2 , T ) plane is mapped out in figure 7. The second order line extends from λ = −λ tri to λ = λ tri . We have omitted the high temperature regime in this figure which contains the strip of first order phase transitions. It will be discussed in more detail in the next subsection.
As is discussed in Appendix A.3.1 for M ≥ 3 there is no line of second order phase transitions in the ( λ, T ) plane and the only second order point at λ = 0.
Moreover, we expect a cascade of phase transitions between phases with the symmetry breaking pattern U(M ) → U(j)×U(M −j) with j = 0, 1, . . . , M , as in the high temperature phase, the system runs through all possible j from j = 0 to j = M when increasing λ. We have corroborated this by numerical minimization of the potential (4.1) for M = 3 and M = 4 (see figure 8) where in both cases we have chosen g 2 = 3 < | g tri 2 |. Again we did not consider the high temperature phase.
4.3.2 Low temperature regime with g tri 2 < | g 2 | < g cr 2 (γ cr 2 < γ 2 < γ tri 2 ) In this regime we have a richer phase diagram which is mapped out in figure 9 using e 1 − e 2 as an order parameter. The most notable feature is that the strip with the cascade of phase transitions is split into two pieces. The strips end in second order points that are visible in the (γ 2 , T ) plane as the two transitions from region I to region II. The two parts join each other at g 2 = g cr 2 . The function g(e) is shown in figure 9 for three different temperatures, T = 0.59 corresponding to the lower part of the strip (green dotted curve), T = 0.7 in between the two strips (blue solid curve) and T = 2.5 corresponding to the upper part of the strip (red dashed curve).
The transition between region III and region IV is first order. Since the curve separating the regions III and IV describes two minima of v(e) coalescing with the minimum at e = 0, one would expect a second order transition, and it may be accidental that the position of the first order transition is located on this curve (it could also be that our numerically accuracy is not sufficient). In the region IV we have three first order transitions as a function of λ while there are only two transitions in the region II which become second order at an intermediate value of T . (c) Figure 9. (a) The phase diagram for M = 2 (a) and M = 3 (b) at µ = 0 with g 1 = 1 and g 2 = 3.75 in the large-N limit. The magnitude of |e 1 − e 2 | is plotted. The strip of first order transitions for high temperature is interrupted roughly between λ = 1 and λ = 2, but is present close to the broken phase around the origin and at high temperatures. In figures (c) we show a log-log-plot of the function g(e) where γ 2 = 3π/(2 g 2 ) ≈ 0.34 for three different temperatures T = 0.59 (green dashed curve), T = 0.7 (blue solid curve), and T = 2.5 (red dashed curve). The remnant of the strip close to the bulk of phase transitions at the origin can be explained by the existence of a "wiggle" of g(e) which briefly dissolves for larger temperature and reappears anew. For smaller g 2 (larger γ 2 ) the high temperature strip of phase transitions is completely separated from the broken phase near the origin, cf. figures 7 and 8.

Low temperature regime with
For these values of g 2 or γ 2 , the system no longer enters region I with increasing temperature so that the strip with the cascade of first order phase transitions is no longer interrupted. We have numerically analyzed this regime for M = 2, 3 and 4 at g 2 = 5 in figure 10. The cascade of phase transitions at high temperature has been visualized in figure 10(d) where we have plotted the actual solutions e k at the global minimum of the potential (4.1). To understand the nature of the two phases above and below this strip, we interpret the e k as the effective masses of the fermions of the theory. As shown in figure 10, the low-T region is characterized by a large value of |e k | implying that the effective masses are heavy. In contrast, in the high-T region above the strip all |e k | drop nearly to zero, making  U(2) × U(1) (namely, two of the three e k coincide). Yet, in a tiny region, shown in figure 11, all e k are mutually distinct and break the symmetry as U(3) → U(1) × U(1) × U(1). In this phase, one of the bosons is very light but the other two are heavy. Indeed when | g 2 | > | g cr 2 | (or |γ 2 | < |γ tri 2 |), we find a different kind of transition in the shape of the function g(e) in the region IV (see insets in figure 6). One of the consequences is the occurrence of exotic phases corresponding to the symmetry breaking patterns U has three positive slopes so that the potential v(e) has three minima, see Appendix A.2.

Nonzero chemical potential
The zero-temperature potential at µ > 0 can be readily found from (2.12) as is as given in (3.2) and Θ(x) is the Heaviside step function. In dimen- (µ − 2| g 2 e k |) 2 (µ + 4| g 2 e k |)Θ(µ − 2| g 2 e k |) . For analytical considerations, we adopt again the notation of (3.6), where the potential becomes and s = γ 1 ( M k=1 e k − λ). The insets in figure 12 show the different shapes of the function g(e) in the (γ 2 , µ) plane. Taking into account the Heaviside Θ function, in a similar way as at µ = 0 and nonzero temperature, the regions are separated by the following three curves: i) The vanishing of the slope at e = 0, g (e = 0) = 2γ 2 + 3µ − 3 = 0, ii) The curve (red curve in figure 12) defined by g(µ) = 0, (5.6) has the explicit solution is given by iii ) The third curve is given by The limit is introduced so that the Heaviside Θ function is equal to one (green curve in figure 12). This equation can be solved explicitly with the solution given by From that expression we obtain the critical point These three curves partition the (γ 2 , µ) plane into four regions which are anew referred to by Roman numerals, see figure 12. The curves meet at the tricritical point where the minima of the potential coalesce. Combining 2γ 2 + 3µ = 3 with 2γ 2 = 3( 1 + µ 2 − µ) we find that the tricritical point is at µ tri = 0 and γ tri 2 = 3/2 (blue point in figure 12). Since the shapes of g(e) at nonzero µ and T = 0 are similar to the shapes of g(e) for µ = 0 and nonzero T , we expect a similar phase diagram where we again can distinguish three regions depending on the value of γ 2 relative to γ cr 2 and γ tri 2 . Especially, we see a cascade of phase transitions between phases of the form U(M ) → U(j) × U(M − j). We have checked this numerically for M = 2 and M = 3 with g 2 = 3 (γ 2 = 0.5236), see figure 13 (a) and (b). As is the case at nonzero T and µ = 0 for γ 2 < γ cr 2 the strip of first order phase transitions is connected. Both inside this strip and in the region around (µ, λ) = (0, 0) we observe a cascade of first order phase transitions.
When increasing the chemical potential at fixed | g 2 | > g cr 2 (or γ 2 < γ cr 2 ) we enter region IV, where v(e) has three minima. This opens the possibility of phases with the symmetry breaking pattern U(M ) → U(j) × U(k) × U(M − j − k). This has been indeed observed by us for M = 3, see figure 14. The region where this kind of symmetry breaking pattern happens is very narrow as it has been the case for finite temperature. Note the similarity of figure 14  we found an exotic phase in which U(4) is broken down to U(2) × U(1) × U(1). For M = 5 we found a phase with the breaking pattern U(5) → U(2) × U(2) × U(1). When γ 2 < γ cr 2 the strip of first order phase transitions is connected to the broken region around the origin of the (λ, µ) plane. For γ cr 2 < γ 2 < γ tri 2 the strip is interrupted exactly as in the finite T and µ = 0 case. The part that is connected to the broken region about the origin disappears for γ 2 < γ tri 2 . The strip divides the U(M )-symmetric part of the phase diagram into two regions. The qualitative distinction between these two regions is clear from the behavior of the {e k } (figure 13). They start with a large value at low µ and then successively drop to a small value as µ increases.
At nonzero chemical potential the region between the tricritical point and γ 2 = 3/2, which gives a second order phase transition for M = 2 at low temperature, is absent. Therefore, we do not expect second order transitions at µ = 0 and T = 0, even for M = 2. Only at the end points of the cascades of phase transitions where all first order transitions run together, we expect second order phase transitions.
Finally, we consider the fermion number density n. Since it is proportional to N it is  In figure 15 we show the µ-dependence of the fermion number density. It jumps at first-order transitions. The plots for λ = 0.6 (M = 2) and λ = 0.8 (M = 3) reveal that the two U(M )-symmetric phases at small µ and large µ are physically quite distinct, as was suggested in figure 13 as well. At small µ, fermions are rather heavy and the number density is quenched to zero. In contrast, at large µ, fermions are nearly massless and their density is enhanced: the large bare mass originating from κψψ in the Lagrangian is dynamically screened by interactions.

Phases at nonzero T and µ
Finally, in this section we investigate the effect of simultaneously nonzero T and µ. From (2.12) the potential in this case is found to be To speed up numerical minimization, we used the Taylor expansion of Li 2 (z) and Li 3 (z) around z = −1 up to 11th order to compute values for −1 ≤ z < 0, and then used functional identities relating Li s (z) to Li s (1/z) [68, 69] to compute values for z < −1.
De novo we use the notation of (3.6) in which the potential reads  The coefficient of the linear term determines the plane (blue dotted plane in figure 16) that separates region III from regions I and IV, which is explicitly given by The intersections of the blue surface with the µ = 0 and T = 0 planes are given by the blue curve in figures 6 and 12, respectively. For M = 2 and λ = 0, the part of this surface between regions I and III allows second order phase transitions. It continues to exist for not too large values of the chemical potential. When γ 2 < γ 2 (T, µ) the phase will have the symmetry breaking pattern U(2) → U(1) × U(1), and for γ 2 > 3/2 the flavor symmetry remains unbroken at low temperature and chemical potential. There are two additional planes that divide the (µ, T, γ 2 ) space. As is the case at zero chemical potential, the first one is given by (green dotted surface in figure 16) g (e) = g (e) = 0. (6.6) On this surface, that separates regions I and II, the extrema of g(e) in region II join so that g(e) in region I becomes monotonous. The second plane is given by the equation (red dotted plane in figure 16) g(e) = g (e) = 0. (6.7) On this plane the minimum of g(e > 0) touches the e-axis so that the corresponding potential will have three minima in region IV. The tricritical points at T tri = 1/2 for µ = 0 becomes now a tricritical curve (see black curve in figure 16 for finite µ, namely when the cubic term in g(e) is also vanishing, which is at There are bounds for the location of this curve, particularly T tri ∈ [0, 1/2], γ tri 2 ∈ [3/2(1 − log 2), 1.5] and µ tri ∈ [0, 0.45] (the latter number is an approximation for the maximum of the right hand side of (6.8)). Whenever µ < µ tri for M = 2 the system experiences a second order phase transition at γ 2 = γ 2 (see eq. (6.5)). Larger values of M remain untouched and all phase transitions are of first order apart from the critical points where first order transition lines end. Also for suitably large λ all second order phase transitions will vanish and what remains are first order transitions.
In region IV at fixed γ 2 the potential will have three minima and we can again expect exotic phases of the form U(M ) → U(j) × U(k) × U(M − j − k). They will certainly appear only for small µ since this has been already the case for either T = 0 or µ = 0.
For suitably large γ 2 , µ and T , we find either a strictly monotonous g(e) (region I) or one which has local maxima and minima symmetrically about the origin in two separate quadrants (region II). The latter signals again the existence of a strip of cascades of phase transitions in the high T and µ region. Yet, the shape of g(e) in region II can be also found for a small region for γ 2 < γ tri 2 which will show itself as a remaining appendix of this strip of phase transitions at the phase region about the origin, cf. figure 9.
For simplicity of exposition we limit our numerical analysis to M = 3. Our main results are summarized in figure 17 (for g 2 = 3) and figure 18 (for g 2 = 5). Figure 17 shows that the phase structure depends on λ in a nontrivial way. At small λ, there is a large region at low T and low µ in which U (3)  boundary of this phase, there is a narrow strip in which U(3) is broken to U(1)×U(1)×U(1) (cf. figure 14). As λ increases, this strip gradually disappears. At λ = 0.4 there are three symmetry-broken phases: they have the same symmetry (U(2) × U(1)) but are separated by first-order phase transitions. As λ increases further, the symmetry gets restored in the low-T low-µ region but remains broken in the cold dense region. At stronger coupling a qualitatively new feature emerges. In figure 18 we observe that the symmetry-broken phase forms a thin annulus, separating the low-T low-µ region from the high-T high-µ region. This annulus never disappears even at very large λ, although it is shifted to higher T and µ gradually. By monitoring the behavior of |e k | we found that the U(3)-symmetric phase below the annulus is characterized by very heavy fermions, while the other U(3)-symmetric phase above the annulus is characterized by massless fermions. Although these two phases cannot be distinguished by symmetries, they host quite different physics.

Conclusions and outlook
In the present article, we investigated various aspects of Dirac fermions with nonstandard quartic interactions in two spatial dimensions. We showed within the mean-field approximation that the model experiences a cascade of phase transitions when the flavor-symmetric parity-breaking mass is varied, in a way quite analogous to the behavior of QCD 3 [44,51]. At nonzero temperature and chemical potential we provided analytical and numerical arguments that show how a complicated phase diagram embellished by exotic symmetry breaking patterns emerges. In particular, we showed (in figures 10, 13 and 18) that, at strong coupling, the low-(µ, T ) phase with heavy bosons is separated from the high-(µ, T ) phase with almost massless bosons by a series of M phase transitions, through which M species of fermions become light one after another. At finite temperature there is a subtlety about symmetry breaking due to enhanced infrared singularities and we gave a speculative comment on this. Summarizing above, our results shed light on previously unnoticed novel dynamics of Dirac fermions in 2 + 1 dimensions and have potential implications for planar gauge theories as well as planar condensed matter systems.
There are several directions in which this work can be extended. First, the present analysis in the large-N limit could be generalized to incorporate finite-N corrections. Fluctuations of bosonic fields can be conveniently included by employing methods such as the functional renormalization group [70]. At finite N , the Jacobian (the squared Vandermonde determinant) associated with the diagonalization of the matrix field can no longer be neglected and will affect ground state properties. Secondly, it would be interesting to see what happens if our assumption g 2 2 > M g 2 1 is relaxed. Thirdly, various topological excitations arise in our model. For instance, in the phases depicted in figures 11 and 14, π 2 (U(3)/U(1) 3 ) = Z × Z, implying there are two kinds of Skyrmions. Fourthly, while we have only considered fermion-anti-fermion condensates, a di-fermion condensate may form at high density [71,72]. The competition of two kinds of condensates may be an interesting subject of research. Finally, it would be challenging but quite important to take into account the possibility of an inhomogeneous condensate that spontaneously breaks translation symmetry. While the existence of such a condensate has been firmly established in some (1 + 1)-dimensional models at finite density [73][74][75], the situation is elusive in higher dimensions [76,77].

A Phase diagram of the effective potential and a toy model
In this Appendix, we outline the general strategy for analyzing the phases for an effective potential of the general structure considered in the main body of the text given by the sum of a confining (γ 1 > 0) harmonic collective potential and confining "single-particle" terms with potential v(x), (A.1) Generically, we assume that v(x) has the shape of a double well potential that increases faster than linear for large |e| (i.e. lim |e|→∞ v(e)/|e| = ∞)). Such potentials show a similar behavior such as the cascade of phase transitions and the kind of symmetry breaking patterns we have found in the physical system. Moreover, the mechanism when and how the system experiences a second order phase transition is very similar for different v(e).
In the last part of this Appendix, we illustrate the general arguments with the properties of a much simpler toy model (indicated by the subscript tm) with the confining potential v tm (e j ; T ) = (e 2 This model is motivated by the analysis in section 4 and its numerical observations. It is essentially a truncation of the expansion of the single particle part of the potential (4.2) to fourth order, which is expected to capture some essential part of the physics in the vicinity of a second-order phase transition. Regardless of its simplified form, this toy model already exhibits generic features for general v(e). Hence, we would like to underline that most conclusions apply for a more general confining potential v(e).

A.1 Saddle point equation and its asymptotic solutions
What has to be studied are the M saddle point equations The extrema of the potential V (e, λ) are determined by the intersections of g(e) with −2s, which also select the possible phases, especially, which symmetry breaking patterns, the system can exhibit as a function of λ and the parameters of the potential. Note that s can have different values for the same values of the parameters. The reason is that the solution of the saddle point equations is not unique. One particular ingredient is the asymptotic behavior of the solutions of (A.3) for large |λ|. The asymptotic super-linear growth of v(e) implies also an asymptotic growth of |g(e)|. Particularly we have three cases to consider where the asymptotic value can be either vanishing (c ± = 0), be finite (0 < c ± < ∞), or diverge (c ± (e) = ∞).
Depending on which case the asymptotic solution for e j becomes unique and takes the form c sign(λ) = 0, , 0 < c sign(λ) < ∞, for all k = 1, . . . , M . The function g −1 is the inverse of g in the asymptotic regime; for instance when g(e) grows like e L , the inverse is essentially e 1/L . In the physical system in the main text we have the situation of a finite c + = c − while the toy model (A.2) leads to c sign(λ) = ∞. An asymptotic behavior with c + and c − in a different class is possible, but we do not consider that in the present work. Regardless which case the asymptotic satisfies, we obtain the same conclusions for the global minimum of the potential. First, all e k are degenerate for suitably large |λ|. Second, the modulus of the auxiliary parameter s in (A.3) also grows asymptotically, and its sign is the opposite of λ and e k . As a physical conclusion we find that for suitably large |λ| we always have a solution with all e k equal which has unbroken flavor symmetry.

A.2 Local extrema of g(e) and implications on the possible phases
The solutions of the saddle point equation can be either minima or maxima for V (e, λ). For a saddle point with g (e k ) > 0 for all e k the potential has certainly a minimum, not necessarily the global one we are looking for. This follows from the fact that the Hessian at the saddle point is given by The Hessian is positive definite if the determinants det(H jk ) j,k=1,...,n = det{2γ 1 + g (e k )δ jk } j,k=1,...,n > 0, for all n ≤ M. (A.7) The term proportional to 2γ 1 is of rank 1 which simplifies the evaluation of the determinant drastically and it is equal to det(H jk ) j,k=1,...,n = 1 + 2γ 1 In general, we may have a solution with L different e k . At most one of the e k may have g (e k ) < 0. The reason is that the term of the Hessian that is proportional to γ 1 is of rank one. A rank one addition can maximally switch one eigenvalue of a matrix from positive to negative and vice versa, regardless how large its prefactor is. Let us label the intersection with g (e k ) < 0 as e M . Then all subdeterminants up to n = M − 1 are positive, and the condition for the positive definiteness of the Hessian matrix is given by the positivity of its determinant, If the solution with L different e k is a global minimum, this would result in the symmetry breaking pattern U(M ) → U(j 1 ) × · · · × U(j l ) with L l=1 j l = M . The unbroken symmetry associated with e k with g (e k ) < 0 can only be a U(1) factor. However, not all of these saddle points are global minima of V (e, λ) which is the hard part of the analysis.
The simplest case is M = 2. Then the only possible flavor symmetry breaking pattern is U(2) → U(1) × U(1), that means a transition from a solution with e 1 = e 2 to a solution with e 1 = e 2 . This transition can only be of second order if the solutions join continuously to the point where the determinant of the Hessian vanishes, i.e. for parameter values with g (e 1 ) = g (e 2 ) = 0 (or at g (e 1 ) = g (e 2 ) = −4γ 1 , but that is not allowed for γ 1 > 0). There is no second order phase transition when g(e) increases monotonically. More generally, the phase transition is first order when the global minimum is not a continuous function of λ.
As is shown for the potential in the main text and for the toy model (A.2) a second order phase transition is realized for M = 2. However, as we will show below, for M > 2 the phase transition for a potential of the form (A.1) is always first order.

A.3 Cascade of phase transitions
In this subsection, we discuss the phases of the general potential (A.1) with v(e) having locally the shape of a confining (not necessarily symmetric) double well. The notion "locally" means that there is region for s bounded by a local maximum and local minimum of g(e) where the saddle point equation (A.3) has only three solution when fixing s. The integral of g(e) in this region looks like a double well potential.
For the solutions of the saddle point equations we use the Ansatz diag(e) = diag(X 1 1 M −j , X 2 1 j ), with X 1 < X 2 without restriction of generality. The two variables X 1 and X 2 must satisfy the equations −2γ 1 ((M −j)X 1 +jX 2 −λ) = g(X 1 ) and −2γ 1 ((M −j)X 1 +jX 2 −λ) = g(X 2 ). (A.10) These equations can be solved for a real j ∈ R. We would like to highlight the difference of j ∈ R and j = 0, . . . , M ; while the former is real and can take an optimal position minimizing the latter can be only an integer and only approximately minimizes the potential. The saddle point equation for j is given by which yields a unique minimum for j in terms of X 1 and X 2 . This follows from the second derivative in j which is always positive when X 1 = X 2 , Note that for X 1 = X 2 the saddle point equation for j is satisfied trivially. The minimizer j will generally not lie on one of the integers j = 0, . . . , M . Yet, the convexity of V j (X 1 , X 2 , λ) in j shows that only those closest to j minimize the potential V j (X 1 , X 2 , λ) with j = 0, . . . , M .
The saddle point equations (A.10) are invariant under as is the saddle point equation (A.12) for j. We therefore must have Since λ is a function of j, the saddle point solutions are functions of j only and increasing λ will increase j. For the discretized version j = 0, . . . , M the solutions X 1 and X 2 get an explicit dependence on λ though it has only a limited impact as j tries to be as close as possible to j. The convexity of the potential as a function of j also tells us that there can be only phase transitions from the phase with flavor symmetry U(j) × U(M − j) to the phase with flavor symmetry U(j + 1) × U(M − j − 1) or U(j − 1) × U(M − j + 1) for j = 1, . . . , M − 1. We still could have a phase transition between a solutions with all e k the same. As we will see below, in the toy model, this happens when T > 1. Also in the physical system at finite temperature and/or at finite chemical potential there is a region where the system may experience such a direct transition from all e k equal and negative to all e k equal and positive, see sections 4, 5 and 6. In summary, the system runs through all phases corresponding to U(M ) → U(j) × U(M − j) from j = 0, . . . , M as the real minimizing set ( j, X 1 ( j), X 2 ( j)) will depend continuously on λ. The kinks and, hence, phase transitions only originate from the discreteness of j rather than the continuous variable j → j implying X 1,2 ( j) → X 1,2 (j, λ).
Let us underscore the following. What is locally required for the discussion above is a potential with two minima and the validity of our assumption that the bipartite Ansatz is valid. However, we could not exclude the possibility that the flavor symmetry is broken according to U(M ) → U(j) × U(M − j − 1) × U(1) when the potential has a minimum with three different e k , two with g (e k ) > 0 and one with g (e k ) < 0. This possibility has to be investigated case by case.

A.3.1 No-go statement for second order phase transitions
The question that remains is whether any of these phase transitions are of second order. The transition from j to j + 1 for j = 1 and j = M necessarily has to be of first order, because if X 1 = X 2 is at a second order phase transition point, we would have a transition from a state with j of the e k at X 2 to a state with all of the e k at X 2 as the positivity condition (A.7) does not allow any other possibility. The only exceptions are the transitions from j = 0 to j = 1 and from j = M − 1 to j = M . We now consider the latter, while the former can be worked out in the same way.
As before, the arguments below apply to a potential v(e) that has locally the form of a double well potential and which is super-linearly increasing for large and small e, this means its derivative g(e) = v (e) has the shape of a wiggle. Because the potential is homogeneous in the x k the derivatives of the potential with respect to x 1 and x 2 also vanish for this Ansatz. Since the determinant of the Hessian vanishes at e 0 , there is at least one direction in which the second order fluctuations vanish. To find a decreasing direction, we thus have to Taylor expand the potential at least to third order In the direction of vanishing second derivative, given by the third order term behaves as For the assumed shape of the potential we have g (e 0 ) > 0 (minimum of g(e) at e 0 ) so that generally the third order term becomes negative. One exception is M = 2. In that case the fourth order term of the expansion is positive. So for M = 2 the point e 1 = e 2 = e 0 can be a global minimum. For M > 2, we always have a decreasing direction excluding the possibility of a second order phase transition. When g (e 0 ) = 0, we need to expand to higher orders. For a local double well shape of v(e) or local wiggle shape of its derivative g(e) the first non-vanishing derivative of v(e) must be even. For L = 3, 5, . . ., we obtain because it must be g (L−1) (e 0 ) > 0 if we consider the transition from j = M − 1 to j = M as g(e) has to have a minimum at e 0 . In summary, we can say that for M ≥ 3 the system always experiences a cascade of first order phase transitions. The only requirement is that the potential has locally the shape of a double well or its derivative g(e) has the shape of a "wiggle" meaning a maximum followed by a minimum and then growing again. This is the situation also for the physical system in the main text. There are surely more complex situations when the potential v(e) has more than two minima in a restricted region so that the saddle point equation (A.3) has more than three real solutions. For instance this happens in the middle temperature regime discussed in section 4. However, generally the mechanism for the phase transition is similar. This function can have only two distinct shapes depending on the temperature T . For T ≥ 1, it is monotonously increasing so that all e k need to be equal, and the flavor symmetry is not broken in this case. For T < 1, the function develops a local minimum and maximum. Therefore, equation (A.3) can exhibit three solutions for a fixed s, two at e (+) > 0 > e (−) with g tm (e (±) ) > 0, and one at e (0) with g tm (e (0) ) < 0. In agreement with the asymptotic analysis in subsection A.1, for sufficiently large or small λ only one solution exists when all e k are the same. In the toy model (A.2), we can have a transition from a broken phase with e 1 = e 2 to a phase with unbroken flavor symmetry with e 1 = e 2 . At a second order transition curve we therefore must have e 1 = e 2 = e 0 while the determinant of the Hessian must vanish. At this point we also must have that g tm (e 0 ) = 0. To determine the possible critical behavior we calculate the Hessian at the saddle point in terms of X = (e 1 + e 2 )/2 and ∆ = (e 1 − e 2 )/2 and simplify it with the difference of the two saddle point equations, The tricritical point, when ∆ changes from ∆ = 0 to ∆ = 0, is located at and λ given by (A.29). If e 1 = e 2 = e 0 is a global minimum, the line of second order phase transitions has to end at T = 1 − γ 1 /6 for γ 1 < 6. In particular, at zero temperature the phase transition has to be of first order for γ 1 < 6. The phase diagram will look like the sketch in figure 20. There is a region for 1 − γ 1 /6 < T < 1 where a second order phase transition happens. For T < 1 − γ 1 /6 we only have a first order phase transitions. To determine if e 1 = e 2 = e 0 is indeed a global minimum, we need explicit expressions for the solutions and substitute them in the potential which will be done in the remainder of this subsection.
We will solve the saddle point equations in terms of X = (e 1 +e 2 )/2 and ∆ = (e 1 −e 2 )/2. From the difference of the two saddle point equation (A.26) we see that we can have two different cases (A.32) The solution ∆ = 0 always exists and is even a solution for general g(e). It corresponds to unbroken flavor symmetry. The non-trivial solution for ∆ obviously represents the case U(2) → U(1) × U(1). The potential in terms of ∆ and X reads as follows 33) which implies that ∆ 2 = 1 − T − 3X 2 is the global minimum whenever this solution is allowed. The resulting saddle point equations for X of the two cases of ∆ differ and are given by for ∆ = 0, and for ∆ = ± 1 − T − 3X 2 1 . The latter one can be simplified by adding the two equations with signs ± which gives (A.36) Using eq. (A.29) for λ at the saddlepoint, this equation can be factorized as The solutions of the quadratic equation are given by One can easily see that the solution X + has the lower potential. Next compare the difference of the potential for X 1 = (1 − T )/3 and X + . With some work one can show that This is negative for T > 1 − γ 1 /6 which shows that the solutions with the second order phase transition is the global minimum. As we have seen before from the analysis of the Hessian, T = 1 − γ 1 /6 is the tricritical temperature. Summarizing, there is a second order phase transition whenever The situation for larger values of M is more complicated. The saddle point equation have 3 M solutions, most of them complex. However, we find a substantial number of real solutions with different values of k e k . A necessary condition for a global minimum is that the Hessian is positive definite, but to uniquely identify the solution we have to substitute it in the potential. For example, for M = 3, excluding permutations, we find three real solutions with a positive definite Hessian for a significant range of parameters. We also have a saddle point with all three e k different but this had never been a minimum for the values of the parameters we have analyzed.
However, as we argued in previous sections, the global minimum of the potential can have at most two different e k . Using this as an Ansatz, this substantially simplifies the saddle point equations. We still have that a second order phase transition can only happen at the minimum of g tm (e) which is at For γ 2 < 3/2 the function g(e) develops a local minimum and maximum. Thus, the derivative of g(e) with respect to e has to vanish at these points. Its symmetry tells us that there is one zero of g (e) for e > 0 and one for e < 0. Indeed, the equation g (e) = 0 can be rewritten as which is the local minimum of g(e). The local maximum lies at −e min . Hence, for a fixed s ∈ [g(e min )/2, g(−e min )/2] we find three solutions for g(e) = −2s. Two solutions of g(e) = −2s, that we denote by e (−) (s) < 0 < e (+) (s), come with a positive slope g (e (±) (s)) > 0. The two solutions e (+) (s) and e (−) (s) satisfy the relation e (−) (s) = −e (+) (−s) due to the symmetry of g(e). Thus, it is enough to state the solution for e (+) (s) = (x + − x −1 + )/2 with with a + and b + as in (B.12), since its limit for s = 0 should be x + = 1. For s ∈ [g(e min )/2, 0], we can use again the symmetry of g(e) = −g(−e), meaning the solution is then e (0) (s) = −e (0) (−s).

C Comment on the role of bosonic fluctuations
In the main text we have explored the pattern of symmetry breaking in the large-N limit. In this limit the fluctuations of the bosonic fields are completely negligible, but this is no longer the case at finite N . Actually the celebrated Coleman-Mermin-Wagner-Hohenberg (CMWH) theorem [78][79][80] stipulates that continuous symmetries cannot be spontaneously broken at nonzero temperature in 2+1 dimensions in the absence of long-range interactions. Thus any symmetry-breaking condensate at T = 0 must disappear as soon as nonzero temperature is turned on. No massless Nambu-Goldstone modes can appear; in fact they acquire nonzero masses, as demonstrated explicitly for O(N )-invariant models in [81][82][83][84][85].
In our model, the ground state has to be E ∝ 1 M everywhere on the phase diagram at T > 0. The second-order phase transition line in figure 7 will be wiped out at finite N and becomes a crossover. As long as N 1, the first-order transitions in figures 7, 8, and 10 may well persist. However, if thermal fluctuations were so strong that the critical point at λ = 0 (present in figures 8) is destroyed, then M first-order transition lines emanating from the T = 0 axis would probably end at M distinct critical points.
From the viewpoint of the CMWH theorem it may seem that there is no point in talking about symmetry breaking for T > 0. However this is not the case. Preceding analyses [81][82][83][84][85] have shown that the mass of the would-be Nambu-Goldstone modes m NG is of order F 2 exp(−cF 2 /T ), where c is an O(1) pure number and F 2 is a square of the "pion decay constant", also known as spin stiffness in the literature of quantum magnets. It is well known that F 2 ∝ N in the large-N limit [86], so we have parametrically m NG ∼ N exp(−N/T ) which gets exponentially small at low temperatures or large N . In experiments or numerical simulations, the correlation length ∼ m −1 NG can easily exceed the system size. The system is then virtually indistinguishable from a genuine symmetry-broken phase. The interactions of the would-be Nambu-Goldstone modes get weaker and weaker as the energy scale goes down, but start to increase at the scale ∼ m NG . As far as physics at length scales m −1 NG is concerned, it is perfectly sensible to adopt a description based on spontaneously broken symmetry. A detailed discussion on the consistency between the CMWH theorem and the utility of low-energy effective theories of Nambu-Goldstone modes can be found in [87].
To go beyond the mean-field analysis of this paper we must employ nonperturbative methods such as Monte Carlo simulations on a lattice, which seems feasible since the statistical weight (2.6) is real and nonnegative for even N .