Perturbative analysis of the Neuberger-Dirac operator in the Schr\"odinger functional

We investigate the spectrum of the free Neuberger-Dirac operator $\Dov$ on the Schr\"odinger functional (SF). We check that the lowest few eigen-values of the Hermitian operator $\Dov^{\dag}\Dov$ in unit of $L^{-2}$ converge to the continuum limit properly. We also perform a one-loop calculation of the SF coupling, and then check the universality and investigate lattice artifacts of the step scaling function. It turns out that the lattice artifacts for the Neuberger-Dirac operator are comparable in those of the clover action.


Introduction
Chiral symmetry plays an important role in the understanding of the strong interaction. A solution to a realization of the exact chiral symmetry on the lattice is proposed by Neuberger [1] as the Neuberger-Dirac operator (overlap). Recently, dynamical overlap lattice QCD simulations started in Ref. [2] and see [3] for an overview of recent progress. Nowadays, thanks to developments of algorithms and powerful current computers, large scale simulations are feasible as shown by the JLQCD collaboration in Ref. [4]. In that course, even after removing systematic errors 1 , finite size effects, cutoff effects and an ambiguity of chiral extrapolation, non-perturbative renormalization becomes an essential element for accurate quantitative predictions.
One elegant solution to this issue is the Schrödinger functional (SF) scheme [5] which is an intermediate scheme connecting the perturbative and hadronic regime. This method was shown to be useful to study the non-perturbative evolution over a wide range for various quantities, the coupling constant [6], the quark masses [7], the structure function [8,9], and the weak matrix elements [10]. By making use of the scaling technique, one can complete a perturbative matching safely at relatively high energy, and then the renormalization group invariant quantities, like the lambda parameter [11], the masses for the light quark [12] and the heavy quark by HQET [13], are determined without worrying about the systematic error of truncation of the perturbative expansion. Although so far the method was mainly used for Wilson type fermions [14], there are several attempts for other fermion formulations, like staggered fermions [15,16,17] and domain wall fermions [18]. Recently, formalisms for the Neuberger-Dirac operator on the SF have been proposed by Taniguchi [19] and Lüscher [20]. The former employs the orbifolding technique and the latter is based on universality considerations. In Ref. [21,22], the latter formulation is examined in the framework of the Gross-Neveu model. Here in this paper, we implement the formulation in QCD and study the spectrum of the free operator and calculate the SF coupling to one-loop order. Furthermore we investigate lattice artifacts of the step scaling function by comparing with the clover fermion action.
The rest of the paper is organized as follows. In Section 2.1, we summarize the definition of the Neuberger-Dirac operator on the SF which is given by Lüscher, and give some practical details about building the operator. Then we show the spectrum of the free Neuberger-Dirac operator in Section 3. In Section 4, we present results for the fermion part of the SF coupling to one-loop order by making use of the Neuberger-Dirac operator. Furthermore we investigate lattice artifacts of the step scaling function in Section 5. Finally, we give some concluding remarks and outlook in Section 6. In Appendices, we give an explicit form of the Wilson-Dirac operator on the SF in a time-momentum space for later use (A), and a discussion about a determination of a boundary coefficient at tree level (B), and we summarize some tables of numerical results (C). 2 Neuberger-Dirac operator on the SF

