Eigenvalue spectrum and scaling dimension of lattice $\mathcal{N} = 4$ supersymmetric Yang-Mills

We investigate the lattice regularization of $\mathcal{N} = 4$ supersymmetric Yang-Mills theory, by stochastically computing the eigenvalue mode number of the fermion operator. This provides important insight into the non-perturbative renormalization group flow of the lattice theory, through the definition of a scale-dependent effective mass anomalous dimension. While this anomalous dimension is expected to vanish in the conformal continuum theory, the finite lattice volume and lattice spacing generically lead to non-zero values, which we use to study the approach to the continuum limit. Our numerical results, comparing multiple lattice volumes, 't Hooft couplings, and numbers of colors, confirm convergence towards the expected continuum result, while quantifying the increasing significance of lattice artifacts at larger couplings.


Introduction
Four-dimensional maximally supersymmetric Yang-Mills theory (N = 4 SYM) is widely studied in theoretical physics. Its many symmetries-including conformal symmetry and an SU(4) R-symmetry in addition to Q = 16 supersymmetries-make it arguably one of the simplest non-trivial quantum field theories in four dimensions, especially in the large-N c planar limit of its SU(N c ) gauge group [1]. This simplicity enabled its role as the conformal field theory of the original AdS/CFT holographic duality [2], provided early insight into S-duality [3], and continues to inform modern analyses of scattering amplitudes [4].
At the same time, the non-triviality of N = 4 SYM makes it important to explore the lattice regularization of the theory. In addition to providing in principle a non-perturbative definition of N = 4 SYM, lattice field theory is also a way to numerically predict its behavior from first principles, even at strong coupling and away from the planar limit. A prominent target for such predictions is the spectrum of conformal scaling dimensions, which depend on the 't Hooft coupling λ = N g 2 YM . Although lattice field theory has been very successfully used to analyze nonsupersymmetric vector-like gauge theories such as quantum chromodynamics (QCD), it has proven more challenging to apply this approach to supersymmetric systems. In large part this is because supersymmetry is explicitly broken by the lattice discretization of space-time. See Refs. [5][6][7] for reviews of these difficulties and the significant progress that has been achieved to overcome them in recent years.
In particular, using ideas borrowed from topological field theory and orbifold constructions, a lattice formulation of N = 4 SYM has been developed which preserves a closed supersymmetry subalgebra at non-zero lattice spacing a > 0 [5][6][7]. Although 15 of the 16 supersymmetries are still broken away from the a → 0 continuum limit, the single preserved supersymmetry significantly simplifies the lattice theory. These simplifications are sufficient to establish that at most a single marginal coupling may need to be tuned to correctly recover the full symmetries of N = 4 SYM in the continuum limit [8,9]. In addition, the moduli space of the lattice theory matches that of continuum N = 4 SYM to all orders in perturbation theory, and the renormalization group (RG) β function vanishes at one loop in lattice perturbation theory [10].
Of course, numerical lattice field theory calculations require both a non-zero lattice spacing that corresponds to a ultraviolet (UV) cutoff scale 1/a, as well as a finite lattice volume (L·a) 4 that introduces an effective infrared (IR) cutoff, explicitly breaking conformal scale invariance. It is also necessary to softly break the single preserved supersymmetry in order to regulate flat directions and make the lattice path integral well defined [7,11]. These facts make it challenging to analyze the approach to the a → 0 continuum limit in the strongly interacting regime where lattice perturbation theory is unreliable. 1 It is therefore essential to carry out detailed numerical studies of the non-perturbative RG properties of the lattice theory.
In this paper we present progress investigating RG properties of lattice N = 4 SYM. Specifically, we compute the eigenvalue mode number of the fermion operator, and use this to estimate the 'mass anomalous dimension' γ * (λ) that would appear in the scaling dimension of the corresponding fermion bilinear. In the continuum theory this anomalous dimension is expected to vanish, γ * = 0, for all values of the 't Hooft coupling λ. Computing γ * (λ) is hence a way to assess the effects of breaking supersymmetry and conformal symmetry in numerical lattice calculations, and verify that the properties of N = 4 SYM are correctly reproduced in the continuum limit. Ref. [12] presented a first preliminary investigation of this topic, and preliminary results from a similar project studying the anomalous dimension of the Konishi operator more recently appeared in Ref. [7]. While Ref. [12] proceeded by numerically computing the fermion operator eigenvalues, that approach quickly becomes inefficient as the lattice volume increases. Here we instead apply stochastic techniques to estimate the mode number.
After reviewing the basic features of lattice N = 4 SYM in the next section, we summarize these stochastic techniques in Section 3, also discussing how the resulting mode number provides information on the anomalous dimension. In Section 4 we consider the free (λ = 0) lattice theory, both to check our methods and to explore discretization artifacts. Our numerical results for both the low-lying eigenvalues and the stochastic mode number are presented in Section 5, and in Section 6 we use these to estimate the anomalous dimension. After looking more closely at the dependence of the results on the gauge group and lattice volume in Section 7, we conclude in Section 8 with some discussion of the next steps for lattice analyses of N = 4 SYM.
2 Lattice formulation of twisted N = 4 SYM As mentioned above, the lattice formulation of N = 4 SYM that we use has its origins in both topologically twisted [13][14][15] and orbifolded [16][17][18][19] approaches, which ultimately produce equivalent constructions [20,21]. Here we will use the twisted language, which organizes the Q = 16 supercharges of the theory into integer-spin representations of a twisted rotation group SO(4) tw ≡ diag [SO(4) euc ⊗ SO(4) R ], where SO(4) euc is the Lorentz group Wick-rotated to euclidean space-time and SO(4) R is a subgroup of the SU(4) R-symmetry. It is then convenient to combine these representations into 1-, 5-and 10-component sets Q, Q a and Q ab = −Q ba , respectively. These transform under the S 5 point-group symmetry of the A * 4 lattice we use to discretize space-time, which consists of five basis vectors symmetrically spanning four dimensions. The twisted-scalar supersymmetry is nilpotent, preserving the subalgebra {Q, Q} = 0 even at non-zero lattice spacing where Q a and Q ab are broken.
The lattice action, just like the continuum theory, is now the sum of the following Q-exact and Q-closed terms [8][9][10][11][22][23][24]: where n indexes the lattice sites and repeated indices are summed. In this paper we will present results in terms of the input lattice 't Hooft coupling λ lat , which differs slightly from the continuum λ [7,23]. The fermion fields η, ψ a and χ ab = −χ ba transform in the same way as the corresponding twisted supercharges, and are respectively associated with the lattice sites, links and oriented plaquettes. The gauge and scalar fields are combined into the five-component complexified gauge links U a and U a , which appear in the finite-difference operators D (+) a and D (−) a [20,21]. These complexified gauge links imply U(N ) = SU(N ) × U(1) gauge invariance, with flat directions in both the SU(N ) and U(1) sectors that need to be regulated in numerical calculations, as mentioned in Section 1. To achieve this, we work with the improved action introduced by Ref. [11], 2 which adds two deformations to Eq. (1). The first is a simple scalar potential with tunable parameter µ, which regulates the SU(N ) flat directions while softly breaking the Q supersymmetry. The second deformation is Q-exact, and replaces the term in S exact , with tunable parameter G. This deformation picks out the U(1) sector through the determinant of the plaquette oriented in the a-b plane, P ab (n), which is an N c × N c matrix at each lattice site n. Expanding the Q transformation produces terms that modify both the fermion operator and the bosonic action [11], as expected for a supersymmetric deformation. Using this improved action, we have generated many ensembles of field configurations using the rational hybrid Monte Carlo (RHMC) algorithm [26] implemented in parallel software that we make publicly available [24]. 3 In addition, we have modified this software to implement the stochastic estimation of the mode number discussed below. For this stochastic computation, it is convenient to rescale some of the fermion field components, which is irrelevant for the path integral since it introduces only a constant prefactor. This rescaling is done only for the measurement of the mode number, not yet in RHMC configuration generation.
The particular rescaling we perform is chosen to put the fermion operator into its most symmetric form, which simplifies analytic considerations and changes the degeneracies of eigenvalues. Reference [10] previously reported on the analytic structure of the lattice theory. In its conventions, which differ slightly from Eq. (1), the fermion operator D has the form where Ψ = (η, ψ a , χ ab ) collects the 16 fermion fields into a vector. We adjust this operator by reducing summations for χ ab to the relevant part over a < b, compensating a factor of 2 by rescaling χ ab = 2 χ ab . The same result up to an overall factor of two is obtained from Eq. (1) by rescaling η → η 2 . The more symmetric fermion operator D defined in this way is equivalent to the original operator in Ref. [18], and a rescaling was also done to discuss symmetries in Ref. [8]. In the free theory, the squared operator D † D is now block diagonal in momentum space, D † D ∼ f (p)I 16N 2 c , implying a 16N 2 c -fold degeneracy of the eigenvalues. The function f (p) on the A * 4 lattice is 2 There is ongoing exploration of alternative lattice actions that address this issue in different ways [25]. 3 github.com/daschaich/susy where p µ are any four of the five linearly dependent lattice momenta [10]. This 16N 2 c -fold degeneracy is lifted in the interacting theory, but for any 't Hooft coupling λ lat ≥ 0 the eigenvalues of the lattice fermion operator (with or without rescaling) occur in +/− pairs, so that the non-negative eigenvalues of the squared operator are always 2-fold degenerate.

