Random matrix approach to three-dimensional QCD with a Chern-Simons term

We propose a random matrix theory for QCD in three dimensions with a Chern-Simons term at level $k$ which spontaneously breaks the flavor symmetry according to U($2N_{\rm f}$) $\to $ U($N_{\rm f}+k$)$\times$U($N_{\rm f}-k$). This random matrix model is obtained by adding a complex part to the action for the $k=0$ random matrix model. We derive the pattern of spontaneous symmetry breaking from the analytical solution of the model. Additionally, we obtain explicit analytical results for the spectral density and the spectral correlation functions for the Dirac operator at finite matrix dimension, that become complex. In the microscopic domain where the matrix size tends to infinity, they are expected to be universal, and give an exact analytical prediction to the spectral properties of the Dirac operator in the presence of a Chern-Simons term. Here, we calculate the microscopic spectral density. It shows exponentially large (complex) oscillations which cancel the phase of the $k=0$ theory.


I. INTRODUCTION
It is believed that in the QCD vacuum the strong interactions of gluons and quarks induce spontaneous breaking of chiral symmetry SU(N f ) R × SU(N f ) L → SU(N f ) V when the number of massless Dirac fermions N f is below the conformal window. The quantum fluctuations of the gauge fields in the broken phase manifest themselves in characteristic spectral fluctuations of the Dirac operator in the microscopic domain [1], which can be exactly reproduced by a zero-dimensional matrix model with the same global symmetries as QCD known as chiral random matrix theory [2][3][4]. We refer to [5][6][7][8][9] for reviews.
In three-dimensional spacetime, whether the symmetry of fermions is dynamically broken or not has remained a matter of debate for decades [10][11][12]. Three-dimensional gauge theories are of distinguished importance in various contexts, ranging from domain walls and surface (boundary) states in four dimensions to quantum Hall effects, graphene, spin liquids and high-temperature superconductivity [13][14][15][16][17]. Part of this rich physics stems from the existence of a Chern-Simons term. In three-dimensional QED (QED 3 ) the interplay of a Chern-Simons term and fermionic symmetry breaking was investigated in [18][19][20][21][22][23][24]. A consensus from these studies is that dynamical symmetry breaking is generally suppressed by a Chern-Simons term, because photons acquire a gauge-invariant mass term which in turn quenches quantum fluctuations. However, recent lattice simulations [25][26][27] report that dynamical mass generation of fermions does not occur in QED 3 even in the absence of a Chern-Simons term (see also [28] for a RG study with the same conclusion).
In contrast, in three-dimensional QCD (QCD 3 ) with an even number 2N f of massless two-component flavors and without a Chern-Simons term, dynamical symmetry breaking 1 is believed to take place through fermion bilinear condensation when N f is below a certain threshold [30][31][32]. The fermion condensate has been observed in quenched lattice simulations [33,34]. A non-chiral matrix model corresponding to (1) is also known [35]; see [36][37][38][39][40][41][42][43][44][45] for further developments. However, not much is known about QCD 3 at nonzero Chern-Simons level k. 2 This is partly due to the sign problem that makes a direct lattice simulation prohibitively hard. Recently, it was argued [47] that there is a finite window of N f in which a novel symmetry breaking occurs. New boson-fermion dualities describing the transition region of (2) were also proposed [47]. While a proof is not available yet, this conjecture passes nontrivial tests such as the matching of symmetries and anomalies, and consistency under mass deformations. Related discussions can be found in [48,49].
In this paper, we propose a new random matrix model that realizes the symmetry breaking scenario (2). 3 This is made possible through a judicious choice of a non-Gaussian weight for matrix elements in which k enters as a parameter. We show that in the large-N limit with N the matrix size, the model reduces to a sigma model with a target space of the complex Grassmannian U(2N f )/[U(N f + k) × U(N f − k)]. When k is varied there occurs a sequence of first order phase transitions that separate phases with different complex Grassmannians. By solving the model we delineate the structure of the Dirac operator spectrum that underlies the exotic symmetry breaking (2). Under the assumption that (2) indeed characterizes the vacuum of QCD 3 with a Chern-Simons term, our approach offers an entirely new way to probe the interplay of strongly coupled fermion dynamics and a topological term within a tractable framework. This paper is organized as follows. In Section II we define the model and derive key properties. The partition function of the model is computed and the phase structure as a function of k and the fermion mass is investigated. Section III is devoted to the spectral functions that are first derived at finite matrix dimension N . As k grows the spectral density evolves from a smooth semicircle to a distorted complex oscillatory form. In addition, we compute the large-N microscopic limit of the spectral density in the quenched and unquenched ensemble. The details of the calculations to get these results are given in Appendix A. Concluding remarks are made in section IV.

II. RANDOM MATRIX MODEL
We introduce three new random matrix models labeled by the Dyson index β in Subsection II A. They are associated with QCD 3 in the presence of a Chern-Simons term with fermions transforming in a complex/pseudoreal/real (β = 2/1/4) representation of the gauge group, respectively. The class β = 2 comprises quarks in the fundamental representation of SU(N c ) with N c ≥ 3; β = 1 includes quarks in the fundamental representation of USp(N c ) (here, N c must be even) 4 ; and β = 4 corresponds to quarks in the adjoint representation of SU(N c ) and in the vector representation of SO(N c ). In Subsection II B we give a discussion on the renormalization of the Chern-Simons term due to the dynamical quarks that are related to the η-invariant.
Each of the above three random matrix models produces a universal non-linear sigma model that is derived in detail for β = 2 in Subsection II C. As in the case of three dimensional QCD, the model experiences a phase transition from one to another sigma model due to the Chern-Simons-like term. In Subsection II D we show how the mechanism works in general, and in Subsection II E we illustrate our findings by studying the two-flavor case for β = 2. Therein, we also present finite results for the partition function at finite matrix dimension.