Definition
A massless Neuberger-Dirac operator on the SF (with size T × L 3 ) [20] is defined by withā = a/(1 + s). The operator follows the modified Ginsberg-Wilson (GW) relation where ∆ B is the exponentially local operator. The A in the kernel operator X of the inverse square root is given by where D w is the massless Wilson fermion on the SF [14]. The tunable parameter s is taken in a range −0.6 ≤ s ≤ 0.6 in the following. The boundary coefficient c represents the strength of the boundary operator P which is given by In Ref. [20], it is shown that the kernel operator X = A † A + caP is bounded from below by the spectral gap of A † A on the infinite lattice if c ≥ 1 holds, and furthermore it is mentioned that is the nearly optimal choice in order to achieve tree level O(a) improvement. We investigate this point in some detail in Appendix B, and we conclude that, for the precision of our calculation here (and maybe for future simulations), this formula is accurate enough. Therefore in the following calculations, apart from the appendix, we alway set c to the value in eq.(7). In this paper, we restrict ourself to the massless case of the Neuberger-Dirac operator.
To carry out perturbation calculations, it usually is convenient to move to momentum space. In the SF setup, however, the Fourier transformation can be done only in the spatial directions, due to a lack of translation invariance for the time direction. Therefore, we work in a time-momentum space, The explicit expression of the free D w in time-momentum space with the Abelian background gauge field in the SU(3) group [6] is shown in Appendix A. Here we show a matrix expression of A in time-momentum space for a fixed spatial momentum p and a color b(= 1, 2, 3), which is a totally 4(T /a − 1) × 4(T /a − 1) matrix. Block elements P ± = (1 ± γ 0 )/2 and g b (p, x 0 ) have Dirac spinor structure, and the latter is given by whereq andq are function of the spatial momentum p, the time x 0 and θ, and they are defined in Appendix A. The dependence of b is caused by the color diagonal background field. The angle θ comes from the generalized periodic boundary condition for the spatial directions, for k = 1, 2, 3. In time-momentum space, the boundary operator in eq.(6) is represented as and it's matrix expression is given by As it is clear from the explicit form, even in the free case we cannot have an analytic form of the Neuberger-Dirac operator on the SF, due to the presence of the background field. Therefore we have to rely on an approximation to the inverse square root, even for perturbative calculations. We will return to this issue how to build the operator in Section 2.3.

Distribution of ǫ
When one approximates the X −1/2 by a polynomial of X following Ref. [23], one needs information about the lower u and upper v bound of X. When the ratio u/v is not too small, one can obtain its approximation with lower degree. Usually, u and v are set to the values of minimal and maximal eigen-value (or norm) of the X, therefore it is important to know the spectrum of X. Here we discuss the spectrum of X at tree level in the presence of the background field.
At fixed momentum p and color b, we evaluate the minimal eigen-value and the norm of the kernel operator and set them as the lower and upper range of the approximation Essentially, ǫ b (p) = u b (p)/v b (p) controls the cost of the computation, since for a given precision it determines the degree (N ) of the approximation polynomial. Its distribution over p and color b in the case of θ = 0 with the non-zero background field are shown in Figure 1 for several lattice sizes L/a = T /a = 6, 12, 24 and s parameters s = 0.0, 0.5. As you can see in the figure, ǫ is distributed in a relatively large range 0.01 ǫ 1, therefore we concluded that it is better to calculate coefficients of the polynomial expansion for each p and b for saving time. The calculation of the coefficients is much cheaper than summing up the polynomial expansion. (You can use common lower and upper bound u comm = min p,b {λ min (X b (p))} and v comm = max p,b ||X b (p)|| for all momentum and color sector, but this is clearly inefficient.)  over the all spatial momenta p and the color b = 1, 2, 3 for the case of θ = 0 with the non-zero background field. We show here for lattice sizes L/a = T /a = 6, 12, 24 (from top to bottom) and s = 0.0 (Left panels) and s = 0.5(Right panels). Due to the permutation symmetry for the spatial momentum, we only have to investigate a momentum configurations p 3 ≤ p 2 ≤ p 1 . We did not correct the weight in the histogram, therefore the distribution shown here is not fully correct one. However the distribution is practically meaningful, since we do not do a computation on a full configurations with size (L/a) 3 , but only on the configurations p 3 ≤ p 2 ≤ p 1 (about (L/a) 3 /6).