The eigenvalue spectrum and stochastic estimation of the mode number
The mode number, which is the integrated eigenvalue density of the fermion operator, allows for a precise estimate of the mass anomalous dimension [27][28][29]. On the lattice the most practical definition is obtained from the spectral density of the massless fermion operator D, Here the eigenvalues λ k should not be confused with the 't Hooft coupling λ. The mode number ν(Ω) is defined to be the number of eigenvalues λ 2 k of the non-negative operator D † D that are smaller than Ω 2 : where ρ is the spectral density of D † D and the second equality follows from the eigenvalue pairing mentioned above. 4 Throughout the paper quantities like the eigenvalues λ k and the scale Ω are provided in lattice units. The anomalous dimension γ * governs the dependence of the mode number on the scale Ω 2 : Additional terms present in the Wilson-fermion case (see Ref. [27]) do not appear here. This makes it possible to define a scale-dependent effective anomalous dimension from any two values of the mode number: where Ω 2 ≡ (Ω 2 1 + Ω 2 2 )/2. In addition to depending on the choice of scales Ω 1 and Ω 2 , the determination of γ eff on any ensemble of lattice field configurations will be affected by lattice artifacts. As we discuss in Section 4, even the free theory with λ lat = 0 only recovers the continuum γ * = 0 after extrapolation to the continuum limit.
If Ω 1 and Ω 2 are close to each other and the lattice is coarse, the results of this naive method are quite unstable and fluctuate significantly. We will show below that fits provide more stable results. An alternative approach to improve stability is to normalize the mode number with respect to the free-theory ν free , using some fixed reference scale Ω 1 , Now that we have seen how an effective anomalous dimension can be extracted from the mode number, we review our stochastic estimation of ν(Ω) using the wellestablished projection method proposed in Ref. [30]. This method is based on a rational approximation of the projection operator P for eigenvalues in the region below a given threshold. The mode number is then where the trace is obtained by stochastic estimation. The projection operator is approximated in terms of the step function h(x) using The parameter Ω * ≈ Ω is adjusted to minimize the error of the approximation-see Ref. [30] for further details. More recently, a different method based on a Chebyshev expansion of the spectral density ρ has been proposed [31,32]. In this Chebyshev expansion method, the spectrum is rescaled to the interval [−1, 1] by defining where λ 2 max and λ 2 min are the maximal and minimal eigenvalues of D † D. We consider the integral of the spectral density ρ M of the rescaled operator M multiplied by the nth term T n of a Chebyshev polynomial of order N p , which we estimate stochastically using N s random Z 4 noise vectors v l : Based on the orthogonality relations for the T n , the spectral density ρ M is now approximated by The spectral density of D † D is then obtained by mapping the interval [−1, 1] back to the original eigenvalue region [λ 2 min , λ 2 max ], and can be integrated analytically to provide the mode number via Eq. (7). The work we present below will focus on measurements of the mode number using the Chebyshev expansion method. We use polynomials of order 5,000 ≤ N p ≤ 10,000 depending on the spectral range [λ 2 min , λ 2 max ], which increases for stronger 't Hooft couplings. We cross-checked these results through the more computationally expensive projection method of Eq. (11), as well as by directly computing the lowlying eigenvalues of D † D using a Davidson-type method provided by the PReconditioned Iterative Multi-Method Eigensolver (PRIMME) library [33]. In addition to checking the stochastic results for small Ω, these direct eigenvalue measurements can provide an alternative estimate of the effective anomalous dimension, from the volume-scaling relation [12] where y k = 2/(1 + γ eff ) and λ 2 k is the average kth eigenvalue of D † D. We will see below that the mode number provides more precise results than the individual eigenvalues, as expected [28].