A. Partition function of our model
The partition functions of our model are defined by where A is a complex hermitian N ×N matrix, A is a real symmetric N ×N matrix, and A is a self-dual N ×N matrix whose elements are real quaternions 5 . The measures dA, dA and dA are the corresponding Lebesgue measures, in particular the products of the differentials of the real independent matrix entries. The square root of the determinants for β = 4, see (5), is exact and may be implemented as a Pfaffian determinant, The masses are gathered in the diagonal matrix M = diag(m 1 , . . . , m 2N f ).
The positive constant α β determines the effective strength of the "Chern-Simons" coupling. Starting with section II E it will be chosen In sections II C and II D the detailed form of α β is irrelevant apart from the convergence requirements of the integrals and, thence, remain unspecified therein. Indeed the integrability of the variable Tr A is guaranteed when α 1 < 1, α 2 < 1 and α 4 < 2, which is satisfied by Eq. (7). In Eq. (5) A is regarded as a complex 2N × 2N matrix, using (1 2 , iσ a ) as the quaternion basis. Note that all matrices are square, reflecting the absence of topological zero modes in 2 + 1 dimensions. The real parameter k corresponds to the Chern-Simons level, and 2N f represents the number of two-component Dirac fermions for β = 1, 2 and of two-component Majorana fermions for β = 4. The case of an odd number of flavors will not be considered in this paper. We expect that the mechanism described below should work similar to the even number of flavors case though their is a significant difference; the Goldstone manifold is disconnected for an odd number of flavors [35,40,41].
The models (3), (4) and (5) differ from the conventional random matrix models for QCD 3 [35,40,41] by the presence of the squared trace term in the exponent. 6 At k = 0, the latter makes the statistical weight complex-valued, just as the Chern-Simons term does in Euclidean QCD 3 causing the infamous "sign problem". It is not problematic for us because we can still solve the matrix models exactly without recourse to numerical simulations. Our motivation to include a squared trace term is that this deformation changes the pattern of flavor symmetry breaking. The microscopic large-N limit [2,3] makes this more lucid. For this purpose, we take N → ∞ and m f → 0 with m f = N m f and k fixed. 7 If N is identified with the volume of space-time, this limit is equivalent to the leading order of the ε-expansion in chiral perturbation theory [1,62], in which the partition function reduces to a non-linear sigma model of static Nambu-Goldstone modes. If k ∈ Z with |k| ≤ N f , one can show for the partition functions (3), (4) and (5) in the microscopic limit reduce to where M ≡ diag( m 1 , . . . , m 2N f ) and The Haar measure of the respective groups are denoted by dµ. This result is derived in the next subsection. Effectively we do not integrate over the whole group but a coset. These cosets are the Goldstone manifolds and reflect the patterns of flavor symmetry breaking in the "chiral" limit given by and N 2 f − k 2 Nambu-Goldstone modes, respectively. When k = 0 they recover the usual symmetry breaking patterns proposed for parity-invariant QED 3 and QCD 3 with no Chern-Simons term [10,12,30,35,40,41,43,45]. This agreement is nontrivial because the partition functions (3), (4) and (5) are different from those in [35,40,41] even at k = 0 due to the squared trace term. It highlights the universality of the large-N limit. When k = 0, the symmetry breaking schemes (12) coincide with the generalizations proposed recently [47] for QCD 3 with a Chern-Simons term at level k.
In (8), (9) and (10) we omitted overall multiplicative factors, which are all proportional to (−1) N k . Therefore choosing even N is mandatory to ensure positivity of the partition function, although the overall normalization of Z does not affect physical expectation values. The number of these configurations decreases according to a2/N . The limit a2/N → 0 (the physical case) becomes a Dirac delta distribution at r = π/2 (black vertical line). Since the histograms with fixed quotient a2/N agree almost perfectly, they are barely distinguishable, we are save to assume that the plots show the limiting large N behavior.

