Scale setting for large-$N$ SUSY Yang-Mills on the lattice

In this paper we study the large $N$ limit of four-dimensional Supersymmetric Yang-Mills on the lattice using twisted reduced models. We have generated configurations with dynamical massive gluinos and various lattice 't Hooft couplings, and verified that the Pfaffian remains positive. We have determined the lattice spacing in terms of various observables obtaining compatible results. Extrapolating to the massless gluino limit we obtain the lattice spacing dependence on the bare couplings for the supersymmetric theory. The observed dependence goes along the expected behaviour predicted by perturbation theory.

1 Introduction This paper deals with N = 1 Supersymmetric Yang-Mills theory in the large-N limit. Hence, it combines two of the ingredients which generally add a simplification to the pure Yang-Mills theory and appear in less traditional approaches to the study of quantum field theories [1]. Our work, however, is based on the more standard Lattice Gauge Theory approach, which has proven well suited to address the non-perturbative dynamics of QCD and other gauge theories. In this context, however, both Supersymmetry and the large N limit pose additional problems to be faced within this approach 1 . Supersymmetry is broken by the lattice formulation and has to be recovered in the continuum limit (For an introduction to Supersymmetry on the lattice the reader can consult ref. [2]). The advantage in our case is that the N = 1 SUSY Yang-Mills theory is just pure Yang-Mills with an additional massless Majorana fermion in the adjoint representation: the gluino. The argument made by Kaplan [3] and Curci and Veneziano [4] is that in the limit in which the mass of the gluino vanishes one recovers SUSY invariance. With this idea in mind several authors have studied the problem in the past, usually for small rank groups as SU (2) and SU (3). Given the important role played by chiral invariance in this context a good deal of work has focused upon overlap and domain-wall fermion lattice formulations [5][6][7][8][9][10][11].
Another well-explored possibility is to use Wilson fermions [4,12,13]. One disadvantage of this last approach is that the positive-definite character of the Pfaffian is not guaranteed and has to be checked or corrected for. There is an extensive study done for SU (2) and SU(3) by a collaboration involving Desy, Münster, Regensburg and Jena Universities [14][15][16][17][18]. We will essentially follow the same discretization and basic methodology. The other ingredient in our work is the study of the large N limit of the theory which, as mentioned at the beginning of the first paragraph, is a further element in matching with results accessible to other methodologies. Although both perturbation theory and strong coupling results [19] simplify in the large N limit, dealing with the increasing number of degrees of freedom in Monte Carlo methods in a theory like this with dynamical fermions represents a formidable challenge. Our approach to the large N limit is based on volume reduction [20], which makes a more efficient use of both spatial and internal degrees of freedom. For the pure gauge part we employ the so-called twisted Eguchi-Kawai model [21][22][23] which has been used successfully to study pure Yang-Mills theory in the large N limit [24][25][26]. As explained in ref. [22] the perturbative proof of the reduction mechanism also holds for theories with fermions in the adjoint representation. This was applied previously to non-supersymmetric theories with various flavours of adjoint Dirac fermions [27]. It was used to study the string tension and other features [28][29][30][31] and most importantly to analyze the possible infrared conformal behaviour of these theories and compute its mass anomalous dimension [32]. Here we will apply the same methodology to the case of one Majorana fermion. There are some inherent limitations of the volume reduction method which prevent us, for the time being, of accessing interesting quantities as the meson and glueball spectrum. However, there are other physical quantities that can be computed and in particular the study of the scale-setting of the theory, a necessary step in connecting lattice quantities to continuum ones. Some partial results were already presented at the 2021 Lattice conference [33,34].
In the next sections we will spell out the details of our lattice approach and computational techniques. Section 2 presents the model used in our study, involving a volumereduced version of Wilson action with twisted boundary conditions and a fermionic part of the action corresponding to Wilson fermions in the adjoint representation. The Majorana character shows up in the appearance of the Pfaffian, which poses potential non-positivity problems. We also explain how for finite N , our gluons and gluinos can be thought of as propagating in a four-dimensional box of size a √ N . This is an important piece of information in interpreting the range of parameters used in our study. The present work is a rather extensive study involving four values of the lattice coupling, three values of N and several values of the gluino mass (hopping parameter). The methodology used in generating the configurations is explained in Section 3. We also present our results regarding the sign of the Pfaffian. Then we proceed to locate the value of the hopping parameters corresponding to massless gluinos in Section 4. These are precisely the points at which our model will represent the N = 1 Supersymmetric Yang-Mills theory. The following two sections are dedicated to the main problem under consideration: determining the value of the lattice spacing in some units for each one of our simulation parameters. The first part (Section 5) explains the underlying philosophy and presents several possible continuum observables that will be used to set the scale. The following one (Section 6) elaborates on the lattice counterparts, the technical problems involved and the results obtained. Redundancy is a proof of the robustness of our results. The closing section is then devoted to extracting the corresponding scales for the N = 1 SUSY Yang-Mills theory, by extrapolating the previous results to the massless gluino limit. A comparison with expectations from perturbation theory in the continuum seems to work remarkably well for the range of values explored in the study. We close with a brief summary of the results.

N = 1 SUSY Yang-Mills at large-N on the lattice
In this section we specify the model which is used to generate our results. The starting point is the basic Wilson discretization of the N = 1 SUSY Yang-Mills action for arbitrary N . For the pure gauge part this is given by Wilson plaquette action: (1 − γ µ )V adj µ (n)δ(n +μ, m) + (1 + γ µ ) V adj µ (n −μ) † δ(n −μ, m) .
where κ a is the adjoint Majorana-fermion hopping parameter and V adj µ is the link variable in the adjoint representation. The charge conjugation matrix C satisfies the following properties

4)
C t = −C. (2.5) ensuring the antisymmetric character of the CD W matrix.
In the large N limit the twisted reduction prescription [22] implies that the model is equivalent to a matrix model obtained as follows (see also refs. [35,36]). One introduces 4 SU(N ) constant matrices Γ µ satisfying with given in terms of the twist tensor n µν , which is an antisymmetric tensor of integers modulo N . The displacement operator by one lattice point in the direction µ is then replaced by the adjoint action by the matrix Γ µ . For the link matrices this is equivalent to the substitution V ν (n +μ) = Γ µ V ν (n)Γ † µ . (2.8) After doing these substitutions, dropping the n label and changing variables to U µ = V µ Γ µ , the pure gauge action converts into the one of the twisted Eguchi-Kawai (TEK) model [21,22]: The same prescription applies to other fields in the adjoint representation such as the gluino fields. There are two ways to represent the adjoint fermion fields. If we represent them as traceless N × N matrices the prescription is similar to that of the links Ψ(n +μ) = Γ µ Ψ(n)Γ † µ . (2.10) If we represent these fields as N 2 − 1 dimensional vectors the prescription is the equivalent one: Ψ(n +μ) = Γ adj µ Ψ(n). (2.11) In this way one reaches the reduced form of the Wilson-Dirac operator: (2.12) -4 -This procedure was used already in ref. [27] in constructing the SU(N ) gauge theory coupled to several flavours of adjoint Dirac fermions. After integration over the fermion fields, we finally arrive to the model used in our simulations, whose partition function is given by where the integration over the group variables U µ is given by the Haar measure, and Pf(M ) denotes the Pfaffian of the antisymmetric matrix M . We still need to specify the choice of the twist tensor n µν . In our case, we have chosen the so-called symmetric twist given by: 14) with √ N and k coprime integers and with θ representing the Heaviside theta-function. The presence of the Pfaffian in eq. (2.13) poses the question of whether this defines a positive definite probability distribution. Indeed, the general properties of the Dirac operator in the continuum ensures its positivity, but some of these properties do not hold for the lattice Wilson-Dirac operator (they do hold for the overlap [6]). This brings in the menace of the so-called sign problem. Since the lack of positivity could only be driven by lattice artifacts we expect it not to be too severe. We take advantage of the studies done by the Desy-Jena-Regensburg-Munster collaboration which deals with the same situation [14][15][16][17]. We can resort to the well-known reweighting method to deal with the sign-problem. This amounts to defining a probability density given by substitution of the Pfaffian by its absolute value: The expectation value of an observable O in the original model is then evaluated as follows The right-hand side involves the ratio of two expectation values with respect to the probability density P (U ). In the next section we describe our methodology to generate configurations according to this distribution. There is one important consideration that has to be made in interpreting the results of our work. It concerns finite N effects. The reduction mechanism applies at infinite N , while the numerical results necessarily are performed at finite although very large N . Thus, it is very important to understand the nature and size of finite N effects. The perturbative proof of reduction shows that planar diagrams in the reduced model are identical to those of the ordinary model on a finite lattice of size ( √ N ) 4 . This is not the only difference between the finite N versions of the ordinary and reduced model, but very often is the dominant effect. Hence, it is useful to regard √ N as the effective lattice size and finite N errors acquire the form of finite volume effects. Our simulations which use N = 169, 289 and 361 correspond to an effective lattice size of 13 4 , 17 4 and 19 4 . From that viewpoint extended observables are more affected by this type of errors as compared to more local ones. For example, for Wilson loops the effects grow with the size of the loop and for smeared or flowed quantities with the corresponding smearing radius. A detailed study of these finite N errors on Wilson loops has been done both perturbatively [37] and non-perturbatively [38]. The finite-N /finite-effective volume identification is crucial in interpreting our choice of parameters and our results. For example, a large value of b implies a small value of the lattice spacing a in physical units. This translates into an effective lattice box of side a √ N and this should be kept larger than the characteristic correlation lengths of the theory. This limits the maximum values of b in much the same way as the lattice size does for the ordinary lattice simulations. In conclusion, it is clear that going to larger values of N is of course desirable but this is limited by the computational effort involved.