How to build the Neuberger-Dirac operator
To build the inverse square root in the Neuberger-Dirac operator, we adopt the Chebyshev polynomial approximation, for the lower u and upper v bound of the X, and T k is the Chebyshev polynomial of degree k. Main tasks to obtain the operator are two-fold: the first is a computation of the coefficients c k and the second is to sum up in eq. (17). Concerning the latter part, we use the Clenshaw summation scheme to maintain precision. For the former part, we examine two methods to compute the c k , the Remez algorithm and the Chebyshev interpolation, in order to check rounding off errors occurring in the perturbative calculation in Section 4. Following Ref. [24], we implement the Remez algorithm to obtain minimax polynomials. On the other hand, we follow the Numerical Recipes for the Chebyshev interpolation. In calculations of the spectrum of the Neuberger-Dirac operator in Section 3, we only use the Chebyshev polynomial method. While for computations of the one-loop coefficient of the SF coupling in Section 4, we employ both methods and show the results to digits which are common in both methods.
When approximating the inverse square root of X, we demand a precision and we alway check consistency when it is built and we observe that errors on the right-hand side are less than 2 × 10 −13 .
3 Spectrum of the free Neuberger-Dirac operator 3.1 Spectrum of D N Since the U in eq.(2) is not a unitary matrix, there is no guarantee that its spectrum is distributed on a unit circle whose origin is (1, 0) as in the case of the infinite volume. However on the SF, as it is shown in Ref. [20], since one can see that Thus the spectrum ofāD N is contained in a unit disk which is enclosed by the unit circle.
We show the actual distribution of the spectrum in the free case in Figure 2 for s = 0.0 and s = 0.5. In the figures, we only show θ = 0 case, but θ = π/5 is also computed and has a similar tendency. In the figures, results for lattice size L/a = T /a = 6 and with the zero (BG= 0) or non-zero (BG= 1) background field are shown. Most eigen-values are localized near the unit circle, and a remnant of the GW relation is observed. Especially, note that p = 0 is most strongly affected by the boundary. When switching on the background gauge field (BG= 1), some degeneracies are lifted, and you can see 'more' points than at zero background field case (BG= 0). The upper panels are for s = 0.0, while the lower ones are for s = 0.5. The value of BG means that 0: zero background gauge field, 1: non-zero background gauge field (choice A). All eigen-values are enclosed by a black circle whose origin is (1, 0) and radius is 1, on which the spectrum of the GW relation operator lie. The blue circle points which come from the p = 0 sector are located far away from the circle, and they are positioned around a center of the circle. We observe that this sector is strongly affected by boundary effects.

Spectrum of
We also investigate the spectrum of the Hermitian operator Nā D N , because we can take the continuum limit and compare the scaling behavior with that of the Wilson-Dirac operator, L 2 D † w D w = (L/a) 2 aD † w aD w [25]. The numerical results of the lowest 10 eigen-values of L 2 D † N D N with nonzero background field are summarized in Table 2 for s = 0.0, 0.5, θ = 0, π/5 and L/a = T /a = 6, 12, 24. The scaling to the continuum limit are plotted in Figure 3 including those of the Wilson-Dirac and the clover action (c sw = 1) for comparison.
In the figure, the lower modes show a good scaling behavior, while higher modes are strongly affected by lattice artifacts.
4 SF coupling to one-loop order 4

.1 Definition and results
We compute the fermion part of the SF coupling [25] (we set L = T as usual) at one-loop order p 1,1 (L/a) for the massless Neuberger-Dirac operator. The one-loop coefficient 2 is given as with a normalization In the actual calculation, we expand the η derivative and use the fact that the determinant is factorized to the individual spatial momentum p and color sector b, Since D N is not the block tri-diagonal in the time and the spinor index, unfortunately we can not use the nice recurrence formula [25] which was used in the case of the Wilson-Dirac fermion. Therefore, we have to evaluate the inverse of D N directly by making use of a solver routine, and multiply with the η derivative of D N to take the trace. The trace in eq.(25), tr, concerns with the spinor and the time indices.
To compute the one-loop SF coupling, we need to take the η derivative of D N . This can be done analytically. To this end, we have to evaluate, where the dotted defines the derivative with respect to η. This summation can be evaluated by another recurrence relation besides the one needed to compose D N itself (the Clenshaw recurrence relation). The additional recurrence relation can be derived from the three terms recurrence formula for the Chebyshev polynomials and its η derivative formula.
2 Formally we should define the coupling in order to properly define the determinant in the continuum limit, but on the lattice the above form is equivalent to eq.(23) due to γ 5 Hermicity.   Table 3 for s parameters, s = 0.0, 0.5 and θ = 0, π/5. In order to estimate rounding off errors, we perform two methods of the approximation to the inverse square root as mentioned in Section 2.3, the minimax and the Chebyshev interpolation. In the table, we show nine significant digits where both approximations agree with each other. Even though we have used double precision arithmetic and been demanding 10 −13 precision for the inverse square root, we lose three to four digits in the summation step of all momentum and color sectors in eq.(25).

