Numerical determination of monopole scaling dimension in parity-invariant three-dimensional non-compact QED

We present a direct Monte-Carlo determination of the scaling dimension of a topological defect operator in the infrared fixed point of a three-dimensional interacting quantum field theory. For this, we compute the free energy to introduce the background gauge field of the $Q=1$ monopole-antimonopole pair in three-dimensional non-compact QED with $N=2,4$ and $12$ flavors of massless two-component fermions, and study its asymptotic logarithmic dependence on the monopole-antimonopole separation. We estimate the scaling dimension in the $N=12$ case to be consistent with the large-$N$ (free fermion) value. We find the deviations from this large-$N$ value for $N=2$ and $4$ are positive but small, implying that the higher order corrections in the large-$N$ expansion become mildly important for $N=2,4$.


I. INTRODUCTION
Conformal field theories in three dimensions, and renormalization group flows from one fixed point to another induced by the introduction of relevant operators at fixed points have been investigated over the last few years. This involves the computation of scaling dimensions, Δ, of operators at the different fixed points. The operators O at a fixed point could be the usual composites of the field variables, and hence trivially local and amenable to the standard Monte Carlo computations of two-point functions of the local operator The operators could also be topological disorder operators [1] which act as sources to topological conserved currents in the theory. Since such operators cannot be written as simple composites of the field variables, studying their scaling dimensions is a challenge, especially on the numerical side. In the case of theories with U(1) global or gauged symmetry, the topological defects are the monopoles, M Q , which create Q units of flux surrounding it [1], and hence serve as the sources of the otherwise trivially conserved U top ð1Þ current, j top μ ¼ ϵ μνρ F νρ =ð4πÞ. Three-dimensional QED, whose gauge group is U(1) as opposed to R, is one such theory where monopole defects can occur. Depending on whether monopoles are energetically allowed or disallowed in the continuum limit, the three-dimensional QED is classified as compact or noncompact, respectively. The presence of two distinct theories, differing simply by the presence or absence of monopoles, offers theoretical and computational possibilities in understanding the emergence of mass gap.
Pure-gauge compact QED in three dimensions is a rare example in which the emergence of mass gap could be understood through the dual superconductor mechanism where in the electric charges experience a linear confining potential due to the presence of plasma of magnetic monopoles [2,3]. Coupling the compact QED 3 (referred to as c-QED 3 ) to many flavors, N, of massless twocomponent fermions (assumed to be even to preserve parity) gives a possibility to counter the emergent mass gap [4]-above certain critical flavor N C c , the theory is expected to be conformal in the infrared, whereas it develops a mass gap below N C c . In a first exploratory study [5] towards finding N C c , a derivative of free energy required to introduce a single monopole was computed on a lattice with open boundary conditions in the N ¼ 8 flavor compact QED 3 , and no convincing evidence for an infinite free energy signaling monopole confinement was found in the continuum limit. Thereby, this suggested a presence of monopole plasma and a consequent mass gap in compact QED 3 for N ≤ 8. Understanding such infrared quantum phases obtained by tuning parameters of the underlying quantum field theory (QFT) is an ongoing field of research (cf. [6]). Similar studies of the critical number of flavors, N NC c , in noncompact QED 3 (referred to as nc-QED 3 ) have continued to be investigated through ab initio lattice simulations [7][8][9][10][11] as well as through other approximation methods [12][13][14][15][16][17]. Unlike noncompact QED 3 , the presence of monopoles in the compact version, even as the continuum limit is approached, is a technical challenge to numerical studies due to the presence of many small eigenvalues of the three-dimension Dirac operator [18]. An indirect feasible approach is to check whether the monopole operator is marginally relevant in the infrared fixed point of the N flavor noncompact QED 3 [4]. Crucial to this inference is that the monopoles in a gauge theory with N massless fermions break UðNÞ global flavor symmetry to UðN=2Þ × UðN=2Þ symmetry [4,19,20]. Such an approach further assumes that (1) both compact and noncompact QED 3 flow to the same infrared fixed point for N > N C c ; (2) N NC c < N C c . At least in the N → ∞ limit, the compact or noncompact action will be subdominant compared to the induced gauge action from the fermion, and hence, the infrared physics should be the same for both nc-and c-QED 3 . The stronger assumption is that this continues to remain so until N ¼ N C c . The second assumption is based more on numerical works [10,11] that strongly indicate that N NC c < 2. This also means that only the dressed, gauge-invariant monopole operators become relevant at N ¼ N C c and other UðNÞ symmetry breaking operators, such as the four-Fermi operators, remain irrelevant [21]. Therefore, a computation of scaling dimensions in nc-QED 3 and a subsequent direct confirmation of N C c in c-QED 3 is well motivated. Monopole operators also play a similar role to understand quantum phase transitions in lattice systems with gauged U(1) symmetry which was recently analyzed computationally in compact QED 3 [22,23], and in the QED 3 -Gross-Neveu model [24]. In [25], analytical progress was made on monopoles in such lattice systems.
A practical method to determine the monopole scaling dimension Δ Q analytically is by coupling the theory with the U(1) symmetry to the classical, scale-and rotationally invariant Dirac monopole background A Q and study the response of the theory. Analytically, one computes the Casimir energy of the theory defined on S 2 with uniform 2πQ flux over it, which by state-operator correspondence is the same as the scaling dimension Δ Q [26][27][28]. Such computations are usually perturbatively done order by order in 1=N (cf. [29]), and currently it is only up to Oð1=NÞ. Nonperturbative conformal bootstrap has also been applied to QED 3 to find the allowed region in the parameter space of scaling dimensions of Q ¼ 1 and 2 monopoles [30]. Complementary to such bootstrap computations, it was demonstrated [31] that a direct way to compute monopole scaling dimensions using lattice computation is to couple such theories to a background field A QQ ðx; τÞ ¼ A Q ðx; x 0 Þ − A Q ðx; x 0 þtτÞ that gives rise to a monopole at x 0 and an antimonopole at x 0 þtτ, which are separated by a distance τ and compute the scaling of the partition function as τ → ∞. It is the aim of this paper to apply this method and compute Δ Q for a Q ¼ 1 monopole in the infrared fixed points in N flavor noncompact QED 3 . In particular, we compute the finite N corrections to the large-N scaling dimension for small enough values of N where a nonperturbative computation becomes inevitable.