Discretization effects for the free lattice theory
Before presenting numerical results obtained from analyzing the available ensembles of lattice N = 4 SYM field configurations, we consider the λ lat = 0 free theory to test both our methods of stochastically estimating the mode number, as well as our extraction of γ eff . By definition, all anomalous dimensions vanish for the free theory, meaning that any non-zero results we obtain will provide information about the lattice artifacts we want to explore. The free theory is simple enough that we are able to analytically compute the mode number on the A * 4 space-time lattice. The free D † D operator has an enhanced symmetry, resulting in larger degeneracies of the eigenvalues, especially when all fields are subject to periodic boundary conditions (BCs) in all four dimensions. In RHMC configuration generation, and in   to lift a fermion zero mode. Even with those antiperiodic BCs, many free-theory eigenvalues are degenerate, leading the mode number to exhibit distinct steps at certain values of Ω. This stepwise behavior is typically quite difficult to capture with a polynomial approximation, but Figure 1 shows that the Chebyshev approach is able to provide reasonable precision.
While the effective anomalous dimension γ eff can be estimated from any two values of the mode number using Eq. (9), this naive approach produces large fluctuations, especially for small Ω where the stepwise behavior of the mode number is most prominent. More stable results are obtained by fitting the mode number according to Eq. (8) over a window [Ω 2 , Ω 2 + ] of length , which we show in Figure 2. This figure compares = 1 fit results for the free-theory γ eff from two different lattice volumes 8 4 and 16 4 , considering both periodic and antiperiodic BCs for the latter. As expected, we obtain more stable results as we approach the continuum limit by increasing the lattice volume. We also see larger oscillations for the case of periodic BCs, which is due to the larger degeneracies of the eigenvalues in this case.
While the IR limit lim Ω→0 γ eff ≈ 0 for the free-theory results shown in Figure 2, as the scale increases up to Ω 7, the effective anomalous dimension tends towards negative values. In Figure 3 we confirm that this trend is due to lattice artifacts, by comparing three dispersion relations for the free operator in momentum space. For the A * 4 lattice with periodic BCs, this is based on the function f (p) in Eq. (5). Considering instead a continuum-like dispersion relation p 2 produces effective anomalous dimension results more consistent with the true γ * = 0, while the hypercubic-lattice dispersion relation 4 sin 2 (p/2) leads to the same qualitative trends in γ eff as the A * 4 lattice used for N = 4 SYM.
From these investigations of the free theory, we can conclude that the Chebyshev method provides a reasonable approximation of the mode number. Even in this free case where the eigenvalue degeneracies are largest, the stepwise behavior of the mode number is accurately resolved. The resulting oscillations in results for the effective anomalous dimension γ eff are manageably small when we fit the mode number over a window of adequate length , and are expected to decrease for the λ lat > 0 to which we now turn.
In addition to checking our methods, we can also use this consideration of the free theory to improve our main numerical analyses. The ensembles of lattice N = 4 SYM field configurations we analyze span a range of 't Hooft couplings 0.25 ≤ λ lat ≤ 2.5 within which some similarities to the free theory persist. Therefore we can employ the method of Eq. (10) as an alternative approach to improve the stability of the results by using the free theory as a reference for the scaling of the mode number.