Coefficients of Symanzik's expansion
From the Symanzik's analysis of the cutoff dependence of Feynman diagrams on the lattice, one expects that the one-loop coefficient has an asymptotic expansion We can reliably extract first several coefficients by making use of the method in Ref. [26].
For the usual renormalization of the coupling constant, B 0 should be 2b 0,1 where b 0,1 is the fermion part of the one-loop coefficient of β-function for N f flavors QCD, We confirmed B 0 = 2b 0,1 = −0.00844343... to three or four significant digits for all cases (θ = 0, π/5 and −0.6 ≤ s ≤ 0.6). When the tree-level O(a) improvement is realized, we expect that B 1 = 0 holds. We check this to 10 −2 or 10 −4 in all cases. This shows that even though we have been using the approximate formula of the boundary coefficient in eq. (7), it works well to achieve the tree-level O(a) improvement to the precision here. In the following analysis we set exact values B 0 = −1/(12π 2 ) and B 1 = 0. A 0 gives an information about a ratio of Λ-parameters, and we show the obtained values in Table 1. By combining the previous results from Ref. [25,27], the values of A 0 can be obtained, and are shown in the second line (numbers with * ) in each s in Table 1 for θ = 0, π/5. We observe excellent agreements within errors for all s and θ parameters we investigated.
To achieve one-loop O(a) improvement, A 1 is needed to determine the coefficient of the fermion part of the boundary counterterm, c (1,1) t [25]. The resulting values are shown in Table 1. No θ dependence on the A 1 is observed beyond errors. The absolute value of |A 1 | = 0.02 − 0.03 of the Neuberger-Dirac operator is roughly factor two smaller than that of Wilson type fermion, |A 1 | = 0.038282(2) [25]. If one imposes an improvement condition [25], one finds that  Table 1: The coefficients of asymptotic expansion for −0.6 ≤ s ≤ 0.6 for θ = 0, π/5. The value of A 0 with * in the last line in each s parameter block are the values from the previous calculations [25,27]. The error for those values should be on the last digit.
For future reference, we provide interpolation formula of the c for −0.6 ≤ s ≤ 0.6.

Lattice artifacts of the step scaling function to one-loop order
In this section, we investigate lattice artifacts of the step scaling function (SSF) [28] σ(2, u), which describes the evolution of the running couplingḡ 2 (L) = u under changes of scale L by a factor 2, The lattice version of the step scaling function is denoted by Σ(2, u, a/L). Perturbative estimate of the lattice artifacts of the step scaling function can be studied by expanding a relative deviation The one-loop deviation, δ 1 (s, a/L), is decomposed into pure gauge and fermion part [25], δ 1 (a/L) = δ 1,0 (a/L) + N f δ 1,1 (a/L).
We are currently only interested in the fermion part. The fermion part of the oneloop deviation δ 1,1 (a/L) in terms of the one-loop coefficient of the SF coupling p 1,1 is given by Depending on the value of the boundary counter term c    Table 4 and plots in Figure 4, where we include those of the Wilson-Dirac and the clover action for comparison [25]. In the case of the clover action, c (1,1) t is set to be the proper value to achieve one-loop O(a) improvement, and for the Wilson fermion it is set that c (1,1) t = 0. We observe that the lattice artifacts for the Neuberger-Dirac operator are comparable to those of the clover action. As in the case of the clover action, the Neuberger-Dirac operator has less lattice artifacts for the case of θ = π/5 than θ = 0.