II. c-QED 3 , nc-QED 3 AND MONOPOLE CORRELATOR IN nc-QED 3
In this section, we consider different versions of QED 3 that one could construct on the lattice. We consider L 3 Euclidean lattices whose physical volume is l 3 , with the lattice spacing being l=L. Let θ μ ðxÞ ∈ R be the lattice gauge fields which are related to the physical gauge fields θ μ ðxÞ ¼ A μ ðxÞl=L. The notation is such that x,y denote integer valued lattice coordinates. The two-component Dirac fermions in all the cases to be considered are coupled to compact gauge links, U μ ðxÞ ¼ e iθ μ ðxÞ , through an UV regulated massless Dirac operator CðUÞ. In this work, CðUÞ is the 1-HYP smeared Wilson-Sheikhoslami-Wohlert Dirac operator C W with the Wilson mass m w tuned to the massless point [10]. In the parity-invariant QED 3 with an even number of flavors, N=2 of two-component fermions are coupled via CðUÞ and the other N=2 via C † ðUÞ. The partition function for QED 3 can be written in general as where W g is the Boltzmann weight from the pure-gauge part. Since the fermionic determinant is invariant under θ μ ðxÞ → θ μ ðxÞ þ 2πn μ ðxÞ for integer values n μ ðxÞ, this part of the action respects the compactness of the U(1) gauge group. Independent of the choice of W g , we can always restrict the above integral over all θ μ ðxÞ to be from −π to π by simply summing W g over all possible n μ ðxÞ for different x and μ. In this way, the underlying gauge group is always U(1) owing to the usage of the compact links U μ ðxÞ in the Dirac operator, and hence magnetic monopoles are well defined in these theories. Depending on the form of W g , one can study QED 3 with or without monopoles as we elaborate on below, and as also discussed in [32].
All gauge actions will be functions of the fluxes on plaquettes where the flux on the plaquette in the ðμ; νÞ plane with one corner at x is F μν ðxÞ ¼ ∇ μ θ ν ðxÞ − ∇ ν θ μ ðxÞ; ∇ μ fðxÞ ¼ fðx þμÞ − fðxÞ: The Boltzmann weight for the noncompact lattice gauge action, does not favor the presence of monopoles in the continuum limit since the flux on each plaquette is peaked around zero when one takes L → ∞ at a fixed l. The compact Wilson gauge action, on the other hand, does not suppress monopoles in the continuum limit since the flux F μν ðxÞ has multiple peaks around 2πN μν ðxÞ with integer values of N μν -monopoles are counted per cube [33] by writing F μν ðxÞ ¼F μν ðxÞ þ 2πN μν ðxÞ whereF μν ðxÞ ∈ ½−π; πÞ and N μν ðxÞ are integers. The monopole charge inside a cube with one corner at x is given by The Villain gauge action [34,35], is also a compact action but has the advantage that the integer part of the flux per plaquette is made explicit. We have introduced new degrees of freedom N μν ðxÞ and one needs to sum over all integer values to define the partition function. This action is expected to be in the same universality class as the compact Wilson gauge action. The Villain action allows for all values of QðxÞ with the only condition that the sum over all x in a finite lattice with periodic boundary conditions will be zero. The only coupling in all cases is l which can be viewed as the dimensionless extent of the lattice and the lattice spacing is a ¼ l L . One can only consider the part of the above Villain action restricted to the sector QðxÞ ¼ 0 for all x. If the manifold is R 3 , then this automatically implies that for integers n μ . On T 3 , as used in Monte Carlo simulations, the condition in Eq. (9) implies QðxÞ ¼ 0 but further restricts the sum of N μν on any ðμνÞ plane to be zero. In particular, this disallows configurations with net constant flux 2πQ, for integer Q, over any of the ðμνÞ plane in the continuum limit. However, such an extra restriction on T 3 cannot be important in the thermodynamic limit since any equal and opposite fluctuations in flux in different parts of the lattice are allowed. For values of N μν of the form in Eq. (9), one can change θ μ ðxÞ → θ μ ðxÞ − 2πn μ ðxÞ and annul the term N μν . Therefore, the Villain path integral restricted to values of N μν of the type in Eq. (9) is the same as the standard noncompact QED 3 path integral defined using Eq. (5). Similarly, one can constrain the integer valued flux N μν to take a particular value N QQ μν defined via The above constraint corresponds to an insertion of flux Q monopole at a lattice site y and a flux Q antimonopole at y 0 , and this cannot be absorbed by a change of variable of the gauge fields. The monopole correlator in nc-QED 3 can simply be defined as the ratio of path integrals subject to the constraint in Eq.
Instead, we find the gauge field background A QQ μ ðxÞ that minimizes and couple the theory to this classical background field in order to define The advantage of using B QQ μν over using 2πN QQ μν is that background field coupling has no effect in pure-gauge theory, and any effect that is observed in Z Q will arise only due to the presence of fermions. This follows from a simple change of variable θ μ ðxÞ → θ μ ðxÞ − A QQ μ ðxÞ that eliminates A QQ μ ðxÞ only in the case of a pure-gauge path integral. As we already noted, N QQ μν ðxÞ cannot be written as a curl, and hence such a change of variable is not possible even in pure-gauge theory. On R 3 , the resulting A QQ is the field for a Dirac monopole-antimonopole pair. The advantage of minimizing Eq. (11) on toroidal lattice is to take care of both the lattice discretization as well as the periodicity correctly. We checked through a full fledged computation in the case of N ¼ 2 QED 3 that the difference in Z Q =Z 0 between the minimum on the torus as defined above and the discretized field of a Dirac monopole-antimonopole pair as defined in [31] is, however, marginal.
Let us denote the lattice distance between the monopole and antimonopole as T ¼ jy − y 0 j, which is related to the physical separation τ ¼ Ta. Then, the "bare" monopoleantimonopole correlation function in lattice units, G ðQÞ B , is the ratio of partition functions with and without the flux Q monopole-antimonopole insertion [36]: Our specific choice for the location of the monopole and antimonopole in Eq. (10) is realized by and zero for all other directions and lattice points ðx 1 ; x 2 ; x 3 Þ. The square tube with nonzero integer flux running between the monopole at y ¼ ð0; 0; 0Þ and y 0 ¼ ð0; 0; TÞ is the Dirac string. Any other configuration for this Dirac string that is simply connected to the above construction is related through appropriately chosen transformations θðxÞ → θðxÞ þ 2πnðxÞ. The details pertaining to the construction of the background field A QQ μ can be found in [37].