Results for eigenvalues and mode number
Our main numerical analyses involve the 18 ensembles of lattice N = 4 SYM field configurations listed in Table 1, which are a subset of a broader collection of ensembles being used to investigate other aspects of the theory [7]. These were generated using the RHMC algorithm and the improved lattice action described in Section 2. For gauge group U(2) we consider L 4 lattice volumes with 10 ≤ L ≤ 16, and lattice 't Hooft couplings 0.25 ≤ λ lat ≤ 2.5. In addition we also analyze more limited data sets for gauge groups U(3) and U(4), which involve significantly larger computational costs.
As part of the process of configuration generation, we use PRIMME [33] to measure the extremal eigenvalues of D † D on every saved configuration, ensuring that the minimum and maximum across the entire ensemble remain within the spectral range where the RHMC rational approximation is accurate. This information is included in Table 1 Table 1: Summary of the ensembles we consider, with gauge group U(N c ) and lattice volume L 4 . For all ensembles G = 0.05 and the fermion fields are subject to antiperiodic BCs in the time direction. We set µ ≈ √ 5λ lat /L to remove the scalar potential in Eq. (2) in the L → ∞ continuum limit with fixed lattice 't Hooft coupling λ lat . For each ensemble we report the extremal eigenvalues λ 2 of D † D, ensuring that they remain within the spectral range where the corresponding RHMC rational approximation is accurate. The listed number of mode number measurements are carried out on thermalized configurations separated by 10 molecular dynamics time units. From measurements of λ 2 min on those same configurations we also estimate autocorrelation times τ in units of measurements. D † D eigenvalue following thermalization/equilibration, estimated using the 'autocorr' module in emcee [34]. This autocorrelation time provides an indication of how many statistically independent samples are available for each ensemble. Some additional information about these ensembles is collected in the Appendix. Following configuration generation and analysis of thermalization, we additionally compute the extremal eigenvalues of the more symmetric rescaled operator D † D on all thermalized configurations, each separated by 10 molecular dynamics time units. We also use all of these thermalized configurations to stochastically estimate the D † D mode number through the Chebyshev expansion method, with less-extensive cross-checks using the projection method, as described in Section 3. Table 1 reports the number of measured configurations for each ensemble.
The low-lying eigenvalues of D † D provide a first look at the mode number and its scaling. Figure 4 shows the general form of these low-lying eigenvalues, which features a pronounced gap between the lowest two pairs and the rest of the spectrum. This gap is rather stable across each RHMC Markov chain, and its relative size decreases for smaller λ lat due to the lowest two pairs moving to larger values. Similar observations were also made in Ref. [12] for a different N = 4 SYM lattice action.
From these eigenvalues we can directly compute the mode number for small Ω, allowing a first cross-check of the its stochastic estimation through the Chebyshev approach. Figure 5 shows a representative example confirming that the Chebyshev method indeed reproduces all the features of the mode number obtained from the eigenvalue spectrum. In particular, the large gap in the eigenvalue spectrum is clearly reflected by the plateau in the mode number.
In Figure 6 we compare the mode number for all six available values of the lattice 't Hooft coupling λ lat for 12 4 lattices with gauge group U(2). We include the analytic result for the free theory, and can see that the large degeneracies of the free-theory eigenvalues are broken even for the smallest λ lat = 0.25 we consider. As the coupling gets stronger, the smallest eigenvalues move towards zero while larger fluctuations    Table 1 and include small errorbands showing the resulting statistical uncertainty. The analytic result for the free theory is also shown. make the mode number a smoother function of Ω. We now turn to the extraction of the effective anomalous dimension from these numerical results for the lattice N = 4 SYM mode number.