Generation of configurations
In this section we will explain the basic characteristics of our data sample as well as the methodology used to generate it. We employ a rather extensive set of configurations generated with the probability distribution explained in the previous section and corresponding to a total of 46 different simulation parameters. This involves three different values of the rank of the group N = 169, 289 and 361, four different values of the inverse lattice coupling b = 0.34, 0.345, 0.35 and 0.36, and four to six values of the gluino mass for each case. This facilitates the extrapolation to vanishing gluino mass, the analysis of the scaling behaviour of the theory and the appropriate control of all sources of errors. The final list of values of the hopping parameter and the number of configurations generated in each case are listed in table 7 in appendix A.
To generate the configurations we use the strategy explained in the previous section based upon making use of the reweighting method to deal with the problems associated with the sign of the Pfaffian. Hence, we split the problem into two parts. First, we generate configurations according to the positive definite distribution P (U ). This is done by making use of the Rational Hybrid Monte Carlo (RHMC) algorithm [39][40][41][42][43], whose specific features for the case at hand will be explained in the next subsection. Then, we will address the problem of the sign of the Pfaffian by analyzing its distribution over our data. This is explained in the following subsection. In fact our study is based on the analysis of the low-lying eigenvalues of the adjoint Wilson-Dirac operator, which is also interesting per se, and serves other purposes as the determination of the SUSY limit, addressed in the following section.
The interested reader can consult additional interesting technical details of our simulation which are described in appendix A.

Rational Hybrid Monte Carlo algorithm
To generate the ensemble of U µ according to the partition function eq. (2.13) via Markov chain Monte Carlo methods, we have to transform the integrand to be suitable for Monte Carlo simulations taking into account that the Pfaffian could be negative. By using the identity where Q W ≡ D W γ 5 , we transform the partition function as follows: and employ the RHMC algorithm [39][40][41][42][43] for the probability weight: Using the pseudo-fermionic integral, the factor det Q 2 W 1/4 is represented by: where φ is the pseudo-fermion field in bi-fundamental form for the adjoint representation and Tr is the trace over the colour index. The matrix R (3.5) Using the Remez algorithm, sets of the coefficients α (p) j=0,...,N R , β (p) j=1,...,N R have been prepared and tabulated for various cases: the power p = 1/8, −1/4, the dynamic range b/a and the order of the approximation N R .
To generate the pseudo-fermion field φ from the probability exp (−S Q ) in eq. (3.4), we set: where η is drawn from a Gaussian distribution exp − Tr η † η at the beginning of the molecular dynamics (MD) evolution in the RHMC algorithm. The traceless condition is imposed as Tr (η) = 0 and Tr (φ) = 0 to remove the unwanted U(1) contribution.
φ, required for the rational polynomial approximation, is evaluated using the multi-shift conjugate gradient algorithm. To apply D W to the pseudo-fermion field φ in bi-fundamental form, we use The coefficients α j=1,...,N R in the RHMC algorithm are determined so as to minimize the metric: where the interval [a, b] is set by the lowest and highest eigenvalues of Q 2 W . At the beginning and the end of each MD trajectory, a highly accurate approximation for the pseudofermionic field η and action S Q is needed in the HMC Metropolis test to have the correct weight, c.f. eq. (3.4). For this purpose a set of coefficients satisfying the condition: ∆ R < tol = 10 −14 , with the smallest N R , is chosen from the table of set-parameters. The smallest and largest eigenvalues of Q 2 W are calculated using a thick restart type Lanczos algorithm that is optimized to the matrix Q 2 W , having the property J Q 2 [44]. For the coefficient set used during the MD evolution, which should be invariant during the period of generation of configurations, we fix a lower accuracy ∆ R < 10 −5 -10 −8 and select a parameter set for (N R , a, b), with a and b tuned to cover the eigenvalue range determined during the thermalization period.
In addition, we employ the generalized multiple-time step MD integrator scheme [45,46] with the fourth order Omelyan-Mryglod-Folk scheme [47,48]. Three different time step sizes are assigned to the gauge action, and to the UV and IR parts of the pseudo-fermion action. The pseudo-fermion action in eq. (3.4) is separated at an order N split ∈ (1, N R ) as , we assign the first portion to the IR part and the last one to the UV part, resulting in a hierarchy in the magnitude of the MD forces [41]. The finest MD time step size is assigned to the gauge action, the next finer one to the UV part of the pseudo-fermion, and the coarsest step size to the IR part, respectively.

