Testing holography using lattice super-Yang--Mills on a 2-torus

We consider maximally supersymmetric SU(N) Yang--Mills theory in Euclidean signature compactified on a flat two-dimensional torus with anti-periodic (`thermal') fermion boundary conditions imposed on one cycle. At large N, holography predicts that this theory describes certain black hole solutions in Type IIA and IIB supergravity, and we use lattice gauge theory to test this. Unlike the one-dimensional quantum mechanics case where there is only the dimensionless temperature to vary, here we emphasize there are two more parameters which determine the shape of the flat torus. While a rectangular Euclidean torus yields a thermal interpretation, allowing for skewed tori modifies the holographic dual black hole predictions and results in another direction to test holography. Our lattice calculations are based on a supersymmetric formulation naturally adapted to a particular skewing. Using this we perform simulations up to N=16 with several lattice spacings for both skewed and rectangular tori. We observe the two expected black hole phases with their predicted behavior, with a transition between them that is consistent with the gravity prediction based on the Gregory--Laflamme transition.


I. INTRODUCTION
Maximally supersymmetric Yang-Mills (SYM) theory in p + 1 dimensions has been conjectured to provide a holographic description of string theories containing Dp-branes. Specifically, this gauge/gravity duality states that (p + 1)dimensional SYM with gauge group SU(N ) is dual to a Type IIA (even p) or Type IIB (odd p) superstring containing N coincident Dp-branes in the 'decoupling' limit [1,2]. The p = 3 case corresponds to superconformal N = 4 SYM in four dimensions and yields the original AdS/CFT correspondence [3]. In this paper we focus on the maximally supersymmetric Yang-Mills in two dimensions at finite temperature, with the spatial circle direction compactified with periodic fermion boundary conditions (BCs) about it.
In this context, at large N and low temperatures, the dual string theory is well described by supergravities whose dynamics are given by certain charged black holes. Two classes of black hole are required to describe these dynamicsthose that wrap the spatial circle (so-called 'homogeneous black strings') and those that are localized on it ('localized black holes') [4][5][6][7][8][9]. Indeed this system of black hole solutions is related by a simple transform to the static uncharged black holes arising in pure gravity in ten dimensions with one spatial dimension wrapped into a circle, i.e., tendimensional Kaluza-Klein theory [8,10] (for a review of black holes in Kaluza-Klein theory see [11]). The two classes have different thermodynamic behaviors, and there is a first-order Gregory-Laflamme [12] phase transition between them in the gravity dual. According to holography, all this should be reproduced by the thermal physics of the SYM. In particular, the phase transition is a deconfinement transition associated to the spatial circle, the magnitude of the spatial Wilson line giving an order parameter. It is thought that this transition extends to high temperatures where an intricate phase structure has been revealed from numerical and analytic treatments [8,13,14].
The remarkably subtle nature of gauge/gravity duality has meant that whilst SYM thus provides a fundamental and microscopic quantum description of certain gravity systems, there still is no 'proof' or derivation of this black hole thermodynamics from (p + 1)-dimensional SYM directly. Indeed even understanding the local structure of the dual ten-dimensional spacetime which emerges from the strongly coupled SYM theory remains a mystery. While there has been some heuristic analytic treatment for general p that hints how certain aspects of black hole thermodynamics can be seen within the SYM theory [15][16][17], and an approximation scheme developed for p = 0 [18][19][20][21], a full derivation showing the SYM reproduces dual black hole behavior remains an important challenge in quantum gravity.
With only limited success from analytic treatment, it is natural to apply lattice field theory, which is well suited to study the thermodynamics of strongly coupled systems (see for example the recent review [22]). Starting with [23,24], several works over the past decade have studied the thermal behavior of the p = 0 SYM quantum mechanics, where again gravity provides a black hole prediction to be tested, and striking agreement has been seen [25][26][27][28][29][30][31][32][33][34][35][36]. However the dual gravity in that setting is simpler than in the p = 1 case we focus on here, where there are different black holes to probe and a gravity phase transition to observe. Less effort has been directed at this two-dimensional case, where the state of the art until recently was simply to provide evidence for the transition at small N ≤ 4 [37]. 1 One of the main goals of this work is to improve the lattice study of this phase transition, working at larger N and smaller lattice spacings. We will also provide large-N tests of the detailed thermal behavior of the two different classes of black holes. (See also the recent conference proceedings [53,54] for other lattice work in this direction.) Conventionally one studies thermal physics in the canonical ensemble by considering the Euclidean theory. This lives on a flat rectangular 2-torus, with the spatial cycle being the circle of the original theory, and the Euclidean time cycle having anti-periodic BCs for fermions and period equal to the inverse temperature. The path integral then plays the role of a thermal partition function. An important point we emphasize in this work is that one may also consider the Euclidean theory on a flat but skewed torus as discussed in [9]. This no longer corresponds to the Lorentzian theory at finite temperature, but taking anti-periodic fermion BCs about Euclidean time, it may be regarded as a generalized thermal ensemble. The key point is that this skewing is easily accommodated in the dual gravity theory, which can be treated in the Euclidean signature, and its behavior is again given in terms of solutions that may be interpreted as generalized black holes.
Studying such skewed flat tori is natural due to the lattice SYM formulation that we employ. Recently, much progress has been made in lattice studies of the p = 3 theory, N = 4 SYM, using a novel construction based on discretization of a topologically twisted form of the continuum N = 4 action. See Ref. [55] for a review of this approach. The chief merit of this new lattice construction is that it preserves a closed subalgebra of the supersymmetries at non-zero lattice spacing. Numerical studies of the four-dimensional theory are in progress [56][57][58][59][60][61][62], but are quite expensive because of the large number of degrees of freedom. In this regard, lower-dimensional theories are more tractable and can be studied extensively at large N with better control over continuum extrapolations. These lattice constructions are based on non-hypercubic Euclidean lattices, which when made periodic are naturally adapted to skewed tori. We dimensionally reduce an N = 4 lattice system to give a discretization of two-dimensional SU(N ) SYM on an A * 2 lattice, preserving four exact supercharges at non-zero lattice spacing. Applying appropriate BCs we then carry out calculations for N ≤ 16, large enough to see dual gravity behavior. Varying the temporal and spatial lattice extent gives the continuum SYM on tori that may be both skewed and rectangular. We confirm that both phases of dual black hole behavior are seen in the appropriate low-temperature regime, and we see consistency between the generalized SYM thermodynamics and that predicted by gravity. We also see a transition between these phases, again compatible with the expectation from gravity, which extends to high temperature as expected.
The plan of the paper is as follows. In Sec. II we review the known predictions for large-N thermal two-dimensional SYM on a spatial circle-i.e., SYM on a flat rectangular Euclidean 2-torus. Then in Sec. III we discuss how this picture generalizes for a flat skewed Euclidean 2-torus. In Sec. IV we present our lattice construction for this skewed continuum theory. Then in Sec. V we discuss our results, focusing on how the various gravity predictions are confirmed. We end the paper with a brief discussion.

II. REVIEW OF THERMAL LARGE-N (1 + 1)-DIMENSIONAL SYM ON A CIRCLE
We now review the predictions for large-N p = 1 SYM, compactified on a circle of size L at temperature T = 1/β, derived in various limits and using input from the dual gravity theory [1,[6][7][8][9]63]. We treat the thermal theory in Euclidean signature, with Euclidean time τ ∼ τ + β, so that it lives on a flat rectangular 2-torus, with side lengths β and L. Fermions have thermal (anti-periodic) BCs on the Euclidean time circle, and are taken periodic on the spatial circle. Starting in Sec. III we will consider the theory on a skewed torus. However, it will be useful to review the rectangular torus case first, as the skewed case will be similar.
The Euclidean action of the theory is Here X I with I = 2, . . . , 9 are the eight spacetime scalars representing the transverse degrees of freedom of the branes. They are N × N hermitian matrices in the adjoint representation of the gauge group. The fermion Ψ and matrices Γ I descend from a dimensional reduction of a ten-dimensional Euclidean Majorana-Weyl spinor, with Ψ also transforming in the adjoint. The dimensionful 't Hooft coupling λ = N g 2 Y M may be used to construct two dimensionless quantities that control the dynamics: r β = β √ λ and r L = L √ λ. We define the dimensionless temperature t = 1/r β . Since we are interested in the large-N 't Hooft limit we wish to consider N → ∞ with r β and r L fixed. The main observables we consider are thermodynamic quantities related to the expectation value of the bosonic action, and also the Wilson loop magnitudes P β and P L (normalized to 1), 2 where which wrap about the Euclidean thermal circle and spatial circle of the two-dimensional Euclidean torus, respectively. For the large-N theory these act as order parameters for phase transitions associated to breaking of the Z N center symmetry of the gauge group. 3 For the thermal circle this is the usual thermal deconfinement transition, with vanishing Polyakov loop P β = 0 at large N indicating the (unbroken) confined phase, and P β = 0 being the (broken) deconfined phase. We will use similar terminology for P L , namely that P L = 0 indicates 'deconfined' spatial behavior while P L = 0 corresponds to 'confined' spatial behavior.
A. High-temperature limit Consider the high-temperature limit of the SYM [8,9]. Then when r 3 β r L we may integrate out Kaluza-Klein modes on the thermal circle and reduce to a bosonic quantum mechanics (BQM) consisting of the zero modes on the thermal circle. Due to the thermal fermion BCs, this is now the bosonic truncation of the p = 0 SYM, as the fermions are projected out in the reduction. Now the 't Hooft coupling λ BQM is related to the original two-dimensional coupling as λ BQM = λ β and the dynamics implies β A ∼ 0 so that P β = 0 indicating thermal deconfinement. Taking the small-volume limit, L 3 λ BQM 1, the dynamics of this BQM (and hence the full SYM) is governed by a bosonic matrix integral of scalar and gauge field zero modes. These dynamics imply that L A ∼ 0, so that the SYM theory in this regime is spatially deconfined with P L = 0. Following Refs. [64,65] the leading behavior of the BQM energy in this regime goes as E BQM 6N 2 /L. The action behaves as S BQM = −E BQM L/3 (see, e.g., [24]). We expect the BQM action to give the SYM bosonic action, S Bos , when the BQM describes it, since the fermions are decoupled in this limit. Hence we expect the SYM bosonic action to go as S Bos −2N 2 . Since this limit applies when we have integrated out both the temporal and spatial Kaluza-Klein modes, reducing to only a bosonic matrix integral, its behavior is common to any high-temperature, small-volume limit, r L , r β 1. For finite volume, L 3 λ BQM ∼ 1, this BQM has an interesting dynamics at large N . This has been studied numerically and analytically in [8,13,14,66] with the conclusion that there is a deconfinement transition around L 3 λ BQM 1.4. However the order of the transition is difficult to determine [8]. Either it is a first-order transition (as most recently found in [66]) or it is a strong second-order transition (as discussed in the earlier [13,14]), in which case there is another very close-by third-order Gross-Witten-Wadia (GWW) [67,68] transition as well.

B. Dual gravity for low temperatures
At large N and low temperatures r β 1, holography predicts a gravity dual given by D1-charged black holes in IIB supergravity [1]. These have a simple solution, with Euclidean string frame metric and dilaton given as where d 1 = 2 6 π 3 and U 2 . There is a 2-form potential yielding N units of D1 charge, and the spatial circle x corresponds to that in the SYM, with x ∼ x + L. Large N is required to suppress string quantum corrections to 2 We define the Wilson loop to be the trace of the Wilson line Pe i β,L A around a closed path. 3 Since we are at finite volume we can only have a phase transition at large N . the supergravity. In order to suppress α corrections we require 1 r β , and to avoid winding mode corrections about the circle we need r β r 2 L . When r β ∼ r 2 L one indeed finds that this solution is unstable to stringy winding modes on the spatial circle x [4][5][6][7][8][9]. This is seen by passing to a second gravity dual by T-dualizing on this circle direction to obtain a solution in IIA supergravity, which in string frame is Now the spatial coordinatex ∼x + L IIA is compact, but due to the T-duality, has period L IIA = (2π) 2 α /L, and there is a 1-form potential supporting D0 charge. The D1 charge of Eq. (3) is given as a distribution of D0 charge smeared homogeneously over the circle. This gravity solution is a good dual for the SYM at large N and 1 r β (to suppress string quantum and α corrections, respectively). To avoid winding mode corrections about the circle we also require r L r β . In particular, for 1 r β this T-dual frame overlaps the regime r L r β r 2 L where the IIB dual exists and describes the physics. It describes the regime where the IIB solution fails and becomes unstable to winding modes, r β ∼ r 2 L , and also covers smaller circle sizes all the way down to the limit r L → 0 where the physics is that of the dimensionally reduced SYM quantum mechanics.
The above solution is homogeneous on the circle-a 'homogeneous black string'. The black hole horizon wraps over the circle direction and has a cylindrical topology R × S 7 . Being related by T-duality it has precisely the same thermodynamics as the IIB solution above. Namely it predicts the thermodynamic behavior for the SYM free energy density f , with t = 1/r β the dimensionless temperature. However what was a winding mode in the original frame is now a classical Gregory-Laflamme (GL) instability in this IIA frame. One finds the above solution is dynamically unstable when where the constant is determined by numerically solving the differential equation that governs the marginal instability mode [8,10]. Thus at smaller circle sizes the above solution remains, but it is not the relevant one for the dynamics, which instead is given in terms of a 'localized black hole' solution. This is inhomogeneous over the circle direction, with the black hole horizon being localized on the circle and having a spherical topology S 8 . From a gravitational perspective the parameter that is varied to move between different solutions is the dimensionless ratio of the size of the horizon as compared to the size of the spatial IIA circle, L IIA . Translating to our SYM variables, this is proportional to r 2 L /r β . These localized black hole solutions are not known analytically. 4 However, recently the challenging numerical construction of these solutions has been performed [69]. Following expectations, Ref. [69] found that for large r 2 L /r β the thermodynamic behavior is dominated by the homogeneous phase. At there is a first-order phase transition to the localized phase, which then dominates the homogeneous one for smaller r 2 L /r β , having lower free energy density. The value of c grav is determined numerically, and we see it is rather close to c GL .
While the analytic form of these localized solutions is not known generally, they do simplify in the limit that the horizon is small compared to the circle size. In SYM variables this is the case for y = r 2 L /r β 1, where the solutions have an approximate behavior [10] 1 N 2 λ f loc = − 2 21 · 3 2 · 5 7 π 14 7 19 Of relevance for us is that the position of the phase transition found in Ref. [69] is such that the leading term in the above approximation (8) agrees very well (to the percent level) with the full numerical solutions over the full range where this phase dominates the thermodynamics. This means that while one generally requires the numerical solutions of [69] to deduce the thermodynamics of a given localized solution, since we are only concerned with this localized branch for r 2 L /r β where it dominates the thermodynamics we may very accurately approximate the thermal behavior for such solutions using (8). 5 In previous lattice investigations [37] this transition was probed using small N = 3 and 4, finding evidence for consistency with the value for c grav in (7). In this work we will improve the lattice study of the phase diagram, employing larger N and smaller lattice spacings.
We will refer to the homogeneous phase as the D1 phase, since in the IIB duality frame it is the D1-brane solution, although we note that it may also be seen as a homogeneous D0-brane solution in the IIA frame. We will refer to the localized phase as the D0 phase, since it may only be seen in gravity in the IIA frame where it is a localized D0-brane black hole.
Since all these gravity solutions are static black holes, their Euclidean time circle is contractible so we expect a deconfined Polyakov loop, P β = 0. 6 In the IIB frame, as the horizon wraps over the spatial circle for the homogeneous black string, this spatial cycle is not contractible in the bulk solution. Hence at large N we expect spatial confinement, P L = 0, when this homogeneous phase describes the thermodynamics (for 1 c grav r β < r 2 L ), with the thermal behavior given by Eq. (5). The homogeneity of the horizon is taken to indicate that the eigenvalues of Pe i L A are uniformly distributed at large N . On the other hand, upon decreasing the circle size r L at fixed r β we have a first-order transition to the localized phase with thermodynamics given by Eq. (8). Due to the localized horizon, the D0-brane charge is compactly supported on the spatial circle, so we expect the eigenvalue distribution for Pe i L A is likewise compactly supported [8,73]. This implies spatial deconfinement, P L = 0. The phase transition curve r 2 L = c grav r β in the gravity regime, r β 1, therefore corresponds to a first-order spatial deconfinement transition associated to P L . We emphasize that we are interested in temperatures and circle sizes where r β and r L ∼ O(N 0 ) in the large-N limit. If we were to take ultralow temperatures r β → ∞ as some sufficiently large positive power of N , then the gravity predictions above would cease to be valid because the gravity would become strongly coupled near the black hole horizons. In particular, for r β ∼ N the theory is thought to enter a conformal phase described by a free orbifold CFT [1,2], which we will not explore in this work. 5 Using only the leading term in (8) and comparing with (5) gives an approximation for the phase transition with cgrav within 2% of the numerically computed value in (7). Including the subleading term improves this to be consistent with the value in (7) within its numerical uncertainty. 6 Recall that when IIB gravity provides a good dual description of the SYM we expect the Wilson loop (normalized as in (2)) about a cycle in this boundary theory to be non-vanishing if that cycle is contractible when extended into the dual bulk (such as for a cycle about Euclidean time when a horizon exists in the bulk) [70][71][72]. Conversely if a cycle is non-contractible in the bulk, we expect the corresponding Wilson loop to vanish. Note this picture does not hold for the spatial cycle after T-dualizing to the IIA frame. Then, instead, the distribution of D0 charge on the spatial circle is thought to determine the eigenvalue distribution of Pe i L A .

C. Summary for SYM on a rectangular torus
For large-N two-dimensional SYM on a rectangular Euclidean 2-torus we have two dimensionless parameters r β and r L . Assuming that r β , r L ∼ O(N 0 ) in the large-N limit we have the following expectations: • The high-temperature, small-volume limit r β , r L 1 is described by the dynamics of scalar and gauge zero modes. We expect P β , P L = 0 and S Bos −2N 2 .
• The high-temperature limit r 3 β r L reduces to BQM. Here we expect P β = 0. For r 3 L < c BQM r β with c BQM 1.4 we expect P L = 0, with a deconfinement transition to P L = 0 for r 3 L > c BQM r β . • The low-temperature limit, r β 1, admits a gravity-dual black hole description, so P β = 0, with the free energy density depending on the ratio r 2 L /r β . For r 2 L = c grav r β , with c grav = 2.45 there is a first-order deconfinement phase transition with respect to P L . The D1 phase, for r 2 L > c grav r β , has free energy density given by Eq. (5) and P L = 0. The D0 phase, for r 2 L < c grav r β , has P L = 0 with the free energy density well approximated by Eq. (8). This is illustrated in Fig. 1.

III. BEHAVIOUR ON A SKEWED TORUS
We now discuss what happens to the Euclidean theory when it is placed on a skewed torus. The motivation is twofold. First we emphasize that skewing the torus provides a new direction to deform the theory and hence a new independent test of holography, given again that gravity dual predictions exist. Second, modern supersymmetric lattice constructions often naturally live on non-hypercubic lattices, which in turn are naturally adapted to giving a continuum Euclidean theory on a skewed torus.
We take the SYM to live on the flat 2-torus generated as a quotient of the two-dimensional plane. Writing the metric as this is generated by the identifications (τ, x) ∼ (τ, x) + β fermions anti-periodic (10) (τ, x) ∼ (τ, x) + L fermions periodic (11) for vectors β ( L) about the thermal (spatial) circles with anti-periodic (periodic) fermion BCs. These vectors defining the cycles have lengths and dot product with |γ| < 1 to ensure a non-degenerate torus. Thus the geometry is determined by three parameters-in addition to the β and L of the rectangular case there is the dimensionless skewing parameter γ. For any value of γ (not just the rectangular case γ = 0) we will be able to compare our numerical SYM results to gravity predictions, and hence test holography.
We have constructed the torus by a quotient of the plane by the vectors L and β. However any SL(2, Z) transformation of these vectors will define the same geometric torus. The anti-periodic BCs about the β circle, and periodic ones about L, complicate this slightly. Consider the transformation generated by the following subgroup of SL(2, Z): where we note that then a ∈ 2Z−1. Then the 2-torus defined by β (anti-periodic fermions) and L (periodic fermions) is the same as that defined by β and L. The fermion BCs restrict us to a subgroup of the full SL action, so that the new β has an odd coefficient multiplying β to maintain anti-periodicity, and likewise L has an even coefficient multiplying β.
In the rectangular case, we have a Lorentzian thermal interpretation of the physics, where we identify β as inverse temperature in a canonical ensemble. In the skewed case this is no longer true. Nonetheless, we may regard this case as a generalized thermal ensemble, with 1/β playing the role of generalized temperature. As for a rectangular torus we again work with dimensionless now supplemented by the dimensionless parameter γ. We denote t = 1/r β the generalized dimensionless temperature, and also define the aspect ratio The 2-volume is then given as Vol T 2 = β 2 α 1 − γ 2 . In practice on the lattice it will be convenient to scan the parameter space by fixing the 'shape' of the torus set by (α, γ) and varying the dimensionless temperature t that controls the size of the torus in units of the SYM coupling λ.
The redundancy in our description of the torus using β and L given in (13) translates into an invariance in this parameterization: a set α, γ, t defined from β and L is equivalent to other parameters α , γ , t similarly defined from β and L . In the usual manner we may describe this using the complex 'modular parameter' τ, given by 7 which encodes the (dimensionless) shape parameters α, γ, and is independent of the torus scale. The parameter transforms under the action (13) as so that m, n, c ∈ Z and a(2m − 1) − 2nc = 1. We term this a restricted modular transformation-it is a usual modular transformation, but restricted to preserve our fermion BCs (anti-periodic on the β cycle and periodic on the L cycle). As for the usual modular invariance of the torus, we may define a fundamental domain for this τ parameter under the action of this restricted transform (13), which gives the set of inequivalent tori (taking into account fermion BCs).
A modular parameter outside this domain can then be mapped back into it using the appropriate (13). The usual modular transformations are generated by τ → τ + 1 and τ → −1/τ, or equivalently τ → τ + 1 and τ → τ τ+1 . We provide some review of this in Appendix A. However for our restricted transform instead we may use the generators τ → τ + 2 and τ → τ τ+1 . The fundamental domain D may then be taken to be These assertions are proved in Appendix A. The lattice construction we use later will give torus geometries with a particular γ (which we will see is γ = −1/2), and we will vary the shape parameter α of the torus and its size t. Some of the torus shapes we study will correspond to τ within the fundamental domain, and others will not. We will generally present results in terms of α and t for this common value of γ, but as we discuss later, for the shapes that fall outside the fundamental domain there is an alternate description with new α , γ , t . 8 An important aim of our lattice calculations is to see the detailed generalized thermodynamic behavior predicted by the dual black holes. However in a skewed setting we are no longer in a strict thermal context and so potentials such as free energy are not strictly meaningful. Instead the relevant quantity is the (logarithm of the) partition function 7 Since 'tau' is conventionally used for both the torus modular parameter and Euclidean time, we attempt to avoid potential confusion by using the symbols τ for the modular parameter and τ for Euclidean time. 8 It is worth emphasizing that the way the modular parameter τ arises is a little different to that in the familiar two-dimensional CFT setting. In the context of Euclidean two-dimensional CFT any 2-torus is Weyl equivalent to a flat torus with modular parameter τ (i.e., one constructed as a quotient of the two-dimensional plane) and unit volume. Hence for any 2-torus the CFT partition function only depends on the geometry through τ. However in our case the SYM is not a CFT-in particular our theory explicitly depends on the size of the torus, and will also depend on the details of its real geometry. We restrict ourselves to only consider 2-tori constructed as quotients of the plane, which are flat but with a skewing parameterized by τ. Hence our SYM depends on τ and the scale (through, say, β) because we have restricted to flat tori. If we also added scalar curvature the situation would be more complicated (unlike for a CFT). However in either case (CFT or our SYM) the dependence is on τ up to modular invariance, so we may always choose τ to be in the fundamental domain. This is simply because the description of the torus has this invariance, and it is not due to any symmetry properties of the field theory.
for the Euclidean action given by Eq. (1), but now on the skewed 2-torus. In a lattice calculation it is hard to determine the value of this integral directly, so instead we focus on expectation values, which are much more convenient to compute. A natural observable is the expectation value of the Euclidean action S. Since the fermionic part of the action is gaussian, its value is simply a constant. Therefore we will focus on measuring the expectation value of the bosonic part of the action, S Bos , given in Eq. (1). We will find it convenient to work with the average bosonic action density, defined from the (renormalized) vev of S Bos in the obvious way, We now show that s Bos is related to a derivative of the partition function ln Z. After scaling the bosonic fields and coordinates as the Euclidean action becomes The 2-torus T 2 is generated by the identifications the former anti-periodic and the latter periodic for the fermions. Thus we have scaled out β. The only explicit β dependence is in the overall coupling, and the integral I[α, γ] depends only on the dimensionless shape parameters α and γ (through the identifications above). By suitably scaling the fermion fields the fermionic action can be chosen to have no β dependence. Then differentiating the partition function with respect to β, keeping α and γ fixed, we obtain Thus computing s Bos as a function of t, α and γ, which may be conveniently done on the lattice, gives the same information as that contained in the partition function.
As in the rectangular case we will be interested in the Wilson loops about the torus temporal and spatial cycles and their magnitudes P β and P L , respectively. Since there are equivalent presentations of the same torus, one can equally consider the loops P β and P L , which correspond to cycles in the original representation that wrap multiples of the cycles generated by β and L. 9 We now proceed to discuss the same limits as in the rectangular case in this skewed geometry. We will find that analogous small circle reductions occur, and that generalized black holes still give gravity predictions for small generalized temperature 1 r β . We first discuss the dimensional reduction of the theory on a skewed torus, and then turn to a discussion of the gravity dual, giving predictions for the observable s Bos .

A. High-temperature limit
We now take the skewing parameter γ to be fixed, and consider the high-temperature limit, finding qualitatively similar behavior to the rectangular torus case we previously discussed. In the small-volume, high-temperature limit 9 The center symmetry transformations of Wilson loops associated to the cycles β and L determine those of β and L . Suppose we define the holonomies W β = Tr Pe i β A and W L = Tr Pe i L A so that center symmetry acts on these as W β → z β W β and W L → z L W L for phases z β,L where z N β,L = 1. Then consider a loop W β associated to β given by the linear combination of cycles β and L in Eq. (13). This will transform now as W β → z c L z 2m−1 β W β and so is determined by the transformations of W β and W L . The same is true for W L . r β , r L 1 precisely the same reduction to the matrix integral occurs. Thus we again expect P β P L = 0, with the bosonic action behaving as S Bos /N 2 −2 at large N . Translating to the bosonic action density we then obtain At finite volume we may again reduce to an effective one-dimensional theory at high temperature. The easiest way to understand this dimensional reduction is to note that skewed tori in the fixed-γ, high-temperature limit become equivalent to nearly rectangular tori (again with small thermal circles) under a suitable transformation (13). Then we may simply dimensionally reduce this equivalent rectangular torus, and pull the result back to the original skewed torus parameterization given by γ.
Thus we consider the high-temperature limit where we fix γ and r L , taking r β → 0 so that α → ∞. We use the transform (13) with c = 0, m = 0 and a = 1, leaving n ∈ Z, to find an equivalent torus that will be approximately rectangular. This maps β, L to This leaves the temperature unchanged, t = t, and relates the shape parameters as In the α → ∞ limit we may choose n appropriately (i.e., taking −αγ/2 n ∈ Z) to obtain an equivalent torus with γ 0 and α α 1 − γ 2 → ∞. This equivalent torus is approximately rectangular and in the high-temperature limit, so from our previous discussion in Sec. II A we may reduce on the thermal circle when We obtain BQM on a circle size L BQM = L with coupling λ BQM = λ/β , so that Thus we have the same relation of the lower-and higher-dimensional coupling as in the rectangular case, but the circle size is related via a skewing-dependent factor. Reducing on the thermal cycle β implies that the Polyakov loop is P β ∼ 1. Given that P β = P β , and this Polyakov loop is trivial, we expect thermal deconfinement, P β = 0. Then P L P L = P BQM , so our previous discussion of BQM implies the deconfinement transition in the limit of Eq. (29) is located at This is associated with a transition from a spatially confined phase with P L = 0 to a deconfined one with P L = 0 as r 3 L /r β is reduced.
B. The low-temperature dual gravity limit -D1 phase Considering the dual IIB gravity we find the same homogeneous D1-charged black hole solution as in Eq. (3), but now we take the two-dimensional torus in the field theory directions to be generated by the identifications anti-periodic fermions Then asymptotically, when U 0 U , the torus spanned by τ and x has our required skewed geometry. We also see that the relation between U 0 and β is exactly the same as for Eq. (3), since the metric is locally the same as in the rectangular case, and the τ circle has the same period β.
We will regard this as a 'generalized black hole' in the sense that for real β and γ = 0 its properties are not related directly to a physical Lorentzian black hole. Nonetheless the Euclidean IIB solution exists. The geometry of this solution differs only globally in the x direction from the rectangular case, and is homogeneous in x. The solution should be a good description of the IIB string theory again for large N and 1 r β . We expect it to become unstable to a winding mode instability for r 2 L ∼ r β . Consider the expectation value of the SYM Euclidean Lagrangian density L E predicted by this solution, which will be homogeneous in both τ and x. Then − ln Z = S = Vol T 2 L E when this black hole is the dominant saddle point of the path integral. Due to the homogeneity in x this Lagrangian density is only a function of β, and has no dependence on the other parameters α and γ determining the shape of the torus. In the rectangular case it simply equals the free energy density of the solution, and thus generally is given by Eq. (5), so L E = −N 2 λ · 2 4 π 5 2 t 3 /3 4 . As before we refer to this homogeneous black string phase as the D1 phase. Hence if this gravitational solution dominates the partition function the theory is in the D1 phase and the bosonic action density is We see explicitly that our observable doesn't depend on the skewing of the torus. As in the rectangular case, ∂/∂τ generates a contractible cycle due to the horizon. Now the Euclidean time cycle of the torus is simply generated by ∂/∂τ . The spatial cycle of the torus is now generated by 1 − γ 2 ∂/∂x + γ∂/∂τ , and since ∂/∂x is not contractible, neither is this spatial torus cycle. Thus when Eq. (3) is the dominant bulk solution, this implies that the Polyakov loop P β = 0, whereas P L = 0. Hence the dual SYM is thermally deconfined, but in a spatially confined phase.
There is an important subtlety in the above discussion: one must be careful whether the solution above (3) does dominate, due to there being other related gravitational dual black holes [9]. In Eq. (32) we have identified with β and L to generate the 2-torus, but as discussed in Ref. [9] we may equally well use any equivalent pair under the transformation (13), since they will give the same flat 2-torus asymptotically. This will yield another inequivalent gravitational dual solution. However, it is the pair of vectors whose corresponding modular parameter τ lies in the fundamental domain that gives the dominant gravitational dual solution. This is understood as follows. Take a τ in the fundamental domain, D, and a temperature t. Then the Euclidean action density is as in Eq. (33). Now transforming this to an equivalent τ outside the fundamental domain results in a new dimensionless temperature t given in terms of t, τ and the transformation. This t is lower than the fundamental-domain temperature t since 10 and from the corollary in Appendix A we have Im(τ )/Im(τ) ≤ 1 for τ ∈ D. Thus we see t ≤ t. Hence the action density of these other gravitational saddle points, which is given by (33) with t → t , is more positive. Since the volume of the torus is preserved under the modular transformation (13), then the Euclidean action is also more positive. Thus these other dual solutions outside the fundamental domain are not the relevant saddle points to determine the partition function behavior.
Hence the subtlety is that the D1-phase prediction is given by Eq. (33) for α, γ, t corresponding to a τ in the fundamental domain. If one has a set of parameters outside the fundamental domain, one must first map them to new parameters α , γ , t in the fundamental domain, and then apply the formula (33) with t → t . 11 In this canonical representation (i.e., τ ∈ D) we will have the prediction P β = 0 and P L = 0 for the Wilson loops. A spatial cycle of the torus is always non-contractible, so we should have P L = 0 in any equivalent representation. However, while P β = 0 for the thermal cycle in the fundamental domain representation, if in the equivalent representation β is a linear combination of both β and L it may correspond to a non-contractible cycle in the gravity dual, so that also P β = 0. We emphasize that this does not imply temporal confinement, since there is some temporal loop (associated to the cycle β ) where P β = 0.
C. The low-temperature dual gravity limit -D0 phase Considering the D1 phase with r β 1 where the gravity is a valid description, fixing γ and reducing the circle size to r 2 L ∼ r β we again expect winding modes on the spatial circle to become important near the horizon, as in Sec. II B. This limit however is straightforward to understand, as it implies r L r β and thus we may play a similar trick as for the small-thermal-circle dimensional reduction in Sec. III A, again mapping to an equivalent almost-rectangular representation. We set a = 1, n = 0 and m = 0 in the transform (13), leaving c ∈ Z. Then so L = L and the equivalence relates the shape parameters as Then by choosing c appropriately (i.e., taking −γ/α c ∈ Z) in the α → 0 limit we obtain an equivalent torus with γ 0 and α = α/ 1 − γ 2 → 0. Thus for fixed non-zero γ and r L r β , so that the modular parameter τ is far outside the fundamental domain, the above transform maps to an approximately rectangular torus with r L = r L and r β = r β 1 − γ 2 . For the rectangular torus we know that the phase transition from the D0 phase to the D1 phase occurs for r 2 L = c grav r β with c grav = 2.45. Mapping this back to our non-rectangular parameterization, we expect a transition at Then for this rectangular torus we may use the approximation (8) for the thermal behavior, replacing t → t . From this we may compute ln Z. Translating back to our original temperature t we may compute the bosonic action density using Eq. (25) to obtain, in our original non-rectangular parameterization, In the rectangular representation when this D0 phase dominates we expect to have P β , P L = 0, so this phase is both spatially and thermally deconfined. Furthermore, since L = L we also have P L = 0. Thus P L remains an order parameter for the transition between the gravity D1 and D0 phases.

D. Summary for SYM on a skewed torus
For large-N SYM on a skewed torus with fixed γ, upon varying r L and r β our expectation is a phase diagram similar to Fig. 1 for the rectangular case. We expect a spatial deconfinement transition line with order parameter P L .
• In the high-temperature, small-volume limit r β , r L 1 we expect P L = 0 and thermal behavior as in Eq. (26).
• For high temperatures r 3 β r L the SYM may be dimensionally reduced to the BQM theory, leading us to expect the phase transitions described in Eq. (31). We will have P L = 0 for r 3 L 1.4r β /(1 − γ 2 ) 3/2 , and P L = 0 otherwise.
• For low temperatures t 1 we expect a IIA or IIB gravity black hole description. The D0 phase, with approximate behavior (38) and P L = 0, dominates for r 2 L < c grav 1 − γ 2 r β . For r 2 L > c grav 1 − γ 2 r β we expect the D1 phase with P L = 0 to dominate with behavior (33), where this formula assumes r L , r β and γ are in the fundamental domain.
We expect all these phases to be thermally deconfined, so assuming the parameters r L , r β and γ are in the fundamental domain then we will have P β = 0. If they are not in the fundamental domain, it is possible in the D1 phase to have P β = 0 even though P β = 0 for a fundamental-domain description of the torus, r L , r β and γ .
In our numerical analyses of the skewed SYM theory it will be convenient to fix α = r L /r β and vary t = 1/r β to scan a 'slice' of the r β × r L plane. For any finite α, at sufficiently high temperature t 1 we will also be in the small-volume regime with P L = 0. As we decrease t, for any finite α we expect to go through a confinement phase transition associated to P L , and for t 1 eventually enter the gravity D1 phase with P L = 0. For large α 1 this will be the confinement transition described by BQM. For small α 1 we expect to enter the gravity regime in the spatially deconfined D0 phase, and encounter the first-order dual Gregory-Laflamme transition to the D1 phase at the lower temperature t = α 2 /(c grav 1 − γ 2 ). These expectations for large, small and intermediate α are illustrated in Fig. 2.

IV. LATTICE FORMULATION
Using ideas borrowed from topological field theory and orbifold constructions it has recently become possible to construct a four-dimensional lattice theory which retains an exact supersymmetry at non-zero lattice spacing and produces N = 4 SYM in the continuum limit. Noting that maximal SYM in any dimension can be thought of as a classical dimensional reduction of the N = 1 SYM in 10 dimensions, it follows that our theory of interest, maximal SYM in two dimensions, can be derived from a dimensional reduction of the four-dimensional N = 4 SYM theory. Thus the approach we take here is to use the existing four-dimensional lattice construction of N = 4 SYM, and reduce this in two directions to obtain a two-dimensional lattice action for our desired two-dimensional SYM.
An interesting subtlety arises, namely that the four-dimensional lattice most naturally has an A * 4 geometry rather than a hypercubic one. When we reduce this lattice action, we obtain a discretization of two-dimensional SYM on a two-dimensional A * 2 (triangular) lattice. Taking the lattice to be periodic, with one direction having thermal (antiperiodic) fermion BCs and the other periodic BCs, we generate two-dimensional SYM on a 2-torus which is skewed, as the A * 2 lattice basis vectors are not orthogonal. However, as emphasized, this should be viewed as a virtue rather than a problem. While there is no direct Lorentzian interpretation of this finite-volume 'generalized' thermal ensemble, as discussed above there are holographic string theory predictions that can be tested, and that is the aim of this paper.
We begin by considering the four-dimensional lattice discretization of topologically twisted N = 4 SYM. We then outline how this is reduced to a two-dimensional system whose continuum limit will be the reduction of N = 4 SYM, giving maximal SYM in two dimensions. Since the reduced lattice has an A * 2 geometry, where we make one lattice direction into the thermal circle of length β and the other into the spatial circle of length L, the continuum limit will be two-dimensional SYM living on a skewed torus. The skewing parameter is then determined from the A * 2 lattice geometry to be γ = −1/2. In Appendix B we provide a detailed discussion of the reduction of a four-dimensional A * 4 lattice theory to the two-dimensional A * 2 lattice for a simpler scalar field theory.

A. Four-dimensional twisted lattice N = 4 SYM
In this section we summarize the important features of this four-dimensional lattice theory before proceeding to its dimensional reduction. The trick to preserving supercharges in a lattice theory is to discretize a topologically twisted formulation of the underlying supersymmetric theory. 12 In the case of N = 4 SYM the twisted construction treats the four-component gauge field and the six massless adjoint scalars of the theory as a five-component complexified gauge field where the roman index 'a' runs from 1 to 5. The four Majorana fermions of the theory are decomposed into χ ab = −χ ba , ψ a and η. The analogous decomposition of the sixteen supercharges provides a twisted-scalar Q corresponding to η, which is nilpotent, Q 2 = 0. The complexified gauge field leads to complexified field strengths where the corresponding complexified covariant derivatives are Using these ingredients we can express the usual N = 4 action as a sum of Q-exact and Q-closed terms, where λ 4 = g 2 N is the usual 't Hooft coupling and we implicitly sum over repeated indices. Here x α are the usual canonical flat-space coordinates with α running from 1 to 4. The action of the scalar supersymmetry charge Q is where d is a bosonic auxiliary field with equation of motion d = D a , D a . Since Q 2 = 0 the Q-exact part of the action is clearly supersymmetric, while Q acting on the Q-closed term vanishes due to a Bianchi identity. The other fifteen supercharges are twisted into a 1-form Q a and antisymmetric 2-form Q ab . This continuum action can be discretized while preserving the single Q supersymmetry as described in Refs. [56,[76][77][78]. This discretization procedure dictates how the continuum fields are placed on the lattice, how derivatives are replaced by lattice difference operators, and even the structure of the underlying lattice itself. Specifically, we must employ the A * 4 lattice whose five basis vectors symmetrically span the four spacetime dimensions. This lattice is a natural generalization of the two-dimensional triangular (A * 2 ) lattice to four dimensions. It possesses five equivalent basis vectors corresponding to the vectors from the center of an equilateral four-simplex out to its five vertices. It has a high S 5 point group symmetry with the dimensions of its low lying irreducible representations matching those of the continuum twisted SO(4) rotation group. Ref. [79] shows that the combination of the Q supersymmetry, lattice gauge invariance and the S 5 global symmetry suffices to ensure that no new relevant operators are generated by quantum corrections. Assuming non-perturbative effects such as instantons preserve the lattice moduli space, only a single marginal coupling may need to be tuned to obtain N = 4 SYM in the continuum limit.
The resultant lattice action takes the form where the lattice difference operators appearing in the above expression are given in Refs. [77,78] and generically take the form of shifted commutators. For example, Remarkably the Q-closed term is still lattice supersymmetric due to the existence of an exact lattice Bianchi identity, Integrating out the auxiliary field d yields The lattice sites in the canonical flat-space coordinates x α of Eq. (42) are arranged as the A * 4 lattice with positions x α = ∆ n ν e α (ν) for n ∈ Z 4 . This discretization is analogous to that discussed explicitly for the scalar theory example in Appendix B. As discussed following Eq. (B4), the resulting continuum action (42) has a coupling related to the lattice coupling as λ 4 = λ lat / √ 5 [56,76]. The presence of an exact lattice supersymmetry allows us to derive an exact expression for the renormalized bosonic action density, which gives the derivative of the partition function with respect to the coupling as in Eq. (25). We find where V denotes the number of lattice sites and S lat B corresponds to the bosonic terms in the lattice action. This definition of the continuum renormalized s Bos has the property that it vanishes as a consequence of the exact lattice supersymmetry, if periodic (non-thermal) BCs are used.
In practice, to stabilize the SU(N ) flat directions of the theory we add to S 0 a soft-supersymmetry-breaking scalar potential with tunable parameter µ. In the dimensionally reduced system this term is particularly important at low temperatures where the flat directions lead to thermal instabilities [80]. This single-trace scalar potential differs from the doubletrace operator used in previous investigations [56][57][58][59][60][61][62] and constrains each eigenvalue of U a U a individually, rather than only the trace as a whole. Exact supersymmetry at µ = 0 ensures that all Q-breaking counterterms vanish as some power of µ. The complexification of the gauge field in Eq. (39) leads to an enlarged U(N ) = SU(N ) ⊗ U(1) gauge invariance. In the continuum the U(1) sector decouples from observables in the SU(N ) sector, but this is not automatic at non-zero lattice spacing [56][57][58]. To regulate additional flat directions in the U(1) sector, we truncate the theory to remove the U(1) modes from U a , making them elements of the group SL(N, C) rather than the algebra gl(N, C). In order to maintain SU(N ) gauge invariance it is necessary to keep the fermions in gl(N, C), explicitly breaking the lattice supersymmetry that would have related U a to ψ a . However, by representing the truncated gauge links as U b = e igaA b , we can argue that the continuum supersymmetry relating A a and ψ a is approximately realized in the large-N limit even at non-zero lattice spacing. This follows from fixing the 't Hooft coupling λ lat = g 2 N as N → ∞, implying g 2 → 0. Then expanding the exponential produces the desired U b = I N + igaA b up to O(ga) corrections that vanish as N → ∞ even at non-zero lattice spacing a. Empirically, when we measure would-be supersymmetric Ward identities we find that they are satisfied up to small (at most percent-level) deviations, and those deviations decrease ∝ 1/N 2 as N increases. (Figure 14 in Appendix C shows some representative results.) Although the constant term in Eq. (49) is no longer exactly 9/2, we expect only comparably small corrections to it and therefore continue to use Eq. (49) to define s Bos . Since our lower-dimensional studies all focus on holographic dualities in the large-N limit, this truncated approach appears viable at least in fewer than four dimensions. Our interest is in two-dimensional maximal SYM on a 2-torus with thermal (anti-periodic) fermion BCs on one cycle and periodic BCs on the other. This two-dimensional maximal SYM is given by the dimensional reduction of the four-dimensional N = 4 theory. Hence to obtain this two-dimensional theory on a 2-torus we simply consider the above four-dimensional lattice discretization of the N = 4 theory, taken on N x × 1 × 1 × N t lattices with periodic BCs in the reduced directions, corresponding to naive dimensional reduction. The gauge fields associated to the reduced directions now transform as site fields ϕ i (n) and are naturally interpreted as the scalar fields arising from dimensional reduction. Their fermionic superpartners, now also site fields, correspond to additional exact lattice supersymmetries [78].
Exactly as for the scalar example discussed in Appendix B, such a reduction results in a continuum theory in two-dimensional flat space where the lattice geometry is that of A * 2 . We take appropriate periodicity conditions on the extended lattice directions, n ∼ n + (0, 0, 0, N t ) with anti-periodic fermions and n ∼ n + (N x , 0, 0, 0) with periodic fermions, thus generating the 2-torus. Since the two-dimensional lattice has A * 2 geometry the 2-torus we generate is not rectangular but skewed, with γ =ẽ 1·ẽ2 |ẽ1||ẽ2| = −1/2. The lengths of the two cycles are then where we see from Appendix B that the two-dimensional lattice spacing is a 2 = ∆ 2/3. Similarly, the two-dimensional continuum gauge coupling is where the second equality from Eq. (B13) considers the coupling as arising from an appropriate Kaluza-Klein reduction of the continuum four-dimensional N = 4 theory. Thus in terms of our dimensionless couplings our lattice action corresponds to two-dimensional N = (8, 8) SYM on a skewed torus with γ = −1/2 and which, being dimensionless, are independent of the scale ∆ as they should be. Noting that in two dimensions the continuum SYM is super-renormalizable, we do not expect any renormalization of the classical geometry of this 2-torus. Finally, we add an additional soft-Q-breaking term to ensure the dimensionally reduced lattice theory correctly reproduces the physics of the continuum theory: This is gauge invariant since ϕ i (n) transform as site fields. In the absence of this term we observe correlated instabilities in the scalar eigenvalues and in Tr [ϕ i ] for low dimensionless temperatures t 1. This is not the correct behavior required by Kaluza-Klein reduction in the continuum. Instead it corresponds to a center-symmetric phase for the reduced dimensions, which could correspond to Eguchi-Kawai reduction at large N in the presence of adjoint fermions. We avoid this center-symmetric phase by using S center to explicitly break the center symmetry. In this work we use c 2 W = µ 2 for low t 1 and c 2 W = 0 for high temperatures t 1, again extrapolating µ 2 → 0 in the former case.

C. Torus geometries
The geometry of our tori is determined by the skewing parameter γ = −1/2 set by our lattice discretization, and by the aspect ratio α = r L /r β = N x /N t , where N x and N t are respectively the numbers of lattice points generating the spatial and temporal cycles. As mentioned earlier, we will typically discuss results specifying the torus with the skewing γ = −1/2, but this parameterization may represent a modular parameter outside the fundamental domain.
Here we review the geometries we will consider and their fundamental parameterization. In particular, while we use a skewed lattice, some of our geometries in fact are those of rectangular tori when mapped to the fundamental domain.
In Table I we list the lattice sizes N x ×N t we numerically analyze, together with their shape parameter α for skewing γ = −1/2. When the corresponding modular parameter τ doesn't fall in the fundamental domain we give a modular 1 Skewed Table I. The lattice geometries we numerically analyze. Our lattice discretization naturally picks γ = −1/2, and by varying the spatial and temporal lattice extents Nx and Nt we generate tori with different aspect ratios α. When these (α, γ) denote a torus with modular parameter τ outside the fundamental domain, we give an appropriate modular transformation (as in Eq. (13)) so that the equivalent (α , γ ) lie within it. We also give the relation between the temperatures t /t. The last column states whether the torus, viewed from the fundamental domain, is skewed or rectangular. transformation to an equivalent representation with shape parameters α , γ , and note whether the fundamental representation is rectangular or skewed. We also give t /t, the ratio between the dimensionless temperature in the new representation to that of the original. In the corresponding Fig. 3 we plot the complex τ parameters for the various tori in the natural representation where γ = −1/2, and in the cases where these lie outside the fundamental domain we draw an equivalent τ contained in it.

V. NUMERICAL RESULTS
We now discuss our numerical results obtained using the lattice formulation described above. Before studying the low-temperature regime relevant for supergravity, we first consider the high-temperature, small-volume limit and then the phase structure of the theory.

A. High-temperature limit
Fixing the shape of the torus, with constant α for γ = −1/2, we vary t → ∞. Following our earlier discussion in Sec. III D, this is the high-temperature, small-volume limit where we expect the theory to be spatially deconfined with P β , P L = 0, and to have bosonic action density (26).
We investigate three different aspect ratios (α = 1, 4 and 6) in the high-temperature regime and plot the bosonic action density in Fig. 4. Qualitative agreement is seen both in the power of t and the α-dependent coefficient, providing -sBos N 2 λ t SU(9), 8nt8 SU(12), 8nt8 SU(9), 24nt6 SU(12), 24nt6 SU(9), 24nt4 SU(12), 24nt4 HT for α = 1 HT for α = 4 HT for α = 6 Figure 4. Bosonic action density versus dimensionless temperature t for three aspect ratios α = 1, 4 and 6 (from top to bottom), considering gauge groups SU(9) and SU (12). The temperature range probed here corresponds to the high-temperature, smallvolume limit, and the prediction (26) for the behavior is given by the dashed curves marked HT. a test of the dimensional reduction that relates the lattice coupling λ lat to the dimensionless continuum parameters r L and r β .
In Fig. 5 we show distributions of the phases of the N eigenvalues of spatial Wilson lines Pe i L A on 24 × 4 lattices (α = 6) at a high temperature t ≈ 11.4, for SU(N ) gauge groups with N = 6, 9 and 12. The phases are measured relative to the average phase of each Wilson line. In order to compute the usual Wilson lines from the complexified gauge links U a of the lattice formulation, we use a polar decomposition U a = H a · U a to separate each link into a positive-semidefinite hermitian matrix H a (containing the scalar fields) and a unitary matrix U a corresponding to the gauge field. To compute the Wilson lines we simply multiply the unitary matrices, We construct P L and P β by taking the trace (normalized to 1), averaging over lattice sites in the temporal and spatial direction (respectively), and then computing the ensemble average of the magnitude. The expectation that P L ∼ 1 implies we should expect a localized distribution of the phases of the spatial Wilson line eigenvalues, which is consistent with the results in Fig. 5. The distributions show little dependence on N , though the N = 6 case has a few outliers with large fluctuations from the average phase. As t decreases we expect a transition with P L → 0, with the eigenvalue distribution spreading over the angular circle and becoming uniform on it. For t 9 we do indeed see the distributions spread out over the full angular period, as we discuss in more detail below.

B. Phase structure of the SYM theory
We have explored the phase structure of the SYM theory by scanning in t = 1/r β for fixed aspect ratio α = r L /r β . From our previous discussion we expect the theory to be thermally deconfined, but to have an interesting phase  structure associated with spatial confinement. Our numerical results for the temporal Wilson loop magnitude P β appear consistent with the theory being thermally deconfined. We now focus on the spatial Wilson line and order parameter P L .
In Fig. 6 we show the jackknife average magnitude of the Wilson line P L vs. r β for α = 4, along with the corresponding susceptibility The results indicate a large-N transition at t c = 4.6(2) separating a spatially deconfined phase with P L = 0 at small r β (high temperatures) from a spatially confined phase at large r β (low temperatures) where P L → 0 as N → ∞. This transition strengthens with larger N , while the general agreement between results from 16 × 4 and 24 × 6 lattices indicates that discretization effects are small. As discussed in Sec. IV C (Table I), the geometry α = 4 for γ = −1/2 is equivalent to a rectangular (γ = 0) torus with α = 2 √ 3 and r β = r β . In Fig. 7 we plot distributions of the spatial Wilson line eigenvalue phases, following the same procedure as described for Fig. 5 while considering a lower temperature t ≈ 3.8 on 24 × 6 lattices (α = 4). Since t < t c we expect to be in a spatially confined phase, with P L → 0 and correspondingly a uniform density of eigenvalue phases on the angular circle. As expected, we do observe these phases spreading out around the angular period in Fig. 7, and the distribution becomes more uniform as N increases. This contrasts with the localized distributions in Fig. 5 for the high-temperature spatially deconfined phase.
Using the Wilson line susceptibility we have mapped the position of the spatial deconfinement phase transition as a function of α for our γ = −1/2. In Fig. 8 we plot our results on the r L -r β plane and compare them with the  Fig. 6). The three symbols mark the ensembles whose Wilson line eigenvalue phase distributions we show in Figs. 5, 7 and 13 (from bottom to top; the point for Fig. 10 lies outside the range of the plot). The solid lines show the expected transitions sketched in Fig. 1: the BQM deconfinement transition at high temperature and the low-temperature Gregory-Laflamme transition. The dashed lines indicate constant aspect ratios 1/2 ≤ α ≤ 4 from top to bottom. expected transitions sketched in Fig. 1. The error bars in this figure are the full widths at half maximum of the susceptibility peaks like those shown in the right panel of Fig. 6, considering gauge group SU(12) and a fixed lattice size for each aspect ratio α. We have checked that alternate determinations of the transition produce consistent results. These include identifying the transition as the r β for which P L = 0.5, motivated by Ref. [8], and using a large-N generalization [81] of the separatrix introduced for SU(3) by Ref. [82].
For α 4 we find the transition occurs at high temperatures r β 1, and nicely agrees with the deconfinement transition behavior predicted by the high-temperature BQM limit we discussed in Sec. III A. Unfortunately, the error bars increase significantly as we approach transitions occurring at lower temperatures (and smaller α) where we expect the dual gravitational prediction (37) to apply. We do not obtain usable susceptibility peaks for α ≤ 1 and cannot conclusively determine the order of the transition for any α. Given their large uncertainties, our results are certainly consistent with the low-temperature behavior predicted by holography, though we are not yet able to test this prediction with great accuracy. There are also systematic uncertainties to be considered. We expect the most significant systematic effects to be those arising from working with gauge group SU(12) and a fixed lattice size for each aspect ratio α rather than extrapolating to the large-N continuum limit. While the mild discretization artifacts and N dependence shown in Fig. 6 give us confidence that these systematic effects will not qualitatively change our results, we are not yet in a position to meaningfully quantify them. Carrying out controlled extrapolations to the limits of N → ∞ and infinitely large lattices is a central goal for future generations of lattice calculations.
To summarize, our numerical results for the phase diagram of the two-dimensional SYM system are consistent with the expectations from holography. We see a phase where the eigenvalues of the spatial Wilson line are uniformly distributed around the unit circle, as expected for a spatially confined phase. This is separated by a transition from a spatially deconfined phase with localized eigenvalue distribution. We now study the thermodynamics of the system at low temperatures in both of these phases, comparing to the predictions from the dual gravity theory.

C. D1-phase thermodynamics
In Fig. 9 we show the bosonic action density versus t for α = 2 lattice sizes 16 × 8 and 24 × 12 with gauge groups SU (12) and SU (16). From the discussion above, the temperature range shown in the figure lies below the spatial deconfinement transition temperature t c 1.2. Hence we expect the system should be spatially confined. In the low-temperature limit t → 0, we expect it to be described by the gravity D1 phase given in Eq. (33). 13 This is corroborated by our lattice results, which lie close to the D1-gravity prediction at low temperatures, t 0.4.
In a manner similar to the well-studied case of supersymmetric quantum mechanics, this system becomes unstable for very low temperatures, and hence our calculations do not extend all the way down to t 1/N . The origin  Figure 9. Bosonic action density versus dimensionless temperature t for 16 × 8 and 24 × 12 lattices (aspect ratio α = 2) with gauge groups SU (12) and SU (16). All points are results of µ 2 → 0 extrapolations. For sufficiently small t our results are consistent with the D1-phase gravity prediction, without significant sensitivity to N or the lattice size. of this instability is well understood and is discussed in Ref. [80]. The scalar potential terms (50) and (54) help stabilize our numerical calculations, but still do not allow access to arbitrarily low temperatures. We then need to extrapolate µ 2 → 0 to remove the soft-Q-breaking effects of these terms (with c 2 W = µ 2 ), which produces the results shown in Fig. 9. Simple linear µ 2 → 0 extrapolations generally have good quality with small χ 2 /d.o.f.; we include a representative example in Appendix C (Fig. 15). Figure 10 shows extended distributions for the Wilson line eigenvalue phases on 24 × 12 lattices at t ≈ 0.33 (with µ 2 ≈ 0.007), which become more uniform as N increases. This behavior supports our conclusion that the system is spatially confined in this region of the phase diagram, corresponding to the gravity D1 phase.
Next, in Fig. 11 we show our corresponding bosonic action density results for α = 1 lattice sizes 8 × 8 and 16 × 16. We compare our results to the D1-phase gravity prediction, and see reasonable agreement for t 0.4. At slightly larger t 0.5 we observe unusual sensitivity to the lattice size and to N . We suspect that the uncertainties on some of these µ 2 → 0 extrapolated results may be underestimated. Figure 8 suggests that for α = 1 the deconfinement transition occurs around t c 0.47, and at the nearby t ≈ 0.46 we observe unusually long autocorrelations for the SU(12) 16 × 16 ensembles with larger µ 2 0.003. Some of the SU(12) 16 × 16 µ 2 → 0 extrapolations around t 0.5 also produce unusually large 2 χ 2 /d.o.f. 4. Future calculations with larger N and larger lattice sizes, as needed to quantify systematic uncertainties and enable controlled extrapolations to the large-N continuum limit, should clarify this behavior. Since the deconfinement transition occurs at quite large t, we do not expect our results for t > t c to be well described by the gravity D0-phase behavior, and indeed they are not. In order to see the gravity D0-phase behavior emerge at low temperature we need to consider smaller α < 1 so that the transition occurs at t c 1.  (16), 8nt8 SU(9), 16nt16 SU(12), 16nt16 D1 gravity prediction Figure 11. Bosonic action density versus dimensionless temperature t for 8 × 8 and 16 × 16 lattices (aspect ratio α = 1) with gauge groups SU(9), SU (12) and SU (16). All points are results of µ 2 → 0 extrapolations and for low temperatures t < 0.4 our results are in reasonable agreement with the D1-phase gravity prediction (solid curve). The points around t 0.5 show unusual sensitivity to the lattice size and to N . This may be related to the nearby transition around tc 0.47 (cf. Fig. 8), which might lead to underestimated uncertainties.  (16) , 8nt16 D0 gravity prediction Figure 12. Bosonic action density versus dimensionless temperature t for 6 × 12 and 8 × 16 lattices (aspect ratio α = 1/2) with gauge groups SU(9), SU (12) and SU (16). All points are results of µ 2 → 0 extrapolations. For sufficiently small t our results are consistent with the D0-phase gravity prediction.

D. D0-phase thermodynamics
The final numerical results we present consider our smallest aspect ratio α = 1/2. Recall from Table I that this lattice geometry is actually equivalent to a rectangular (γ = 0) torus with α = 1/ √ 3, so that r β = √ 3r β /2. In Fig. 12 we plot the bosonic action density for N = 9, 12 and 16 from 6 × 12 and 8 × 16 lattices. Again the points shown are results of µ 2 → 0 extrapolations. From Fig. 8 we expect the system to be spatially deconfined for the low temperature range 0.25 t 0.5 shown here. Eventually at very low temperatures, presumably around t 0.12, it should confine, but we are not yet able to probe such a low temperature regime. The dashed curve is the low-temperature gravity prediction from the (spatially deconfined) D0 phase, Eq. (38), which is indeed consistent with the data for t 0.35. Figure 13 shows intermediate distributions for the Wilson line eigenvalue phases on 8 × 16 lattices at t ≈ 0.46 (with µ 2 ≈ 0.004), which become more compact as N increases. This behavior supports our conclusion that the system is spatially deconfined in this region of the phase diagram, consistent with the dual gravity approaching the D0 phase in the large-N limit over this temperature range.

VI. CONCLUSIONS
We have studied two-dimensional SYM with maximal supersymmetry compactified on a flat but skewed torus in which an anti-periodic boundary condition is imposed on the fermion fields wrapping one of the cycles. The theory contains three dimensionless parameters: r L , r β = 1/t and the skewing angle cos θ = γ. From the holographic conjecture, at low 'generalized temperature' t 1 this theory should give a description of a dual gravitational system containing various types of black holes arising in Type IIA and IIB supergravity. The phase diagram of the gravitational system is expected to contain a region where homogeneous D1 (black string) solutions dominate and another in which localized D0 (black hole) solutions dominate. The critical line separating these two regions in the dual gravitational system is conjectured to be dual to a first-order deconfinement transition with the spatial Wilson loop magnitude P L serving as an order parameter.
We use lattice gauge theory to explore and test this holographic conjecture using a recently constructed lattice action based on a formalism that maintains an exact supersymmetry at non-zero lattice spacing. The construction singles out a particular skewing angle γ = −1/2, which allows us to test holography both for the usual rectangular tori and also-for the first time-for skewed tori as well.
We have mapped out the phase diagram of the SYM system and indeed find a line of transitions separating a spatially confined phase from a deconfined one. The parametric form of this phase boundary agrees with the results from the gravity dual. Furthermore, the action density computed in either phase is consistent at low temperatures with the corresponding black hole thermodynamics. Our results can be seen as the first step in checking the predictions of gauge/gravity duality in a situation that is distinct from SYM quantum mechanics and has a more subtle phase structure. However we expect it will be a considerable technical challenge to reach for this two-dimensional theory the degree of precision that is now state of the art in the quantum mechanical case, featuring fully quantified systematic uncertainties and controlled extrapolations to the large-N continuum limit.
where a, b, c, d ∈ Z and ad − bc = 1, which corresponds to an element a b c d ∈ SL(2, Z). This group is generated by S, T where S(z) = −1/z and T (z) = z + 1. The fundamental domain for this action on τ is In fact we may also take the generators of the group to be R and T , where R(z) = z/(z + 1), since S = T −1 RT −1 .
In this work we are interested in the subset G of the modular group G std that leaves our fermion boundary conditions invariant, namely the above with a ∈ 2Z, b ∈ 2Z − 1 and c ∈ Z. It is easy to see that G is a subgroup of G std . The fundamental domain for the action of this new group, G, on H is It is generated by the transformations R(z) and U (z) = T 2 (z) = z + 2, with R as above and U corresponding to the SL(2, Z) matrices 1 2 0 1 . We now prove these statements using a simple adaptation of Serre's arguments concerning the domain and generators of the usual modular group [83]. Following Serre we firstly consider the subgroup G of G generated by R and U , and later show this is in fact the group G.
Proposition: For any z ∈ H there exists some g ∈ G such that gz ∈ D.
Proof. Let z ∈ H and g be an element of the group G (which is a subgroup of G). We note that and for integers c, d this implies that there exists a g which maximizes Im(gz). Taking such a g, then we may choose an integer n so that z = U n gz has real part between ±1. In fact z ∈ D as we may see by considering Rz and R −1 z . For any w ∈ H we have Thus if z ∈ D, so that either |z + 1| < 1 or |z − 1| < 1, then g = RU n g and g = R −1 U n g are also elements of G , but either Im(g z) > Im(gz) or Im(g z) > Im(gz). However Im(gz) was assumed to be maximal, hence we conclude that z ∈ D.
Appendix B: Scalar example: an A * 4 lattice and its dimensional reduction Here we consider a scalar field theory in four dimensions discretized on an A * 4 lattice. We use this to illustrate explicitly the reduction of such a lattice theory to a two-dimensional theory on an A * 2 lattice. We take a discretization analogous to that we use for N = 4 SYM, having the lattice action with the lattice derivative taken as The lattice variable φ lives at lattice sites n ∈ Z 4 and we have included an interaction term and coupling, normalized in analogy with a gauge coupling. Here the vectors µ (a) have componentsμ α (ν) = δ α ν for α, ν = 1, . . . , 4 andμ α (5) = (−1, −1, −1, −1). The kinetic term is differenced symmetrically in these five directions.
Consider continuum coordinates y µ = ∆ n µ with ∆ the scale setting the lattice size (proportional to the lattice spacing). Taking the continuum limit ∆ → 0 and assuming a suitably smooth scalar field φ we may expand (B3) Then, using n∈Z 4 f (n) d 4 y f (y) for a (suitably smooth) function f in the ∆ → 0 limit, and defining a continuum field Φ = 1 ∆ φ, the lattice action has the continuum limit S lat → S 4-cont , where Here the components of the metric are g µν = δ µν − 1 5 , and |g| = det g µν = 1 5 . The continuum coupling λ 4 is related to the lattice coupling by λ 4 = λ lat / √ 5. We may change to canonical flat-space coordinates Then and the lattice sites are located at x α = ∆ n µ e α (µ) for n ∈ Z 4 . Recognizing these e (µ) as basis vectors for the A * 4 lattice, and noting that Eq. (B1) includes a difference also in the direction e (5) = − 4 µ=1 e (µ) , we see our original lattice theory is indeed defined on an A * 4 lattice. While from the explicit coordinate presentation above it isn't obvious, this lattice is maximally symmetric as we should expect from Eq. (B1), e (a) · e (b) = δ ab − 1 5 a, b ∈ 1, . . . , 5, and the lattice spacing a 4 along all five directions e (a) is a 4 = ∆ 4/5. Now suppose we are interested in reducing this four-dimensional lattice theory to two dimensions. We take the lattice variables to be independent of the e (3,4) lattice directions, and restrict the lattice sum n∈Z 4 only to a twodimensional slice of the original four-dimensional lattice, n = (n 1 , n 2 , 0, 0) with (n 1 , n 2 ) ∈ Z 2 . The lattice action becomes S (red) lat = 1 λ lat (n 1 ,n 2 )∈Z 2 2 a=1 φ n + µ (a) − φ(n) 2 + φ n − µ (1) − µ (2) − φ(n) 2 + φ(n) 4 (B8) As above, taking continuum coordinates y i = ∆ n i with i = 1, 2, and the continuum limit ∆ → 0, for a suitably smooth function f we have (n 1 ,n 2 )∈Z 2 f (n) 1 ∆ 2 d 2 y f (y). The lattice action has continuum limit with the metric components h ij = δ ij − 1 3 so that |h| = 1 3 and the two-dimensional continuum coupling λ 2 = λ lat /(∆ 2 √ 3). We may again move to canonical flat-space coordinatesx m = y iẽm (i) with m = 1, 2, wherẽ Then and the two-dimensional lattice sites are now located atx m = ∆ n iẽm (i) for n ∈ Z 2 . 14 Now definingẽ (3) = −ẽ (1) −ẽ (2) , we see the reduced lattice action, Eq. (B9), is defined using differences generated by ẽ (1) ,ẽ (2) ,ẽ (3) . We recognize this as an A * 2 lattice action, noting that e (a) ·ẽ (b) = δ ab − 1 3 a, b ∈ 1, . . . , 3, and the lattice spacing a 2 along these three directionsẽ (a) is a 2 = ∆ 2/3. Finally we note that the two-dimensional continuum action (B11) is simply the dimensional reduction of the fourdimensional action (B6). More precisely, if we consider the four-dimensional continuum action (B4) and take the y 3 and y 4 coordinates to be periodic with period ∆, so that y 3,4 ∼ y 3,4 + ∆, then for ∆ → 0 we may Kaluza-Klein reduce on these directions. The zero modes determine the reduced two-dimensional action, which after integrating over the y 3,4 directions gives precisely the two-dimensional action (B9), with the two-and four-dimensional couplings related by The factor 5/3 arises as the Kaluza-Klein reduction is over small periodic directions that are not orthogonal to each other or to the extended y i directions.

Appendix C: Numerical details
We use the standard rational hybrid Monte Carlo (RHMC) algorithm [84] implemented in the publicly available parallel software described by Ref. [57]. 15 In the course of this work we have improved this software to enable the SU(N ) truncation of the gauge links discussed near the end of Sec. IV A, as well as to add the scalar potential terms in Eqs. (50) and (54). These additions, along with related improvements to the large-N performance of the code and other advances, will soon be presented in another publication [85].
The results presented in the body of this paper involve eight aspect ratios α = r L /r β = 1/2, 1, 3/2, 2, 8/3, 4, 6 and 8, investigated for up to five SU(N ) gauge groups with N = 3, 6, 9, 12 and 16. In order to ensure that the soft-Q-breaking scalar potential terms (50) and (54) introduce only small effects, we require that µ 2 , c 2 W λ lat . Specifically, as λ lat varies we fix the ratio µ 2 /λ lat = 0.01, 0.02, 0.03 or 0.04, with either c 2 W = 0 or c 2 W = µ 2 . In Fig. 14 we show violations of a Q Ward identity (discussed in detail by Refs. [56,61]), which are small (percent-level at most) and decrease roughly ∝ 1/N 2 with fixed µ 2 /λ lat = c 2 W /λ lat = 0.01. This establishes that both the scalar potential and the SU(N ) truncation of the gauge links have insignificant numerical effects for sufficiently large N .  Fig. 11. The RHMC algorithm treats the factor of e −S in the partition function (19) as a Boltzmann weight, requiring that the Euclidean action S be real and non-negative. However, gaussian integration over the fermion fields of N = (8,8) SYM produces a pfaffian that is potentially complex, where D is the fermion operator. Therefore all our numerical calculations 'quench' the phase e iφ → 1 [57].  Table II. Tests of pfaffian phase fluctuations for SU(3) ensembles, considering three aspect ratios α and a range of temperatures 0.38 ≤ t ≤ 4.56. In all cases 1 − cos φ 1 corresponds to very small fluctuations in the phase itself, e iφ pq ≈ 1 so that phase reweighting has no practical effect.
Previous lattice studies of N = (2, 2) and N = (8, 8) SYM theories in two dimensions found e iφ pq ≈ 1 even at nonzero lattice spacing, with deviations from unity vanishing rapidly upon approach to the continuum limit [47,48,86,87]. Since we use a different lattice action than those considered previously, for a few SU(3) ensembles we have checked that this remains true in our current work. Table II collects the results of these tests, considering 1 − cos φ to ensure that positive and negative fluctuations in φ cannot cancel out. In all cases we find 1 − cos φ 1, corresponding to e iφ pq close enough to unity that reweighting has no practical effect.