B. Chern-Simons term and η-invariant
The parameter k is not the only source of the Chern-Simons-level as we will show for the Dyson index β = 2. For M → 0 the phase of the fermion determinant also contributes by the η invariant η(A) = N j=1 sign(λ j ), see [63][64][65][66], as follows This can be combined with the imaginary part of the first term in the exponent of Eq. (3), which can be written as This approximation can be seen by considering the integral (3) for k = 0 which gives a quenched approximation for Tr A. The integral can be rewritten as For the physically interesting limit of , which is much larger than the average level spacing of O(1/N ). These fluctuations are collective, meaning that all eigenvalues of A move up and down with x in the same way. Therefore, starting with a configuration with an equal number of eigenvalues on the left and right of zero, a fluctuation where k eigenvalues move to the right changes Tr A by δ Tr A = N k∆λ (16) and ∆λ = π/N the level spacing so that The sum of the sign of the eigenvalues changes by This results in the ratio which is the desired relation (14). This relation (14) has also been checked numerically, see Figure 1, where the distribution of the ratio r = N j=1 λ j / N j=1 sign(λ j ) has been generated with Monte Carlo simulations. Summarizing, the phase of the fermion determinant renormalizes the bare Chern-Simons level as Therefore, the action occurs as in three-dimensional QCD, see Eq. (1.4) in [47] 8 . We will see in the ensuing discussion that the pattern of chiral symmetry breaking does not depend on α β , as long as the integrals are convergent, which is certainly the case for any value of α β < N/(N + 4N f /β) < 1. The slightly smaller bound than 1 avoids that the integral over T rA does not diverge (this can be seen after splitting A into its trace and a traceless part).
C. Derivation of the sigma model at large N The derivations of the partitions functions (8), (9) and (10) are similar, and we, therefore, outline only the β = 2 class here. The procedure follows standard steps [2,67]. First we linearize the squared trace term at the expense of a new Gaussian integral over an auxiliary variable x and, afterwards, shift A → A − x1 N to eliminate the linear term in A. This makes it clear that the partition function of our model is nothing but a reweighted integral of the ordinary Gaussian matrix model 9 Upon rewriting the determinant in terms of Grassmann variables, one can easily integrate out A, The quartic term in the fermions can be brought into bilinear form by means of the Hubbard-Stratonovich transformation. For this purpose we introduce an auxiliary Hermitian 2N f × 2N f matrices H. This allows to integrate out the Grassmann variables leading to the result After shifting H → H + ix1 2N f via analytic deformation of the contours, we perform the x-integral for which we need the stricter bound α 2 < N/(N + 2N f ). Then, we obtain This is the finite N result that is still exact without any approximations. When employing the choice (7) the parameter α 2 simplifies, i.e., α 2 = 1. In the following we keep α 2 fixed. Let us consider the large-N limit with M = N M fixed. The integral is dominated by saddle point manifolds and fluctuations around them. When diagonalizing H = U † ΛU with a real diagonal matrix Λ = diag(λ 1 , . . . , λ N f ), the saddle-point equation in the chiral limit M = 0 reads In general, there are multiple real solutions to this equation. However, we look for the minimum of the real part of S Λ which is achieved for Λ = Λ k , cf. Eq. (11), when k ∈ Z with |k| ≤ N f . Indeed, the lower bound to the real part is which is saturated only by Λ k and permutation of its diagonal elements. The second inequality follows from the fact that λ 2 i /2 − log |λ i | is a concave function with its two minimums at λ i = ±1. The fluctuations about Λ k give an overall constant, that comprises the sign (−1) N k from det N Λ k , and can be incorporated exactly. As a result we obtain to first order in M , where we exploited that Λ −1 k = Λ k because k is an integer, cf. (11). This result hold regardless of the value of α 2 as long as α 2 is fixed. The latter implies that 1 − α 2 is of the order 1/N . The proper normalization of this partition function is computed in Appendix A 2.
The result (30) realizes the symmetry breaking pattern (12) for the β = 2 class. The integration for U is effectively over the coset ]. This modified symmetry breaking pattern is evidently enforced by the squared trace term in (26). The idea of using squared trace terms to constrain the symmetry realization is similar to the squared trace deformation in Polyakov-loop models [72][73][74]. The symmetry breaking patterns for β = 1 and 4 in (12) are realized by the same mechanism.
Spectral sum rules for the eigenvalues of the Dirac operator can be derived by matching the quark mass expansion of the effective finite-volume partition function (30) with that of the partition function in QCD 3 . They have already been obtained in eq. (5.17) of [75] for general k from a mathematical perspective. Adapting [75] to our convention, we obtain for a few low order sum-rules. Here {iζ n } are the Dirac eigenvalues rescaled by the average level spacing with ζ n ∈ R, while · · · k denotes the average with respect to the QCD 3 action with a Chern-Simons term at level k. These sumrules are a direct consequence of the pattern of spontaneous symmetry breaking and are independent of the specific details of the random matrix model.