Eigenvalue spectrum and the sign of the Pfaffian
As mentioned in the previous section and the preceding subsection our strategy has been to deal with the effect of the sign of the Pfaffian sign(Pf(CD W )) in eq. (3.2) by means of the reweighting method. In this context the expectation value of an observable O is evaluated with where . . . P (U ) denotes the expectation value with the non-negative weight P (U ) which is statistically evaluated with the ensemble generated with the RHMC algorithm. To compute the sign of the Pffafian for our configurations we make use of the results of refs. [18,[49][50][51][52], according to which, the sign can be determined by counting the number of negative real eigenvalues of D W . To perform this search, we carried out a complete determination of all the complex eigenvalues of D W near the origin of the complex plane, using the ARPACK library [53]. The shift invert mode with the shift parameter σ = −0.1, which computes the extreme eigenvalues of (D W − σ) −1 , is used so that the eigenvalues having negative real part can be captured. Figure 1 shows, for all configurations, the 100 complex eigenvalue of D W closest to z = −0.1 in the complex plane. Given the symmetry of (CD W ) t = −(CD W ), each point represents a two-folded eigenvalue.
For heavy gluino masses we are not expecting flips of sign of the Pfaffian. Hence, our analysis focused on the lightest adjoint fermion masses at each b and N . We did not observe any negative-real eigenvalues in the spectrum. This absence was also checked for several heavier fermion masses. Therefore we conclude that the sign of the Pfaffian is always positive for the model parameters we employed, simplifying the reweighting method by validating the distribution using the absolute value of the Pfaffian. Similar results showing that the negative sign of the Pfaffian are rare even with moderately light adjoint masses has been observed in lattice SUSY models [16,54].

SUSY restoration limit
As mentioned in the introduction, N = 1 supersymmetry is broken by the lattice discretization and only appears as an emerging feature as one approaches the continuum limit and properly tunes the gluino mass to zero [3,4]. From a practical point of view, continuum and chiral (massless-gluino) limits are delicate to perform on the lattice. Nevertheless, while extrapolations to the continuum limit are quite standard in every lattice simulation, the methodology regarding how to explore the limit in which the gluino mass vanishes requires a special consideration, particularly when, as in our case, Wilson-Dirac fermions are employed to simulate the gluino. With Wilson fermions chiral symmetry is explicitly broken by an O(a) mass term. The consequence of this is that the bare mass of the fermion acquires an additive renormalization and a tuning of the hopping parameter κ a is required to reach the chiral limit. The way this tuning is done in QCD takes advantage of the fact that pions are pseudo-Goldstone bosons for the spontaneous breaking of chiral symmetry, and their mass is proportional to the square root of the renormalized quark mass. On the other hand, in N = 1 SUSY Yang-Mills, the chiral-symmetry breaking pattern is different, and no Goldstone excitation appears in the spectrum, one has instead an excitation analogous to the QCD η . Several methods have been proposed in the literature to attain in this case the limit where SUSY is restored. One of them is to define a renormalized gluino mass using the supersymmetric Ward identities regularized on the lattice and tune this mass parameter to zero [15]. Another possibility is to use the connected part of the adjoint-η correlator, this leads to a non-singlet adjoint-pion correlator that can be seen as composed of the gluino and a quenched Majorana fermion. On the basis of partially-quenched chiral perturbation theory [17], this adjoint-pion is a pseudo-Goldstone boson for the spontaneous breaking of the "partially quenched" chiral-symmetry and can be tuned to vanishing mass to restore supersymmetry [14,16].
In this work, the strategy used to determine the point of SUSY restoration consists of two very different methods. The first one makes use of the analysis of the eigenvalues of the Wilson-Dirac matrix, whose calculation has already been tackled in section 3.2. Let us call λ 2 min the minimum eigenvalue of Q 2 W . At infinite volume, we expect this quantity to vanish in the chiral limit. In our case this corresponds to it vanishing at large N , as already observed in the context of Yang-Mills theories with N f = 2 Dirac flavors [32]. Hence, we can use its dependence on the fermion mass to define a proper massless gluino limit. The second one is in the same spirit of the aforementioned adjoint-pion mass tuning. Although in the target supersymmetric theory the gluino field is described by a single flavour, we make use of an additional "quenched" flavor which allows us to define a pseudoscalar P =Ψγ 5 Ψ and an axial current A µ =Ψγ µ γ 5 Ψ. We can therefore define an analog of the PCAC mass as where |π stands for a generic state with the quantum numbers of the adjoint-pion. On the basis of partially-quenched chiral perturbation theory, we expect this mass to be directly proportional to the renormalized gluino mass and we tune it to zero to find the SUSY  Table 1: Values of the adjoint critical hopping parameter κ c a determined from the vanishing, as a function of 1/(2κ a ), of |λ min | or the m pcac mass. restoration limit.
4.1 Determination from the eigenvalues of D W γ 5 As already mentioned, the first method that we have used to determine the point of SUSY restoration is based on the eigenvalues of the Dirac matrix. In the continuum limit the minimum eigenvalue of the massive Dirac operator at infinite volume goes to zero linearly with the fermion mass. On the lattice one expects a similar behaviour for the operator D W /κ a which tends to the continuum one up to a renormalization factor. As explained earlier the lowest lying spectra of Q W = D W γ 5 has been determined for all our configurations. This includes |λ min |, whose average value is listed in table 8 for all our datasets. Notice that this quantity has very small errors, making it a perfect observable for the determination of the zero mass point. Indeed, the aforementioned continuum behaviour implies that |λ min | /κ a should depend linearly on 1/(2κ a ) and vanish at κ a = κ c a , signalling the point of vanishing gluino mass. Since finite N corrections amount to finite volume effects, this behaviour is only expected to hold at infinite N . Building on this analogy we would expect a dependence of the following form: where the function δ(1/N 2 ) should vanish at large N . A correction of this type was observed earlier [32] for the large N reduced model coupled to two flavours of adjoint fermions. In our case the formula describes our data very well as shown in figure 2. The dashed lines gives our best fit with δ(1/N 2 ) = B/N 2 and the dotted line the extrapolation to infinite N . The resulting critical hopping parameters κ c a are presented in table 1. In the next subsection we will present an alternative method for determining the critical hopping values.

Determination from the PCAC mass
An alternative way of extracting the point of SUSY restoration is based on the vanishing of the non-singlet adjoint-pion correlator and the associated PCAC mass. Contrary to what happens in standard lattice simulations, and due to the form in which it is computed in the reduced setup, the signal-to-noise ratio in the adjoint-pion correlation function deteriorates significantly at large time [55]. For this reason, it turns out to be more convenient to do the  tuning in terms of the PCAC mass which can be determined with very good precision, and this is the choice we have used in this work. As we will see below, the results obtained in this way turn out to be consistent with the less precise determination from the adjoint-pion mass.
In this and other sections we will be using correlators of meson operators made of fundamental quarks or adjoint gluinos. The methodology is simple as it amounts to allowing the fermions to travel in a lattice of varying sizes, including infinite, in the background field of the reduced model gauge field. This has been presented in ref. [55] and used extensively in ref. [26]. For other representations, like the adjoint, the philosophy is the same (details will be given elsewhere ref. [56]). Some specific features about the implementation in this work are contained in appendix D.
On the lattice, the PCAC mass can be determined from a discretized version of eq. (4.1). For that, we consider an ultralocal version of the axial vector current and an interpolating operator for the pion with optimized projection onto the pseudoscalar channel ground state, obtained as detailed in appendix D. The PCAC mass is then determined by fitting to a constant the ratio of correlation functions: where C(x 0 ; A, B) stands for the two-point correlation function of the operators corresponding to channels A and B.  provided in table 2. We include in the table results for our datasets with N = 289 and 361. The general agreement, within errors, of the results for the two different values of N gives an indication of the small influence of finite N effects on this quantity. In the chiral limit the PCAC mass is expected to be proportional to the fermion mass. This implies that close to this limit m pcac should depend linearly on 1/κ a and vanish at the critical value. Indeed, in figure 3 we display our results, for the datasets with larger value of N , together with a fit to a linear function of 1/(2κ a ). We see that the straight line provides a very good qualitative description of the data. Given the small errors, an additional quadratic term is needed in some cases to get a fit with χ 2 per degree of freedom of the order of one, if we include also the largest masses. In any case, for the purpose of determining the point of vanishing of the PCAC mass both fits give identical results.
Our final results for the critical value of κ a , corresponding to the largest value of N simulated in each case, are collected in table 1. The values obtained are in remarkable agreement with the determination based on the eigenvalues of the Dirac operator. In the rest of this work, the vanishing of the PCAC mass will be used as the criteria to tune all other quantities to the zero gluino mass limit.
Finally, to check for the consistency of our approach, we have as well determined the adjoint-pion mass and analyzed its dependence on the hopping parameter κ a . Results are collected in table 2. Figure 4 displays m 2 π vs 1/(2κ a ). The coloured lines shown in the plot correspond to linear fits where, for each value of the bare coupling b, the critical value of the hopping parameter is fixed to the one determined from the PCAC mass; only the slope is left as a free parameter. Linearity works very well in all cases, giving χ 2 per degree of freedom of the order of one, except on the smallest lattice with bare coupling b = 0.36, where a deviation is observed for the lighter masses (the two lightest ones are excluded from    the fit). This difference can be due to finite size effects, which come out more pronounced for the pion mass than for the PCAC mass. For this reason it is also preferable to use the latter in order to determine the point of SUSY restoration.

Scale setting: choosing the observables
In this section we will present our results on the scale setting for our lattice model. By extrapolation to the massless gluino limit we will later achieve our main goal of determining the scale for the N = 1 Supersymmetric gauge theory. Scale setting is a standard and rather crucial step one has to face in lattice gauge theories. Different procedures start by defining what observable is going to be used as a unit of energy. On the lattice this observable will appear as a dimensionless quantity given by the product of the continuum quantity times the lattice spacing. Hence, by measuring this particular lattice observable we are actually measuring the lattice spacing in units of the inverse of that physical quantity. It is obvious that there are infinitely many ways to fix the scale corresponding to the infinitely many observables that one can choose. In theories with a massive spectrum one could simply take one of the masses as unit. For QCD or Yang-Mills theory the lowest glueball mass, the vector meson mass or the square root of the string tension are typical units of this kind. However, in practice making the right choice of a unit is a crucial one. A unit must be of the appropriate size to the objects to be measured, it should be precise, easy to compute and as insensitive as possible to typical statistical or systematic errors. From that point of view the natural type of units mentioned earlier might not be the optimal, since they involve taking asymptotic limits. Other possible family of choices involve dimensionless monotonous functions of an energy scale. The unit of energy is taken to be the value at which the function takes a particular numerical value. To this class belong some of the most common scales used nowadays by the lattice community such as the Sommer scale [57], based on the quark-antiquark potential, and those relying on the use of the gradient flow such as t 0 [58] or w 0 [59]. For a review of these and other common methods to fix the scale, the interested reader can consult ref. [60].
For theories like pure Yang-Mills or N = 1 SUSY Yang-Mills with a single relevant operator the ratio of two scales is just a constant, given by the ratio of the corresponding units. On the lattice this is only true close enough to the critical point and the phenomenon is called scaling.
In this work we will study three different methods to fix the scale. This will allow us to measure the lattice spacing in each of the three units, and in addition to test whether our data shows scaling at the values of b at which we are simulating. Obviously, scaling is expected to work better for the larger values of b. In making our choices of units we have taken into account three factors. First of all the precision of the corresponding observable. Second its insensitivity to the gluino mass. And third the insensitivity to finite N effects. The latter, as mentioned in Section 2, are equivalent to finite volume effects, and represent a very important source of systematic errors.
Here we will mention briefly what the methods are about and in the next three subsections develop in more detail what are the corresponding observables. The technical details of its implementation on the lattice as well as the final results will be covered in the following section.
The first method is a modification of gradient flow based methods devised to reduce finite volume effects in the determinations of t 0 and w 0 , equal in spirit to the one introduced for periodic boundary conditions in [61], but making use in this case of twisted boundary conditions that have well known advantages from the point of view of perturbation theory. The second, a variant of which has been already used in large N pure Yang-Mills theory [24], is analogous to the Sommer scale but involves Wilson loops of fixed aspect ratio. The advantage of this is that the extrapolation to asymptotically large time needed to compute the quark-antiquark potential is no longer required. Finally, the third method involves using the spectra of mesons made of quarks in the fundamental representation to set the scale. In the large N limit fundamental fermion loops are suppressed and the quenched approximation is exact. The computation of the meson spectrum can therefore be determined at no additional cost on our set of dynamical adjoint fermion configurations. These meson masses depend on an additional scale, which is the fundamental quark mass, but in the zero mass limit this dependence disappears and the mass of the lightest non-singlet vector meson becomes a natural scale for SUSY Yang-Mills.
We emphasize that the methods presented in the following subsections can also be used at finite N and for other gauge theories and are particularly well suited when finite size effects are an important source of concern.

Gradient flow observables
The gradient flow [58,[62][63][64] is a smoothing technique that has received much attention in recent years. The idea is to replace the original gauge fields A µ (x) by a set of flow timedependent fields obtained by solving the gradient flow equations: ) and leading to an effective smearing of the gauge fields over a length scale √ 8t. The flow time, having physical dimensions of a length squared, induces natural candidates for scale setting; one just has to find a dimensionless, flow time dependent, quantity and determine the flow time at which it equals a particular pre-fixed value. The most common choices used for this purpose are based on the quantity: where E(t) stands for the flowed energy density: Particular examples are the ones obtained by solving either of the two following implicit equations: where s = 0.1 corresponds to the standard choices in the literature denoted by t 0 [58] and All these considerations take place in infinite volume. However, numerical simulations are implemented on a L-site lattice of finite physical size l = aL, a fact that leads to finite size effects in the determination of the scale whenever the flow radius extends over a significant fraction of the box size. Our proposal, analogous to the one in ref. [61], is to take an alternative version of the observable Φ(t) with the correct infinite volume limit but with reduced finite size effects. The idea makes use of the fact that in infinite volume one can define a renormalized gradient flow (GF) coupling constant proportional to Φ(t) [58]: with a flow-time independent proportionality constant K(N ) fixed by imposing that the GF coupling equals the bare one at leading order in perturbation theory. On a finite torus of size l with twisted boundary conditions, this proportionality factor, denoted from now on by N (c(t), N ), acquires a volume and flow time dependence parameterized in terms of the dimensionless variable c(t) ≡ 8t/N l 2 , for further details see appendix B. We then define the finite volume quantity: which tends to the gradient flow coupling in the infinite volume limit, and, by analogy, a modified flowed energy density:Φ As we will see in subsection 6.1, finite volume (finite N in our set-up) effects are considerably reduced with this choice. Our scale definitions are then obtained by replacing Φ byΦ in eqs. (5.3) and (5.4), which at infinite volume coincide with the standard ones. Let us finally mention that the twisted finite volume normalization factor, N (c(t), N ), has been computed and used before for a series of step scaling studies [25,65,66]; explicit formulas are given in appendix B, where we also give some more details on the perturbative expansion and the relevant references. This procedure also accounts for correcting lattice artefacts at leading order if the normalization factor is computed instead in lattice perturbation theory. In subsection 6.1 we will discuss the specific lattice implementation of this method and the methodology applied to extract the scale from the implicit eqs. (5.3) and (5.4).

Wilson loops
We move now onto the discussion of the second alternative to fix the scale, based on Wilson loops. As already mentioned, a particular version of this procedure was used to set the scale in large N pure Yang-Mills theory [24]. Standard physical quantities such as the string tension or glueball masses are difficult observables to measure as they involve taking asymptotic limits. This makes them more affected by corrections and systematic uncertainties. It is in this spirit that Sommer scale [57], which is based on the quarkantiquark potential but not at infinite separation, is superior. However, the Sommer scale still involves the study of loops that are asymptotically long in time. Our proposal is based on fixed aspect ratio Wilson loops and hence it involves no limit.
Ultimately all gluon observables are based on Wilson loop expectation values. However, the Wilson loops themselves are UV divergent quantities and thus not suitable observables.
Our proposal is based upon the following observable function where W(r, r ) are expectation values for rectangular r × r Wilson loops. The double derivative gets rid of perimeter and corner singularities and one gets a well-defined continuum quantity depending on two scales. We can reduce it to a single scale by fixing the aspect ratio r /r. Different choices give different definitions. Notice that Sommer scale involves the limit r r . Here we will restrict ourselves to the opposite limit given by symmetric loops r = r (a restriction taken after the derivative is evaluated).
Indeed we claim that F (r, r) is a very interesting physical quantity for SU(N ) Yang-Mills theory, which has been computed in ref. [24] both for finite and infinite N . Obviously the string tension is given by However, it is better to fix the scale in a different way. For that purpose one notices that F (r, r) has dimensions of energy square, so we consider the dimensionless observable G(r) = r 2 F (r, r). A physical scaler(f 0 ) can be defined,á-la Sommer, as follows where f 0 is some numerical value that can be chosen arbitrarily. This is essentially the method used in ref. [24] to fix the scale with f 0 = 1.65. Here we will consider a variant of the method that makes use of the gradient flow and results advantageous from the point of view of the lattice implementation. One can use the gradient flow to construct flow-time dependent Wilson loops W t (r, r ) and use them to define a function F t (r, r), in analogy to eq. (5.8). Our new proposal is to use the flowed functions at non-zero flow time to set the scale. In order to define a function of a single energy scale in the problem, we fix the flow smearing radius to be proportional to the loop size: √ 8t = 8t(r, s) ≡ sr. This amounts to defining Wilson loops having fat edges with thickness proportional to the size of the loop. This choice defines a dimensionless function of a single length r G(r; s) = r 2 Ft (r,s) (r, r). (5.11) Applying now the, by now well-known, prescription we can define a physical scaler(f 0 , s) as followsĜ Any choice of f 0 and s defines a different scale, but they should all be proportional to each other. The concrete choice of s and f 0 is dictated by practical reasons of accessibility and insensitivity of the corresponding lattice observable to different error sources. From that viewpoint the choice used previouslyr(1.65, 0) is not the most adequate here.

Fundamental meson spectrum
In the large N limit, fermions in the fundamental representation play a completely different role than adjoint ones. If the limit is takená-la 't Hooft, sending N to infinity while keeping fixed the number of fermion flavours N f , fundamental-quark loops are suppressed and the quenched approximation is exact. As already mentioned, these fundamental quarks give rise to a meson spectrum, that can be used to determine an additional scale for the N = 1 SUSY Yang-Mills theory at large N . Even though the fundamental spectrum introduces an additional scale, the fundamental fermion mass, this scale is removed in the chiral limit (not to be confused with the massless gluino limit described in previous sections). It is in this particular limit when the masses of non-Goldstone fundamental mesons can be considered a natural scale of the SUSY theory. For the purpose of this paper we have selected the lightest vector meson state to determine this additional scale. In subsection 6.3 we will present our results for the determination of the lattice spacing using this method.

Scale setting: lattice determination of the scale
In this section, we are presenting the results we obtained by implementing on the lattice the three different scale setting methods described in the previous one.

Setting the scale with the flow
We will firstly present our results for determining the scale using the improved version of the flow discussed in subsection 5.1. We start by addressing the lattice implementation and by analyzing its efficiency in reducing finite size effects; as we will see the method works remarkably well inside a properly defined scaling window. The final analysis leading to the determination of the lattice spacing a(b, κ a ) in units of the flow scale is presented in subsection 6.1.2.

Methodology
So far our discussion of the flow observables in subsection 5.1 has been in the continuum. A generalization to the lattice is rather straightforward, one has to select a discretization of the energy density and of the flow equations. We adopt the clover version of the field strength which, in our one-site reduced lattice, leads to a discretized version of eq. (5.2) given by: (6.1) As for the flow, the flow time is discretized in units of the lattice spacing as: t = T a 2 , where we will from now on use capital letters to denote lattice, dimensionless quantities. We employ the so-called Wilson flow [58] and integrate the discretized flow equations by using a 3rd order Runge-Kutta integrator with constant time interval ∆T = 0.03125. We have checked that changing the time step does not produce any sizeable difference in the  In terms ofÊ(T ), the naive dimensionless flowed energy density can be estimated on the lattice from: The lattice equivalent of the normalization factor N (c(t), N ), corresponding to a concrete discretized definition of the energy density, is also easily determined by a tree level calculation in lattice perturbation theory. For our choice of the clover discretization and for the one-point lattice, an explicit expression has been computed in ref. [25] and is provided in Appendix B. One advantage of using the lattice determined instead of the continuum norm is that one corrects finite lattice artefacts on top of finite size effects at tree level. With all this, our final formula for the discretized version of the flow is given by: with N L (x, N ) given by eq. (B.10).
In order to illustrate the kind of improvement attained with the use of eq. (6.3) we present in figure 5a the dependence on T of the naive expression for three different values of N , corresponding to N = 169, 289 and 361. As can be observed, finite N corrections are a sizable source of systematics. The plot of figure 5b displays instead the flowed energy density obtained using the improved observableΦ. In this case, the window in which the three curves collapse to a single one extends over a much larger range. Generically, this window is set by the ratio c(t) = 8T /N , determining the fraction of the box occupied by the smearing radius. An empirical observation is that the correction is efficient over a scaling window given by: where the lower end of the interval is set to have under control the remaining lattice artefacts.

Results
Before presenting our results, there is an additional consideration that has to be made.
To have a precise determination of the scale it is important that the value of T s ≡ t s /a appearing in the lattice counterparts of eqs. (5.3) and (5.4) falls well within the scaling window given by eq. (6.4) and can therefore be reached by interpolation. With our set of parameters this is best attained by choosing s = 0.05 instead of the commonly used value of s = 0.1; the corresponding lattice scale will be denoted by T 1 from now on. However, on our smaller lattices, particularly those corresponding to b = 0.36, even this value leads to a scale that can only be reached by extrapolation. We address this issue by simultaneously fitting our data, for all values of b and hopping parameter κ a , to a single universal curve depending on T /T 1 (b, N, κ a ). Our combined data covers a window that runs from the perturbative small flow-time region up to values around T 0 /T 1 . In playing this game we are neglecting violations of universality that may come from lattice artefacts, remnant finite size effects or a dependence of the flow on the gluino mass. As we will see below, these assumptions are well satisfied by our results.
In appendix C we give full details of the form of the universal fitting functional used to describe the flow, here let us just indicate that we have selected it so as to describe appropriately the perturbative domain for small values of T /T 1 . For that one can use that the infinite volume flow defines a renormalized coupling constant and the small T /T 1 regime is therefore well described by using the two-loop perturbative β-function. As we will see below, this two-loop function describes very well the time dependence of the flow in a remarkably large window of flow times.
To serve the purpose of illustrating how well universality holds, we display in figure 6 the dependence ofΦ(T ) on T /T 1 for N = 289. The different data displayed in the plot correspond to different values of b and κ a restricted to the scaling window eq. (6.4) with γ = 0.28. The values of T 1 have been obtained from a universal fit to the data as the one described by eq. (C.3) with three parameters of the β-function, in addition to the two universal ones corresponding to N = 1 Supersymmetric Yang-Mills, and gives a χ 2 per degree of freedom of 1.2 (we obtain χ 2 /#dof = 3.1 for the datasets with N = 361). As becomes evident from the plot, the advantage of the joint fit is that it allows to constrain the time dependence of the flow in a region of scales much larger than the actual fitting window and permits to determine T 1 even when it falls outside it. Finally, we also display for comparison the prediction of two loop perturbation theory at infinite N , given by the black line in the plot, which, as mentioned previously, describes quite well our results in a large window of flow times.  Let us now present our results. In order to have an additional check on the effective reduction of finite N effects, we have fitted separately the N = 289 and 361 data sets, obtaining compatible results within errors. Our final values for T 1 are given in table 5, the first quoted error is statistical and the second systematic. The latter is determined so as to cover various determinations of the scale corresponding to different fitting functional and ranges as detailed in appendix C.
We end this section by discussing the relation of the scale t 1 to those more common in the literature, as t 0 or w 0 . We have already mentioned that in most of our simulations t 0 falls out of the scaling window. Nevertheless, and under the assumption of scaling, we can use the universal fitting functional describing the flow to obtain a determination of the ratio R = T 0 /T 1 . The result is R = 1.624(50) and R = 1.631 (70) for N = 361 and 289 respectively. This is in perfect agreement with the ratios obtained at the few cases where we can determine T 0 directly by interpolation. The error quoted in all cases covers for the   systematics in the fitting functional and fitting ranges following the same procedure used to determine t 1 .
Finally, we have also determined the scales w 1 and w 0 derived by solving the implicit equation eq. (5.4) with s = 0.05 and s = 0.1 respectively. The strategy to determine these scales is very similar to the one used for t 1 . We rely on the universality of the flow and fit tdφ(t)/dt as a function of t/a 2 simultaneously for all our datasets within the scaling window corresponding to γ = 0.28. In this case we use a degree-two polynomial fit, with the systematic error obtained by varying the fitting range. The resulting scales can be compared to t 1 . Excluding the sets at b = 0.36, which have very large systematic errors and give nevertheless results consistent within errors, and restricting to the cases with m pcac √ 8t 1 > 0.3, a fit of the dimensionless ratio of scales to a constant gives: w 1 / √ 8t 1 = 0.4535 (49) and w 0 / √ 8t 1 = 0.586(10) with χ 2 per degree of freedom equal to 0.24 and 0.15 respectively.
We collect our final results for the ratio of √ 8t 0 , w 0 and w 1 to √ 8t 1 in table 3. These ratios can be used to convert all the results given in this section, in particular the values of the lattice spacing as a function of the bare coupling and gluino mass, to the other more standard units used in the literature.

Setting the scale with Wilson loops (Creutz ratios)
In this subsection we focus on the second class of observables chosen to set the scale. These are the logarithm of the Wilson loop expectation values or rather its derivative F (r, l). As usual we have to look for lattice counterparts of these observables. These turn to be very well-known lattice quantities, the Creutz ratios, defined as follows: Notice that the first term is universal, being independent of the lattice bare coupling, once r is measured in units of an implicitly definedr. A practical problem which appears when implementing this method is that Creutz ratios χ(R, R ) are very noisy quantities for R and R large. To solve this problem one can use, as mentioned in the previous section, the gradient flowed (or APE smeared) equivalents of these observables: W t (r , r), F t (r, r) and G t (r). The lattice counterparts employ the same methodology as in the previous subsection allowing the definition of a lattice equivalent toĜ(r; s) in eq. (5.11):Ĝ L (R; s). In this case, the flow time in lattice units is fixed in terms of R as follows: √ 8T = sR. Although, this function is well-defined in the limit s → 0, performing this extrapolation as in ref. [24] is unnecessary. For an effective statistical error reduction it is enough to keep s larger than 0.1 and smaller than 1. Notice that intuitively our lattice observables are square Creutz ratios obtained from fat links of thickness proportional to its edge length.
Having defined the lattice observables to be used, we will now describe the process leading to the determination of the lattice spacing in units ofr(f 0 , s), c.f. eq. (5.12). This involves a collection of technical steps that we list below.
• For each configuration of our simulation parameters (b, κ a , N ) we evolve the lattice link variables using the same discretized flow that was explained earlier.
• At each flow time T we compute the Wilson loops and Creutz ratios for R = R =1.5, 2.5, 3.5, 4.5 and 5.5, and determine from them the corresponding value ofĜ L (R; s) and its error by averaging over configurations having the same simulation parameters.
• We investigated the N dependence of these values for those cases in which we have at least two values of N . As expected the sensitivity depends on the quantity R/ √ N giving the ratio of the loop size to the effective lattice size. This is indeed what happens when computing these quantities in perturbation theory at finite volume. To keep the finite volume correction smaller than 1-2%, one should set R/ √ N ≤ 0.25. In pure Yang-Mills theory we could reach N = 841 and N = 1369 which allowed us to go up to R max = 7.5 and R max = 9.5 respectively. Here our largest value of R max = 5.5 gives ratios of 0.29, 0.32 and 0.49 for N = 361, 289 and 169 respectively, which are larger than 0.25. Indeed, forĜ L (5.5; s) the difference between the value at N = 289 and N = 361 can reach up to 10%. Thus, in order to process the results we first extrapolate to N = ∞ using the N = 289 and 361 data and assuming a 1/N 2 dependence as predicted by perturbation theory for large enough N . In practice this limits our determination of the scale to the cases in which there is N = 361 data available.
• To extract the value of a/r(f 0 , s) one should deal with two additional issues. The first is that the values of r obtained on the lattice are multiples of the lattice spacing.
Determining the valuer at whichĜ(r; s) = f 0 must be done by interpolation. The second is thatĜ L differs from the continuum functionĜ by terms of order a 2 as described in eq. (6.7). One can deal with both issues simultaneously by a method  that gives a more robust determination of the scale. It involves a simultaneous fit to all our data points with 2.5 ≤ R ≤ 5.5 by a function where α, γ and δ only depend on the value of s and f 0 , and 1/R = a(b, κ a )/r(f 0 , s) expresses the lattice spacing inr(f 0 , s) units. Thus, the fitted data contains 4 values of R for each of the 14 total simulation parameters, and the fit parameters are the 14 values ofR and the three additional parameters α, γ and δ. The rationale behind the parameterization is given by the flowed equivalent to eq. (6.7). The universal function Ft (r;s) (r, r) is well described by a second order polynomial in (r(f 0 , s)/r) 2 forced to be equal to f 0 forr(f 0 , s)/r = 1. This parameterization is inspired by the results at s = 0 in which α is given by the string tension σr 2 (f 0 , s) and γ by a Lüscher-type term. The parameter δ is introduced to account for the a 2 correction appearing in eq. (6.7).
The procedure can be performed for various values of s and f 0 and the results should be compatible up to a change in the unit. In particular we chose two values of f 0 (0.65 and 1) and two values of s (0.5 and 0.65) to check consistency. The results for different f 0 are perfectly compatible since they involve fitting the same data points. The data just predicts the ratior(1, 0.65)/r(0.65, 0.65) = 1.611 andr(1, 0.50)/r(0.65, 0.50) = 2.045. On the other hand a change in s involves data at different flow times so that the comparison serves to check independence of this choice. If we fit the ratio of scales to a constant we get perfect compatibility with a constant value ofr(1, 0.65)/r(1, 0.50) = 1.120 (6).
Finally, we will express our lattice spacing in units ofr(1, 0.65) which are the ones affected by smaller errors. The results are given in table 5. Notice that the final values come from global fit to the data which assumes scaling. Hence, the errors do include a part associated to the amount of scaling violation present in our data. A visual determination of how well our data satisfies scaling can be obtained by plotting the best fit to the continuum functionĜ(r; 0.65). This is given in figure 7. Together with the function we plot all our data points after subtraction of the lattice artefact δ term. The horizontal errors come from the errors in the determination of the scalesR. The overall agreement is very good.

Setting the scale with fundamental meson spectrum
In this section we will explain the most basic details of the lattice implementation of the determination of the meson spectrum of fundamental quarks in this theory. More technical aspects will be collected in appendix D.
Fundamental quarks are quenched if the number of flavours over the number of colours tends to zero so that meson masses can be considered observables of the gauge-gluino theory. From that point of view the methodology applied to the determination of the fundamentalmeson spectroscopy does not differ from the one applied for the large N pure gauge twisted reduced models in previous publications [26,55]. We refer the reader to these publications for derivations and a more detailed explanation. In any case, as mentioned previously, the main idea is quite simple. We let the fermions propagate in an extended lattice on the background of the reduced model gauge field. This background field is periodic up to a twist and repeats itself after a √ N periods in each direction. Effectively it amounts to having gauge fields defined on a periodic lattice of size ( √ N ) 4 . The fermion fields can live in a much bigger lattice including infinite. For fundamental quark fields it is natural to impose that they propagate on a lattice of the same size ( √ N ) 4 , which considerably simplifies the formulas. For practical reasons in computing correlators in time it is useful to duplicate the lattice temporal length L 0 = 2 √ N . To achieve our final goal we have first to compute the correlation functions of bilinear quark operators. Here we restrict to the pseudoscalar, axial and vector channels. The next step is to determine the mass of the lightest state having the corresponding quantum numbers. This is done for various values of the fundamental quark mass. Then we extrapolate those masses to the fundamental chiral limit. The pion should behave as a Goldstone boson and the PCAC relation should hold as we approach this limit. The procedure is repeated for all our gauge couplings and gluino masses, and the resulting vector meson mass in that chiral limit is precisely the observable that we use to set the scale of the theory. All these steps will be spelled out and the results presented in the following susubsections.

Methodology
The meson correlators in time of fundamental quarks are expectation values where O A and O B are gauge invariant bilinear quark operators projected to zero spatial momentum. As explained previously, for the purpose of this paper we restrict ourselves to pseudoscalar and vector-meson operators, as for example: n Ψ(n 0 , x) γ 5 Ψ(n 0 , n) and n Ψ(n 0 , n) γ µ Ψ(n 0 , n). However, there are infinitely many operators with the same quantum numbers and this is an essential advantage that is used by us and other researchers. In practice we will be considering spatially non-local operators obtained by applying two well-known algorithms: the Wuppertal smearing of fermion bilinears [67][68][69] and the three-dimensional APE smearing [70]. The explicit expressions have been used in previous papers and are recalled in appendix D. The advantage of using these operators is that their coupling to the lowest mass state in each channel is enhanced with respect to excited states. In any case, it is useful to consider linear combinations of these operators designed to optimize the coupling to the ground state in each channel relative to other low mass states. This variational procedure is by now quite standard and goes by the name GEVP [71,72]. As an example of the outcome of the method, we display in figure 8 the correlators of the optimal operators in the pseudoscalar and vector channels for the case b = 0.35, κ a = 0.1825, κ f = 0.1525, with κ f the fundamental hopping parameter. The time-dependence of these correlators shows a clear exponential decay from which the mass of the π-meson (pseudoscalar) and ρ-meson (vector) can be extracted. The numbers obtained and the χ 2 of the fits to an exponential are also displayed on the figure. We also show the signal that allows one to extract the fundamental PCAC mass m f pcac , defined as in eq. (4.3).
Summing up, the procedure to obtain the final results is the one described in the example. It amounts to determining the optimal operator and then to fit the corresponding correlator in a certain interval to an exponential (rather to a hyperbolic cosine function, taking into account the periodic nature of the temporal direction). There are of course some specific details which reflect the selection of the operator and the choice of fitting interval that affect the final numerical value for the mass. In the end, variations of these types which are of similar statistical significance are accounted for as a systematic error of the determination. The specific details of the GEVP method that we have used will be collected in appendix D. Here we will comment briefly on the choice of the fitting interval. The main points to be taken into account in this selection are finite volume effects, lattice artefacts and the contamination of excited states. Typically, finite-size effects are more relevant close to the chiral limit, as the pion mass goes to zero and its Compton wavelength becomes comparable to the effective volume. To avoid too severe effects our selected values of the hopping parameter should stay sufficiently far from the chiral limit. In particular in our data the lightest cases still had m π a √ N ∼ 3.4. Nonetheless, finite-size effects may still reflect in the appearance of a constant term in the correlator arising from the propagation of quarks wrapping around the finite extent of the lattice [73]. Although this effect disappears in the large volume (large N ) limit, we have observed that in some cases the addition of a small constant to the hyperbolic cosine was required to obtain a good fit. Finally, we should comment about the lower limit of the fitting interval. A smaller time separation implies a smaller relative error at the expense of a larger contamination of excited states. A balance is then necessary, our choice has been to fix the lower limit of the interval in physical units setting it equal to R √ 8t 1 . We have taken in all cases R = 0.95 except for b = 0.36 for which, to increase the signal-to-noise ratio, a value R = 0.8 was preferred. As for the upper limit of the interval, it was set in all cases so as to have at least 6 lattice points in the fitting range.  The values of m f pcac , m π and m ρ for fermions in the fundamental representation and for each (b, κ a , κ f ) value at N = 289 are given. The fitting procedure is the one described in sec. 6.3.1, from which we obtained χ 2 per degrees of freedom in general smaller than one.  Figure 9: Extrapolation of am f pcac to 0. In the labels we report the corresponding value of κ (c) f extracted from the fit and the corresponding χ 2 per degrees of freedom. Errors are calculated using standard jackknife techniques.

Results
The results we obtained for the PCAC mass, the pion mass and the vector meson mass of fundamental fermions are reported in lattice units in table 4. We use these results to explore the (fundamental) chiral limit of the theory and extract an alternative scale to the ones presented in subsections 6.1 and 6.2.
The first check we do is to analyze the dependence of the PCAC mass on the fundamental hopping parameter κ f . In order to determine the critical hopping parameter where the fundamental fermions become massless we follow the same strategy as for gluinos, i.e. we analyze the dependence of the PCAC mass on κ f and determine the critical hopping parameter as the point where the fundamental PCAC mass vanishes. In figure 9 we plot am f pcac as a function of 1/(2κ f ). Performing separate linear fits for each value of (b, κ a ), we  Table 5: Values of the inverse lattice spacing in units of the three different scales determined in this work. For flow related scales, the first error is statistical, the second is systematic.
extracted in each case the critical value of the hopping parameter κ (c) f . As signalled by the χ 2 per degree of freedom reported in the legend of the plot, the observed linear dependence is very good and confirms what one would expect from chiral symmetry restoration in the limit of vanishing quark masses for Wilson fermions.
We finally come to the determination of the scale based on the vector meson mass. To determine the chiral extrapolation of this quantity we perform a global fit of am ρ as a linear function of am f pcac for each value of (b, κ a ), by imposing a common slope and extracting am χ ρ (b, κ a ) as the intercept at vanishing PCAC mass. We obtain a χ 2 per degree of freedom of 0.56.
The final values of the inverse lattice spacing in units of m χ ρ obtained in this way are reported in our summary table 5. As we will discuss in the next section, the ratio between this scale and √ 8t 1 turns out to be quite close to one. 8t 1 . The red bands represent the average weighted over the errors, while their width represent the statistical uncertainty over the average ratio. When needed, points corresponding to the same κ a have been slightly shifted in the x-axes to avoid overlapping.

Scale comparison
We dedicate this last part of the section to the comparison of the results obtained with the three different scale setting methods which are summarized in table 5. In each method we have determined the lattice spacing in terms of a different unit: a/ √ 8t 1 , a/r and am χ ρ . Although each of these quantities changes considerably when we change b (the inverse lattice 't Hooft coupling) and the hopping parameter κ a (related to the gluino mass), scaling dictates that the ratio should stay constant and be given by the ratio of the corresponding units of energy. In figure 10 we display the two independent ratios √ 8t 1 /r(1, 0.65) and systematic ones. The results are compatible with being a constant within errors. From the best fit we estimate that the conversion factor between the two units √ 8t 1 andr(1, 0.65) is 0.8708 (54), while between √ 8t 1 and 1/m χ ρ it is 1.144 (23). A more visual impression of how scaling works can be seen in figure 11 where the three lattice spacing determinations of the inverse of the lattice spacing are displayed side by side after applying the conversion factors determined earlier. All datasets for whichr(1, 0.65) was available are displayed. The figure shows how the three scales change considerably within all the datasets following the same trend in a consistent way. One can also notice the relative size of the errors of the three determinations of the lattice spacing. It looks as if the scale based on the Creutz ratios is the most precise, but the N dependence had to be corrected for and systematic errors might be underestimated. Furthermore, our flow-based scale is the one that has been determined for all our datasets and hence, it will be used for the determination of the scale in the Supersymmetric limit to be presented in the next section.
Concerning the rho mass, although also consistent, it is much less precise than the other two. Indeed, the scaling behaviour also holds for the meson masses built from massive fundamental quarks. This can be seen in figure 12 in which the rho mass in 1/ √ 8t 1 units is plotted against the fundamental PCAC mass in the same units. Although with large errors all the data points seem to follow the same trend.  After having completed the determination of the scale with good precision for our massive gluino values, we will here attempt achieving our main goal of determining the scale for the supersymmetric theory. This will be done by extrapolation of the scale to the masslessgluino limit. Given that we have used different scale setting observables and units which are mutually compatible, we will here focus on the flow unit √ 8t 1 since it covers all our simulation points and is both precise and relatively insensitive to finite N corrections. Hence, we extrapolate the lattice spacing a in those units to the limit of vanishing m pcac always expressed in physical units. The resulting plot is shown in figure 13. In the plot, different markers were used to display the points belonging to N = 289 and N = 361. By looking at the plot, it is visible by eye that the points corresponding to different markers are compatible within errors, showing that we were able to control finite-volume effects. The points are well fitted by a straight line and no higher-order polynomial terms are necessary to perform the extrapolation, which is remarkable, taking into account the wide range of m pcac values covered.
The values of the scale extrapolated to the massless-gluino limit are reported in table 6   For comparison with other authors we also convert to √ 8t 0 and w 0 units using the conversion factor determined in the previous section table 3. The last column displays the plaquette expectation value extrapolated to the massless gluino limit.
in physical units.
Our results provide the value of the lattice spacing as a function of the inverse lattice 't Hooft coupling λ. This is precisely the dependence that follows from the β-function of the theory as given by the following formula: where Λ is just an integration constant. Thus, by assuming some functional form for the β-function, integrating it and comparing it with our data values we can determine the parameters appearing in this functional form. This β-function is the one describing the bare coupling and depends on the discretization procedure, and on the definition of this bare coupling. However, perturbation theory predicts the leading behaviour of the function up to next to leading order. For Yang-Mills theory coupled to N f flavours of adjoint Dirac fermions this gives with the coefficients b 0 and b 1 given by There is one particular scheme in which the β-function is known to all orders. This is the so-called Novikov-Shifman-Vainshtein-Zakharov (NSVZ) β-function [74]: This functional form is particularly well-suited for performing the integration of the inverse of the β-function which is what is needed for eq. (7.1). One gets Thus, although we are certainly not in NSVZ scheme, this formula incorporates nicely the universal part of the β-function and provides a suitable parametrization which allows adding extra terms proportional to the coupling and to higher powers of it. Now let us apply these ideas to our data. It is well-known that for lattice QCD the naive lattice coupling 1 /b is not particularly well-suited for comparison with the perturbative predictions of scaling. Different authors have proposed improved couplings that do behave much better in this respect. We will then consider these same definitions for the Supersymmetric Yang-Mills theory. One possible choice is the one given in ref. [75]: where P (b) represents the average value of the plaquette extrapolated to the massless-gluino limit, whose values are reported in  Figure 14: Dependence of the logarithm of the lattice spacing as a function of the improved coupling λ I defined in eq. (7.7) (left panel), and λ I defined in (7.8) (right panel). Points have been fitted with eq. (7.6) leaving only Λ free (solid orange line), and with the same analytical form leaving also N f as a free parameter (dashed green line).
fitted parameter is ( √ 8t 1 Λ) −1 = 45.2 (7). If we modify the fit to include a higher order term in the β-function the fit gives a worse χ 2 /#dof and the additional parameter comes out compatible with zero. We conclude that our data does not have enough sensitivity to determine modifications to the NSVZ β-function. However, our data does have sensitivity to the leading coefficients of the β-function. A two-parameter fit leaving N f and Λ free having χ 2 /#dof = 1.38, gives N f = 0.31(20) (24π 2 b 0 = 9.76 (80)) to be compared with the perturbative result N f = 0.5 (24π 2 b 0 = 11 − 4N f = 9).
One can repeat the procedure with another choice of improved coupling constant, like the one proposed in refs. [76,77]: The one-parameter fit to the NSVZ β-function gives χ 2 /#dof = 1.42 with a best fit ( √ 8t 1 Λ) −1 = 472 (7). Again a two-parameter fit leaving also N f free gives N f = 0.30(22) (24π 2 b 0 = 9.78(87)), completely compatible with that obtained for the other improved coupling. In figure 14 we display side-by-side the logarithm of the lattice spacing as a function of both improved couplings together with the lines corresponding the fits described before.
In summary, we emphasize that the behaviour of the scale in the range explored in our study is certainly not far and even compatible with the dependence predicted by perturbation theory. The data also shows the tendency expected by the addition of an adjoint Majorana fermion since N f comes out larger than the value 0 corresponding to the pure gauge theory and not incompatible with the value 0.5, expected for the SUSY Yang-Mills theory at asymptotic small values of the coupling. This implies that the leading perturbative coefficient agrees within a 10% with the expected value b 0 = 9/(24π 2 ).

Summary
The purpose of this paper was to study the large N limit of the N = 1 SUSY Yang-Mills theory on the lattice. More precisely our goal was the determination of the dependence of the lattice spacing on the bare coupling constant, a necessary first step in connecting lattice to continuum results. The paper includes dealing with many aspects of a technical nature, although some of them have been moved to appendices to facilitate readability. Also with this idea in mind we consider appropriate to add this short final summary to point to our main results.
We have used redundancy to ensure the robustness of our final conclusions. Hence, we have computed the lattice spacing using three different and rather independent observables. First we did this for all our different simulation runs having various values of the lattice coupling and the gluino mass. The results are collected in table 5. Since each observable expresses the lattice spacing in its associated unit, compatibility demands that the values are proportional to each other with a proportionality constant given by the ratio of units. This is actually happening as can be deduced from the table or more visually from figure 11. The proportionality is in principle only exact in the continuum limit, but our results show that this holds for our lattice sample within errors.
The last step is to extrapolate the results to the massless gluino limit giving the desired value of the lattice spacing for each bare lattice coupling (table 6). Remarkably, the resulting dependence of the beta function for two common improved lattice couplings comes close to the prediction of perturbation theory ( figure 14). This also shows that the range of bare couplings that we have explored is within the precocious scaling regime. hp190004) and Subsystem B of ITO system (Kyushu U.). We acknowledge the use of the Hydra cluster at IFT and HPC resources at CESGA (Supercomputing Centre of Galicia).

A Simulation Parameters and Statistics
In this appendix we collect extra information regarding the simulation in addition to the one described in section 3. In particular, the full set of model parameters and the number of configurations generated for them are collected in table 7. The hyper-parameters of the RHMC algorithm and the statistics for each ensemble are also included in the table. Each configuration is separated with five trajectories with a trajectory length of τ = 1. The number of the MD time steps in the generalized multiple-time step integrator is denoted by N G for the gauge action, and N UV (N IR ) for the UV (IR) part of the pseudo-fermion action, respectively. Table 8 shows the HMC acceptance rate, the expectation value of the plaquette, and the expectation value |λ min |, with λ 2 min the minimum eigenvalue of The choice of hopping parameters have been performed having in mind the ultimate goal of taking the massless limit. Thus, our κ a values approach the critical value from below. On the opposite edge, for small hopping parameter the theory becomes equivalent to the TEK model. This model exhibits a first order phase transition at a certain b c ∼ 0.346 for (N, k) = (169, 5) and below this point one enters a strong coupling region of the model. Indeed, this phase transition extends into a first order phase transition line in the current model in the b-κ a plane. Hence, we first surveyed the phase diagram for the weak coupling region from which the proper continuum limit can be approached, and the model parameters in the table are chosen to stay in the weak coupling region.
The RHMC algorithm program was primarily written in Fortran 2003/2008 language and uses OpenMP thread parallelization. To speed up the generation of configurations for the cases with N = 361, the pseudo-fermion force computation in the MD evolution is offloaded to the corresponding GPU kernel written in the CUDA language. To further speed up the configuration generation we made several replica with the same parameter set but with different random number sequence so that several computational resources are concurrently available on multi-node systems. We discarded the first 200-500 trajectories in each replica for the thermalization, where the seed configurations are chosen in such a way that the system remains in the weak coupling region avoiding the phase transition.    N , b and κ a we list the HMC acceptance rate, the expectation value of the plaquette P , and |λ min | , with λ 2 min the lowest eigenvalue of Q 2 W = (D W γ 5 ) 2 .

B Finite size effects in the gradient flow
In this appendix we summarize the basic ingredients for the derivation of the improved expression of the flow presented in eq. (6.3) and provide also explicit formulas for computing the normalization constant N L (x, N ) on the lattice. The perturbative expansion of the infinite volume flow in terms of the MS 't Hooft coupling at scale µ = 1/ √ 8t is given by [58]: where: where N stands for the number of colours. This expression is used to define the gradient flow (GF) coupling constant [58]: showing that c 1 stands for the finite one-loop renormalization constant that relates the Λ parameters in the two schemes. One can derive an analogous expression for the finite volume flow. In our case, on a box of size l 4 with twisted boundary conditions as the ones used in this work, one obtains at second order in the coupling [78]: where finite volume effects depend on the dimensionless variable: and are encoded in the functions N (c(t), N ) and C(c(t), N ). The former has a simple expression in terms of Jacobi θ 3 functions and reads: with: In the large volume limit, taken by sending l to infinity at fixed flow time, i.e. by sending c(t) to zero, this normalization factor tends to the infinite volume one, c.f. eq. (B.2). The function C(c(t), N ) goes to zero in this limit and the infinite volume expansion of the flow presented above is recovered. Let us also mention that, on account of volume independence, the same expressions can be obtained for SU (∞) by taking the c(t) to zero limit in a different way, i.e. by sending N to infinity at fixed torus size l (in the particular case of the one-point lattice, l = a). In that limit, the explicit N dependence of the normalization K(N ) factor disappears. This exercise indicates a simple way to correct the flow at leading order in the coupling. We introduce the quantityλ which, at leading order in the coupling, has the correct perturbative expansion and we define a modified flowed energy density given by: This removes finite size effects at tree level. At second order in λ, the remnant volume dependence comes from the c(t) dependence of the function C (c(t), N ), which has not yet been computed in perturbation theory for the case envisaged in this work (with a dynamical Majorana fermion and the symmetric twist). As discussed in section 6.1.2, within the numerical accuracy we have achieved, this correction is very small for values of c(t) 0.3. On the lattice, an analogous correction, which in addition takes into account lattice artefacts, is derived by computing the tree level factor N (c(t), N ) in lattice perturbation theory. For our choice of the clover discretization of the energy density and for the one-point lattice [25]: where the lattice momentum is given by q µ = 2 sin(q µ /2), with q µ taking values: for m µ = 0, · · · , √ N − 1. The prime in the sum excludes the zero momentum mode, with m µ = 0, ∀µ. With all this, our final formula for the discretized version of the flow is given by eq. (6.3).

C Fitting strategy to determine the flow-based scale
In this appendix we discuss the fitting strategy employed to determine the scales based on the gradient flow and present a detailed account of the extraction of T 1 .
As mentioned in section 6.1, our results show that one can use a universal fitting function to describe the flow-time dependence of the energy density, allowing to extract the scale even in the cases where T 1 can only be reached by extrapolation. The parameterization of the flow-time dependence we have used, to be described below, relies on the connection between the infinite volume flow and the gradient flow renormalized coupling constant λ gf , c.f. eq. (5.5) 3 .
Starting from the renormalization group equation defining the β function and integrating it between two different reference scales t s and t one arrives at: The left hand side of this equation can be easily integrated using a general parameterization of the β-function inspired by the Novikov-Shifman-Vainshtein-Zakharov (NSVZ) β-function [79]: with coefficients chosen to reproduce the universal expansion of the β-function to second order in λ, c.f. a 0 is set to b 1 /b 0 , where b 0 and b 1 stand for the first two universal coefficients of the β-function of the N = 1 SUSY Yang-Mills theory, c.f. eqs. (7.3), (7.4) with N f = 1/2. After integration, one arrives at the following identity: One can use this identity to fit the numerical results for the flow, leaving as free parameters the non-universal coefficients of the β-function and the t s scale. In the continuum, at infinite volume, and for one masless Majorana fermion, this relation is universal and dictated by the non-perturbative β-function of the N = 1 SUSY Yang-Mills theory in the gradient flow scheme. Departures from universality arise due to lattice artefacts, finite volume effects and the fermion mass-dependence of the flow. Our results, presented in section 6.1, indicate that violations of universality remain small, within the numerical accuracy we have been able to attain, as long as one stays within the scaling window given by eq. (6.4).
To implement the fitting procedure and determine the lattice scale T 1 (b, κ a , N ), we set λ(t s ) = 0.05/K(∞) and minimize the χ 2 function defined as:  (8)(144) 5.323 (7) Table 9: Values of the scale T 1 determined from the flow data with N = 289 and N = 361. Results labelled as T 1 come from a fit to eq. (C.3) with 3 free parameters of the β-function on top of the two universal ones, as described in the text. Values labelled as T 0 1 are directly determined by interpolation in the cases in which T 1 falls within the scaling window given by eq. (6.4) with γ = 0.28. where: in a procedure analogous to the one used to parameterize the step scaling function in ref. [80]. Our final determination of T 1 is obtained by fitting simultaneously in this way all our simulations at fixed value of N , all values of b and all values of the hopping parameter κ a . For the final fit, the β-function has been parameterized with 3 free coefficients, in addition to the two universal ones. Fits were performed restricting the data for different values of b and κ a to the corresponding scaling window T ∈ [1.25, γ 2 N/8], with γ = 0.28, as shown in figure 6.
Our final results are given in table 9. We give separately results of the scale determined from N = 289 and N = 361 simulations, whose compatibility serves as a test of the absence of finite size effects. In addition to the values obtained from the universal fitting function, we also provide the values T 0 1 obtained by interpolation in those cases in which the scale falls well within the scaling window. We have assigned a systematic error to the final result that covers for the difference between them as well as for other types of fits and fitting ranges. Among those, we have included fits to eq. (C.3) with only one and two free parameters and also a second-degree polynomial fit performed in a reduced fitting window corresponding to γ = 0.22. In addition, separate fits have also been performed, including joint fits to all datasets at fixed value of the bare coupling b. The final systematic error quoted in the table remains in general below a 3% relative error, going up to about 6% for the, in physical units, smaller lattices.

D Meson correlators in the reduced model
In this appendix we give technical details about computing meson correlators in both the fundamental and adjoint representation and their subsequent use to obtain the masses of the lowest state with the corresponding quantum numbers. The underlying ideas and derivations are given in the text and references supplied. For the case of fundamental representation fermions there have been extensive studies done previously [26,55] using the same techniques applied in the paper. For that reason we will mostly focus on the adjoint representation formulas.
The meson correlation functions are expectation values of products of fermion bilinears separated in time and averaged over space. The correlation function is computed as a Fourier transform in time where the trace is over spin, space and colour degrees of freedom. In this work we take the time period to be 2 √ N so that the temporal momentum takes values q 0 = πm/ √ N with integer m. The symbols A and B specify the spin-parity quantum numbers of the operator and the indices i and j run over a family of operators with the same quantum numbers. The symbol D −1 W (0) denotes the inverse of the Wilson-Dirac operator of the corresponding fermion (fundamental or adjoint). For the reduced model this inverse greatly simplifies. In the case of the adjoint the operator is simply the one appearing in eq. (2.12). The inversion is performed using the BiCGStab algorithm or the Conjugate Gradient algorithm whenever the former does not converge.
One can define D W (p) as the corresponding operator with the substitution If |Λ p |, the number of elements in Λ p , is large, this might imply many inversions. In practice, what we do is to perform this average stochastically: For each configuration we generate p randomly and use it to perform the inversion of D W (p). For the case of fundamental fermions things are very similar except for the expression of the Wilson-Dirac operator. We refer the reader to the literature for further details [26,55]. Now we have to specify the selection of operators O

(i)
A used to obtain the masses. This is based in applying Wuppertal smearing [67][68][69] to the bilinear operator combined with three-dimensional APE smearing to the link variables [70]. The particular application to the reduced model in the adjoint representation is explained below.

Wuppertal fermion smearing in adjoint representation
This amounts to replacing the fermion bilinear associated to an element of the Clifford algebra γ A as follows Ψγ A Ψ −→ ΨO In the case of adjoint fermion we only used up to 6 operators, while for fundamental ones we used up to 9. The symbolŪ k appearing in eq. (D.5) corresponds to the 10 times APE-3d smeared link to be explained below.

APE-3d link smearing
This smearing procedure is an iterative one which in our case maps three spatial SU(N ) matrices onto new ones: where P is an operator projecting on SU(N ). The starting point of the iteration are the link matrices U For the case of fundamental meson correlators we applied the same method using f = 0.15.

GEVP methodology
We expect the signal to decay in time as an (infinite) sum of exponentials, corresponding to the contribution of the ground state plus other heavier excited states. Given the typical hierarchy in the mass spectrum, we expect the excited states to decay faster than the ground state, whose signal dominates for large-enough time separation. We want to extract the mass of the ground state from an exponential fit in a region where excited states give no systematic contribution. In order to do so, one has to maximize the projection onto the ground state to have a single-exponential decay setting in early on. Fermion smearing in eq. (D.4) should provide an improvement in this regard since it increases the overlap onto the ground state wave function. Beyond that, we make use of a variational approach, in the same philosophy of ref. [26]. Using the definition of the correlator matrix given in eq. (D.1), where the index i and j run over the smearing levels, we applied the so-called GEVP method. Given two timeslices τ 0 and τ 1 (with τ 0 < τ 1 ) (not to be confused with gradient flow scales) we solve numerically the generalized eigenvalue problem (GEVP) where v (n) and λ (n) are the eigenvectors and the eigenvalues, respectively for a given choice of τ 0 and τ 1 (in order to lighten the notation, we also omitted the indices A and B that appeared in eq. (D.1)). In this work we used τ 0 = a and τ 1 = 2a. Among the basis composed by the eigenvectors v (n) , we choose the one whose corresponding eigenvalue is the biggest among the others. Let us denote this maximum eigenvector with v max , which we use to define an optimal operator by rotating the original correlator-matrix C opt (n 0 , τ 1 , τ 0 ) = v max i * C ij (n 0 , τ 1 , τ 0 )v max j . (D.9) The ground state mass is extracted from the exponential decay at large time of this correlator.