Estimates for the anomalous dimension
From Eq. (8) we see that the effective anomalous dimension appears in the slope 2/(1 + γ eff ) of the mode number vs. Ω 2 on double-logarithmic axes, as shown in Figure 6. The general behavior shown in this figure therefore already reveals the main features of our results for γ eff and its dependence on λ lat . As we saw for the free theory in Section 4, the stepwise behavior of the mode number at very small Ω obstructs precise extraction of γ eff . For larger Ω the slope decreases with increasing λ lat , which implies a larger γ eff . In addition, for stronger couplings this region of larger γ eff extends towards smaller Ω 2 , though the non-linearities in the results imply that γ eff decreases for small Ω and may become consistent with the IR limit lim Ω→0 γ eff ≈ 0 observed for the free theory and expected for the interacting continuum theory. We now confirm these main features through more quantitative analyses enabled by our precise data for mode number, focusing on the L = 16 ensembles with gauge group U(2). As a first investigation we use Eq. (8) to fit the stochastic Chebyshev mode number data across the complete window [0, Ω 2 ], obtaining the results shown in Figure 7. We perform correlated fits and omit from the figure results from fits  . We again include free-theory results (for = 1, cf. Figure 2), and now omit results from any fits that produce a correlated with standard fit errors larger than 0.1. Within uncertainties the results for λ lat ≤ 2 clearly converge to the expected anomalous dimension of zero in the IR. The stepwise behavior of the mode number in the IR makes this the most difficult region to resolve, as we illustrate by including in Figure 7 results from fitting the analytic free-theory mode number. While the strongest coupling λ lat = 2.5 still shows a trend of γ eff approaching zero in the IR, the current data do not suffice to establish this concretely. Either (or both) larger lattice volumes or more refined analyses are needed for this coupling.
To begin exploring alternative analyses, in Figure 8 we present results obtained by fitting the stochastic Chebyshev mode number data across windows [Ω 2 , Ω 2 + ] of fixed size . This is the procedure we presented in Section 4, and as explained there even the free theory shows significant deviations from zero. We have tested several possible window lengths , excluding those that are so small they produce large oscillations in γ eff , as well as those that produce large correlated χ 2 /d.o.f. > 10. In Figure 8 we again include results from the corresponding fits of the analytic free-theory mode number with = 1, which were previously shown in Figure 2. Compared to Figure 7 this analysis produces much smaller uncertainties, with γ eff ≈ 0 in the IR for λ lat ≤ 1.5. There are also clear trends towards zero for λ lat = 2 and 2.5, but those results remain non-vanishing down to the smallest Ω 2 we can access with this approach on 16 4 lattices.
Finally, in Figure 9 we show results from a third method, which uses Eq. (10) to improve stability by normalizing the stochastic Chebyshev mode number data with respect to the free theory. We choose the reference scale Ω 2 1 = 8 to be beyond of the range considered in the figure without becoming too large. In this approach oscillations from the stepwise behavior of the mode number are clearly visible for small Ω 2 , but the main features discussed above remain the same: γ eff increases for stronger couplings λ lat , while approaching zero in the IR. Again, the 16 4 lattice volume doesn't suffice to completely resolve the convergence to zero for the strongest couplings we consider.