D. Phase transitions at non-integer k
As the symmetry breaking pattern changes when k is shifted with unit increment, there must be a phase transition at non-integer values of k for |k| < N f . 10 To determine the locus of phase transitions, we have to solve the 2N f coupled saddle point equations (28).
In the first step we consider as the (2N f + 1)'st variable. Then, the equations (28) for λ 1 , . . . , λ 2N f decouple and can be solved in terms of the auxiliary variable λ 0 , yielding with L n = ±1 a sign which is not fixed yet. The sum of these signs is denoted by 2k L = 2N f j=1 L j , which plays the role of 2k for a non-integer k. To obtain the solution for λ 0 , we sum over all n in (33) and find This equation has a unique real solution, because the right hand side plus 2λ 0 /α 2 is strictly monotonically increasing. The unique solution is Summarizing, equations (33) and (35) yield all 2 2N f saddle points. The solutions only depend on the still free This quantity has to be minimized in the integer k L = −N f , . . . , N f . When k is an integer with |k| ≤ N f , (36) has a unique minimum at k L = k when λ 0 = 0. Note that this minimum is completely independent of N f and α 2 . Thus the discussion is valid for any number of flavors. The N f and α 2 dependence enters the game when we studying the phase transition point with a real-valued k. Then we have to compare the actions S ( k ) Λ and S ( k ) Λ with the floor function . and the ceiling function . yielding the largest integer smaller than or equal to k and the smallest integer larger than or equal to k, respectively. The phase transition then happens when This is a transcendental equation in k which always has one solution in each interval [j, j obviously exhibits a kink in the parameter k L . Therefore these phase transitions are of first order.
Phase transition points for general N f are located symmetrically on the positive and negative sides of the k axis, and hence it is sufficient to look for solutions to (37) for k > 0. Table I is a summary for 1 ≤ N f ≤ 5 with the Dyson index β = 2 and the "Chern-Simons" coupling α 2 chosen as in (7). Only in the limit N f → ∞ do the half integers k = n + 1 2 (n ∈ Z) become the phase transition points. For a finite number of flavors we get corrections, which can be computed via a large-N f expansion of (37) and the solution k = (n + 1/2) ∞ j=0 c j /(N f + 1/α 2 ) j . Assuming n = O(1), the corrected transition points are for n ∈ Z. The residue O (N f + 1/α 2 ) −6 may depend on n, as well, though it seems to be a very weak dependence.  I. Location of first-order phase transitions in each interval of k for the Dyson index β = 2 and the choice α2 = 1, in particular we have chosen Eq. (7) for the "Chern-Simons" coupling α2. The critical point in the interval [N f , N f + 1] is obtained by solving (37) for N f + ε with ε → 0.