Conclusion and outlook
In this paper, we have explored the free Neuberger-Dirac operator on the SF. We investigated the spectrum of the operator, and then we confirmed that the spectrum of D N is enclosed by the unit circle and the spectrum of D † N D N has the expected scaling behavior (1/L 2 ) and the correct continuum limit. We also performed the oneloop computation of the SF coupling by making use of the operator. We confirmed the universality, and the fermion part of the O(a) boundary counterterm at oneloop order, c   artifacts for the Neuberger-Dirac operator is almost the same as that of the clover action. Thus, we may expect small lattice artifacts for the non-perturbative SSF of the Neuberger-Dirac operator, as in the case of the clover action [11]. In Appendix B, we demonstrate that the choice of the boundary coefficient c = 1 + s given in [20] at tree level is almost optimal. This formula may be precise enough for actual simulations.
We exclusively considered the massless case. Although the massive case can be explored, we leave it as a future task. Comparison of scaling behavior with the other formalism [19] and a consistency check is also interesting. Before starting non-perturbative computations, we have to compute some improvement coefficients to one or two loop order. Furthermore perturbative calculations of the renormalization factors are in a to-do list. By combining techniques in Ref. [29], a two-loop calculation of the SF coupling including the Neuberger-Dirac operator as a fermion part may also be feasible.
Apart from the one-loop computations, next target would be a computation of the renormalization constant of the flavor singlet scalar density Z S non-perturbatively, since the bare quark condensate in two flavor QCD was already computed by JLQCD [4] 3 . Due to the chiral symmetry, Z S is identical to the renormalization constant for the non-singlet flavor pseudo scalar density, Z P = Z S [30,31]. Actually, the non-perturbative renormalization group running of Z P is already known in Ref. [12] for the SF scheme. A missing piece to obtain the renormalization group invariant quark condensate is a low energy matching factor, Z P (g 0 , µ = 1/L max ) in the SF scheme for the overlap fermion. This is an urgent and possibly doable task in the near future.

A Free Wilson-Dirac operator on the SF in time-momentum representation with non-zero background field
In the presence of the background gauge field [6], the free part of the Wilson fermion action in time-momentum space has a form [25] with the boundary conditions The massless part of the Wilson-Dirac operator on the SF is given as where indices b, c refer to color and h b (p; x 0 ) is given by and The spatial component of the momentum p is given by for n k = 0, · · · , L/a − 1. The boundary phases φ b and φ ′ b are given by Switching off the phases φ b = φ ′ b = 0 correspond to the zero background field. One can consider the effect of the non-zero background field a shift for the spatial momentum. The θ also affect as a constant shift for the momentum in a global manner, on the other hands, the background field provides a time dependent (local) shift.

B Determination of boundary coefficient c at the tree level
The boundary coefficient c in the kernel of the overlap operator on the SF is expanded in terms of the coupling constant g 2 0 , In this appendix, we determine the tree coefficient c (0) which depends on s parameter. Ref. [20] gives the formula in eq.(7), and we will examine it carefully here. We consider the SF with T = 2L and θ = 0 in the presence of the non-zero background field. The massless Neuberger-Dirac operator is assumed also in this appendix. The basic correlation functions [32] we use are given by where boundary fields [20] are given by At the tree level, f A (x 0 ) f P (x 0 ) are given as with Γ = γ 0 γ 5 , γ 5 for f   We extract the order a coefficient A 1 from the Symanzik's expansion for the ratio f (0) We estimate the error of A 1 by making use of the method in [26]. We determine c * (0) such that A 1 (c * (0) ) = 0 for the range −0.6 ≤ s ≤ 0.6 (improvement condition). In Figure 5, we plot c * (0) (s) as a function of s. By fitting the data points with a functional form c * (0) (s) = 1 + k 1 s + k 2 s 2 + k 3 s 3 + k 4 s 4 + k 5 s 5 , This curve is also shown in Figure 5. For larger s, a discrepancy between the above formula and eq.(7) can be seen, and their difference is maximally 10% in the range −0.6 ≤ s ≤ 0.6. As a consistency check, by making use of the value of c in eq.(60,61), we compute the ratio in eq.(59) with θ = π/5, and then A 1 = 0 is confirmed up to 10 −4 in the range of s.  [20] is also shown as solid red. The error bar of the points are too small to see in this scale.
In order to measure an effect of the difference on a physical quantity, we compare the one-loop coefficient p 1,1 with different values of c from the different formulae of c (eq.(7) and eq.(60)) at s = 0.5. It turns out that a difference in the p 1,1 is less than one percent on the lattice size L/a = 4, ..., 48. Furthermore, the resulting Symanzik's coefficients of p 1,1 in eq. (27) do not change within errors. Therefore, we concluded that, to the precision in our calculation, the formula c = 1 + s is accurate enough to achieve the tree level O(a) improvement.      Table 3: The one loop coefficient p 1,1 (L/a) for s = 0.0, 0.5 with θ = 0, π/5. The last digits may be affected by rounding off errors.  for s = −0.5, 0.0, 0.5 for θ = 0(upper) and θ = π/5(lower). Results for the Wilson and the clover action [25] are also included for comparison.