III. METHOD AND SIMULATION DETAILS
In a Monte Carlo simulation, it is only possible to compute ensemble averages and not the partition function itself. A brute force way to implement the correlator in Eq. (13) is to compute the average where hÁ Á Ái 0 is the ensemble average with respect to Z 0 for N flavor theory. The problem with such an approach is the absence of overlap between the configurations sampled by Z 0 and Z Q . In order to avoid this overlap problem, we couple QED 3 to the background gauge field ζA 11 through a generalization of Eq. (12) to noninteger values of Q, where ζ is a tunable auxiliary variable [38]. Consistent with the previously introduced notation, the resulting partition function is Z ζ . From this, we can compute the lattice free energy F ðQÞ B ðτ; l; aÞ to introduce the monopole-antimonopole pair separated by physical distance τ in an l 3 torus at finite lattice spacing a as where Thus, WðζÞ can be computed in the Monte Carlo simulation of Z ζ through the measurement of F μν ðxÞ − ζB ð1Þ μν ðxÞ on the gauge fields that are sampled. In this paper, we will only study Q ¼ 1 monopoles and we drop labels for Q henceforth.
A way to determine the correlator in Eq. (13) is to compute G B ðτ; l; aÞ at different large values of τ in an l 3 box. At each fixed τ, one should first convert the lattice correlator to a renormalized physical one, then take the continuum limit L → ∞ at a fixed l, followed by the infinite volume limit l → ∞. Finally, one can consider the asymptotic τ → ∞ limit to study its τ −2Δ scaling. However, such a method is not practical since it requires computations of multiple values of τ per Monte Carlo sample point in the parameter space, and further introduces unwanted systematic errors from the l → ∞ extrapolations at fixed τ. As was demonstrated in the case of monopole correlators [31], a better method is to make use of the scaling of correlators near the infrared fixed point. That is, one expects the scaling G B ðτ; l; aÞ ¼ a 2d G R ðτ; lÞ; The conversion factor a 2d takes the bare correlator to the renormalized correlator of the naive dimension d monopole operator. 1 The subtle issues with this will be addressed in the next section. In addition, the first expression is only true up to finite a, or equivalently finite 1=L, corrections. Assuming we have obtained the renormalized correlator, the second expression exhibits its scaling near the infrared fixed point. We do not have to make any further assumption about the form of Gðτ=lÞ if we fix τ=l ¼ ρ as l is varied. Here, we take ρ ¼ 1=4. Equivalently, the free energy to introduce a monopole-antimonopole pair separated by the distance τ ¼ ρl would be up to higher-order corrections in 1=l. Since we keep ρ fixed in this paper, we keep its dependence implicit. It will 1 The symbol d should not be confused with the Euclidean space-time dimension, which is always 3 in this paper. be useful to consider the free energy per two-component flavor as In the limit of N → ∞, both f 0 0 ðρÞ and Δ=N have welldefined limits. In 1=N expansion, one finds with k < 0. The large-N value Δ ∞ was computed using the free fermion coupled to the monopole background since it was argued that the fluctuations in dynamical gauge fields are suppressed by 1= ffiffiffiffi N p . Such an analysis gave Δ ∞ ¼ 0.265 [1,28,31]. For the Q ¼ 1 monopole we consider here, the leading correction was computed to be k ¼ −0.0383 [4].
In the current work, we studied N ¼ 2, 4 and 12 flavors of fermions-the idea being that we can use N ¼ 12 to check for consistency with large-N expectations, and use N ¼ 2, 4 to study the effect of smaller N. We sampled configurations from Z ζ using 50 000 trajectories of hybrid Monte Carlo simulation. For each value of l, L and N, we simulated 24 different equally spaced values of ζ from 0 to 1. At each ζ, we computed WðζÞ using jackknife analysis to take care of autocorrelation, and performed the numerical integration in Eq. (16) after smoothly interpolating the 24 data points for WðζÞ. We used different values of l ranging from l ¼ 1 to l ¼ 250 at each fixed value of L. To estimate the continuum limit of the l dependence of the free energy, we used L 3 lattices with L ¼ 16, 20, 24 and 28. In Fig. 1, we show WðζÞ as determined for N ¼ 2 at four different values of l on 20 3 lattice as a sample. The area under each of those curves gave the bare free energy F B ¼ − log G B .