E. Partition function at finite and large N
In this section, we evaluate the partition function at finite N and use this result to derive its large-N limit. This discussion serves two purposes. First and foremost, the finite N results enable us to study the approach to the thermodynamic limit. Second, it provides an independent consistency check of the large N result (30).
The partition function (21) is a one-parameter integral over a GUE partition function with 2N f flavors. The latter will be rewritten using the identity [76,77], where the integration is over hermitian N × N matrices, and we employ the notation ∆ L (.) for the Vandermonde determinant with the exceptional case ∆ 1 ≡ 1. The kernel is given by with the monic polynomials where H (x) = e x 2 (−∂) e −x 2 are the Hermite polynomials. The polynomials P (N ) (s) satisfy the orthogonality relation and the normalization constant is given by In our case we have 2N f flavors, each with its own mass m f , so that we have no natural splitting into two sets of masses. We choose λ j = x + im j and µ j = x + i m j ≡ x + im N f +j for j = 1, · · · , N f . Applying (39) to the partition function (21), we obtain, up to irrelevant normalization, where we used the translation symmetry and homogeneity of the and likewise for m. The sign (−1) N f results from pulling a factor i out of each determinant in (21). In the simplest case k = 0, due to the the Gaussian factor, the variable x scales as O(1/ √ N ) while the masses m and m are of order O(1/N ). When exploiting the asymptotic form, Thereupon, the integral over x factorizes, and we are left with the simpler expression, For a parity-invariant mass m = −m it coincides with the partition function obtained previously [44,45]. One can also obtain the result (47) from the low-energy limit of the partition function given in (8) for k = 0. For this purpose we exploit the parameterization [78,79] and the Lagrangian takes the form The integrals over U 1 and U 2 are each Harish-Chandra-Itzykson-Zuber integrals [80,81], and similar for m. The remaining integrals in over Φ can be performed with the Andréief identity [82] yielding (47). Let us examine the simplest case of two flavors ( after substituting x → 2/N x − im. This integral can be carried out with the help of a formula in [83,Sec. 7.374, eq. 9], with the result FIG. 2. The free energy for two flavors with equal mass in the chiral limit for varying N for zero mass (left) and different masses for N = 100 (right). The value at k = 0 is subtracted for each N . The spacing between curves in the right plot has been adjusted by hand for better visibility.
This form is suited for fast numerical evaluation. Since Z(m; k) = Z(−m; −k), one may assume m ≥ 0 without loss of generality. The k dependence of the partition function can now be investigated via (54). Figure 2 (left) displays the evolution of the free energy with N in the microscopic limit. Evidently there are quite strong deviations from the large N limit even for quite moderate matrix sizes like N = 50 which usually yields close to perfect agreement for the microscopic level density. As N grows, there appear kinks that get sharper. They represent first order transitions in the thermodynamical limit. The region around the origin is the flavor symmetry broken phase U(2) → U(1) × U(1) with massless Nambu-Goldstone modes. The other two regions are the symmetry-restored phases. The right plot of Figure 2 illustrates the shift of of the phase transition points as a function of the masses that are of order O(1) instead of O(1/N ), and hence they are outside the microscopic domain. The two kinks move towards negative k (positive k) for m > 0 (m < 0), respectively. Figure 3 (left plot) shows the mass dependence of the free energy at k = 0, again at an N independent mass, i.e., N m 1. The two pronounced kinks at m ≈ ±0.34 indicate a strong first-order transition that corresponds to the passage of a kink over the origin in the right plot of Fig. 2. The middle region corresponds to a symmetrybroken phase with massless Nambu-Goldstone bosons 11 whereas the outer regions are symmetric gapped phases. Such phase transitions at nonzero masses were argued to exist in [47] and our matrix model serves as a toy model for this phenomenon.  Finally, in Fig. 3 (right) we show the k dependence of the quark-antiquark condensate defined as for three different values of the mass, m 1 = m 2 = 0.0, m 1 = m 2 = 0.25 an m 1 = m 2 = 0.75. Hence, also in Fig. 3 we do not show the microscopic limit of the chiral condensate. The masses are of order O(1) and not O(1/N ). In the microscopic limit the two kinks are at exactly the same position as in Fig. 2 (left). The condensate has still a discontinuity at the values of k for which the free energy has a kink for masses which are of order one. Yet, for increasing mass the kinks move to infinity, approximately as ∼ 3 2 m, and the quenched result is recovered for m → ∞. For masses with opposite sign, m 1 = −m 2 , the k-dependence of the free energy remains symmetric about zero, but the kinks move away from zero, again like ∼ 3 2 m for large m, so that for m → ∞ the quenched limit is recovered.

III. SPECTRAL CORRELATION FUNCTIONS
This section is mostly devoted to the level density, the quark-antiquark condensate and their microscopic large-N limit, though in Subsection III A we also study all k-point correlation functions at finite N . In this subsection we obtain an exact expression for the spectral density at finite N . To illustrate what happens when N is increasing at fixed Chern-Simons coupling k, we consider the quenched level density (N f = 0) at finite N in Subsection III B. This result is amenable to a saddle point approximation which allows us to obtain the microscopic limit of the spectral density (see section III C) defined as For more flavors the microscopic limit can be derived much more easily from an expression where the spectral density is given by a ratio of an N f + 2 flavor partition function and an N f flavor partition function as is discussed in section II E. Starting from this result we obtain in section III D the microscopic spectral density for N f flavors. In particular, we will work out the one-flavor case in detail.
For N f = 1 the partition function is given by the sum of three saddle points (see section II D) Because of the reweighted structure of the partition function, the spectral density is given by with the level density corresponding to the saddle point k L given by Since the partition functions behave as where f k L is the free energy of the k L 'th saddle point, we have that in the large-N limit, the spectral density is dominated by one saddle point (unless k is exactly at the phase transition point). To derive the results of subsection III D, we make use of results obtained in Appendix A.
A. General Nf at finite N The matrix model approach gives us a way to investigate spectral fluctuations of the Dirac operator in QCD 3 . The simplest way to compute the n-point correlation functions R(λ 1 , . . . , λ n ; M ; k) of the matrix A in the ensemble (3) would be to make use of the reweighted structure (21). If R(λ 1 , . . . , λ n ; M ; x) is the n-point spectral correlation function of A for the GUE ensemble Z N,N f (M − ix) with 2N f flavors, we simply have .
The spectral correlation functions corresponding to the partition function Z N,N f (M − ix) have been computed in [35,36] for the massless case, and in [37,38,44,45] for nonzero masses with pairwise opposite signs. By substituting the results from [37,38,44,45] into R we obtain the spectral functions for our matrix model. In this paper we only explicitly work out the one-point function (spectral density) for arbitrary 2N f masses. Recalling the shift A → A − x1 N , we get from (22) The remaining integral is nothing but the partition function of GUE with 2N f + 2 flavors. Thus, it can be expressed as an integral over a hermitian where we have defined m 2N f +1 = m 2N f +2 = −iλ. This matrix integral can be computed with the help of (39). In doing so, the degeneracy of m 2N f +1 and m 2N f +2 must be lifted slightly to avoid an apparent singularity in (39). Labeling and combining (39), (44), (61) and (63), the spectral density is finally obtained as where we employ the notation The overall normalization of R is fixed by Hence, Re R(λ; M ; k) is an even function of λ while Im R(λ; M ; k) is odd. Another representation of the level density, which is convenient for the derivation of the microscopic limit, is .
(68) This result is obtained by shifting back to A → A + x1 N and B → B + x1 N −1 and then integrating over x in the denominator as well as in the numerator. The normalization follows from integration over λ and combining this integral with the B-integral to the A-integral in the denominator. We recall the value of α 2 = N/(N + 2N f + 1), cf. (7). Indeed this result could be directly derived from the partition function (3), by setting one of the eigenvalues of A equal to λ.
When additionally shifting B → B + λ 1 N −1 with λ = λ/(2N f + 2) in Eq. (68), we reduce the level density to a quotient of two almost identical partition functions, B. Quenched limit at finite N In the quenched limit N f = 0, equation (65) reduces to after shifting x → x − λ. Using the orthogonality relations for the Hermite polynomials, the normalization can be easily evaluated. The x integral may also be performed with the help of the formula [83, Sec. 7.374, eq. 9], leading to the result Since H 2N −2k is an even function, R(λ; k) = R(−λ; −k), and one can assume k ≥ 0 without loss of generality.
As can be seen from Figure 4 for the quenched spectral density at k = 0, in contrast to the standard GUE, the oscillatory structure of the spectral density due to peaks of individual eigenvalues is not present even for small N . This feature was also seen in other one-parameter-reweighted ensembles [70,71]. We expect that this feature will carry over to three-dimensional QCD as well. This figure also shows that the large N -limit, given by the semi-circle ρ(x) = 1 − (x/2) 2 /π, is already well-approximated for N = 100. When increasing k for a fixed matrix dimension N , say N = 20, the spectral "density" becomes complex-valued. We illustrate this in Figure 5 where the real and imaginary parts of the level density R(λ; k) are shown. At nonzero k > 0, the semi-circle undergoes a dramatic deformation of its shape. First, small oscillations at the two edges appear. They even change the sign of the spectral density for small regions regardless of how small k is. The amplitude of these oscillations grows with k. While keeping k below a threshold k c , see the ensuing subsections, the oscillations die out around the origin and we can expect a well-defined microscopic limit. Yet, when increasing k beyond k c , the oscillations intensify and move into the bulk of the spectrum, see Figure 6, such that even there the spectral density does not remain strictly positive. The amplitudes of the oscillations grow rapidly with k, even though the normalization condition dλ R(λ; k) = N is strictly satisfied. A similar oscillation of the spectral density was also observed in matrix models for QCD at nonzero chemical potential [8,[84][85][86] and for QCD with nonzero theta angle [87,88].
The question is how this oscillatory behavior carries over to the large-N limit while keeping k fixed. Three things may happen. Either the oscillations do not reach the origin; then we expect the universal results from GUE with a possible reweighting since the level spacing is changing. Second, the oscillations reach the origin but are not strong enough to make the microscopic limit ill-defined, in particular the amplitude does not grow with the matrix dimension N . And third, the oscillations become so dominant that the microscopic limit is not well-defined at the origin. The latter will usually happen at about k c ≈ N f + 1/2 as we will see below.  7. r(N, x) and its large-N approximation (76) for N = 20.

C. Quenched microscopic large-N limit
The next task is to evaluate the microscopic limit of the quenched density (70) where we evaluate its large-N limit at fixedλ ≡ λN . Incorporating the normalization factor, we write (70) as where and which is normalized as ∞ −∞ dx r(N, x) = 1. While r(N, x) is similar to that of the Wigner-Dyson ensemble, the N dependence is slightly different. However, for large N it also approaches a semi-circle. To obtain more quantitative results we use the uniform asymptotic expansion of the Hermite polynomials [89] valid for |x| < √ 2N . For large N , we thus obtain In Figure 7 we compare (74) and (76) for N = 20. The agreement is excellent except near the edge of the semi-circle.
Returning to the quenched density, it is comprised of two pieces associated with the two terms in (76), where R (a) (λ; k) corresponds to the semi-circle part and R (b) (λ; k) to the oscillatory part. The prefactor 1/N results from the scale on which we want to zoom in about the origin.
To evaluate the first contribution at large N we note The saddle point can be approximated as x = 2ik since λ = λ/N , so The average over the oscillatory part can be evaluated as We have dropped the term exp[−N λ 2 /2] since λ ∼ 1/N in the microscopic large-N limit. Furthermore we introduced the function For large N the integral can be evaluated with the saddle point method. Solving yields the solutions To check this, we want to point out that which forbids a saddle point of f + (x) for k > 0 and of f − (x) for k < 0.
The saddle point expansion leads to The explicit form of the function f −sign(k) (x c ) is given by To derive this intermediate result we used the identity Let us collect all results, so that the microscopic level density reads This indicates that the amplitude of the oscillation is controlled by 1 + k 2 + 2arcsinh[(k 2 − 1)/2|k|], which is shown in Figure 8. The function changes sign at k = ±k c with Thus, for N 1, the amplitude of the oscillations grows exponentially for |k| > k c , but it dies out for |k| < k c . In the limit k → 0 we smoothly recover the well-known microscopic spectral density of GUE, ρ mic ( λ) = 1/π, see [35]. It is intriguing to note that the value of k c in Eq. (90) coincides with the phase transition point obtained from (37) in the limit N f → +0. This observation is not surprising when considering an alternative derivation given in the Appendix A 3.
In Figure 9 we numerically compare (89) with the exact density (71) for various k. In all cases they show good agreement in the region λ ∼ 1/N despite the relatively small matrix dimension N = 20. When increasing N for |k| > k c the oscillations become dominant; the amplitude grows exponentially, and a microscopic limit does not exist.
The quark-antiquark condensate in the quenched case can be readily calculated since the microscopic spectral density is R( λ; k) = √ 1 + k 2 /π. Hence, the quark-antiquark condensate is equal to This result only holds for |k| < k c .

D. Unquenched microscopic level density
To derive the microscopic level density with dynamical quarks, we start from Eq. (69) where we have to compute two kinds of partition functions. For this purpose, we first need to approximate the prefactor in Eq. (69) which in the microscopic limit simplifies to We recall that in this limit λN ≡ λ and m f N ≡ m f with λ and m f fixed in the limit N → ∞. The two partition functions in the numerator and denominator of (69) are computed in detail in the Appendix A with the aid of random matrix methods. Let us briefly revisit the quenched density. To obtain this quantity we combine Eqs. (A3), (A30) and (92), the latter two for N f = 0, in (69). Since the partition function Z N f =0 = Z q is always in the trivial phase k L = 0 for |k| < k c ≈ 0.527697 and M = −i λ/21 2 , we obtain Since λ ± = k ± √ 1 + k 2 in the present case (see eq. (A12) for the definition for arbitrary N f ), we can simplify this expression to This is the leading order term in Eq. (89). The oscillatory part which becomes dominant for |k| > k c can be obtained from Eq. (A30) in the limit k L → 1 instead of k L = 0. The limit is important since some terms seem to diverge; nevertheless they cancel with other terms that vanish so that one needs to employ l'Hospital's rule several times. In the case of dynamical quarks, we collect the terms in Eqs. (A18), (A30) and (92) and find Now, we combine the saddle point solutions λ ± into and their sum, and we perform the two Harish-Chandra-Itzykson-Zuber integrals [80,81] in (95) which yields the simplification Thence, the density for k L = 0 is essentially the one with no Chern-Simons term. Interestingly, at the phase transition points the sign of the level density jumps. Especially in the case of k being an integer, the result reduces to Here, we want to underline as before that we assume that |k| is smaller than a critical value k c > N f + 1/2. Above this value, the microscopic level density is governed by the oscillations as in the unquenched system, and a microscopic limit does not exist. The corresponding oscillatory part can be obtained by choosing k L = N f in Eq. (A18) and k L = N f + 1 in Eq. (A30). The amplitude will again grow exponentially with N . The critical value k c for a fixed number of flavors 2N f can be obtained from (37) by comparing k = N f and k = N f + 1.
Finally we want to present the result for two flavors (2N f = 2). For the phase k L = 0, the level density is equal to with ω = 1 + k 2 /4, cf. Eq. (98), and ρ 2N f =2 mic (k = 0) being the microscopic level density without a Chern-Simons term. Interestingly, the whole spectrum is only rescaled by the factor ω, in particular the mean level spacing is not anymore π but π/ω. We recall that |k| has to be smaller than the critical value of 0.50551 . . ., see Table I. When m 1 = −m 2 = m the density is evidently real It is shown in the left plot of Figure 10 for two different quark masses. When taking the masses to infinity m → ∞ the quarks decouple and we recover the quenched case ρ mic ( λ) = 1/π as expected. The corresponding quark-antiquark condensates for m = m 1 = m 2 and m = m 1 = − m 2 are in this phase equal to and respectively. Since only the normalization constant in front of the level density changes with k inside this phase (note the unfolding with ω), we have chosen k = 1. We recall that the imaginary part is an odd function around the origin while the real part is an even function.
In the nontrivial phase k L = 1 (the case k L = −1 is very similar), the level density has the form ]. This spectral density has always a non-trivial imaginary part even when we set m 1 = −m 2 = m in which case it simplifies to From these results we can read off several things. First of all, the level density exhibits complex oscillations in the phase k L = 1, cf., middle and right plot in Figure 10, which also holds for any |k L | > 0 phase when considering even more flavors. Additionally, the amplitude of these oscillations grows exponentially in the quark masses. Therefore, we cannot expect that the quenched limit exists in this phase; especially the reduction of the number of flavors does not work anymore. Finally, the change from one phase to another, say k L → k L + 1, drastically changes the microscopic spectral density. There is not a smooth transition of the microscopic level densities in the various phases and it completely breaks down when |k| crosses a critical k c , see discussion above.
The quark-antiquark condensate as a function of k readily follows for from the partition function which for 2N f has three components The explicit expressions are given in Eq. (A18). In the thermodynamic limit only one of the three components contributes to the condensate depending on the value of k. For m 1 = m 2 the integrand in Eq. (A19) does not depend on U resulting in a pure exponential mass dependence. In the microscopic limit, we thus find a mass independent chiral condensate given by with k c = 0.50551 (see Table I). Up to 1/N corrections, this result (black curve in Fig. 11) is in agreement with the k dependence of the condensate for close to massless quarks obtained from the exact partition function (54)  in Fig. 11) which coincides with the result from the microscopic partition function (A18) (red curve in Fig. 11). The small discrepancy between the last two curves and the analytical result is due to 1/N corrections -taking the quark masses closer to zero does not change the curves.