Gauge group and volume dependence
In Section 6 we focused on results for the effective anomalous dimension γ eff for the L = 16 ensembles with gauge group U(2) listed in Table 1. We now investigate the dependence of our results on the lattice volume L 4 and the gauge group U(N c ). In Figure 10 we fix the lattice 't Hooft coupling λ lat = 0.5 and observe reasonably consistent behavior in the mode number for all available L and N c , with similar slopes on double-logarithmic axes implying similar γ eff for sufficiently large Ω 2 . The low-lying eigenvalues clearly depend on the volume and gauge group, as expected, but away from the small-Ω 2 region the only outlier is the U(2) ensemble with L = 12, which may be related to the choice of µ for this ensemble. Figure 11 shows the same consistency for the stronger couplings λ lat = 1 and 1.5 where we have multiple L and N c to compare.
As shown by Eq. (17), we can use the volume scaling of individual eigenvalues to obtain alternative estimates of the effective anomalous dimension. In part because we only directly compute the 2 × 6 lowest D † D eigenvalue pairs, we expect this approach to provide only very rough estimates. Indeed, the results y k = 2/(1+γ eff ) ≈  Table 2 are significantly different from the expected y k = 2. Similar results can be seen from the data in Ref. [12] for the lowest eigenmodes, while the higher modes are more consistent with y k = 2. While we expect the lowest eigenmodes to be those most affected by finite-volume effects, analyzing the free-theory eigenvalues in this way produces reasonable agreement with the expected scaling, implying that the results in Table 2 cannot be directly attributed to lattice artifacts.
We can obtain a more complete picture of the volume scaling by considering the stochastic Chebyshev mode number data, which we show in Figure 12 for fixed λ lat = 0.5 and N c = 2. In the left panel of this figure, we plot the normalized mode number against L 2 Ω 2 so that the approximate volume independence across intermediate scales corresponds to the expected γ eff ≈ 0. However, the lowest modes clearly depart from this scaling, and by plotting these same data vs. L 4 Ω 2 in the right panel we can confirm that they prefer the y ≈ 4 shown in Table 2. It is possible that this behavior could be caused by the interactions in the theory inducing fermion (near-)zero modes despite the deformation in Eq. (3) and the antiperiodic BCs we use. The sampling of near-zero modes by the RHMC algorithm would be suppressed due to the large forces that would arise in the molecular dynamics evolution used to update the field configurations. It may be that the interplay between the interactions in the theory vs. this algorithmic effect could be responsible for both the λ lat dependence of the minimum eigenvalues in Table 1 as well as their unexpected volume scaling in Table 2 and Figure 12.
We complete our discussion with a short comment on the N c dependence of our results. In the free theory the degeneracy of the lowest eigenvalues scales with N 2 c , so we might expect the stepwise patterns in the mode number at small Ω to increase ∼ N 2 c in the interacting theory as well. As shown in Figure 13, we do not see such behavior in our stochastic Chebyshev mode number data. For λ lat ≤ 2, all the lattice volumes and gauge groups we analyze exhibit a cluster of 4 lowest modes followed by a second cluster of 12 additional modes, well separated from the rest of the spectrum. While the size of the lowest eigenmodes λ 2 min scales approximately proportional to 1/N c , this empirical observation might be affected by the interplay between interactions vs. algorithmic details discussed above.

