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

EEigenvalue spectrum and scaling dimensionof lattice N = 4 supersymmetric Yang–Mills Georg Bergner ∗1 and David Schaich †21 University of Jena, Institute for Theoretical Physics,Max-Wien-Platz 1, D-07743 Jena, Germany Department of Mathematical Sciences, University of Liverpool,Liverpool L69 7ZL, United Kingdom

Abstract

We investigate the lattice regularization of N = 4 supersymmetric Yang–Mills theory, by stochastically computing the eigenvalue mode number of thefermion operator. This provides important insight into the non-perturbativerenormalization group ﬂow of the lattice theory, through the deﬁnition of ascale-dependent eﬀective mass anomalous dimension. While this anomalousdimension is expected to vanish in the conformal continuum theory, the ﬁnitelattice volume and lattice spacing generically lead to non-zero values, whichwe use to study the approach to the continuum limit. Our numerical results,comparing multiple lattice volumes, ’t Hooft couplings, and numbers of colors,conﬁrm convergence towards the expected continuum result, while quantifyingthe increasing signiﬁcance of lattice artifacts at larger couplings. Four-dimensional maximally supersymmetric Yang–Mills theory ( N = 4 SYM) iswidely studied in theoretical physics. Its many symmetries—including conformalsymmetry and an SU(4) R-symmetry in addition to Q = 16 supersymmetries—makeit arguably one of the simplest non-trivial quantum ﬁeld theories in four dimensions,especially in the large- N c planar limit of its SU( N c ) gauge group [1]. This simplicityenabled its role as the conformal ﬁeld theory of the original AdS/CFT holographicduality [2], provided early insight into S-duality [3], and continues to inform modernanalyses of scattering amplitudes [4].At the same time, the non-triviality of N = 4 SYM makes it important toexplore the lattice regularization of the theory. In addition to providing in principlea non-perturbative deﬁnition of N = 4 SYM, lattice ﬁeld theory is also a way to ∗ [email protected] † [email protected] a r X i v : . [ h e p - l a t ] F e b umerically predict its behavior from ﬁrst principles, even at strong coupling andaway from the planar limit. A prominent target for such predictions is the spectrumof conformal scaling dimensions, which depend on the ’t Hooft coupling λ = N g .Although lattice ﬁeld theory has been very successfully used to analyze non-supersymmetric 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 dis-cretization of space-time. See Refs. [5–7] for reviews of these diﬃculties and thesigniﬁcant progress that has been achieved to overcome them in recent years.In particular, using ideas borrowed from topological ﬁeld theory and orbifoldconstructions, a lattice formulation of N = 4 SYM has been developed which pre-serves a closed supersymmetry subalgebra at non-zero lattice spacing a > a → N = 4SYM in the continuum limit [8,9]. In addition, the moduli space of the lattice theorymatches that of continuum N = 4 SYM to all orders in perturbation theory, and therenormalization group (RG) β function vanishes at one loop in lattice perturbationtheory [10].Of course, numerical lattice ﬁeld theory calculations require both a non-zerolattice spacing that corresponds to a ultraviolet (UV) cutoﬀ scale 1 /a , as well as aﬁnite lattice volume ( L · a ) that introduces an eﬀective infrared (IR) cutoﬀ, explicitlybreaking conformal scale invariance. It is also necessary to softly break the singlepreserved supersymmetry in order to regulate ﬂat directions and make the latticepath integral well deﬁned [7, 11]. These facts make it challenging to analyze theapproach to the a → It is therefore essential to carry out detailednumerical studies of the non-perturbative RG properties of the lattice theory.In this paper we present progress investigating RG properties of lattice N = 4SYM. Speciﬁcally, we compute the eigenvalue mode number of the fermion operator,and use this to estimate the ‘mass anomalous dimension’ γ ∗ ( λ ) that would appearin the scaling dimension of the corresponding fermion bilinear. In the continuumtheory this anomalous dimension is expected to vanish, γ ∗ = 0, for all values ofthe ’t Hooft coupling λ . Computing γ ∗ ( λ ) is hence a way to assess the eﬀects ofbreaking supersymmetry and conformal symmetry in numerical lattice calculations,and verify that the properties of N = 4 SYM are correctly reproduced in the con-tinuum limit. Ref. [12] presented a ﬁrst preliminary investigation of this topic, andpreliminary results from a similar project studying the anomalous dimension of theKonishi operator more recently appeared in Ref. [7]. While Ref. [12] proceeded bynumerically computing the fermion operator eigenvalues, that approach quickly be-comes ineﬃcient as the lattice volume increases. Here we instead apply stochastic The vanishing β function in the continuum limit makes this problem much more diﬃcult thanthe case of lattice QCD, where asymptotic freedom can guarantee weak coupling at the UV scaleof the lattice spacing, even when the lattice volume is large enough to access strong coupling atlong distances. N = 4 SYM in the next section, wesummarize these stochastic techniques in Section 3, also discussing how the resultingmode number provides information on the anomalous dimension. In Section 4 weconsider the free ( λ = 0) lattice theory, both to check our methods and to explorediscretization artifacts. Our numerical results for both the low-lying eigenvaluesand the stochastic mode number are presented in Section 5, and in Section 6 weuse these to estimate the anomalous dimension. After looking more closely at thedependence of the results on the gauge group and lattice volume in Section 7, weconclude in Section 8 with some discussion of the next steps for lattice analyses of N = 4 SYM. N = 4 SYM

