Scale-invariance and scale-breaking in parity-invariant three-dimensional QCD

We present a numerical study of three-dimensional two-color QCD with $N=0, 2, 4, 8$ and 12 flavors of massless two-component fermions using parity-preserving improved Wilson-Dirac fermions. A finite volume analysis provides strong evidence for the presence of $Sp(N)$ symmetry-breaking bilinear condensate when $N \le 2$ and its absence for $N \ge 8$. A weaker evidence for the bilinear condensate is shown for $N=4$. We estimate the critical number of flavors below which scale-invariance is broken by the bilinear condensate to be between $N=4$ and 6.


I. INTRODUCTION
Three-dimensional gauge theories coupled to an even number of two-component massless fermions can be regularized to form a parity-invariant theory. The parity-invariant fermion action for N flavors of two-component fermions coupled to SU (N c ) gauge-field A µ is where C(A) is the two-component Dirac operator, with an ultraviolet regularization being imposed implicitly. In the continuum, with σ µ being the three Pauli matrices. One can think of the above parity-invariant theory of N flavors of two-component fermions to be equivalent to a theory of N 2 flavors of four-component fermions ψ i with a Hermitian Dirac operator D: with the following identifications Since C † = −C in the continuum, the theory has a global U (N ) flavor symmetry. However, the N c = 2 theory is special, and there is a larger Sp(N ) global symmetry following from the property σ 2 τ 2 Cσ 2 τ 2 = C t , where σ 2 and τ 2 are Pauli matrices in spin and color space respectively [1]. This is also related to the fact that the operator D can be made real symmetric in a suitable basis [2].
Consider such a theory on an Euclidean periodic l 3 ph torus. Physics depends on the dimensionless size, = l ph g 2 ph , where g 2 ph is the physical coupling constant. Since g 2 ph has the dimension of mass, these theories are super-renormalizable, and the continuum limit in a lattice regularization at a fixed can be obtained by setting the lattice coupling constant (same as the lattice spacing) to g 2 lat = L on a periodic L 3 lattice and taking L → ∞. The physics of this theory will smoothly cross-over from a non-interacting theory at small to a strongly interacting theory as → ∞; this strongly interacting theory could either be scale-invariant or scale-breaking. Scalebreaking is expected to produce a parity-preserving non-zero fermion bilinear condensate, that for a generic SU (N c ) gauge theory breaks the U (N ) flavor symmetry to U N 2 × U N 2 [3,4], but in the case of SU (2) gauge theory breaks the Sp(N ) symmetry to Sp N 2 × Sp N 2 [1].
The possibility of the generation of any parity-breaking bilinear condensate has been argued against in [4].
A numerical study of the Abelian U (1) gauge theory using Wilson fermions showed no evidence for a bilinear condensate for any even value of N > 0 [5]. Since the Wilson fermions realize the U (N ) flavor symmetry only in the continuum limit, overlap fermions were used to study the theory with the symmetry present even away from the continuum. This enabled us to study the scale-invariant properties of the theory in further detail [6]. Assuming the N = 2 theory to be scale-invariant, a strong infra-red duality [7,8] predicts an enhanced O (4) symmetry and this was also verified using overlap fermions [9]. Since the Abelian theory is scale invariant for all even values of N , we turn to SU (N c ) non-Abelian theory in this paper in order to study a transition from a scale invariant theory to one that breaks scale invariance as one changes the number of flavors. It is worth noting that three dimensional SU (2) gauge theories with even number of massless fermions appear in the study of spin liquids [7,[10][11][12].
The non-Abelian theory in the 't Hooft limit (number of colors going to infinity at a fixed number of flavors) has been shown to have a non-zero bilinear condensate [13]. Dagotto, Kocić and Kogut [14] numerically investigated the SU (2) theory with staggered fermions on small lattices (mainly on L = 8). They tried to deduce the presence or absence of the condensate in the massless limit by computing it at two different fermion masses and then extrapolating it to zero.
The four-dimensional SU (N c ) gauge theory coupled to massless fermions becomes infra-red free and loses asymptotic freedom above a certain number of fermion flavors. In 4− dimensions, -expansion calculation suggests the infra-red free behavior develops into a non-trivial conformal [15,16]. Also, an earlier analysis using the Schwinger-Dyson equations [17] using 1/N -expansion in three dimensions, similar to the one for the Abelian theory [18], suggests that the theory is scale invariant if N > 256 Nc . However, a study of the flow of four-fermion operators in the -expansion [16] suggests a new IR fixed point, different from the one at large-N , might be present if N < 11N c + 12 . It is tempting to identify the lower bounds on N as a critical N * below which scale invariance is broken but one has to consider the possibility that it is a point that separated one class of IR fixed points from another. There is indication of such a scenario in QED 3 where the IR fixed point at N = 2 has an O(4) symmetry unlike the IR fixed point at N = ∞. Therefore, a first principle numerical estimate of an N * that is below any of the bounds obtained from other calculations ( with their own underlying assumptions) has interesting implications.
We are, therefore, motivated to perform a careful numerical study of the N c = 2 theory with N flavors of massless two-component fermions with the sole aim of obtaining the critical number of fermions, N * , such that the theory develops a scale for N < N * and is scale-invariant for N ≥ N * . Since we are interested in numerically studying several different values of N and our aim is only to locate the critical number of flavors, we use improved Wilson fermions [5] in order to reduce the computational cost. A careful study of the transition from above N * to below N * will require overlap fermions [6] or domain wall fermions [19,20]. The type of lattice fermions one uses is irrelevant for the continuum physics; however, the main difference arises at finite lattice spacings where the propagator of two-component Wilson fermions is not anti-Hermitian but the propagator of two-component overlap fermions is anti-Hermitian. We reserve the computation of using U (N ) symmetry preserving lattice fermions for the future, which will give us more information on the fermion bilinear correlators and serve the purpose of cross-checking the continuum results in this study.

LESS DIRAC OPERATOR
The basic philosophy of the finite volume analysis in this paper is the same as in [5]. The eigenvalues of the Hermitian Dirac operator D in Eq. (4) occur in positive-negative pairs ±λ given by the equation The eigenvalues λ are gauge invariant and discrete at finite . Therefore, we study the ordered discrete spectrum of eigenvalues of √ C † C: obtained as observables in the N flavor theory of massless fermions. Henceforth, we will use λ i to denote the expectation value of the i-th lowest eigenvalue of √ C † C over the gauge-field path integral. The asymptotic behavior of λ i ( ) as → ∞ falls into three types: Type 1: Free field behavior will result in λ i ( ) ∝ 1 with the proportionality constants determined by the mean value of the gauge field which is equivalent to the induced boundary conditions on the fermions. Such a behavior is certainly expected in small due to the asymptotic freedom in the theory. Type 2: A complete level repulsion between the eigenvalues will result in λ i ( ) ∝ 1 3 with the proportionality constants determined by one value of the bilinear condensate Σ that breaks scale-invariance.
Type 3: In an interacting scale-invariant theory, the eigenvalues which have the naive dimension of mass, will scale as λ i ( ) ∝ 1 1+γm with 0 < γ m < 2 being the mass anomalous dimension. Here, we are assuming there are no further restrictions in three dimensions on the possible values γ m can take between 0 and 2.
Since the theory at very small will be non-interacting, we expect a smooth cross-over from the 1 behavior at small to one of the above cases as → ∞ 1 . As an abuse of notation, we will use γ m = 2 to mean the scale-breaking type (2) behavior even though the exponent is no longer an anomalous dimension in this case.
We have summarized the numerical details pertaining to the lattice simulation and the extraction of the low-lying spectrum in Appendix A. Having chosen a numerical approach to study the theory as a function of , we have to face its limitations. The lattice spacing is given by L and in spite of improving the lattice action we have to face the fact that as gets large we have to make L also large to reduce lattice spacing effects. Furthermore, if L > 2.3 we will be in a lattice theory with strong coupling and there is a bulk cross-over [22,23] that separates an unphysical strong coupling phase from the continuum phase that we want to study. Using the computing resources available to us, we were able to go up to L = 28 in a theory with N ≥ 2, and up to L = 32 in the quenched limit (N = 0). If we now ask for acceptable levels of lattice spacing effects and also require that we are in the continuum phase of the theory we are led to study the theory for ≤ 17 and this is what we have performed here.
At this point, it is appropriate to make a few remarks about the analysis performed earlier in [14]. The lattice coupling constant in [14] is defined as β = 4 present, one could put their observations on a firmer footing by following a similar approach 1 One could have taken a different approach and kept the physical extent in one of the directions fixed at a value 1/T while taking → ∞ in other directions, in which case we would be studying the theory at finite temperature T , and there might be singular behavior around some T c . We are not taking that approach here (c.f. [21] for such an approach in three-dimensional quenched SU(3) theory).
by simulating a wide range of fermion masses m at different lattice volumes in order to take a controlled thermodynamic as well as the massless limits. In this case, one should also follow an unbiased approach by assuming a Σ(m) ∼ m 2−γ + O(m) mass dependence of the condensate at infinite volume [24], in order to allow for the massless theory to be scale-invariant. However, we use the finite-size scaling of eigenvalues of the massless Dirac operator to determine the presence or absence of condensate along the lines of our previous studies of QED 3 . This method is also advantageous because we avoid the presence of two scales and m at once in the problem.
If the theory has a bilinear condensate, the behavior of the eigenvalues as a function of will smoothly cross-over from the type (1) to type (2) as we increase . The type (2) asymptotic behavior will be given by where Σ is the value of the bilinear condensate per fermion flavor at = ∞ and z i 's are universal numbers given by an appropriate random matrix model [1,25]. If 1/ √ Σ sets the typical spontaneously generated length-scale in the system, then this cross-over to the asymptotic behavior happens only for box sizes ∼ O 1/ √ Σ , which renders the measurement of small values of Σ computationally difficult. As we get closer to the critical number of flavors N * from below, the value of the bilinear condensate would get smaller, making it more difficult to decide if the theory has a non-zero bilinear condensate or not. If the theory does not have a bilinear condensate the behavior of λ i ( ) will smoothly cross-over from type (1) to type (3) as is increased and the asymptotic behavior will set in early in if we are well above N * since we expect γ m to approach zero (free field behavior) as N → ∞.
With the above picture in mind as we change the number of flavors, we will start by defining where the estimate of λ i is made by averaging over the eigenvalues measured in different gauge configurations sampled by Monte Carlo. We use z i simply to scale the right-hand side and we make no assumption about the presence of a non-zero bilinear condensate; for the type (3) conformal behavior of λ, Σ( ) as defined above will approach 0 as −2+γm . We fit the data for where the three fit parameters, a i k (γ m ), k = 0, 1, 2 depend on the choice of γ m and i. Note that our choice implies that a i 0 (γ m ) cannot be zero implying that the term in front of the parenthesis  is the asymptotic behavior. By studying the χ 2 per degree of freedom (χ 2 /DOF) of the fits as a function of γ m , we will be able to find the value that best fits the data for a given i. We expect the best value of γ m , where the χ 2 is minimized, to be independent of i. As per our discussion in the previous paragraph, we expect this approach to work with relative ease for value of N away from N * . Contrary to the form used in Eq. (7), we can also use where we are assuming the possibility for a non-zero condensate Σ, with the finite volume correction that is Taylor expandable in 1/ . This form does not allow for an anomalous dimension, which is indeed the case when there is a condensate. Since we will have results for λ i ( ) in a finite range of , both Eq. (9) and Eq. (10) should result in the same physical conclusion well away from N * . But, we expect conflicts between these two forms closer to N * . In both the ansätze, we could have included more orders of 1/ corrections. But empirically, we find that 1/ and 1/ 2 corrections are enough to describe our data well. Therefore, our conclusions have to be interpreted in a Bayesian sense -given the priors that only 1/ and 1/ 2 finite volume corrections are important in the data we have, and assuming this continues to be the case at even larger where we do not have the data, we ask for the probable values of the anomalous dimension or the condensate.
For N < N * , we will be able to further substantiate our results via a comparison with the appropriate random matrix models. Random matrix models appropriate for describing low-lying eigenvalues in a three dimensional gauge theory coupled to massless fermions that generates a non-zero bilinear condensate can be found in [25] where the fermion operator C for each fermion flavor is realized as a random anti-Hermitian matrix. Under parity, C → C † and therefore these random matrix models will be parity-invariant for even number of flavors since the Haar measure of a random Hermitian matrix is parity-invariant. Since our gauge group is and it is assumed that M is taken to infinity. Ordering of the eigenvalues of C is according to the absolute value of C since this matches with the definition in a parity invariant theory.
We opted to numerically simulate the random matrix model, and the universal numbers, z i , appearing in Eq. (7) are the averages of the eigenvalues so ordered. We have listed their values in Table I.

III. RESULTS
We present the results of our numerical analysis in this section. We have given the list of our simulation points ( , L, N ) along with the corresponding averages of the first four eigenvalues of the massless Hermitian Wilson-Dirac in Tables-III,IV  We fit Σ i ( ) as a function of using Eq. (9) with different values of the exponent γ m . In Figure   2, we show the χ 2 /DOF for these fits as a function of γ m for different N , as labelled on top of each panel. The different curves in the panels correspond to χ 2 /DOF of the fits to the four different Σ i . In order to easily interpret the plot, a rule of thumb is that the fit using a value of γ m is good if its χ 2 /DOF is about 1, while it being 2 or above is indicative of the fit describing the data poorly. In the different panels, the limit γ m → 2 points to a theory with a bilinear condensate and the limit γ m → 0 points to a free field behavior. even larger where we do not have the data. If one assumes a well-behaved 1/ expansion with successively smaller higher order terms, with no cross-over from one type of leading behavior to another, one would favor the smaller of the allowed values of γ m , which is around 0.4 to 0.5 for both N = 8 and 12. At all flavors, the allowed range of γ m as deduced from Σ 1 is broader than as allowed by other Σ i . The most likely cause is that the behavior of the lowest eigenvalue is affected the most by the need to fine tune the Wilson mass to realize massless fermions on the lattice as explained in Appendix A. This is enhanced at N = 4 as is evident from the flat behavior of the χ 2 for i = 1 at N = 4. For N = 4, the range of allowed γ m as deduced from the other i are also broad and includes γ m = 2, and hence we are unable to rule out scale-breaking in this case. In the following subsections, we analyze the different N separately.
A. N = 0 Using the four values of z i for N = 0 in Table I, we obtained Σ i ( ) from Eq. (8). In Figure   3, we have plotted Σ i ( ) as a function of 1/ . One can see that the asymptotic behavior has set in for ≥ 14 with a finite value of the bilinear condensate in the → ∞ limit. In the same plot, we have also shown the error bands of fits to the data using Eq. (10). We find the values of Σ i extrapolated to infinite for i = 1, 2, 3, 4 to be 0.0165 (8), 0.0154(7), 0.0147(5) and 0.0143(4) respectively. It is reassuring that the values of Σ i for different i converge to about the same value within errors in the → ∞ limit. By taking the average over all the four values of Σ i ( = ∞), we quote the value of the bilinear condensate per color degree of freedom for with the error being purely statistical. We take the spread of values in the four different Σ i to be a measure of the systematic errors in the various extrapolations, and we conservatively estimate this systematic error to be about 0.0011. In the N = 0 quenched theory, the fermions are used merely as a probe with no back-reaction on the gauge fields. Since the pure SU (2) theory does have a scale, it is not a surprise to find a bilinear condensate in this case. We remind the reader that the above value of condensate is dimensionless and in units of the coupling g 2 ph . When the above value of condensate per color is measured in units of the 't Hooft coupling, N c g 2 ph , it is roughly a factor of 2 smaller than the corresponding value in the 't Hooft limit in [13]. It is also interesting to note that a linear extrapolation in 1 Nc of the N c = 1 value in [9] and the N c = 2 value obtained here of the quenched condensates in units of 't Hooft coupling is consistent with the value in the 't Hooft limit.

B. N = 2
We use the four values of z i for N = 2 in Table I to obtain Σ i ( ) from Eq. (8). In Figure   4, we show the dependence of Σ i . Unlike the N = 0 theory, we do not see an asymptotic plateauing of the condensate even at the largest we were able to simulate. We have shown the fits of the data to Eq. (10) by the solid lines in the same plot for the different values of i.
From these extrapolations, we estimate the condensate at infinite volume for i = 1, 2, 3, 4 to be 0.0030(4), 0.0039(4), 0.0042(4), 0.0042(4) respectively. We find the condensates estimated from different eigenvalues to converge to significantly non-zero values at = ∞. As in the case of N = 0, these extrapolated values of Σ i from i = 2, 3 and 4 agree within errors. But, the extrapolated central value of Σ 1 is 30% lower than the others, but still significantly larger than zero. We think this is due to the difficulty in tuning the Wilson mass to obtain exactly massless fermions. One can rectify this in a future study with overlap fermions. Taking the average over the estimates of Σ from all the four low-lying eigenvalues, the condensate per color degree of For N = 4, we again use the four corresponding values of z i in Table I to obtain Σ i ( ) from Eq. (8). From the third panel of Figure 2 which shows the χ 2 /DOF as a function of γ m for N = 4, we find a wide range of values of γ m , including γ m = 2, that fits the data well. Therefore, we analyze the data first assuming the presence of non-vanishing condensate, in which case we force γ m = 2 and extrapolate the condensate to infinite volume using Eq. (10). This procedure is shown on a linear scale on the left panel of Figure 5. This leads to the values 0.0015(6), 0.0026(6), 0.0030(5), 0.0032(5) for Σ 1 , Σ 2 , Σ 3 and Σ 4 respectively in the infinite volume limit.
We stress that this analysis involves a prior assumption that γ m = 2, and with this assumption  we find a possible non-zero condensate. However, on the right panel of Figure 5, we again show Σ i ( ) as a function of 1/ , but on a log-log plot; a power-law behavior will be seen as a straight line on this scale with the slope being the scaling exponent. Now we can see a possible scaling behavior setting in for > 13. Assuming a perfect −2+γm behavior for > 13, we find γ m ≈ 0.60 if the theory is scale-invariant. Thus, the two different analyses leads to two different conclusions. A more careful analysis, using even larger values of than we have used, could place this theory on either side of the critical value N * .

D. N = 8 and 12
Having found an indication of transition from scale-broken to conformal phase at N ≈ 4, we studied N = 8 and 12 to see if we find strong evidence for a scale-invariant behavior. In

F. A flow from UV to IR
As explained in the beginning of this section, we used Eq. (9) to describe the data and also to find the value of γ m that best describes the asymptotic finite-size scaling behavior of the low-lying eigenvalues. Once we have fit the data using Eq. (9), we can define an dependent γ m as which will flow from γ (i) = 0 in the UV limit → 0, to γ (i) = γ m in the IR limit → ∞ for all i.
In Figure 8

IV. CONCLUSIONS
We have performed a numerical analysis of three dimensional SU (2) gauge theory coupled to an even number of massless fermions in such a way that parity is preserved. We used Wilson fermions instead of overlap fermions to reduce the computational cost. The price to pay was the absence of the full Sp(N ) flavor symmetry away from the continuum limit but this did not prevent us from extracting the critical number of fermion flavors.
We studied the finite volume behavior of the low lying Dirac spectrum in a periodic finite physical volume of size 3 using two different forms, namely, Eq. (9) and Eq. (10). The first one has the anomalous dimension, γ m , as a parameter that we varied to find the best fit. We found that the data clearly favored γ m = 2 for N = 2 and N = 4 and the data clearly favored γ We therefore conclude that the critical number of flavors is somewhere between N = 4 or N = 6 and we are not able to exclude N = 4 using the analysis presented in this paper. Having narrowed down the critical value of the number of flavors to one of two integers, the next step is to study the N = 4 and N = 6 theories using overlap [6] or domain-wall fermions [19,20].
Since the U (N ) flavor symmetry as well as the larger Sp(N ) global symmetry will be exact in the lattice theory, the behavior of all low lying eigenvalues can be used in the numerical analysis with equal confidence. Furthermore, one can also study the propagator of scalar and vector mesons and extract the behavior of their masses in finite volume.
The analysis presented in this paper clearly shows that one can use SU (2) gauge theories with massless fermions in three dimensions to numerically study the transition from scale invariant behavior to one that generates a scale using continuum finite volume analysis. It will also be interesting to test the predictions for symmetry breaking in [27] Acknowledgments All computations in this paper were made on the JLAB computing clusters under a class In our lattice simulation, both the gauge field as well as the N flavors of fermions are dynamical. That is, the Boltzmann weight e −S in our simulation using the Hybrid Monte Carlo (HMC) technique [28] uses S = S f + S g , where S f and S g are the fermion and gauge action respectively. The lattice variables are the SU(2) gauge-links; a gauge link U µ (n) is an SU (2) matrix that represents the parallel transporter from site n = (n 1 , n 2 , n 3 ) to site (n +μ), with 1 ≤ n i ≤ L. We impose periodic boundary condition for both the gauge field as well as the fermions in all three directions of the lattice; for SU(2) theory this boundary condition is sufficient since both 1 as well as −1 are part of the gauge group. We use the standard single plaquette gauge action, namely, Tr P µν (n), where P µν (n) is the parallel transporter around a plaquette in µν-plane at lattice site n: We reduced lattice spacing effect in fermion action as well as fermionic observables by smoothening the gauge field that enters the Dirac operator using the technique of gauge-link smearing.
For this, we used 1-level improved stout links [29], denoted by V µ (n): We used an optimum value s = 0.65 where the value of the smeared plaquette was maximized.
For the regularized two-component Dirac operator in Eq. (4), we used the two-component Wilson-Dirac operator C W given by where σ µ are the Pauli matrices, and M P is the Wilson mass which needs to be fine-tuned to non-zero values at finite lattice spacings in order study massless continuum fermions. The above two-component Wilson-Dirac operator satisfies C t W (n, m) = σ 2 τ 2 C W (n, m)σ 2 τ 2 , and hence there is a Sp(1) symmetry associated with a single flavor of two-component Wilson fermion even at finite lattice spacing. We incorporated the resulting fermion determinant, det C W C † W N/2 , from the N/2 parity-invariant pairs of two-component Wilson-Dirac fermions in our HMC simulation using using N/2 pseudo-fermion random vectors as explained in [5]. In the HMC, we were able to marginally optimize the molecular dynamics stepsize by using the Omelyan symplectic integrator [30] as well as by tuning the stepsize at run-time such that the acceptance is above 80%.  As discussed earlier, the eigenvalues of the four-component Hermitian Wilson-Dirac operator, D W , appear in ± pairs. The operator C W is not anti-Hermitian but becomes essentially one upon tuning M P to achieve massless fermions. As such, the Sp(1) × Sp(1) symmetry at finite lattice spacing becomes the full Sp(2) symmetry only in the continuum limit. But, the positive eigenvalues of D W can be used to study the presence or absence of a bilinear condensate as discussed in the context of U (N ) flavor symmetry [5]. In order to realize massless fermions, we tuned the value of M P at each simulation point to that value where the lowest eigenvalue of C W C † W is minimized when measured over a small ensemble of thermalized configurations at that simulation point. Since these tuned M P are required for any future computation at larger and L, and also for normalizing the eigenvalues of the overlap operator which makes use of the Wilson-Dirac kernel [6], we parametrize the tuned value of M P ( ) for different L and N and tabulate these parameters in Table II.
We used Ritz algorithm [31] to compute the four low-lying eigenvalues of C † W C W . We