Conclusion
We have presented initial non-perturbative investigations of the RG properties of N = 4 SYM regularized on a space-time lattice, as part of a broader program of numerical investigations of this theory [7]. Employing ensembles of lattice field configurations generated using the RHMC algorithm with an improved lattice action, we stochastically estimated the Chebyshev expansion of the mode number of the fermion operator D † D, and analyzed these data to obtain an effective anomalous dimension γ eff that is expected to vanish in the conformal continuum theory, γ * (λ) = 0. These RG properties are quite challenging to study in discrete lattice space-time, due to the necessary breaking of conformal invariance and 15 of the 16 supersymmetries despite the preservation of a closed supersymmetry subalgebra by our formulation of lattice N = 4 SYM. Our work reported here provides new information about the resulting lattice artifacts and the recovery of N = 4 SYM in the continuum limit. In addition to our main non-perturbative numerical analyses, we also considered the free theory on the A * 4 lattice. This allowed us to check our stochastic estimation of the mode number, to test our extraction of the effective anomalous dimension, and to explore the lattice artifacts that lead to non-zero γ eff even in this case. We carried out further validation of our main Chebyshev results by checking them against direct measurements of the low-lying eigenvalues as well as a more computationally expensive stochastic projection method. All three of these analyses are provided in our public parallel software. We compared three strategies for extracting the effective anomalous dimension from the Chebyshev mode number, observing the same general features in each of the corresponding Figs. 7-9. These show γ eff increasing for stronger lattice 't Hooft couplings λ lat , while still approaching zero in the IR, with the convergence to zero not completely resolved for the strongest λ lat = 2.5 we consider, which may need to be analyzed on larger lattice volumes.
Finally, we considered the dependence of our results on the lattice volume L 4 and the gauge group U(N c ). While we observed the expected L 2 scaling of our stochastic mode number results in an intermediate range of scales, the lowest eigenvalues instead scaled like L 4 , which we speculated could be connected to the RHMC algorithm used to generate lattice field configurations. The multiplicities of those low-lying eigenvalues also don't display the expected dependence on N 2 c , while their size scales approximately proportional to 1/N c , which may also be affected by algorithmic details.
Overall, while our numerical results indicate that lattice artifacts are increasingly significant at larger couplings λ lat , the convergence towards the expected γ * (λ) = 0 in the IR provides reassurance that the correct superconformal continuum limit remains accessible from 16 4 lattice volumes for λ lat ≤ 2. For larger λ lat ≥ 2.5 it seems larger lattices will be needed in order to be confident that the continuum limit is under control. This is useful input for other ongoing studies of lattice N = 4 SYM, investigating inter alia the static potential and the Konishi operator scaling dimension [7]. Similar work can also be considered for alternative N = 4 SYM lattice actions currently being explored [25]. Our results also highlight features of the lowest-lying eigenmodes that are not yet clearly understood, and merit further consideration.  Table 3: Additional information about the ensembles summarized in Table 1. These violations of the three Q Ward identities discussed in the text contain complementary information about the approach to the continuum limit. As expected [11], the violations decrease as λ lat decreases, as L increases, and as N c increases.
symmetries of N = 4 SYM. In particular, Ref. [8] shows how the recovery of all 16 supersymmetries of the theory results from the restoration of Q combined with a set of discrete R-symmetries, subgroups of the full SU(4) R .
In this Appendix we supplement Table 1 by reporting numerical results for the violations of three Q Ward identities for each of the 18 ensembles we consider. Here we briefly describe the three Ward identities under consideration, which are discussed in detail in Ref. [11]. Each can be expressed as the vacuum expectation value of the supersymmetry transformation of a suitable local operator, QO . Such a local operator already appears in the Q-exact part of the lattice action shown in Eq. (1). Because the fermion action is gaussian, this Ward identity fixes the bosonic action per lattice site to be s B = 9N 2 c /2, and we can therefore define as a normalized measure of its violation.
Another local operator was pointed out by Ref. [23]: Here we introduce the shorthand "F " for the second term that involves the ηψ a fermion bilinear and "D" for the first term that depends on the equations of motion for the bosonic auxiliary field d, which are affected by the deformation in Eq. (3). Again we define a normalized measure of the violations of this Ward identity, estimating the fermion bilinear stochastically using random gaussian noise vectors. Finally, the presence of Eq. (3) in the improved lattice action makes Q Tr η = Tr d non-trivial. (The finite-difference term in Eq. (19) vanishes identically upon averaging over the lattice volume.) Defining det P to be the average of the plaquette determinant over all orientations and lattice sites, our final Ward identity violations are simply W det ≡ | Re det P − 1|, which is sensitive only to the U(1) sector of U(N ) = SU(N ) × U(1). In Table 3 we collect numerical results for these three Q Ward identity violations from the 18 lattice N = 4 SYM ensembles summarized in Table 1. These results provide complementary information about the approach to the continuum limit where Q is restored along with all the other symmetries of the theory. Although we include normalization factors in Eq. (18) and Eq. (20), these Ward identity violations can (and clearly do) all have different overall scales. All that matters is that they vanish in the continuum limit, and Ref. [11] found (considering λ lat ≤ 2) that the improved action we use in this work produces effective O(a) improvement in these continuum extrapolations. Here we will be content to note that the Ward identity violations in Table 3 all systematically decrease as L increases towards the continuum limit-and also as λ lat decreases or as N c increases.