As mentioned above, the lattice formulation of N = 4 SYM that we use has its ori-gins in both topologically twisted [13–15] and orbifolded [16–19] approaches, whichultimately produce equivalent constructions [20,21]. Here we will use the twisted lan-guage, which organizes the Q = 16 supercharges of the theory into integer-spin rep-resentations of a twisted rotation group SO(4) tw ≡ diag [SO(4) euc ⊗ SO(4) R ], whereSO(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 theserepresentations into 1-, 5- and 10-component sets Q , Q a and Q ab = −Q ba , respect-ively. These transform under the S point-group symmetry of the A ∗ lattice weuse to discretize space-time, which consists of ﬁve basis vectors symmetrically span-ning four dimensions. The twisted-scalar supersymmetry is nilpotent, preservingthe subalgebra {Q , Q} = 0 even at non-zero lattice spacing where Q a and Q ab arebroken.The lattice action, just like the continuum theory, is now the sum of the following Q -exact and Q -closed terms [8–11, 22–24]: S exact = N λ lat X n Tr (cid:20) Q (cid:18) χ ab ( n ) D (+) a U b ( n ) + η ( n ) D ( − ) a U a ( n ) − η ( n ) d ( n ) (cid:19)(cid:21) S closed = − N λ lat X n Tr (cid:20) ε abcde χ de ( n + b µ a + b µ b + b µ c ) D ( − ) c χ ab ( n ) (cid:21) , (1)where n indexes the lattice sites and repeated indices are summed. In this paperwe will present results in terms of the input lattice ’t Hooft coupling λ lat , whichdiﬀers slightly from the continuum λ [7, 23]. The fermion ﬁelds η , ψ a and χ ab = − χ ba transform in the same way as the corresponding twisted supercharges, andare respectively associated with the lattice sites, links and oriented plaquettes. Thegauge and scalar ﬁelds are combined into the ﬁve-component complexiﬁed gaugelinks U a and U a , which appear in the ﬁnite-diﬀerence operators D (+) a and D ( − ) a [20,21].These complexiﬁed gauge links imply U( N ) = SU( N ) × U(1) gauge invariance,with ﬂat directions in both the SU( N ) and U(1) sectors that need to be regulated innumerical calculations, as mentioned in Section 1. To achieve this, we work with the3mproved action introduced by Ref. [11], which adds two deformations to Eq. (1).The ﬁrst is a simple scalar potential with tunable parameter µ , S scalar = N λ lat µ X n X a (cid:18) N Tr h U a ( n ) U a ( n ) i − (cid:19) , (2)which regulates the SU( N ) ﬂat directions while softly breaking the Q supersym-metry. The second deformation is Q -exact, and replaces the term Q (cid:18) η ( n ) D ( − ) a U a ( n ) (cid:19) −→ Q η ( n ) D ( − ) a U a ( n ) + G X a = b (det P ab ( n ) − I N c (3)in S exact , with tunable parameter G . This deformation picks out the U(1) sectorthrough the determinant of the plaquette oriented in the a – b plane, P ab ( n ), which isan N c × N c matrix at each lattice site n . Expanding the Q transformation producesterms that modify both the fermion operator and the bosonic action [11], as expectedfor a supersymmetric deformation.Using this improved action, we have generated many ensembles of ﬁeld conﬁgur-ations using the rational hybrid Monte Carlo (RHMC) algorithm [26] implementedin parallel software that we make publicly available [24]. In addition, we havemodiﬁed this software to implement the stochastic estimation of the mode numberdiscussed below. For this stochastic computation, it is convenient to rescale someof the fermion ﬁeld components, which is irrelevant for the path integral since it in-troduces only a constant prefactor. This rescaling is done only for the measurementof the mode number, not yet in RHMC conﬁguration generation.The particular rescaling we perform is chosen to put the fermion operator intoits most symmetric form, which simpliﬁes analytic considerations and changes thedegeneracies of eigenvalues. Reference [10] previously reported on the analytic struc-ture of the lattice theory. In its conventions, which diﬀer slightly from Eq. (1), thefermion operator f D has the formΨ T f D Ψ = e χ ab D (+)[ a ψ b ] + η D † ( − ) a ψ a + 12 ε abcde e χ ab D † ( − ) c e χ de , (4)where Ψ = ( η, ψ a , χ ab ) collects the 16 fermion ﬁelds into a vector. We adjust thisoperator by reducing summations for e χ ab to the relevant part over a < b , compens-ating a factor of 2 by rescaling χ ab = 2 e χ ab . The same result up to an overall factorof two is obtained from Eq. (1) by rescaling η → η .The more symmetric fermion operator D deﬁned in this way is equivalent to theoriginal operator in Ref. [18], and a rescaling was also done to discuss symmetriesin Ref. [8]. In the free theory, the squared operator D † D is now block diagonalin momentum space, D † D ∼ f ( p ) I N c , implying a 16 N c -fold degeneracy of theeigenvalues. The function f ( p ) on the A ∗ lattice is f ( p ) = 4 X µ =1 sin ( p µ /

2) + 4 sin − X µ =1 p µ / , (5) There is ongoing exploration of alternative lattice actions that address this issue in diﬀerentways [25]. github.com/daschaich/susy p µ are any four of the ﬁve linearly dependent lattice momenta [10]. This16 N c -fold degeneracy is lifted in the interacting theory, but for any ’t Hooft coupling λ lat ≥ / − pairs, so that the non-negative eigenvalues of the squared operatorare always 2-fold degenerate. The mode number, which is the integrated eigenvalue density of the fermion oper-ator, allows for a precise estimate of the mass anomalous dimension [27–29]. Onthe lattice the most practical deﬁnition is obtained from the spectral density of themassless fermion operator D , ρ ( ω ) = 1 V X k h δ ( ω − λ k ) i . (6)Here the eigenvalues λ k should not be confused with the ’t Hooft coupling λ . Themode number ν (Ω) is deﬁned to be the number of eigenvalues λ k of the non-negativeoperator D † D that are smaller than Ω : ν (Ω) = Z Ω b ρ ( ω ) dω = 2 Z Ω0 ρ ( ω ) dω , (7)where b ρ is the spectral density of D † D and the second equality follows from theeigenvalue pairing mentioned above. Throughout the paper quantities like theeigenvalues λ k and the scale Ω are provided in lattice units.The anomalous dimension γ ∗ governs the dependence of the mode number onthe scale Ω : ν (Ω) ∝ (Ω ) / (1+ γ ∗ ) . (8)Additional terms present in the Wilson-fermion case (see Ref. [27]) do not appearhere. This makes it possible to deﬁne a scale-dependent eﬀective anomalous dimen-sion from any two values of the mode number: γ eﬀ (Ω) = 2 log(Ω ) − log(Ω )log( ν (Ω )) − log( ν (Ω )) − , (9)where Ω ≡ (Ω + Ω ) /

2. In addition to depending on the choice of scales Ω andΩ , the determination of γ eﬀ on any ensemble of lattice ﬁeld conﬁgurations will beaﬀected 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 continuumlimit.If Ω and Ω are close to each other and the lattice is coarse, the results of thisnaive method are quite unstable and ﬂuctuate signiﬁcantly. We will show below Lattice QCD experts may expect the upper limit of integration in Eq. (7) to involve Λ = p Ω − m R , with m R a renormalised fermion mass [27]. In this work m R = 0 and Λ = Ω. ν free , using someﬁxed reference scale Ω , γ eﬀ (Ω) = log( ν free (Ω)) − log( ν free (Ω ))log( ν (Ω)) − log( ν (Ω )) − . (10)Now that we have seen how an eﬀective anomalous dimension can be extractedfrom the mode number, we review our stochastic estimation of ν (Ω) using the well-established projection method proposed in Ref. [30]. This method is based on arational approximation of the projection operator P for eigenvalues in the regionbelow a given threshold. The mode number is then ν (Ω) = h Tr P (Ω) i , (11)where the trace is obtained by stochastic estimation. The projection operator isapproximated in terms of the step function h ( x ) using P (Ω) ≈ h ( X ) X = 1 − ∗ D † D + Ω ∗ . (12)The parameter Ω ∗ ≈ Ω is adjusted to minimize the error of the approximation—seeRef. [30] for further details.More recently, a diﬀerent method based on a Chebyshev expansion of the spectraldensity ρ has been proposed [31, 32]. In this Chebyshev expansion method, thespectrum is rescaled to the interval [ − ,

1] by deﬁning M = 2 D † D − λ − λ λ − λ , (13)where λ and λ are the maximal and minimal eigenvalues of D † D . We considerthe integral of the spectral density ρ M of the rescaled operator M multiplied by the n th term T n of a Chebyshev polynomial of order N p , c n = Z − ρ M ( x ) T n ( x ) dx ≤ n ≤ N p , (14)which we estimate stochastically using N s random Z noise vectors v l : c n ≈ N s N s X l =1 h v l | T n ( M ) | v l i . (15)Based on the orthogonality relations for the T n , the spectral density ρ M is nowapproximated by ρ M ( x ) ≈ π √ − x N p X n =0 (2 − δ n ) c n T n ( x ) . (16)The spectral density of D † D is then obtained by mapping the interval [ − ,

1] backto the original eigenvalue region [ λ , λ ], and can be integrated analytically toprovide the mode number via Eq. (7). 6 . . . .

81 0 2 4 6 8 10 12 14 16 18 ˜ ν / ( N c V ) Ω Chebyshev approximationanalytic result − . . . . . . . . . . .

01 0 0 . . . ˜ ν / ( N c V ) Ω Chebyshev approximationanalytic result

Figure 1: Free-theory mode number analytic result and Chebyshev approximationfor lattice volume V = 8 with N c = 2 and antiperiodic BCs, the latter using N p = 1000 and N s = 10 estimators. The right plot zooms in on Ω ≤ . N c V .The work we present below will focus on measurements of the mode number usingthe Chebyshev expansion method. We use polynomials of order 5 , ≤ N p ≤ , λ , λ ], which increases for stronger ’t Hooftcouplings. We cross-checked these results through the more computationally ex-pensive projection method of Eq. (11), as well as by directly computing the low-lying eigenvalues of D † D using a Davidson-type method provided by the PRecon-ditioned Iterative Multi-Method Eigensolver (PRIMME) library [33]. In additionto checking the stochastic results for small Ω, these direct eigenvalue measurementscan provide an alternative estimate of the eﬀective anomalous dimension, from thevolume-scaling relation [12] D λ k E ∝ L − y k , (17)where y k = 2 / (1 + γ eﬀ ) and h λ k i is the average k th eigenvalue of D † D . We willsee below that the mode number provides more precise results than the individualeigenvalues, as expected [28]. Before presenting numerical results obtained from analyzing the available ensemblesof lattice N = 4 SYM ﬁeld conﬁgurations, we consider the λ lat = 0 free theoryto test both our methods of stochastically estimating the mode number, as well asour extraction of γ eﬀ . By deﬁnition, all anomalous dimensions vanish for the freetheory, meaning that any non-zero results we obtain will provide information aboutthe lattice artifacts we want to explore. The free theory is simple enough that weare able to analytically compute the mode number on the A ∗ space-time lattice.The free D † D operator has an enhanced symmetry, resulting in larger degen-eracies of the eigenvalues, especially when all ﬁelds are subject to periodic boundaryconditions (BCs) in all four dimensions. In RHMC conﬁguration generation, and in7 γ e ﬀ Ω L = 16 periodic L = 16 antiperiodic L = 8 antiperiodic Figure 2: Eﬀective anomalous dimension γ eﬀ obtained from using Eq. (8) to ﬁt theanalytic mode number of the free lattice D † D operator across the window [Ω , Ω +1].Two diﬀerent A ∗ lattice volumes 8 and 16 are compared, the latter considering bothperiodic and antiperiodic fermion BCs in the time direction.Figure 1, antiperiodic BCs are applied to the fermion ﬁelds in the time direction,to lift a fermion zero mode. Even with those antiperiodic BCs, many free-theoryeigenvalues are degenerate, leading the mode number to exhibit distinct steps atcertain values of Ω. This stepwise behavior is typically quite diﬃcult to capturewith a polynomial approximation, but Figure 1 shows that the Chebyshev approachis able to provide reasonable precision.While the eﬀective anomalous dimension γ eﬀ can be estimated from any twovalues of the mode number using Eq. (9), this naive approach produces large ﬂuc-tuations, especially for small Ω where the stepwise behavior of the mode numberis most prominent. More stable results are obtained by ﬁtting the mode numberaccording to Eq. (8) over a window [Ω , Ω + ‘ ] of length ‘ , which we show in Fig-ure 2. This ﬁgure compares ‘ = 1 ﬁt results for the free-theory γ eﬀ from two diﬀerentlattice volumes 8 and 16 , considering both periodic and antiperiodic BCs for thelatter. As expected, we obtain more stable results as we approach the continuumlimit by increasing the lattice volume. We also see larger oscillations for the case ofperiodic BCs, which is due to the larger degeneracies of the eigenvalues in this case.While the IR limit lim Ω → γ eﬀ ≈ (cid:46)

7, the eﬀective anomalous dimension tends towardsnegative values. In Figure 3 we conﬁrm that this trend is due to lattice artifacts, bycomparing three dispersion relations for the free operator in momentum space. Forthe A ∗ lattice with periodic BCs, this is based on the function f ( p ) in Eq. (5). Con-sidering instead a continuum-like dispersion relation p produces eﬀective anomalousdimension results more consistent with the true γ ∗ = 0, while the hypercubic-latticedispersion relation 4 sin ( p/

2) leads to the same qualitative trends in γ eﬀ as the A ∗ lattice used for N = 4 SYM.From these investigations of the free theory, we can conclude that the Chebyshev8 − . . .

52 0 2 4 6 8 10 12 14 γ e ﬀ Ω L =8 L =16 (a) f ( p ) − . . . . γ e ﬀ Ω L =8 L =16 (b) p − . − . . . . . γ e ﬀ Ω L =8 L =16 (c) 4 sin ( p/ Figure 3: Free-theory results for the eﬀective anomalous dimension γ eﬀ obtained fromusing Eq. (8) to ﬁt the analytic mode number over the window [Ω , Ω +1]. Each plotcompares L lattice volumes with L = 8 and 16, using periodic BCs. The left plotshows the A ∗ lattice result using f ( p ) from Eq. (5). The center plot considers insteada naive continuum-like discretization of the free operator in momentum space, whilethe right plot corresponds to the hypercubic-lattice dispersion relation 4 sin ( p/ γ eﬀ are manageably small when we ﬁt the mode number overa window of adequate length ‘ , and are expected to decrease for the λ lat > N = 4SYM ﬁeld conﬁgurations we analyze span a range of ’t Hooft couplings 0 . ≤ λ lat ≤ . Our main numerical analyses involve the 18 ensembles of lattice N = 4 SYM ﬁeldconﬁgurations listed in Table 1, which are a subset of a broader collection of en-sembles being used to investigate other aspects of the theory [7]. These were gen-erated using the RHMC algorithm and the improved lattice action described inSection 2. For gauge group U(2) we consider L lattice volumes with 10 ≤ L ≤ . ≤ λ lat ≤ .

5. In addition we also analyze morelimited data sets for gauge groups U(3) and U(4), which involve signiﬁcantly largercomputational costs.As part of the process of conﬁguration generation, we use PRIMME [33] to meas-ure the extremal eigenvalues of f D † f D on every saved conﬁguration, ensuring thatthe minimum and maximum across the entire ensemble remain within the spectralrange where the RHMC rational approximation is accurate. This information isincluded in Table 1, where we also report the autocorrelation time of the smallest9 c L λ lat µ e λ e λ Spect. range Meas. τ ( e λ )2 10 0.5 0.16 9 · −

25 [1 · − ,

45] 700 4.12 12 0.25 0.095 6 · −

23 [5 · − , · −

25 [1 · − , · −

29 [1 · − ,

50] 240 2.61.5 0.23 2 · −

33 [5 · − , · −

36 [5 · − , · −

43 [1 · − , · −

25 [1 · − ,

45] 700 10.32 16 0.25 0.07 2 · −

22 [1 · − ,

50] 180 11.30.5 0.1 2 · −

25 [1 · − ,

45] 180 8.01.0 0.14 1 · −

29 [1 · − ,

50] 410 4.21.5 0.17 6 · −

33 [1 · − ,

50] 230 4.82.0 0.2 4 · −

38 [1 · − ,

50] 250 7.52.5 0.22 2 · −

44 [1 · − ,

50] 340 6.13 12 0.5 0.15 3 · −

24 [5 · − , · −

31 [5 · − , · −

45 [5 · − , · −

23 [5 · − , N c ) and lat-tice volume L . For all ensembles G = 0 .

05 and the fermion ﬁelds are subject toantiperiodic BCs in the time direction. We set µ ≈ √ λ lat /L to remove the scalarpotential in Eq. (2) in the L → ∞ continuum limit with ﬁxed lattice ’t Hooft coup-ling λ lat . For each ensemble we report the extremal eigenvalues e λ of f D † f D , ensuringthat they remain within the spectral range where the corresponding RHMC rationalapproximation is accurate. The listed number of mode number measurements arecarried out on thermalized conﬁgurations separated by 10 molecular dynamics timeunits. From measurements of e λ on those same conﬁgurations we also estimateautocorrelation times τ in units of measurements.10 . . . . . . e i g e n v a l u e s MC λ lat = 0 . λ lat = 1 . λ lat = 1 . Figure 4: Monte Carlo history of the 2 × D † D eigenvalue pairs for gaugegroup U(2) on 12 lattices with diﬀerent values of the ’t Hooft coupling λ lat . Theeigenvalues are measured on thermalized conﬁgurations separated by 10 ‘MC’ unitson the horizontal axis. f D † f D eigenvalue following thermalization/equilibration, estimated using the ‘auto-corr’ module in emcee [34]. This autocorrelation time provides an indication ofhow many statistically independent samples are available for each ensemble. Someadditional information about these ensembles is collected in the Appendix.Following conﬁguration generation and analysis of thermalization, we addition-ally compute the extremal eigenvalues of the more symmetric rescaled operator D † D on all thermalized conﬁgurations, each separated by 10 molecular dynamicstime units. We also use all of these thermalized conﬁgurations to stochasticallyestimate the D † D mode number through the Chebyshev expansion method, withless-extensive cross-checks using the projection method, as described in Section 3.Table 1 reports the number of measured conﬁgurations for each ensemble.The low-lying eigenvalues of D † D provide a ﬁrst look at the mode number andits scaling. Figure 4 shows the general form of these low-lying eigenvalues, whichfeatures 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 sizedecreases for smaller λ lat due to the lowest two pairs moving to larger values. Similarobservations were also made in Ref. [12] for a diﬀerent N = 4 SYM lattice action.From these eigenvalues we can directly compute the mode number for small Ω,allowing a ﬁrst cross-check of the its stochastic estimation through the Chebyshevapproach. Figure 5 shows a representative example conﬁrming that the Chebyshevmethod indeed reproduces all the features of the mode number obtained from theeigenvalue spectrum. In particular, the large gap in the eigenvalue spectrum isclearly reﬂected 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 lattices with gauge group U(2). We include the analyticresult for the free theory, and can see that the large degeneracies of the free-theoryeigenvalues are broken even for the smallest λ lat = 0 .

25 we consider. As the couplinggets stronger, the smallest eigenvalues move towards zero while larger ﬂuctuations11 × − × − × − × − × − × − × − × − × − .

005 0 .

01 0 .

015 0 .

02 0 .

025 0 .

03 0 .

035 0 .

04 0 . ˜ ν / ( N c V ) Ω eigenvaluesChebyshev expansion Figure 5: Normalized mode number obtained by two diﬀerent methods on a V = 12 lattice for gauge group U(2) at λ lat = 1, averaging over all 240 measured conﬁgura-tions. One method is the direct computation of the 2 × D † D eigenvalue pairsshown in Figure 4. The other method is a stochastic estimation of the Chebyshevexpansion with N p = 5000. .

001 0 .

01 0 . ˜ ν / ( N c V ) Ω λ lat = 0 . λ lat = 0 . λ lat = 1 . λ lat = 1 . λ lat = 2 . λ lat = 2 . free theory Figure 6: Normalized mode number obtained for gauge group U(2) on 12 latticesfor all six values of the coupling λ lat , on double-logarithmic axes. We average overall the stochastic Chebyshev measurements speciﬁed in Table 1 and include smallerrorbands showing the resulting statistical uncertainty. The analytic result for thefree theory is also shown. 12 . − . − . . . . . . . . . γ e ﬀ Ω λ lat = 0 . λ lat = 0 . λ lat = 1 . λ lat = 1 . λ lat = 2 . λ lat = 2 . λ lat = 0 . Figure 7: Eﬀective anomalous dimension γ eﬀ obtained from using Eq. (8) to ﬁt thestochastic Chebyshev mode number across the window [0 , Ω ], for 16 lattices withgauge group U(2). Results from ﬁtting the analytic free-theory mode number areincluded for comparison. The errorbars are dominated by the standard ﬁt error,and we omit results from ﬁts that produce errors larger than 0 . N = 4SYM mode number. From Eq. (8) we see that the eﬀective anomalous dimension appears in the slope2 / (1 + γ eﬀ ) of the mode number vs. Ω on double-logarithmic axes, as shown inFigure 6. The general behavior shown in this ﬁgure therefore already reveals themain features of our results for γ eﬀ and its dependence on λ lat . As we saw forthe free theory in Section 4, the stepwise behavior of the mode number at verysmall Ω obstructs precise extraction of γ eﬀ . For larger Ω the slope decreases withincreasing λ lat , which implies a larger γ eﬀ . In addition, for stronger couplings thisregion of larger γ eﬀ extends towards smaller Ω , though the non-linearities in theresults imply that γ eﬀ decreases for small Ω and may become consistent with the IRlimit lim Ω → γ eﬀ ≈ L = 16 ensembles with gaugegroup U(2). As a ﬁrst investigation we use Eq. (8) to ﬁt the stochastic Chebyshevmode number data across the complete window [0 , Ω ], obtaining the results shownin Figure 7. We perform correlated ﬁts and omit from the ﬁgure results from ﬁts13 . − . − . . . . . γ e ﬀ Ω λ lat = 0 . λ lat = 0 . λ lat = 1 . λ lat = 1 . λ lat = 2 . λ lat = 2 . free theory Figure 8: Eﬀective anomalous dimension γ eﬀ obtained from using Eq. (8) to ﬁt thestochastic Chebyshev mode number across windows [Ω , Ω + ‘ ] with Ω ≥ . . ≤ ‘ ≤ ensemble with gauge group U(2). We againinclude free-theory results (for ‘ = 1, cf. Figure 2), and now omit results from anyﬁts that produce a correlated χ / d.o.f. > .

1. Within uncertainties the results for λ lat ≤ λ lat = 2 . γ eﬀ approaching zero in the IR, the current data do not suﬃce to establishthis concretely. Either (or both) larger lattice volumes or more reﬁned analyses areneeded for this coupling.To begin exploring alternative analyses, in Figure 8 we present results obtainedby ﬁtting the stochastic Chebyshev mode number data across windows [Ω , Ω + ‘ ] ofﬁxed size ‘ . This is the procedure we presented in Section 4, and as explained thereeven the free theory shows signiﬁcant deviations from zero. We have tested severalpossible window lengths ‘ , excluding those that are so small they produce largeoscillations in γ eﬀ , as well as those that produce large correlated χ / d.o.f. > ‘ = 1, which were previously shown in Figure 2.Compared to Figure 7 this analysis produces much smaller uncertainties, with γ eﬀ ≈ λ lat ≤ .

5. There are also clear trends towards zero for λ lat = 2 and2 .

5, but those results remain non-vanishing down to the smallest Ω we can accesswith this approach on 16 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 data14 . . . . . γ e ﬀ Ω λ lat = 0 . λ lat = 0 . λ lat = 1 . λ lat = 1 . λ lat = 2 . λ lat = 2 . Figure 9: Eﬀective anomalous dimension γ eﬀ obtained from Eq. (10) with referencescale Ω = 8, for all six 16 ensembles with gauge group U(2).with respect to the free theory. We choose the reference scale Ω = 8 to be beyondof the range considered in the ﬁgure without becoming too large. In this approachoscillations from the stepwise behavior of the mode number are clearly visible forsmall Ω , but the main features discussed above remain the same: γ eﬀ increasesfor stronger couplings λ lat , while approaching zero in the IR. Again, the 16 latticevolume doesn’t suﬃce to completely resolve the convergence to zero for the strongestcouplings we consider. In Section 6 we focused on results for the eﬀective anomalous dimension γ eﬀ for the L = 16 ensembles with gauge group U(2) listed in Table 1. We now investigatethe dependence of our results on the lattice volume L and the gauge group U( N c ).In Figure 10 we ﬁx the lattice ’t Hooft coupling λ lat = 0 . L and N c , with similarslopes on double-logarithmic axes implying similar γ eﬀ for suﬃciently large Ω . Thelow-lying eigenvalues clearly depend on the volume and gauge group, as expected,but away from the small-Ω 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 sameconsistency for the stronger couplings λ lat = 1 and 1 . L and N c to compare.As shown by Eq. (17), we can use the volume scaling of individual eigenvalues toobtain alternative estimates of the eﬀective anomalous dimension. In part becausewe only directly compute the 2 × D † D eigenvalue pairs, we expect thisapproach to provide only very rough estimates. Indeed, the results y k = 2 / (1+ γ eﬀ ) ≈ e-081e-071e-061e-050.00010.0010.010.1 .

01 0 . ˜ ν / ( N c V ) Ω N c = 2 , L = 10 N c = 2 , L = 12 N c = 2 , L = 14 N c = 2 , L = 16 N c = 3 , L = 12 N c = 4 , L = 12 Figure 10: Normalized mode number for diﬀerent volumes L and gauge groupsU( N c ) at a ﬁxed lattice ’t Hooft coupling λ lat = 0 . .

001 0 .

01 0 . ˜ ν / ( N c V ) Ω N c = 2 , L = 12 N c = 2 , L = 16 N c = 3 , L = 12 (a) λ lat = 1 .

001 0 .

01 0 . ˜ ν / ( N c V ) Ω N c = 2 , L = 12 N c = 2 , L = 16 N c = 3 , L = 12 (b) λ lat = 1 . Figure 11: Normalized mode number for diﬀerent volumes L and gauge groupsU( N c ) at couplings λ lat = 1 (left) and λ lat = 1 . λ lat L y y y k = 2 / (1 + γ eﬀ ) obtained from Eq. (17) for U(2) en-sembles using the corresponding L = 16 results for reference. The two estimatesconsider the average second and eleventh eigenvalues, h λ i and h λ i .16 . . ˜ ν / ( N c ) L Ω L = 10 L = 12 L = 14 L = 16 (a) L scaling . . ˜ ν / ( N c ) L Ω L = 10 L = 12 L = 14 L = 16 (b) L scaling Figure 12: Scaling of the N c = 2 normalized mode number with the lattice volume L at λ lat = 0 .

5. For an intermediate range of Ω the results are consistent with L scaling, whereas the lowest modes scale like L .4 for two representative k = 2 and 11 in Table 2 are signiﬁcantly diﬀerent from theexpected y k = 2. Similar results can be seen from the data in Ref. [12] for thelowest eigenmodes, while the higher modes are more consistent with y k = 2. Whilewe expect the lowest eigenmodes to be those most aﬀected by ﬁnite-volume eﬀects,analyzing the free-theory eigenvalues in this way produces reasonable agreementwith the expected scaling, implying that the results in Table 2 cannot be directlyattributed to lattice artifacts.We can obtain a more complete picture of the volume scaling by consideringthe stochastic Chebyshev mode number data, which we show in Figure 12 for ﬁxed λ lat = 0 . N c = 2. In the left panel of this ﬁgure, we plot the normalizedmode number against L Ω so that the approximate volume independence acrossintermediate scales corresponds to the expected γ eﬀ ≈

0. However, the lowest modesclearly depart from this scaling, and by plotting these same data vs. L Ω in theright panel we can conﬁrm that they prefer the y ≈ λ lat dependence of the minimum eigenvalues in Table 1 as well as theirunexpected volume scaling in Table 2 and Figure 12.We complete our discussion with a short comment on the N c dependence of ourresults. In the free theory the degeneracy of the lowest eigenvalues scales with N c ,so we might expect the stepwise patterns in the mode number at small Ω to increase ∼ N c in the interacting theory as well. As shown in Figure 13, we do not see suchbehavior in our stochastic Chebyshev mode number data. For λ lat ≤

2, all thelattice volumes and gauge groups we analyze exhibit a cluster of 4 lowest modesfollowed by a second cluster of 12 additional modes, well separated from the rest ofthe spectrum. While the size of the lowest eigenmodes λ scales approximately17 .

02 0 .

04 0 .

06 0 .

08 0 . ˜ ν Ω λ lat = 2 . λ lat = 2 . λ lat = 1 . λ lat = 1 . λ lat = 0 . λ lat = 0 . N c = 4 , L = 12 , λ lat = 0 . N c = 3 , L = 12 , λ lat = 0 . N c = 2 , L = 12 , λ lat = 0 . Figure 13: Total mode number in the small-eigenvalue regime (Ω ≤ .

1) for diﬀerentvolumes L , gauge groups U( N c ) and couplings λ lat . If not otherwise indicated, theresults are for N c = 2 and L = 16. The lowest modes form clusters of 4 + 12eigenvalues independent of L , N c and λ lat ≤ /N c , this empirical observation might be aﬀected by the interplaybetween interactions vs. algorithmic details discussed above. 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 programof numerical investigations of this theory [7]. Employing ensembles of lattice ﬁeldconﬁgurations generated using the RHMC algorithm with an improved lattice ac-tion, we stochastically estimated the Chebyshev expansion of the mode number ofthe fermion operator D † D , and analyzed these data to obtain an eﬀective anom-alous dimension γ eﬀ that is expected to vanish in the conformal continuum theory, γ ∗ ( λ ) = 0. These RG properties are quite challenging to study in discrete latticespace-time, due to the necessary breaking of conformal invariance and 15 of the16 supersymmetries despite the preservation of a closed supersymmetry subalgebraby our formulation of lattice N = 4 SYM. Our work reported here provides newinformation about the resulting lattice artifacts and the recovery of N = 4 SYM inthe continuum limit.In addition to our main non-perturbative numerical analyses, we also consideredthe free theory on the A ∗ lattice. This allowed us to check our stochastic estimationof the mode number, to test our extraction of the eﬀective anomalous dimension, andto explore the lattice artifacts that lead to non-zero γ eﬀ even in this case. We carriedout further validation of our main Chebyshev results by checking them against18irect measurements of the low-lying eigenvalues as well as a more computationallyexpensive stochastic projection method. All three of these analyses are provided inour public parallel software. We compared three strategies for extracting the eﬀectiveanomalous dimension from the Chebyshev mode number, observing the same generalfeatures in each of the corresponding Figs. 7–9. These show γ eﬀ increasing forstronger lattice ’t Hooft couplings λ lat , while still approaching zero in the IR, withthe convergence to zero not completely resolved for the strongest λ lat = 2 . L and the gauge group U( N c ). While we observed the expected L scaling of ourstochastic mode number results in an intermediate range of scales, the lowest ei-genvalues instead scaled like L , which we speculated could be connected to theRHMC algorithm used to generate lattice ﬁeld conﬁgurations. The multiplicities ofthose low-lying eigenvalues also don’t display the expected dependence on N c , whiletheir size scales approximately proportional to 1 /N c , which may also be aﬀected byalgorithmic details.Overall, while our numerical results indicate that lattice artifacts are increasinglysigniﬁcant at larger couplings λ lat , the convergence towards the expected γ ∗ ( λ ) = 0in the IR provides reassurance that the correct superconformal continuum limitremains accessible from 16 lattice volumes for λ lat ≤

2. For larger λ lat ≥ . N = 4SYM, investigating inter alia the static potential and the Konishi operator scalingdimension [7]. Similar work can also be considered for alternative N = 4 SYMlattice actions currently being explored [25]. Our results also highlight features ofthe lowest-lying eigenmodes that are not yet clearly understood, and merit furtherconsideration. Acknowledgements

We thank Simon Catterall and Joel Giedt for ongoing collaboration on lattice N = 4SYM, and particularly appreciate helpful conversations with Joel concerning themode number and mass anomalous dimension. GB acknowledges support fromthe Heisenberg Programme of the Deutsche Forschungsgemeinschaft (DFG) GrantNo. BE 5942/3-1. DS was supported by UK Research and Innovation Future LeaderFellowship MR/S015418/1 and STFC grant ST/T000988/1. Numerical calculationswere carried out at the University of Liverpool, the University of Bern, and onUSQCD facilities at Fermilab funded by the US Department of Energy. Appendix: Additional information on ensembles

In prior work [11], the approach to the continuum limit of lattice N = 4 SYM wasmainly monitored by measuring violations of Ward identities for the twisted-scalarsupersymmetry Q . These violations are introduced by the soft breaking of Q dueto Eq. (2), and must vanish in the continuum limit in order to recover the full19 c L λ lat

Eq. (18) Eq. (20) Eq. (21)2 10 0.5 0.00051(5) 0.014(1) 0.00323(9)2 12 0.25 0.0003(1) 0.005(2) 0.0005(2)0.5 0.00039(3) 0.0085(7) 0.0020(2)1.0 0.00098(7) 0.0197(9) 0.0086(4)1.5 0.00168(7) 0.025(1) 0.0148(6)2.0 0.00226(9) 0.0284(8) 0.0204(6)2.5 0.00336(6) 0.042(1) 0.0351(8)2 14 0.5 0.00022(3) 0.0067(6) 0.00149(6)2 16 0.25 0.00016(6) 0.0024(8) 0.0004(1)0.5 0.00014(4) 0.0054(7) 0.0013(2)1.0 0.00045(3) 0.0106(5) 0.0046(2)1.5 0.00084(6) 0.0152(7) 0.0088(4)2.0 0.00139(4) 0.0203(5) 0.0145(3)2.5 0.00181(3) 0.0239(4) 0.0204(3)3 12 0.5 0.00023(6) 0.007(1) 0.0015(3)1.0 0.00044(5) 0.0089(7) 0.0038(3)1.5 0.00100(5) 0.0108(7) 0.0063(5)4 12 0.5 0.00015(4) 0.0023(8) 0.0006(2)Table 3: Additional information about the ensembles summarized in Table 1. Theseviolations of the three Q Ward identities discussed in the text contain complement-ary information about the approach to the continuum limit. As expected [11], theviolations 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 16supersymmetries of the theory results from the restoration of Q combined with aset of discrete R-symmetries, subgroups of the full SU(4) R .In this Appendix we supplement Table 1 by reporting numerical results for theviolations of three Q Ward identities for each of the 18 ensembles we consider.Here we brieﬂy describe the three Ward identities under consideration, which arediscussed in detail in Ref. [11]. Each can be expressed as the vacuum expectationvalue of the supersymmetry transformation of a suitable local operator, hQOi . Sucha local operator already appears in the Q -exact part of the lattice action shown inEq. (1). Because the fermion action is gaussian, this Ward identity ﬁxes the bosonicaction per lattice site to be s B = 9 N c /

2, and we can therefore deﬁne W s B ≡ |h s B i − . N c | . N c (18)as a normalized measure of its violation.20nother local operator was pointed out by Ref. [23]: Q Tr " η X a U a U a = Tr d X a U a U a − Tr η X a ψ a U a ≡ D − F. Here we introduce the shorthand “ F ” for the second term that involves the ηψ a fermion bilinear and “ D ” for the ﬁrst term that depends on the equations of motionfor the bosonic auxiliary ﬁeld d , d = D ( − ) a U a + G X a = b (det P ab − I N c , (19)which are aﬀected by the deformation in Eq. (3). Again we deﬁne a normalizedmeasure of the violations of this Ward identity, W Bilin ≡ (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)* D − F √ D + F +(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) , (20)estimating the fermion bilinear stochastically using random gaussian noise vectors.Finally, the presence of Eq. (3) in the improved lattice action makes hQ Tr η i = h Tr d i non-trivial. (The ﬁnite-diﬀerence term in Eq. (19) vanishes identically uponaveraging over the lattice volume.) Deﬁning det P to be the average of the plaquettedeterminant over all orientations and lattice sites, our ﬁnal Ward identity violationsare simply W det ≡ | h Re det

Pi − | , (21)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 viola-tions from the 18 lattice N = 4 SYM ensembles summarized in Table 1. Theseresults provide complementary information about the approach to the continuumlimit where Q is restored along with all the other symmetries of the theory. Althoughwe include normalization factors in Eq. (18) and Eq. (20), these Ward identity vi-olations can (and clearly do) all have diﬀerent overall scales. All that matters isthat they vanish in the continuum limit, and Ref. [11] found (considering λ lat ≤ O ( a ) improvementin these continuum extrapolations. Here we will be content to note that the Wardidentity violations in Table 3 all systematically decrease as L increases towards thecontinuum limit—and also as λ lat decreases or as N c increases. References [1] N. Arkani-Hamed, F. Cachazo, and J. Kaplan, “What is the SimplestQuantum Field Theory?,”

JHEP (2010) 016, arXiv:0808.1446 .[2] J. M. Maldacena, “The Large- N limit of superconformal ﬁeld theories andsupergravity,” Adv. Theor. Math. Phys. (1998) 231–252, hep-th/9711200 .[3] H. Osborn, “Topological Charges for N = 4 Supersymmetric Gauge Theoriesand Monopoles of Spin 1,” Phys. Lett.

B83 (1979) 321–326.214] H. Elvang and Y.-t. Huang,

Scattering Amplitudes in Gauge Theory andGravity . Cambridge University Press, 2015.[5] S. Catterall, D. B. Kaplan, and M. Ünsal, “Exact lattice supersymmetry,”

Phys. Rept. (2009) 71–130, arXiv:0903.4881 .[6] G. Bergner and S. Catterall, “Supersymmetry on the lattice,”

Int. J. Mod.Phys. A (2016) 1643005, arXiv:1603.04478 .[7] D. Schaich, “Progress and prospects of lattice supersymmetry,” Proc. Sci.

LATTICE2018 (2019) 005, arXiv:1810.09282 .[8] S. Catterall, J. Giedt, and A. Joseph, “Twisted supersymmetries in lattice N = 4 super-Yang–Mills theory,” JHEP (2013) 166, arXiv:1306.3891 .[9] S. Catterall and J. Giedt, “Real space renormalization group for twistedlattice N = 4 super-Yang–Mills,” JHEP (2014) 050, arXiv:1408.7067 .[10] S. Catterall, E. Dzienkowski, J. Giedt, A. Joseph, and R. Wells, “Perturbativerenormalization of lattice N = 4 super-Yang–Mills theory,” JHEP (2011) 074, arXiv:1102.1725 .[11] S. Catterall and D. Schaich, “Lifting ﬂat directions in lattice supersymmetry,”

JHEP (2015) 057, arXiv:1505.03135 .[12] D. J. Weir, S. Catterall, and D. Mehta, “Eigenvalue spectrum of lattice N = 4super-Yang–Mills,” Proc. Sci.

LATTICE2013 (2014) 093, arXiv:1311.3676 .[13] F. Sugino, “A Lattice formulation of super-Yang–Mills theories with exactsupersymmetry,”

JHEP (2004) 015, hep-lat/0311021 .[14] F. Sugino, “Super-Yang–Mills theories on the two-dimensional lattice withexact supersymmetry,”

JHEP (2004) 067, hep-lat/0401017 .[15] S. Catterall, “A Geometrical approach to N = 2 super-Yang–Mills theory onthe two dimensional lattice,” JHEP (2004) 006, hep-lat/0410052 .[16] A. G. Cohen, D. B. Kaplan, E. Katz, and M. Ünsal, “Supersymmetry on aEuclidean space-time lattice. 1. A target theory with four supercharges,”

JHEP (2003) 024, hep-lat/0302017 .[17] A. G. Cohen, D. B. Kaplan, E. Katz, and M. Ünsal, “Supersymmetry on aEuclidean space-time lattice. 2. Target theories with eight supercharges,”

JHEP (2003) 031, hep-lat/0307012 .[18] D. B. Kaplan and M. Ünsal, “A Euclidean lattice construction ofsupersymmetric Yang–Mills theories with sixteen supercharges,”

JHEP (2005) 042, hep-lat/0503039 .[19] M. Ünsal, “Twisted supersymmetric gauge theories and orbifold lattices,”

JHEP (2006) 089, hep-th/0603046 .2220] S. Catterall, “From Twisted Supersymmetry to Orbifold Lattices,”

JHEP (2008) 048, arXiv:0712.2532 .[21] P. H. Damgaard and S. Matsuura, “Geometry of Orbifolded SupersymmetricLattice Gauge Theories,”

Phys. Lett.

B661 (2008) 52–56, arXiv:0801.2936 .[22] S. Catterall, P. H. Damgaard, T. Degrand, R. Galvez, and D. Mehta, “PhaseStructure of Lattice N = 4 Super-Yang–Mills,” JHEP (2012) 072, arXiv:1209.5285 .[23] S. Catterall, D. Schaich, P. H. Damgaard, T. DeGrand, and J. Giedt, “ N = 4supersymmetry on a space-time lattice,” Phys. Rev.

D90 (2014) 065013, arXiv:1405.0644 .[24] D. Schaich and T. DeGrand, “Parallel software for lattice N = 4supersymmetric Yang–Mills theory,” Comput. Phys. Commun. (2015)200–212, arXiv:1410.6971 .[25] S. Catterall, J. Giedt, and G. C. Toga, “Lattice N = 4 super-Yang–Mills atstrong coupling,” JHEP (2020) 140, arXiv:2009.07334 .[26] M. A. Clark and A. D. Kennedy, “Accelerating Dynamical-FermionComputations Using the Rational Hybrid Monte Carlo Algorithm withMultiple Pseudofermion Fields,”

Phys. Rev. Lett. (2007) 051601, hep-lat/0608015 .[27] A. Patella, “A precise determination of the ¯ ψψ anomalous dimension inconformal gauge theories,” Phys. Rev. D (2012) 025006, arXiv:1204.4432 .[28] A. Cheng, A. Hasenfratz, G. Petropoulos, and D. Schaich, “Determining themass anomalous dimension through the eigenmodes of Dirac operator,” Proc.Sci.

LATTICE2013 (2014) 088, arXiv:1311.1287 .[29] Z. Fodor, K. Holland, J. Kuti, D. Nógrádi, and C. H. Wong, “The chiralcondensate from the Dirac spectrum in BSM gauge theories,”

Proc. Sci.

LATTICE2013 (2014) 089, arXiv:1402.6029 .[30] L. Giusti, C. Hoelbling, M. Luscher, and H. Wittig, “Numerical techniques forlattice QCD in the epsilon regime,”

Comput. Phys. Commun. (2003)31–51, hep-lat/0212012 .[31] Z. Fodor, K. Holland, J. Kuti, S. Mondal, D. Nogradi, and C. H. Wong, “Newapproach to the Dirac spectral density in lattice gauge theory applications,”

Proc. Sci.

LATTICE2015 (2016) 310, arXiv:1605.08091 .[32] G. Bergner, P. Giudice, G. Münster, I. Montvay, and S. Piemonte, “Spectrumand mass anomalous dimension of SU(2) adjoint QCD with two Diracﬂavors,”

Phys. Rev. D (2017) 034504, arXiv:1610.01576 .2333] A. Stathopoulos and J. R. McCombs, “PRIMME: preconditioned iterativemultimethod eigensolver–methods and software description,” ACM Trans.Math. Softw. (2010) 21.[34] D. Foreman-Mackey, D. W. Hogg, D. Lang, and J. Goodman, “emcee: TheMCMC Hammer,” PASP (2013) 306–312, arXiv:1202.3665arXiv:1202.3665