IV. CONCLUSION AND OUTLOOK
We have constructed a random matrix theory for QCD in three dimensions (QCD 3 ) with a Chern-Simons term of level k that reproduces the pattern of spontaneous symmetry breaking according to U as proposed recently by Komargodski and Seiberg [47]. This random matrix model is an extension of the random matrix for QCD 3 without a Chern-Simons term (k = 0). The Chern-Simons term of the random matrix model is in some aspects different in character from the Chern-Simons term of QCD in 3 dimensions but agrees in other aspects. In particular, the level k is not quantized, which is not surprising since the random matrix model does not have a local gauge invariance. However, the effect of the Chern-Simons term on the eigenvalues is similarit adds a phase proportional to k to the phase of the fermion determinant. It is remarkable that, in all cases we know of, random matrix theories with global symmetries of QCD-like theories reproduce their pattern of spontaneous symmetry breaking and break the symmetry in such a way that the corresponding condensate has the maximum global symmetry, see Ref. [90]. The present work shows that a complex action can violate this feature, even in the case of random matrix theory.
What we have learned from earlier work on random matrix theory with a complex action is that the imaginary part of the action can move the phase boundaries of the phase quenched theory. For QCD at nonzero chemical potential the phase of the fermion determinant moves the critical chemical potential of half the pion mass to 1/3 of the baryon mass. For QCD at nonzero theta angle, the chiral condensate does not change sign when one of the quark masses does not change sign. Keeping this in mind, it is not unexpected that the phase due to the Chern-Simons term can change the phase of the theory: the imaginary part of the action nullifies the leading phase so that the subleading phase becomes dominant. At k = 0, the phase with the standard pattern of chiral symmetry breaking is canceled, so that phases with asymmetric breaking of spontaneous symmetry breaking becomes dominant.
Appendix A: Derivation of some partition functions In this appendix we work out the explicit computation of the two partition functions in Eq. (69). The one in the denominator has to be dealt separately for the quenched (Subsection A 1) and unquenched (Subsection A 2) ensemble while this distinction is not relevant for the partition function in the numerator, which is evaluated in Subsection A 3.