IV. RESULTS
First, we show the dependence of the lattice free energy per flavor, f B ðl; LÞ ¼ N −1 F B ðρl; l; l=LÞ, for introducing an monopole-antimonopole pair at a distance τ ¼ l=4 from each other on the box size l in Fig. 2 for different fixed values of L. The plots from left to right are for N ¼ 2, 4 and 12 flavors of two-component fermions, respectively. This dependence as computed using L ¼ 16, 20, 24 and 28 are shown as different colored symbols. As expected, the bare lattice free energy from different L do not fall on a universal curve since the lattice spacing a ¼ l=L keeps changing as l is varied at fixed L. In fact, as it stands the result seems unphysical-the free energy decreases with increasing l at fixed L. Therefore, we have to first convert the lattice correlator G B to the correlator in physical units, G R by determining d in Eq. (18). Since QED 3 is superrenormalizable, d for a local operator would simply be its naive dimension (e.g., the flavor triplet vector operator O V In other words, the background field method singles out the scaling operator of lowest naive dimension only in the large-L limit at the Gaussian FP. (B) The effective action for the background field is apparently nonlocal, and numerically showing that it can be renormalized by a simple a 2d factor is nontrivial. As discussed above, for the background field coupling at finite L at l ¼ 0 (or a ¼ 0), one can only obtain an Ldependent effective value of the scaling dimension, dðLÞ, which we expect to approach the free fermion value d in the limit of very large L that are not feasible in the computation presented here. The value of dðLÞ relevant in the range of L studied here can be determined numerically via an L −2dðLÞ fit to G B ðρLa; La; aÞ in the limit of a → 0 over a small range of L. With this L-dependent value of d determined at the Gaussian fixed point, we can best approximate the renormalized correlator at other nonzero l and a by using which automatically ensures that the correlator G R ðρl; lÞ has no L dependence for l ≈ 0. For larger l, the residual L dependence in G R ðρl; lÞ is a lattice artifact, which can be removed by L → ∞ continuum extrapolations. To actually study QED 3 in the strict l → 0 limit at finite L, one needs to integrate the fermion determinant over the still unsuppressed constant modes of the gauge fields in all three directions of the torus. This is nontrivial to implement, and hence we consider the result from the fixed, small a ¼ 1 7 as an approximation of the strict l ¼ 0 results. In Fig. 3, we have shown the logðLÞ dependence of the free energy F B ðρLa; La; aÞ, simply denoted as F B ðLÞ, at this fixed a ¼ 1 7 over the range of L used in this paper. This corresponds to changing l ¼ 112 at L ¼ 16 to l ¼ 196 at L ¼ 28 in each of the panels in Fig. 2. It was possible to fit the data over the range of L from 16 to 28 using f 0 þ 2d logðLÞ and thereby obtain the value of dðLÞ over this range of L. The fits are shown as the black straight lines in the three panels. We obtained the slope 2dðLÞ as 2.34(22), 3.81(21) and 11.17(48) for N ¼ 2, 4 and 12, respectively, which correspond to dðLÞ=N of 0.585(56), 0.476(53) and 0.465 (20), respectively. Since the value of G B at l ¼ 0 depends on the distribution of constant gauge fields in the three directions that are allowed at l ¼ 0 for any finite N, the value of dðLÞ=N for intermediate L can depend on N. Instead of the above method, where l is varied at fixed L thereby forcing us to construct G R from G B , we could have instead studied G B at fixed lattice spacing. However, achieving larger physical volumes at fixed small lattice spacing would become numerically prohibitive.
In Fig. 4, we show the resulting free energy per twocomponent flavor, f R ðlÞ ¼ N −1 F R ðρl; lÞ. This was obtained by adding 2dðLÞ logð l L Þ to F B ðρl; l; l=LÞ and then computing the resulting renormalized free energy per two-component fermion. The three panels are for the three different values of N. The data from different L, made distinct by the colored symbols, now fall on near universal curves. This data collapse is quite nontrivial and supports the assumption that we have defined the correlator of an operator that has a local description in the continuum. Contrary to the behavior of the bare lattice free energy, the renormalized free energy starts increasing with l as physically expected since one does not expect monopoles to be spontaneously created in noncompact QED 3 . For all N, including N ¼ 2, the dependence of f R ðlÞ shows no evidence of a linear l-dependent piece corresponding to an exponential fall, G R ðlÞ ∼ expð−μρlÞ, with a mass μ that could set the scale for a scale-broken theory. Assuming QED 3 with N ¼ 2, 4, 12 flow to infrared fixed points as l → ∞, the asymptotic values of slope in this linear-log plot would give the values of 2ΔðNÞ. As can be seen, the slope changes with l and extracting the value of Δ will require extrapolations. Therefore, first we focus on modelindependent inferences from the data.
From the three panels in Fig. 4, where we have kept the range of l and f R in the plots to be the same, we find the free energy per flavor f R ðlÞ for N ¼ 12 shows a weaker dependence on l compared to N ¼ 2, 4. To make this quantitative, in Fig. 5, we show the difference between the free energy per flavor in N and N 0 flavor theories. 2 In the infrared, we expect such a difference to be as l → ∞. For N 0 , we choose the largest value, N 0 ¼ 12, that we have. We have shown δðl; 2; 12Þ in the left panel, and the difference δðl; 4; 12Þ in the right panel as functions of logðlÞ. At leading order in 1=N, this difference vanishes.
Instead, we find that both δðl; 2; 12Þ and δðl; 4; 12Þ increase with l making it quite evident that ΔðN ¼ 2Þ=2 and ΔðN ¼ 4Þ=4 are larger than ΔðN ¼ 12Þ=12. This effect is arising purely due to the finite value of N. Quite surprisingly, δðl; N; 12Þ shows a logarithmic dependence on l over the entire range of l used. Perhaps this is due to the finite l corrections to the infrared scaling get approximately canceled between f R at N and N 0 . Therefore, we performed a combined fit to the data for δðl; N; 12Þ at all l from different L using an ansatz, with a 0 , a 1 , b 0 and b 1 as fit parameters. The value of b 0 will then give us an estimate of the difference ΔðNÞ=N − Δð12Þ=12 in the continuum limit L → ∞. We find such a fit ansatz to describe the data for both N ¼ 2 and 4 well with χ 2 =dof < 2 with 81 data points in each fit. We find the best fit parameters ½a 0 ; a 1 ; b 0 ; b 1 for N ¼ In the expansion in 1=N around the large-N fixed point, the leading contribution to the difference is ΔðNÞ 0383. Thus, to order 1=N, the value of the scaling dimension decreases at finite N from the large-N value, and the slope of δðl; N; N 0 Þ as a function of logðlÞ should be negative at this order. Our numerical result for N ¼ 2 and 4, on the other hand, suggests the opposite behavior for the values of N which are Oð1Þ. This implies that higher-order terms in 1=N that are of opposite sign cannot be ignored, or perhaps a breakdown of the 1=N expansion. However, it is true that the corrections themselves are also small. Now, it remains to be shown that for N ¼ 12, its scaling dimension Δð12Þ is consistent with the large-N expectation, as one might naively expect the large-N expansion to hold here. Unlike the above conclusion about the correction to large-N behavior, extraction of Δð12Þ requires modeling and extrapolations since the free energy does not exhibit a pure logðlÞ dependence over the entire range of l used in this computation. In Fig. 6, we focus only on values of l > 64. By fitting a simple ansatz to the N ¼ 12 data for f R ðlÞ over the larger values of l from different L, we estimate the value of Δð12Þ=12 as the best fit value of β 0 . The best fit logðlÞ dependence for L ¼ 16, 20, 24 and 28 are shown along with the data in Fig. 6, which are described by ½α 0 ;α 1 ;β 0 ;β 1 ¼½−0.42ð18Þ;−19ð4Þ;0.261ð19Þ; 2.22ð41Þ. Though the data seem to be well described by such an ansatz, the χ 2 =dof is about 3 due to the much smaller errors in the data for N ¼ 12. The black dashed line is the estimated logðlÞ dependence in the L → ∞ limit. We estimate the slope of this continuum dependence as β 0 ¼ 0.26ð2Þ. For comparison, the value of Δð12Þ=12 from large N up to leading order in 1=N is 0.262. Our estimated value of Δð12Þ is consistent with this value, and thereby lends further support for our numerical work. This implies that the monopole scaling dimension for N ¼ 12 is estimated to be  3.24 (24). This is consistent with N ¼ 12 being the critical flavor where the Q ¼ 1 monopole operator becomes just marginally relevant in the infrared fixed point.

V. CONCLUSION AND DISCUSSION
In this paper, we presented an ab initio lattice computation of the monopole correlator in N ¼ 2, 4 and 12 flavor massless QED 3 by using the background field method. To avoid the overlap problem which would make the computation of ratio of partition functions with and without a monopole-antimonopole background field, we slowly increased the value of monopole flux from 0 to integer Q. One of the noteworthy results in this paper is the feasibility of this approach itself, seen via the good signal to noise ratio of the monopole free energy (which is the negative logarithm of the monopole correlator). This encourages the application of this method to other QFTs where monopole operators can be defined. We demonstrated empirically that the monopole correlator behaves like a local operator and can be simply "renormalized" by the factor a 2d at lattice spacing a, where d is the naive monopole dimension as obtained on L 3 lattice in the limit a → 0. The key numerical result for the free energy to introduce a monopole-antimonopole pair in N ¼ 2, 4 and 12 flavor QED 3 is shown in Fig. 4. Using these data, we demonstrated that the scaling dimension for N ¼ 12 QED 3 is consistent with the large-N expectation and that N ¼ 12 is consistent with being the critical flavor where the monopole scaling dimension takes the marginally relevant value of 3. By computing the differences in monopole free energy between N ¼ 2 and 12, and between N ¼ 4 and 12, we found evidence for the deviation of the scaling dimension in N ¼ 2, 4 theories from the N ¼ 12 theory to be positive. This is in contradiction to the 1=N analysis up to Oð1=NÞ. It remains to be seen if this tension can be resolved by inclusion of higher-order corrections in 1=N in the analytical expressions, or points to a breakdown of the 1=N framework itself where in the fixed point for smaller N belongs to a different family than the one in large N. In conclusion, the results in the paper along with slightly different analytical results from large-N analysis support the direct computation of N flavor compact QED 3 around N ≈ 12, which however requires algorithmic development to deal with large number of near-zero Wilson-Dirac modes.
In the paper, we did not demonstrate in the lattice regularization framework that monopoles carry flavor quantum number and break the UðNÞ flavor symmetry to UðN=2Þ × UðN=2Þ symmetry. Demonstrating this is not important to the computation presented in this paper, but central to the UðNÞ flavor symmetry breaking in N c flavor compact QED 3 . In the continuum, one shows [1,19] this by noting that the ground state of the Hamiltonian of massless fermion on S 2 with constant flux 2πQ has QN zero modes. Thus, the gauge-invariant, CP-invariant vacua are obtained by filling the QN zero modes with QN=2 fermions picked amongst N=2 flavors and QN=2 antiparticles picked amongst another N=2 flavors. Thus, the vacua with monopole background transform under the irreducible representations of the UðNÞ flavor symmetry. With the lattice regularized fermions on a spherical monopole background, we do not have a similar derivation to study the flavor structure of the vacuum. A difficulty is defining the lattice fermion on S 2 . Therefore, to gain an understanding of a similar flavor symmetry breaking mechanism on the lattice, we consider a background field on T 3 corresponding to constant flux 2πQ on all T 2 spatial slices, as an analogue of constant flux background on S 2 × R for the spherical monopole. It was shown [39] that the two-dimensional transfer matrix for the two-component Wilson-Dirac operator, C W , has 2L 3 þ Q eigenvalues greater than 1 and 2L 3 − Q eigenvalues less than 1 for m w > 0; a consequence of gauge field topology in two dimensions. Thus, the vacuum has total Q fermions. Similarly, for C † W with m w < 0, the vacuum has a total of Q antifermions. With N=2 fermions coupled via C W and another N=2 fermions coupled via C † W , the vacuum has N=2 fermions from N=2 flavors and N=2 antifermions from the other N=2 flavors. This is very similar to the structure seen in a spherical monopole background. The flavor symmetry breaking on this particular background is determined by the need to preserve parity: a choice of N=2 flavors coupling to C W and the other ones to C † W . We expect a similar mechanism to be true for a massless Wilson fermion on a monopole background.