Schr\"odinger functional boundary conditions and improvement for N>3

The standard method to calculate non-perturbatively the evolution of the running coupling of a SU(N) gauge theory is based on the Schr\"odinger functional (SF). In this paper we construct a family of boundary fields for general values of N which enter the standard definition of the SF coupling. We provide spatial boundary conditions for fermions in several representations which reduce the condition number of the squared Dirac operator. In addition, we calculate the O(a) improvement coefficients for N>3 needed to remove boundary cutoff effects from the gauge action. After this, residual cutoff effects on the step scaling function are shown to be very small even when considering non-fundamental representations. We also calculate the ratio of Lambda parameters between the MS-bar and SF schemes.


I. INTRODUCTION
Asymptotically free theories, such as gauge theories coupled to fermionic matter fields 1 , are characterized by having a coupling which becomes small at short distances.
This property enables reliable perturbative calculations of physical quantities at large energies. A dimensionful scale is dynamically generated through the process of dimensional transmutation. Typically, this scale is associated in perturbation theory with the Λ parameter, i.e. a multiplicative constant of the integrated beta function.
The non-perturbative evolution of the running coupling in different gauge field theories from the low energy sector to the high energy regime has been the central goal of many studies. The standard approach is the use of a finite size scaling technique based on the Schrödinger functional (SF), in which the size of the system is associated to the renormalization scale [1]. This method was successfully used to calculate the scale evolution of the coupling in the SU(2) [2] and SU(3) [3] Yang-Mills theories and in QCD [4]. Motivated by ideas of physics beyond the standard model (BSM), in the last decade this method has also been applied to study the SU(4) pure gauge theory [5] and several theories containing matter transforming under higher dimensional representations of the gauge group or a large number of fermions in the fundamental representation [6][7][8][9][10][11][12][13][14].
However, for lattices accessible in typical numerical simulations SF schemes are affected by lattice artifacts arising from the bulk and from the boundaries of the lattice.
These can be removed, following Symanzik's improvement program, by adding the corresponding counterterms to the action at the bulk and the boundaries. Symanzik's program was successfully carried out in [1, [15][16][17][18], where the improvement coefficients necessary to remove O(a) effects from the coupling were calculated in perturbation theory.
For theories beyond QCD the situation is still inconclusive. A program for the nonperturbative study of SU(N) gauge theories in the large N limit [19] started in the last decade driven by interest from string theory. As part of that program, in [5] the ratio σ between the lambda parameter in the MS scheme and the string tension σ was calculated for the SU(4) theory aiming to obtain extrapolations of the N dependence of the Λ MS / √ σ in the large N limit. There, the dominant systematic errors are due to the lattice artifacts present by using an unimproved action. For the case of theories with non-fundamental fermions, although the O(a) improvement coefficients are known, the remaining higher order cutoff effects have been reported to be very large if the standard setups, which work fine for QCD, are naively exported.
In the last few years a new renormalized coupling based on the gradient flow (GF) has been proposed for step scaling studies [20][21][22][23][24]. Compared to the original SF coupling based on a background field (see bellow), the gradient flow coupling has the advantage that considerably smaller statistics are required for obtaining a similar accuracy.
However, there are some situations where the original SF coupling is superior compared to the gradient flow. First of all, it has been observed that while the GF coupling works better at large physical volumes, at small volumes the SF coupling fares better than the gradient flow [25]. Also, in the pure gauge theory, relevant for the large N limit, the generation of configurations is so fast compared to the measurement of the gradient flow that the reduced accuracy can be overcome with increased statistics. In addition, in BSM lattice studies one is often interested in the existence of a nontrivial infra-red fixed point.
The value of the coupling constant g at the fixed point is a renormalization scheme dependent quantity and it differs between Schrödinger functional and gradient flow schemes.
Therefore, it is possible that in a specific scheme the coupling is too strong at the fixed point or it is on the wrong side of a bulk phase transition. This is true even if the fixed point is visible in other schemes. The only study, we are aware of, that compares these two methods with the same action found the gradient flow coupling to be about twice the Schrödinger functional coupling [26]. Moreover, due to the property of continuum reduction [27], at large N it is possible to do simulations at small lattice volumes where the SF coupling is known to perform well.
This work completes the Schrödiger functional framework to study the phase diagram of strongly interacting gauge theories [28] with any N or representation. In the paper we generalize the boundary conditions for the gauge fields in the SF to obtain a family of schemes useful for arbitrary N with a good signal to noise ratio in lattice simulations. Another appealing property of the present family of schemes is that, together with an appropriate choice of spatial boundary conditions for the fermions, they lead to a setup for which higher order cutoff effects due to fermions are very small even for non-fundamental representations. Preliminary results of this work have been published in [29].
The paper is organized as follows: in section II we recall some concepts concerning the Schrödinger functional and collect a set of formulas useful for the remaining discussion.
In section III the generalized boundary conditions are provided. The calculation of the improvement coefficients is presented in section IV, where we also discuss the effect that the fermionic spatial boundary conditions have on the residual higher order cutoff effects.
The matching of the Λ parameters to the MS scheme is done in section V. We conclude in section VI.

II. SCHRÖDINGER FUNCTIONAL
In this section we briefly recall the ideas introduced in [1, 15,30] and collect the expressions necessary for the subsequent discussion. We refer the interested reader to the original articles for further detail.
The Schrödinger functional is the euclidean propagation amplitude between a field configuration at time 0 and another field configuration at time T, which has a path integral representation given by with Dirichlet boundary conditions specified for the gauge fields U and the fermion fields ψ andψ.
In the present work, we are interested in the O(a) improved Wilson action The pure gauge part is the standard SU(N) Wilson gauge action The spatial components of the gauge fields at the temporal boundaries (t = 0 and t = T) satisfy nonhomogeneous Dirichlet boundary conditions The boundary gauge fields W k and W k can be parametrized as where C k (η) and C k (η) are taken to be homogeneous, abelian and spatially constant [1], and they depend on a dimensionless parameter η. A specific form for these boundary matrices is derived in section III for gauge group SU(N) with arbitrary N. In the spatial directions the gauge filds are taken to be periodic U µ (t, x) = U µ (t, x + Lk).
The weight w(P) = 1 except for the spatial plaquettes at the boundaries for which w(P) = 1 2 . Due to the particular choice of boundary conditions for the gauge fields, the spatial boundary plaquettes give only a constant contribution to the action and can be ignored. It is well known that within SF schemes, the mere presence of temporal boundaries constitutes an extra source of lattice artifacts. Removal of these effects has first been studied in [1, 16,17], where it was shown that the O(a) lattice artifacts coming from the boundaries can be canceled by tuning the weight w(p) = c t (g 0 ) for the temporal plaquettes The fermionic part of Eq. (2) is the standard Wilson fermion action with the clover term where D WD is the improved Wilson-Dirac operator The operator F µν (x) is the symmetrized lattice field strength tensor, σ µν = i 2 γ µ , γ ν and the operators D µ and D * µ are the covariant forward and backward derivatives which are defined in Eq. (A8). The improvement coefficient c sw can be determined perturbatively [16,31] and non-perturbatively [32,33]. To the lowest order in perturbation theory c sw = 1 [18]. The removal of O(a) effects arising from the interplay between fermions and the SF boundaries requires the addition of another dimension 4 counterterm at the boundaries.
Since this does not contribute to the observables studied further in this work at the present order in perturbation theory, we ignore it from now on and refer the reader to the original literature [17] for further details.
The fermionic fields satisfy the following boundary conditions where P ± = 1 2 (1 ± γ 0 ). The boundary conditions in the spatial directions are periodic up to a phase [30]: The phase is usually chosen so that the smallest eigenvalue of the squared Dirac operator is large [30]. In this situation, the condition number (i.e. the ratio between the highest an lowest eigenvalues) is small, which improves the speed of the known inversion algorithms. However, the value of θ also has an effect on the convergence of the 1-loop perturbative coupling to its continuum limit. This is discussed in subsection IV C.
The boundary conditions for the gauge fields in Eqs. (4) and (5) with the effective action Γ = − ln Z. The normalization constant is defined so that g 2 = g 2 0 to the lowest order of perturbation theory. One of the central quantities in numerical simulations is the step scaling function This is required for reconstructing non-perturbatively the scale evolution of the running coupling. In presence of a lattice regulator, the deviations of the lattice counterpart of the step scaling function Σ(u, L/a) from the continuum σ(u) can be used to monitor the size of cutoff effects (see subsection IV C).

A. 1-loop expansion
The renormalized coupling Eq. (11) is suitable for both perturbative and non-perturbative  [34][35][36]. In the present work we extend the previous calculations to arbitrary N. Although the main strategy of the calculation follows closely previous works, some care has to be taken to generalize those ideas without complicating the calculation. In the present subsection we collect some formulas necessary for the subsequent discussion and leave all technical details on the calculation to Appendices A and B.
After going through the gauge fixing procedure [1], the effective action is expanded to 1-loop as Here Γ 0 is the classical action. The 1-loop term Γ 1 in the effective action can be written as where ∆ 0 , ∆ 1 and ∆ 2 are the quadratic ghost, gluonic and fermionic operators respectively.
The explicit forms of the operators ∆ i are given in Appendix A.
According to Eq. (14), the 1-loop coefficient receives an independent contribution from ghost, gauge, and fermionic fields with The gauge and fermionic contributions to Eq. (18) can be calculated independently.
The gauge part is given by The calculation of h 0 (L/a) and h 1 (L/a) has been described in great detail for SU(2) in [1] and the calculation has been done for N = 3 in [3]. In Appendix A we give the generalization of the calculations to N ≥ 3.
The calculation of the fermionic part p 1,1 (L/a) is straight forward to generalize to any boundary fields and to any representation of the gauge group. One just needs to replace the link variables in the Wilson Dirac operator Eq. (7) with their counterparts in the desired representation. Thus we will refer the interested reader to the original paper [15].
The continuum and lattice step scaling functions are given to first order in perturbation theory by with σ 1 = 2b 0 ln(2). The 1-loop coefficient b 0 of the beta function is given in an arbitrary representation by where the color group invariants are defined as in the representation R of SU(N) 3 .
Similarly as in Eq. (18), the step scaling functions can be separated into a gauge and a fermionic part, This allows us to study separately cutoff effects due to gauge and fermion fields independently. 3 The values of the invariants are given by

III. BOUNDARY FIELDS FOR N > 2
In this section we present a generalization of the boundary fields for N > 2. The selection of the boundary fields is only limited by the requirement that there is a unique and stable classical solution to system. In practice, this limits us to Abelian boundary fields W k and W k which can be written as in Eq. (5), where Since W k has to be an SU(N)-matrix the vectors Now the classical solution, i.e. the background field, can be written as where It is shown in [1] that the solution V µ (x) is absolutely stable if the vectors φ and φ satisfy Eq. (26) and These conditions define a fundamental domain, which is an irregular (N − 1)-simplex and has vertices at points To define a renormalized coupling we can choose any two different points inside the fundamental domain to set up the boundary fields. A different choice leads to a different renormalization scheme, which can be matched to each other using perturbation theory (see section V). However, there are practical considerations in selecting the boundary fields, namely the signal to noise ratio in the Monte Carlo simulations and the size of higher order lattice artifacts. Our choice is based on the attempt to maximize the signal to noise ratio as in practice the minimization of the higher order lattice artifacts often leads to a low signal, which neglects the gains of a better continuum extrapolation.
To obtain a maximal signal strength we have two competing requirements. We need to twist the gauge fields as much as possible while staying away from the boundaries of the fundamental domain. This is because the coupling is proportional to the twist and because closeness of the instability of the classical solution increases noise. According to these considerations we choose φ to be in the middle of a line connecting X 1 and the centeroid of the fundamental domain To determine φ we find a transformation which is a map from the fundamental domain to itself and mirrors the vertices. First we define a simple map R i, j (φ) that reflects the points in the fundamental domain with respect to a (N − 2) dimensional hyperplane. The hyperplane R i,j (φ) goes through vertices X k , k i, j and intersects the line connecting X i and X j at the middle. For N > 3 the function R i,j (φ) is not in general a mapping from the where is a mapping from the fundamental domain to itself and written in components it has a simple form To define the coupling we choose a one parameter curve of background fields φ + t(η).
We select it in a way that the results are equivalent to those of the SU(3) theory given in [3] 4 , i.e. we select t(η) so that it changes sign under the mapping M(φ) and points towards the boundaries of the fundamental domain See Fig. 1 for illustration of the fundamental domain and boundary conditions for SU(4) and Table I for the boundary values for N = 3, 4, 5.
In the lattice computations it is advisable not to set t(η) beforehand, but to measure a complete N − 1 dimensional basis which can be used to construct a generic curve. Each curve corresponds to a different renormalization scheme.

A. Fermionic spatial boundary conditions
Recalling the spatial boundary conditions for the fermion fields Eq. (10), we still have to choose a particular value for the angles θ k . For simplicity, we consider the same angle in all spatial directions θ = θ k , k = 1, 2, 3.
We then fix θ, following the criteria introduced in [15], so that the minimum eigenvalue λ min of the fermion operator ∆ 2 is as large as possible. This leads to a small condition number which optimizes the speed of the numerical inversion of the operator.
The values of θ leading to a maximum λ min depend on the background field and also on the fermion representation being considered. For the fundamental representation, the profile of smallest eigenvalues λ min as a function of θ is shown in Fig. 2 for the different gauge groups considered in this work excluding the case of SU (2).
Although the maximum of λ min is achieved at different values of θ for every gauge group considered, the choice θ = π/2 is always close to the maximum and hence leads to a small condition number. For homogeneity in the definition of a renormalization scheme in the subsequent calculations, we will fix θ = π/2 for all values of N. As we will show in Section IV C, this choice of θ together with the family of background fields defined in this work will lead to a setup for which higher order cutoff effects are highly  is taken considering the fundamental representation, this value also leads to reasonably small condition numbers for the symmetric and antisymmetric representations. For the adjoint representation the smallest condition number is obtained for θ = 0. However, we decide to stick to the choice θ = π/2 also in this case since it leads to a situation were higher order lattice artifacts are highly reduced 5 .
For the case of SU(2), we choose θ = 0 for the fundamental representation but leave θ = π/2 for the symmetric/adjoint.   5 We consider the reduction of higher order cutoff effects, which have been shown to be very large [34][35][36], to be of higher importance.
where s 0,i = 2b 0,i and s 1,i = 0 after setting c SW to its tree level value. The boundary improvement coefficients c (1,i) t are determined by demanding linear cutoff effects to be absent in Eq. (37), which is achieved by fixing c (1,i) t = r 1,i /2. The continuum coefficients r 0,i are needed when matching the Λ parameter to other schemes.
In order to extract the coefficients r n,i as accurately as possible we first evaluate p 1,0 (L/a) and p 1,1 (L/a) adapting the strategies in [1,15] to general N (see appendices A and B for details on the calculation). Once the series of data for p 1,i (L/a) is produced, the coefficients r n,i can be extracted using a suitable fitting procedure.
In the pure gauge case, we calculated

Results
In Table II we where C 2 (R) is the quadratic Casimir operator in the representation R.   The data and the fit are shown in Fig. 3. To check the consistency of our results, we have also performed fits adding additional terms to Eq. (38). The coefficients of the additional terms are zero within statistical errors as shown in Table IV. For completeness, we also include here the value of the fermionic part c (1,1) t for an arbitrary group and representation This was calculated for the fundamental representation in [15] and later extended to other represesentations in [35]. In the present work we have been able to reproduce the value of c (1,1) t with similar accuracy, which is a further check on the correctness of the whole calculation.

C. Residual cutoff effects
The determination of the gauge and fermion contributions to c (1) t removes O(a) lattice artifacts coming from the boundaries to 1-loop in perturbation theory. However, cutoff    It is remarkable that the value θ = π/2, established in section IV A to obtain a condition number as small as possible, also leads to a situation where higher order cutoff effects are highly suppressed.
A very similar picture is observed when considering any of the 2-index representations.
Cutoff effects for the adjoint, symmetric and anti-symmetric representations are shown respectively in panels (b), (c) and (d) of Fig. 5. The smallness of the residual lattice artifacts is at first glance surprising, since they have previously been reported to be very large if particular care is not taken in the choice of BF [34][35][36]. The magnitude of δ 1,1 for the 2index representations strongly depends on the angle θ in a very similar way as it is shown 7 Note that the value θ = π/5 chosen in [15] leads to very reduced cutoff effects for their choice of BF. This would not be the case if this value was used with our BF. in Fig. 6 for the fundamental representation. It is then possible to tune θ to minimize cutoff effects without the need of modifying the BF [35,36]. What is remarkable of the family of background fields proposed in this work is that for the fundamental, symmetric and antisymmetric representations, values of θ which lead to small condition numbers also lead to small higher order lattice artifacts in the step scaling function. This is not true for the adjoint representation since, as discussed in section IV A, the condition number is minimized for θ = 0. It is also remarkable that cutoff effects for all the representations considered are minimized for the same value θ = π/2. The Λ parameter is a renormalization group invariant and scheme dependent quantity given by (in a generic scheme X) It is a dimensionfull scale dynamically generated by the theory.
In subsection IV B 2 we have performed the computation of the SF coupling 8 g SF in the Schrödinger Functional scheme to one loop order in perturbation theory, i.e. we have 8 In the following we write explicitly a subindex with the coupling to indicate the scheme. calculated the renormalized coupling as an expansion in terms of the bare coupling g 0 g 2 SF (L) = g 2 0 + p 1 (L/a)g 4 0 + O(g 6 0 ), (42) where, after doing a continuum extrapolation, To be able to compare the results at different values of N, we are interested in a relation between α SF = g 2 SF /4π and some scheme where N is only a parameter. For that we choose the usual MS scheme, defined at infinite volume and at high energies. The relation between the running coupling in the two schemes can be written as an expansion where s is a scale parameter and c i (s) are the coefficients relating the couplings in the two schemes at each order in perturbation theory.
The relation between the Λ parameter in the SF and MS scheme is given by where c 1 (1) is the coefficient of the 1-loop relation (44). Note that Eq. (45) is an exact relation even though it depends on the 1-loop coefficient relating the couplings in two different schemes.
For determining the coefficient c 1 (s) in Eq. (44), we first use the known relation between α MS in the MS and the bare coupling α 0 [17,30,33,40], which at 1-loop is given by The 1-loop coefficient d 1 (s) is given for generic N and fermionic representation R by where The coefficient k 1 of the gauge part is taken from [17,30,33] and reads The coefficient K 1 is a representation independent function of the tree level coefficient c (0) SW given by It was calculated in [41][42][43] for the fundamental representation and extended to arbitrary representations in [40].
with r 0,i being the continuum coefficients in the series (37).
Knowing this, the relation between Λ parameters in Eq. (45) can be given as a function of the parameters N, N f and T R and of the coefficients r 0,i Finally, in Table V we collect the values of the ratio of Λ parameters for the schemes studied in this work and for the pure gauge theory (N f = 0) and for 2 flavors of fundamental fermions 9 . Ratios of lambda parameters for 2 index representations can be recovered using Eq. (52) and the corresponding coefficients from Table (III).

VI. CONCLUSIONS
We have studied the Schrödinger functional boundary conditions and the perturbative The fermionic twisting angle θ is also studied and we found out that the value θ = π/2 is a good compromise between the simulation speed and the minimization of the O(a 2 ) lattice artifacts in the perturbative 1-loop lattice step scaling function. Greek indices α, β, γ, . . . are always summed over unless otherwise stated in the formula.
The operators we are interested in are defined as There is no summation over µ in the r.h.s. of the Eq. (A2). The star product in Eq. (A2) which maps an N × N matrix M and an SU(N) matrix X to an SU(N) matrix is defined as In Eq. (A3) the operator D WD is the same as in Eq. (7) with c sw = 1.
The first step is to find a suitable basis for the SU(N) generators. This is a basis that is invariant under the star product defined in Eq. (A4). In practice we want to find generators X a that satisfy cosh G 0k X a = χ c a X a , sinh G 0k X a = χ s a X a , with arbitrary coefficients χ c a and χ s a . The hyperbolic sine and cosine of the non-zero elements of the field strength tensor are A basis that satisfies Eq. (A5) for the non-diagonal generators are the ladder operators defined as where n and m are the matrix indices and a( j, k) is the color index. The properties of a( j, k) are given in Table VI. The generators X a are normalized as Tr X a X b = − 1 2 δ a,b . The diagonal generators can be chosen in any way that satisfies Eq. (A5). The boundary conditions generate a background field V µ (x) defined in Eqs. (27), (28) and (29) which enters the covariant derivatives and through them to the operators ∆ s . The next step is to calculate the covariant derivatives with the background field V µ (x) when q(x) = q a (x)X a is proportional to a generator X a .
The covariant derivatives can then be written in a general form where and t is the time component of the four vector x = (t, x). The operators ∆ 0 and ∆ 1 can now be decomposed to color subspaces according to the basis selected. The operators are also invariant under spatial translations and thus the determinants can be written as where p = 2πn/L, is the three momentum.
Next we will show how one can calculate the determinant of operators ∆ s . In [1] it has been shown that for an operator ∆ that satisfies for matrices A, B and C and an eigenvalue equation there exists a matrix M(ξ) such that The determinant of ∆ is then given by We will next use these properties of the ∆ s operators.
Since the operator ∆ 0 is invariant under spatial translations and constant diagonal gauge transformations its eigenfunctions are of the form Operating with ∆ 0 on ω a (x) we get with A = C = −1 and Clearly the operator ∆ 0 is similar to the operator in Eq. (A14) and thus the strategy shown can be used. Using the Eq. (A15) with ξ = 0 i.e.
we get a recursion relation for ψ a (t) with initial values ψ a (0) = 0, ψ a (1) = 1 which is According to Eq. (A17) the determinant is then We will then move on to the more challenging case of ∆ 1 . The eigenfunctions of the operator ∆ 1 have the general form where normalization 10 matrix R a µν (t) is a diagonal 4 × 4 matrix with R a 00 (t) = −i, R a kk (t) = e i(p k + f a (t))/2 , k = 1, 2, 3.
Again we can operate with ∆ 1 on the eigenfunction Eq. (A25) which yields where the matrices A a µν (t), B a µν (t) and C a µν (t) are B a 00 (t) = 2λ 0 + 3 k=1 s a k (t) χ c a s a k (t) − iχ s a c a k (t) , B a kl (t) = (λ 0 − 1)s a k (t)s a l (t) + δ k,l (2χ c a + 3 n=1 (s a n (t)) 2 ), B a 0k (t) = χ c a s a k (t) − iχ s a c a k (t) − λ 0 s a k (t), B a k0 (t) = B a 0k (t), We have used the following short handed notation c a k (t) = 2 cos (p k + f a (t))/2 , (A29) s a k (t) = 2 cos (p k + f a (t))/2 , (A30) 10 Adding R a µν ensures that the matrices A a µν (t), B a µν (t) and C a µν (t) in the recursion relation are real.
The operator ∆ 1 is also similar to the case in Eq. (A14) and the same strategy can again be exploited. Additionally the boundary conditions of ψ a µ (t) in Eq. (A27) are With this we can first calculate the determinant of ∆ 1 in the more general case where the boundary conditions are given by Eq. (A33). Setting ξ = 0 in Eq. (A15) we get With the help of these equations we find F a µν (t) which has the property where are the first nonzero components of ψ a µ (t). The matrix F a µν (t) is B a ρσ (1)F a σν (1) + C a ρσ (1)P σν , where the projection operator P µν is With F a µν (t) we will be able to construct a matrix M a µν that couples v a µ from Eq. (A37) and the boundary condition Eq. (A33) at t = L This matrix M a µν turns out to be M a µν = F a µν (L) − P µρ F a ρν (L − 1), (A41) and the determinant of ∆ 1 in this subspace according to Eq. (A17) is det ∆ 1 | (p,a) = det M a µν λ L 0 (N a ) 3(L−1) .
We can then move on to the case of ∆ 1 where a > N 2 − N i.e. for diagonal generators X a and when p = 0. In this case the boundary conditions are given by Eq. (A32) and ψ a 0 (t) and ψ a k (t) components decouple since the matrices A a µν (t), B a µν (t) and C a µν (t) are diagonal where prime indicates partial derivative w.r.t. the parameter η, M (−1) µν is the inverse matrix of M µν and the normalization κ is defined in Eq. (12).

Appendix B: Chosen basis for the diagonal generators and the values of the coefficients which depend on the background field
In Appendix A we showed how the 1-loop coupling can be calculated for a generic background field and basis of generators. In here we will specify the basis that we have selected as well as the values of the coefficients χ c a , χ s a and f a (t). We have chosen a basis given by where The coefficients f a (t) are 11 for a ≤ N 2 − N and f a = 0 for a > N 2 − N.