Quenched A integral
We first consider the quenched partition function and linearize the squared term by introducing an auxiliary x-integral, Next, we shift A → A − x/N 1 N and integrate over A, Finally, we perform the integration over x and arrive at The next quantity we consider is the unquenched partition function As in the previous section, we first linearize the squared term with the help of an x-integral and then introduce a Gaussian integral over a complex Grassmann valued N × 2N f matrix V , Here, we used the anticommuting property of Grassmann variables. The integral over A can again be performed after In the next stage, we integrate over x and find Now we are ready to apply the bosonization formula [91][92][93] and replace V † V by N U with U ∈ U(2N f ). The scaling factor N is chosen for convenience of the saddle point analysis. Thus, we have where M = N M is fixed in the microscopic limit. Next we diagonalize the matrix U = U † zU with z a diagonal matrix of complex phases, Let z k be one saddle point solution of the saddle point equation and maximizing the integrand. The relation to Λ k , discussed in section II D, is z k = Λ −1 k which indeed yields the saddle point equation (28) after plugging this relation into (A11). Therefore, we already know that there are points satisfying these two conditions, where k L is either the integer above k or below k depending on the phase the system is in. In particular we can choose z k = diag(λ −1 We underline that these solutions are always real because k 2 − k 2 L ≥ −|k + k L |/2 ≥ − |k| ≥ −N f . Moreover we have always λ + > 0 and λ − < 0.
The contribution from the fluctuations about the saddle point can be obtained from the expansion where the phase pre-factors reflect the direction of the original contour. Then, the measure transforms as follows Exploiting the identity λ + λ − = −1, one can explicitly write In the next step, we expand the determinant and the exponent where we employed the saddle point equation. The mass dependent term is independent of the massive modes δz ± . As it should be, the exponents proportional to √ N cancel each other so that we are left with the δz ± integrals that can be computed as follows . (A17) Now we are ready to put everything together and apply Sterling's formula. We eventually arrive at In the case when k is an integer, meaning k = k L and λ ± = ±1, this result simplifies drastically to Note that the prefactors that only depend on N and N f could have been absorbed in the definition of the partition function. We would like to emphasize that the partition function for integer as well as non-integer k is real, as can be already seen from its definition, and even positive when omitting the factor (−1) (N f +k)N . The latter can be readily achieved by choosing an even matrix dimension N .

The B integral
In this section we evaluate the integral in the large-N limit. The computation proceeds along the same lines as for Z N f . However, we cannot easily carry over the entire result for N → N − 1 and N f → N f + 1 since the standard deviations do not change in the same way.
Again, we introduce the auxiliary real variable x to linearize the squared trace term and the complex Grassmann valued matrix V , albeit now it has the dimension (N − 1) × (2N f + 2). Collecting the masses in the diagonal matrix M = diag(m 1 , . . . , m 2N f +2 ) + iλ 1 2N f +2 , we find the integral after the integration over B, which is the counterpart of (A7). Next, we apply the bosonization formula [91][92][93] and integrate over the variable x so that we obtain with M = N M . When diagonalizing U = U † zU with z a 2N f + 2 dimensional diagonal matrix of complex phases we obtain We are ready for a saddle point analysis of thez integral whose saddle point equation is By plugging in the choice z k = diag(λ −1 with exactly the same λ ± and k L as in section A 2, one can easily verify that this is a solution. The real part of the two corresponding actions is apart from an N also the same, especially Thus, the phase transition points from the original A integral and from the original B integral are the same as they should. If they were not equal, the spectral density would be either exponentially small or exponentially large in N . Yet, Y N f , with two more flavors than Z N f , has one additional phase transition compared to Z N f at about k c ≈ ±(N f + 1/2). This suggest the presence of an additional phase transition in the spectral density as compared to the original partition function. Indeed, we have found in section III C that the limit of the microscopic level density does not exist when |k| is larger than a critical value k c . Hence we have to stay always in the regime where the two phases of the partition functions Y N f and Z N f agree. When this is not the case the oscillation will become dominant and a microscopic limit does not exist, see the discussions in sections III C and III D. The expansion works exactly the same as before, i.e., z = z k + i diag(δ z + , −δ z − )/ √ N with δ z + ∈ R N f +1+k L and δ z − ∈ R N f +1−k L . Thence, we have for the measure for the determinant and for the exponent Tr δ z 2 + + Tr δ z 2 − + (Tr δ z + − Tr δ z − ) 2 4(N f + 1) , where we again have used the saddle point equation to simplify the result. The integral over the Gaussian fluctuation δ z ± about the saddle point is given by (λ 2 + + 1) Tr δ z 2 + + (λ 2 − + 1) Tr δ z 2 − + (Tr δ z + − Tr δ z − ) 2 4(N f + 1) .

(A29)
The degeneracy of the saddle points is (2N f + 2)!/[(N f + 1 + k L )!(N f + 1 − k L )!]. Combining all contributions we arrive at the main result of this subsection, and for integer k, implying k = k L and λ ± = ±1, it reduces to