1 Introduction

Supersymmetry plays several prominent roles in modern theoretical physics. It is a valuable tool to improve our understanding of quantum field theory (QFT), an ingredient in many new physics models, and even a means to study quantum gravity via holographic duality. Lattice field theory provides a non-perturbative regularization for QFTs, which has had enormous success as a means to analyze QCD and similar vector-like gauge theories. It is therefore natural to explore how lattice field theory can be applied to investigate supersymmetric QFTs, especially in strongly coupled regimes.

In this short review I briefly summarize the recent progress and near-future prospects of lattice studies of supersymmetric systems. This is an update and expansion of Ref. [1], incorporating pedagogical material presented at the 2021 online program “Nonperturbative and Numerical Approaches to Quantum Gravity, String Theory and Holography”, run by the International Centre for Theoretical Sciences in Bengaluru. To connect to the subject of holography featured by this Special Topics issue, I focus on four-dimensional gauge theories and their dimensional reductions to \(d < 4\). (See Refs. [2, 3] for reviews of theories without gauge invariance, such as Wess–Zumino models and sigma models. More recent work in this area includes Refs. [4,5,6,7,8,9].)

Lattice supersymmetry has been investigated for more than four decades [10], and previously reviewed by Refs. [1,2,3, 11,12,13,14], among others. Unfortunately, progress in this area has been slower than for QCD-like theories, primarily because the lattice regularization of QFTs breaks supersymmetry. This occurs in three main ways. First, the super-Poincaré algebra includes the anti-commutation relation \(\left\{ Q_{\alpha }, {\overline{Q}} _{{\dot{\alpha }} }\right\} = 2\sigma _{\alpha {\dot{\alpha }} }^{\mu } P_{\mu }\) that connects the spinorial generators of supersymmetry transformations, \(Q_{\alpha }\) and \({\overline{Q}} _{{\dot{\alpha }} }\), to the generator of infinitesimal space-time translations, \(P_{\mu }\). Lattice regularization formulates the QFT of interest in a discrete space-time where no such infinitesimal translations exist, implying broken supersymmetry.

Second, bosonic and fermionic fields are typically discretized differently on the lattice. In the context of supersymmetric gauge theories, standard discretizations associate the gauginos \(\lambda _{\alpha }(n)\) with lattice sites while the gauge connections are associated with links \(U_{\mu }(n)\) between nearest-neighbor sites. That is, under a lattice gauge transformation \(\lambda _{\alpha }(n) \rightarrow G(n) \lambda _{\alpha }(n) G^{\dag } (n)\) while \(U_{\mu }(n) \rightarrow G(n) U_{\mu }(n) G^{\dag } (n + a\widehat{\mu } )\), where ‘a’ is the lattice spacing. Away from the \(a \rightarrow 0\) continuum limit, these differences prevent supersymmetry transformations from correctly interchanging superpartners. Although scalar fields also tend to be associated with lattice sites, their discretization typically omits the features (e.g., a Wilson term or staggering) required to address the famous fermion doubling problem, resulting in a similar breaking of supersymmetry.

Finally, the Leibniz rule \(\partial \left[ \phi \eta \right] = \left[ \partial \phi \right] \eta + \phi \partial \eta\) plays an important role in ensuring supersymmetry, but is violated by standard lattice finite-difference operators [10]. In discrete space-time, ‘no-go theorems’ presented by Refs. [15, 16] establish that only non-local derivative and product operators can obey the Leibniz rule and hence fully preserve supersymmetry. This implies a trade-off between locality and supersymmetry, which continues to be explored for simple systems such as supersymmetric quantum mechanics, where a lattice field product operator obeying a ‘cyclic Leibniz rule’ suffices to preserve partial supersymmetry and establish non-renormalization [17,18,19,20]. A different non-local ‘star product’ is able to satisfy the Leibniz rule, but in such a way that the lattice spacing no longer acts as a regulator [21]. Even restricted to simple systems without gauge invariance, and mostly in 0+1 dimensions, these constructions already become very complicated, and the remainder of this review will focus on approaches that preserve locality at the expense of broken supersymmetry.

The breaking of supersymmetry in lattice calculations has the consequence that quantum effects generate supersymmetry-violating operators. Of particular concern are relevant supersymmetry-violating operators, for which counterterms have to be fine-tuned in order to recover the supersymmetric QFT of interest in the \(a \rightarrow 0\) continuum limit that corresponds to removing the UV cutoff \(a^{-1}\). Considering the case of four space-time dimensions, many relevant operators can appear for theories that involve scalar fields. Such theories include supersymmetric QCD (superQCD) with scalar squarks, as well as gauge theories with ‘\({\mathcal {N}} > 1\)’ extended supersymmetry, which include scalar fields in the gauge supermultiplet. The mass terms of these scalars introduce fine-tuning problems similar to that of the standard model Higgs, and additional relevant supersymmetry-violating operators can arise from the fermion (quark and gaugino) mass terms, Yukawa couplings, and quartic (four-scalar) terms. Careful counting typically finds \({\mathcal {O}} (10)\) relevant operators for lattice discretizations of supersymmetric theories with scalar fields [11, 22, 23]. Simultaneously fine-tuning counterterms for all of these operators in numerical lattice calculations appears impractical, to say the least.

In order to make lattice supersymmetry practical, the amount of fine-tuning needs to be reduced. The following three sections briefly review three different ways to achieve this for lattice studies of supersymmetric Yang–Mills (SYM) theories. First, the next section discusses dimensional reductions of SYM theories to fewer than four space-time dimensions, which has been the focus of a great deal of recent work. Returning to four dimensions, Sect. 3 considers the special case of minimal (\({\mathcal {N}} = 1\)) SYM, which is significantly simplified by the absence of scalar fields. Maximal (\({\mathcal {N}} = 4\)) SYM is another special case in four dimensions, for which fine-tuning can be vastly reduced by preserving a closed subalgebra of the supersymmetries, as discussed in Sect. 4. We conclude in Sect. 5 by highlighting some prominent challenges to be faced by future lattice studies of supersymmetric QFTs, including superQCD and sign problems that can arise in several contexts.

2 Dimensionally reduced SYM theories

Working in a smaller number of space-time dimensions, \(d < 4\), can make numerical analyses much more tractable. In addition to the smaller number of degrees of freedom corresponding to \(L^d\) lattices, lower-dimensional theories tend to be super-renormalizable and in many cases a one-loop counterterm calculation suffices to restore supersymmetry in the continuum limit [24,25,26]. Focusing here on dimensional reductions of SYM theories, we will label systems by the number Q of supersymmetry generators, or ‘supercharges’, that they have. Starting in four dimensions, \({\mathcal {N}} = 1\), 2 or 4 SYM corresponds to \(Q = 4\), 8 or 16, respectively. These theories can also be considered dimensional reductions of minimal SYM in respectively \(D_* = 4\), 6 or 10 dimensions. For \(d \le 4\), these theories involve a d-component gauge field, fermionic fields with Q total components, and \(D_* - d\) real scalar fields. All fields are massless and transform in the adjoint representation of the gauge group, here taken to be either SU(N) or \(\text {U(}N\text {)} = \text {SU(}N\text {)} \otimes \text {U(}1\text {)}\).

2.1 0 + 1 Dimensions

Dimensional reduction all the way to (0+1)-dimensional ‘SYM quantum mechanics’ (QM) has been the subject of many numerical studies over the past fifteen years, starting with Refs. [27, 28] investigating the \(Q = 4\) case and also including Refs [29,30,31,32,33,34,35,36,37,38,39,40,41,42,43,44,45,46,47,48,49,50,51,52,53]. These SYM QM systems consist of balanced collections of interacting bosonic and fermionic \(N\!\times \! N\) matrices evolving in (euclidean) time at a single spatial point. They are simple enough that lattice regularization may not even be required to analyze them—Refs. [27, 29,30,31, 33, 36, 37] instead employ a gauge-fixed Monte Carlo approach with a hard momentum cutoff, Ref. [9] exactly diagonalizes truncated hamiltonians, and Refs. [7, 8, 50, 54, 55] explore prospects for quantum computing. Another aspect of this simplicity is the proposal that \(Q = 16\) SYM QM can be ‘ungauged’ to produce a scalar–fermion system with SU(N) global symmetry, with both the gauged and ungauged models flowing to the same theory in the IR [48, 53, 56].

Even though SYM QM systems are much simpler to study on the lattice than their four-dimensional SYM counterparts, they remain computationally non-trivial. Let’s consider this in the context of the maximally supersymmetric \(Q = 16\) case, which has attracted particular interest due to its connections to string theory [57]—especially the conjecture by Ref. [58] that the large-N limit of this system describes the strong-coupling (‘M-theory’) limit of type-IIA string theory in light-front coordinates. Another contribution to this Special Topics issue will review this subject in more detail, and earlier reviews from string theory perspectives include Refs. [59, 60]. At finite temperature, this holographic conjecture relates the bosonic action of the deconfined SYM QM system to the internal energy of a dual compactified black hole geometry in eleven-dimensional M-theory.

This quantity is straightforward to compute through numerical Monte Carlo analyses [29, 31, 32, 34, 37, 41, 42, 44, 45], with Ref. [45] representing the state of the art that improves upon earlier results by carrying out controlled extrapolations to the large-N continuum limit. In addition to the role of large N in holography, the absence of any spatial volume means that in these studies the thermodynamic limit itself corresponds to extrapolating \(N \rightarrow \infty\). These controlled extrapolations enable more robust comparisons to dual gravitational predictions, with numerical results convincingly approaching the leading-order gravitational prediction from classical supergravity at low temperatures, providing non-perturbative first-principles evidence that the holographic duality conjecture is correct. In addition, deviations between the lattice results and leading-order supergravity at higher temperatures can be considered a prediction of higher-order quantum gravitational effects that are enormously difficult to calculate analytically.

The main computational challenge of these SYM QM investigations comes from the large numbers of colors N and lattice sizes L that are needed to control the large-N continuum extrapolations. In this context, with a fixed dimensionless temperature \(\widehat{T} \equiv T / \lambda ^{1 / 3}\), the continuum limit corresponds to \(L \rightarrow \infty\). Here the \(~\widehat{\cdot }~\) decoration highlights dimensionless ratios that can be considered consistently in both the lattice and continuum theories—note the ’t Hooft coupling \(\lambda = g_{\text {YM}}^2 N\) has dimension \([\lambda ] = 4 - d\) in d dimensions. Ref. [45] employs \(16 \le N \le 32\) and lattice sizes up to \(L = 64\), with the \({\mathcal {O}} \!\left( N^3\right)\) cost scaling of matrix–matrix multiplication dominating over the \(\sim\) \(L^{5d / 4}\) cost scaling of the rational hybrid Monte Carlo (RHMC) algorithm, and requiring large-scale supercomputing.

Large values of N are also motivated by a thermal instability associated with the non-compact quantum moduli space of \(Q = 16\) SYM QM [34]. At low temperatures the system is able to run away along these flat directions, which is interpreted holographically as D0-brane radiation from the dual black hole. A formal solution is to stabilize the desired vacuum by adding a supersymmetry-breaking scalar potential to the lattice action, which then needs to be removed in the course of extrapolating to the continuum limit [34, 61]. To avoid this complication, Ref. [45] argues that in practice it can be possible to carry out Monte Carlo sampling around a metastable vacuum so long as N is large enough to reduce the tunneling rate to the true run-away vacuum. The necessary value of N increases as the temperature decreases.

An alternative is to consider a deformation of \(Q = 16\) SYM QM introduced by Berenstein, Maldacena and Nastase (BMN) [62], which lifts the flat directions mentioned above while preserving all 16 supercharges. This deformation serves as a supersymmetric regulator that doesn’t need to be removed in the continuum limit. It introduces non-zero masses for the 9 scalars and 16 fermions of the theory, explicitly breaking the SO(9) global symmetry associated with the compactified spatial dimensions of \(d = 10\) minimal SYM, \(\text {SO(}9\text {)} \rightarrow \text {SO(}6\text {)} \!\times \! \text {SO(}3\text {)}\). The deformation depends on a dimensionful mass parameter \(\mu\), which can be combined with the ’t Hooft coupling to define a dimensionless coupling \(\widehat{g} \equiv \lambda / \mu ^3\) (not to be confused with the dimensionful Yang–Mills gauge coupling \(g_{\text {YM}}^2 = \lambda / N\)).

This BMN model has been studied numerically by Refs. [35, 38, 49,50,51,52,53], with particular focus on its finite-temperature confinement transition. The critical temperature \(\widehat{T} _c\) of this transition can be predicted both by perturbative calculations in the weak-coupling regime \(\widehat{g} \ll 1\) [63,64,65] as well as by dual supergravity calculations for strong couplings \(\widehat{g} \rightarrow \infty\) with \(N \rightarrow \infty\) and \(\widehat{T} \ll 1\) [66]. The goals of ongoing lattice calculations include both reproducing these limits as well as non-perturbatively connecting them by mapping out the intermediate regime where perturbative and holographic approaches are unreliable. Good progress has been made by Refs. [49, 51, 52], considering different lattice discretizations, numbers of colors N, lattice sizes L, and couplings \(\widehat{g}\) that span several orders of magnitude.

2.2 1 + 1 Dimensions

Moving from (0+1)-dimensional quantum mechanics to QFTs in two dimensions introduces many additional phenomena to explore, while remaining significantly more tractable than numerical studies in four dimensions. Research on supersymmetric lattice gauge theories in \(2 \le d \le 4\) dimensions features both conceptual work constructing clever lattice formulations that minimize fine-tuning in principle, as well as numerical work that exploits these constructions to carry out practical calculations. The most widely applied reformulations are based on approaches known as topological twisting [67,68,69] and orbifolded dimensional deconstruction [70,71,72,73], which are thoroughly reviewed by Ref. [2]. While the concepts and terminology of these approaches differ, in the end they actually produce equivalent constructions [74, 75].

Here we will only briefly summarize the twisting approach. Restricting our considerations to flat space-time, twisting is just a change of variables that organizes linear combinations of the supercharges into completely antisymmetric p-forms \({\mathcal {Q}}\), \({\mathcal {Q}} _{\mu }\), \({\mathcal {Q}} _{\mu \nu }\), etc., with the feature that any twisted-scalar supercharge is nilpotent, \({\mathcal {Q}} ^2 = 0\). This is only possible if there are at least \(2^d\) supercharges in d dimensions, and provides at most \(\lfloor Q / 2^d \rfloor\) nilpotent \({\mathcal {Q}}\). The p-forms transform with integer ‘spin’ under the so-called twisted rotation group \(\text {SO(}d\text {)} _{\text {tw}} \equiv \text{ diag }\left[ \text {SO(}d\text {)} _{\text {euc}} \otimes \text {SO(}d\text {)} _R\right]\), where \(\text {SO(}d\text {)} _{\text {euc}}\) is the Wick-rotated Lorentz group and \(\text {SO(}d\text {)} _R\) comes from the global R-symmetry. This procedure clearly provides a closed supersymmetry subalgebra \(\left\{ {\mathcal {Q}} , {\mathcal {Q}} \right\} = 0\) that can be preserved at non-zero lattice spacing, and ultimately leads to a \({\mathcal {Q}}\)-invariant lattice action with no need of the Leibniz rule.

Consistently with the discussion in Sect. 1, the other twisted supercharges \({\mathcal {Q}} _{\mu }\), \({\mathcal {Q}} _{\mu \nu }\), \(\cdots\), are all broken by the lattice discretization. However, the preservation of a subset of the supersymmetries and a closed subalgebra of the full super-Poincaré algebra vastly reduces the fine-tuning required to recover the other, broken supersymmetries in the continuum limit. We will comment on this in more detail for four-dimensional \({\mathcal {N}} = 4\) SYM in Sect. 4, but it is useful here to consider an analogy with the more basic Poincaré symmetries of translations and rotations. While these continuous symmetries are also broken in discrete lattice space-times, since the early days of lattice field theory it has been appreciated that the discrete symmetries preserved by (e.g.) hypercubic lattices guarantee their recovery in the continuum limit, with no need for fine-tuning. In a similar way, the preserved twisted-scalar supersymmetries \({\mathcal {Q}}\) guarantee the recovery of the full set of supersymmetries in the continuum limit, with little to no fine-tuning.

For some theories there are multiple, inequivalent ways that twisted lattice systems can be constructed. One approach [74, 76,77,78] combines the gauge and scalar fields into a complexified gauge field. This results in \(\text {U(}N\text {)} = \text {SU(}N\text {)} \otimes \text {U(}1\text {)}\) gauge invariance and non-compact lattice gauge links \(\left\{ {\mathcal {U}} , \overline{{\mathcal {U}}} \right\}\) with a flat measure. Although the \(\text {U(}1\text {)}\) sector decouples in the continuum, at non-zero lattice spacing it can introduce unwanted artifacts, especially at strong couplings. Refs. [78,79,80,81] explore various ways these artifacts could be suppressed for four-dimensional \({\mathcal {N}} = 4\) SYM, which we will revisit in Sect. 4. Numerical studies using this formulation for \(Q = 16\) SYM in two dimensions include Refs. [79, 82,83,84].

A different approach [67, 68, 85,86,87,88,89] works with compact gauge links and gauge group SU(N), at the cost of imposing an admissibility condition to resolve a huge degeneracy of vacua. While Ref. [90] proposes formulations of \(Q = 4\) and \(Q = 8\) SYM in two dimensions that avoid the need to impose admissibility conditions, this issue becomes more problematic in higher dimensions [2, 91]. Numerical studies using this formulation in two dimensions include Refs. [61, 92,93,94,95,96,97,98,99,100] investigating \(Q = 4\) SYM and Refs. [101, 102] investigating \(Q = 16\) SYM.

As for the BMN model discussed in Sect. 2.1, a prominent physics target is to map out the non-perturbative phase diagrams of two-dimensional SYM theories [79, 103, 104]. These systems are now formulated on an \(r_L \!\times \! r_{\beta }\) torus, where \(r_{\beta } = 1 / \widehat{T}\) is the inverse dimensionless temperature while \(r_L = L \sqrt{\lambda }\) is the dimensionless length of the spatial cycle. At high temperatures corresponding to small \(r_{\beta }\), the fermions pick up a large thermal mass, reducing the system to bosonic quantum mechanics (BQM) in one dimension. This is illustrated for \(Q = 16\) SYM by the sketch in the left panel of Fig. 1, which highlights the ‘spatial deconfinement’ transition expected as \(r_L\) decreases in the large-N limit. (The system is always thermally deconfined.) In the BQM limit, this transition appears to be first order—see Ref. [105] and references therein.

A similar first-order transition is predicted by holography in the large-N, low-temperature regime corresponding to large \(r_{\beta }\). Holographically, the large-\(r_L\) spatially confined phase is conjectured to be dual to a homogeneous black string with a horizon wrapping around the spatial cycle, while the small-\(r_L\) spatially deconfined phase corresponds to a localized black hole. Challenging numerical supergravity calculations are required to construct and analyze these dual black hole and black string geometries [106].

Fig. 1
figure 1

A sketch (left) of the expected phase diagram for two-dimensional \(Q = 16\) SYM, compared to numerical results (right), both adapted from Ref. [79]. The numerical calculations use a twisted formulation with gauge group SU(12) and aspect ratios \(\alpha = r_L / r_{\beta } = N_s / N_t\) ranging from 8 to 3/2

Full size image

The right panel of Fig. 1 presents lattice results for the \(Q = 16\) SYM spatial deconfinement transition from Ref. [79], again aiming to reproduce the expected high- and low-temperature limits while non-perturbatively mapping out intermediate temperatures. To this end the calculations vary \(r_L\) with a fixed aspect ratio \(\alpha = r_L / r_{\beta } = N_s / N_t\) defined by the \(N_s\!\times \! N_t\) lattice size, monitoring the spatial Wilson line \(\text{ Tr }\left[ \prod _{x_i} {\mathcal {U}} _x(x_i, t) \right]\) as the order parameter for this transition. As shown by the dashed lines in the plot, larger aspect ratios (up to \(\alpha = 8\) for a \(32\!\times \! 4\) lattice) probe the spatial deconfinement transition at higher temperatures, matching the BQM expectations quite well. As the aspect ratio decreases (down to \(\alpha = 3 / 2\) for an \(18\!\times \! 12\) lattice), the lower-temperature numerical results are consistent with holography, albeit with rapidly increasing uncertainties. At these low temperatures, a supersymmetry-breaking scalar potential is added to the lattice action to lift flat directions, and then extrapolated to zero, as discussed in Sect. 2.1. Ref. [79] also calculates the internal energies of the dual black hole and black string, finding consistency with holographic expectations in both phases, again with large uncertainties. Although the work compares multiple lattice volumes and SU(N) gauge groups up to \(N = 16\), room for improvement remains in terms of carrying out controlled extrapolations to the continuum and thermodynamic (i.e., large-N) limits. In particular, larger values of N should help to reduce uncertainties and access even lower temperatures.

Beyond investigations of phase diagrams, two-dimensional SYM theories also possess rich zero-temperature dynamics that are important to explore non-perturbatively. For example, Ref. [107] analyzes the ‘meson’ spectrum of the \(Q = 4\) lattice theory, using a straightforward Wilson-fermion discretization rather than a twisted construction, and observing a massless supermultiplet predicted by Refs. [108, 109]. This work also checks for spontaneous supersymmetry breaking (SSB), which Ref. [110] suggests might occur for this theory. While SSB is also being explored for supersymmetric QM in 0 + 1 dimensions [5,6,7, 111,112,113], this can be a truly dynamical process in two dimensions [114, 115], rather than being determined by the superpotential. We will revisit SSB in the context of the sign problem in Sect. 5.2, for now simply noting that Ref. [107] saw no evidence of SSB for \(Q = 4\) SYM, consistent with other lattice studies using twisted formulations [93, 94, 97, 116].

2.3 2+1 Dimensions

Three-dimensional SYM remains a prominent frontier for lattice studies, with many compelling physics targets and more modest computational costs compared to \(d = 4\). Supersymmetry-preserving twisted formulations (discussed by Refs. [26, 117, 118]) are even more important than in two dimensions, and are used by all numerical calculations so far [119, 120]. Figure 2 presents some results from these works, investigating the \(Q = 16\) SYM bosonic action that holography relates to the internal energy of the dual black brane geometry in supergravity. Here the calculations use \(L^3\) lattice volumes corresponding to an aspect ratio \(\alpha = 1\), and are kept within the spatially confined phase dual to a homogeneous black D2-brane. The results in the left panel of Fig. 2 approach the corresponding supergravity prediction for sufficiently low dimensionless temperatures. While these calculations consider gauge group U(N) with only a single \(N = 8\), one notable advance are the extrapolations to the continuum limit shown in the right panel of Fig. 2, which were not attempted in the \(d = 2\) study that produced Fig. 1. As in lower dimensions, with a fixed dimensionless temperature and \(\alpha = 1\), the continuum limit corresponds to extrapolating \(L \rightarrow \infty\) while the thermodynamic limit would be provided by extrapolating \(N \rightarrow \infty\).

Fig. 2
figure 2

Bosonic action density for three-dimensional \(Q = 16\) SYM with gauge group U(8) and \(L^3\) lattice volumes, from Ref. [120]. Left: Results for \(L = 8\), 12 and 16 vs. the dimensionless temperature \(\widehat{T}\), compared with the high-temperature expectation \(\propto\) \(\widehat{T} ^3\) and the low-temperature dual-supergravity prediction \(\propto\) \(\widehat{T} ^{10 / 3}\). Right: Linear \(L^2 \rightarrow \infty\) continuum extrapolations for the six lowest temperatures

Full size image

Ongoing work is now building on these results to explore the phase transition between this ‘D2 phase’ and the spatially deconfined ‘D0 phase’ dual to a localized black hole geometry. This transition corresponds to considering an \(r_L \!\times \! r_L \!\times \! r_{\beta }\) three-torus, with other behavior possible if the two spatial cycles are allowed to have different sizes [121]. Ref. [120] also presents initial investigations of \(Q = 8\) SYM in three dimensions, for which new code is being developed within the publicly available SUSY LATTICE parallel software package for twisted lattice supersymmetry [77, 122].

3 Minimal \({\mathcal {N}} = 1\) SYM in four dimensions

In Sect. 1 we emphasized that scalar fields are responsible for most of the relevant supersymmetry-violating operators that require fine-tuning in lattice calculations, including the scalar mass terms, Yukawa couplings, and quartic operators. This implies that upon reaching \(d = 4\) dimensions, the most promising supersymmetric gauge theory to consider would be \({\mathcal {N}} = 1\) SYM, the only one with no scalar fields. \({\mathcal {N}} = 1\) SYM involves only the SU(N) gauge field and its superpartner gaugino—a massless Majorana fermion transforming in the adjoint rep of SU(N). The mass term for this gaugino is the only relevant operator that may need to be fine-tuned in order to obtain the correct continuum limit [123, 124]. In fact, even this single fine-tuning can be avoided by using overlap or domain-wall lattice fermion formulations that obey the Ginsparg–Wilson relation and preserve chiral symmetry at non-zero lattice spacing. Even though the axial anomaly breaks the classical \(\text {U(}1\text {)}\) R-symmetry of \({\mathcal {N}} = 1\) SYM to its \(Z_{2N}\) subgroup, preserving this discrete global symmetry still protects the gaugino mass from additive renormalization. The spontaneous breaking of \(Z_{2N} \rightarrow Z_2\) by the formation of a gaugino condensate, \(\left\langle \lambda \lambda \right\rangle \ne 0\), is also relatively straightforward to investigate with this approach [125,126,127,128].

The downside to overlap and domain-wall fermions is that they are computationally expensive, and after Refs. [125,126,127] there has been very little work applying them to \({\mathcal {N}} = 1\) SYM during the past decade [128]. Instead, recent lattice research into this theory has employed improved Wilson fermions, fine-tuning the gaugino mass to recover both chiral symmetry and supersymmetry in the continuum limit. For example, the DESY–Münster–Regensburg–Jena Collaboration made significant progress using clover-improved Wilson fermions [129,130,131,132,133,134,135,136,137,138]. A second group is exploring a SYM analogue of the twisted-mass fermion action [139], with the aim of improving the formation of composite supermultiplets at non-zero gaugino masses and lattice spacings, and thereby gaining better control over the chiral and continuum extrapolations.

With limited connections to holography, and a spatial lattice volume to provide the thermodynamic limit, there is little motivation for extremely large values of N, and most work considers only gauge groups SU(2) [129,130,131, 135, 138] and SU(3) [132,133,134, 136, 137, 139]. A long-term goal of these studies has been to observe the composite states forming degenerate multiplets in the supersymmetric continuum chiral limit. One such multiplet is expected to contain a scalar particle, a pseudoscalar particle, and a fermionic ‘gluino–glue’ particle. Degeneracy in these channels, in the chiral–continuum limit, was convincingly observed by Ref. [136], overcoming numerical challenges that include fermion-line-disconnected contributions to all physical two-point functions, as well as mixing between glueballs and meson-like singlet states.

An additional challenge is carrying out the chiral extrapolations in the presence of an unprotected gaugino mass. This is done by taking the \(m_{\pi }^2 \rightarrow 0\) limit for an ‘adjoint pion’ defined in partially quenched chiral perturbation theory [140], which is measured from just the connected part of the correlator for the \(\eta '\)-like ‘gluinoball’. Supersymmetric Ward identities provide an alternative means to determine the chiral limit [133, 137]. At a non-zero lattice spacing, any differences between these two determinations can be considered a measure of the supersymmetry-breaking discretization artifacts. Ref. [137] finds that these vanish \(\propto\) \(a^2\), as expected for clover-improved Wilson fermions, supporting the restoration of susy in the chiral continuum limit.

Of course, many other lattice \({\mathcal {N}} = 1\) SYM investigations are valuable to carry out in addition to calculations of the gaugino condensate, composite spectrum, and Ward identities. These include explorations of the finite-temperature phase diagram, which features both a chiral transition related to gaugino condensation as well as a confinement transition related to spontaneous center symmetry breaking. Refs. [129, 141] find that these two transitions occur at roughly the same critical temperature, at least for gauge group SU(2), which was not known a priori. In addition, Refs. [130, 134] investigate the phase diagram on \({\mathbb {R}} ^3 \!\times \! S^1\) with a small radius for the compactified temporal direction, comparing thermal and periodic boundary conditions (BCs) for the gauginos. This work finds evidence that periodic BCs allow the confined, chirally broken phase to persist for weak couplings where analytic semi-classical methods [142] may be reliable.

In these investigations, the gradient flow has helped enable the precise measurement of the gaugino condensate with clover-improved Wilson fermions. The gradient flow is also widely used to set the scale and extrapolate to the continuum limit in lattice studies of \({\mathcal {N}} = 1\) SYM [129, 131,132,133, 135,136,137], including the recent Ref. [143], which uses twisted Eguchi–Kawai volume reduction to study the theory on a ‘lattice’ that consists of a single site. Another potential application of the gradient flow is to define renormalized supercurrents and help guide fine-tuning, by constructing a flow that is consistent with supersymmetry in Wess–Zumino gauge [144,145,146]. This complements the ongoing use of lattice perturbation theory to analyze these supercurrents, and other operators [147,148,149]. Finally, given the progress in algorithms and computing hardware, it is also compelling to continue exploring the use of overlap [128] or domain-wall fermions to investigate \({\mathcal {N}} = 1\) SYM.

4 Maximal \({\mathcal {N}} = 4\) SYM in four dimensions

\({\mathcal {N}} = 4\) SYM, with \(Q = 16\) supercharges in four dimensions, turns out to be another special supersymmetric gauge theory that is promising to consider on the lattice. This is serendipitous, given the large role that \({\mathcal {N}} = 4\) SYM plays in theoretical physics thanks to its many supersymmetries, large SU(4)\(_R\) symmetry and conformal symmetry—with additional simplifications in the large-N planar limit. Among many other important applications, it is the conformal field theory of the original AdS/CFT holographic duality [150], and provided early insight into S-duality [151]. This provides many compelling targets for non-perturbative lattice investigations of \({\mathcal {N}} = 4\) SYM to pursue, complementing the many analytic approaches that have already been brought to bear.

Because \({\mathcal {N}} = 4\) SYM features six massless real scalar fields in addition to the gauge field and four massless Majorana fermions, a naive lattice discretization would be problematic. At least 8 fine-tunings would be required if all the supersymmetries were broken [23]. Fortunately, \({\mathcal {N}} = 4\) SYM is the only \(d = 4\) theory with \(Q = 2^d\) large enough to apply the (equivalent) twisted and orbifolded constructions introduced in Sect. 2.2. This allows a single ‘twisted-scalar’ supercharge \({\mathcal {Q}}\) to be preserved, dramatically improving the situation.

Twisted lattice \({\mathcal {N}} = 4\) SYM [2, 72,73,74,75] combines the bosonic fields into five-component complexified gauge links \(\left\{ {\mathcal {U}} _a, \overline{{\mathcal {U}}} _a\right\}\). In order for these five links to symmetrically span four dimensions, the discretization of space-time needs to employ the \(A_4^*\) lattice, which features a large \(S_5\) point-group symmetry. The twisted p-form fermions \(\eta\), \(\psi _a\) and \(\chi _{ab}\) are respectively identified with the sites, links and oriented plaquettes of this \(A_4^*\) lattice. A single fine-tuning of a marginal operator may be required to recover the continuum twisted \(\text {SO(}4\text {)} _{\text {tw}}\) from the discrete \(S_5\) symmetry. The combination of \({\mathcal {Q}}\) and \(\text {SO(}4\text {)} _{\text {tw}}\) then ensures the restoration of the 15 supersymmetries broken by the lattice discretization [23, 152, 153]. Most numerical calculations so far don’t explore this potential fine-tuning, instead fixing the corresponding coefficient to its classical value.

As described in Sect. 2.2, this twisting procedure results in \(\text {U(}N\text {)} = \text {SU(}N\text {)} \otimes \text {U(}1\text {)}\) gauge invariance, and numerical calculations need to regulate flat directions in both the SU(N) and \(\text {U(}1\text {)}\) sectors. Simple supersymmetry-breaking scalar potentials like those used in lower dimensions (and removed in the course of extrapolating to the continuum limit) only affect the SU(N) sector. Regulating the \(\text {U(}1\text {)}\) sector is more challenging, in large part because the corresponding artifacts appear much more severe. Initial studies [76] with a second supersymmetry-breaking potential—involving the determinant of the plaquette as the simplest gauge-invariant quantity sensitive to the \(\text {U(}1\text {)}\) sector—exhibited far larger discretization artifacts than an improved action [78] that instead incorporates this plaquette determinant into a \({\mathcal {Q}}\)-invariant modification of the moduli equations. More recent work has abandoned the \(\text {U(}1\text {)}\) gauge invariance entirely [81], justified by the decoupling of the \(\text {U(}1\text {)}\) sector in the continuum limit. Despite the large \(\text {U(}1\text {)}\) effects indicating non-decoupling in previous lattice studies, Ref. [81] was able to observe apparently reasonable results out to unprecedentedly strong lattice ’t Hooft couplings \(\lambda _{\text {lat}} \le 30\). Moving forward, detailed comparisons with the gauge-invariant improved action would be worthwhile to clarify the role of discretization artifacts and gauge invariance.

Fig. 3
figure 3

Results from ongoing four-dimensional lattice \({\mathcal {N}} = 4\) SYM calculations. Left: The static potential Coulomb coefficient \(C(\lambda _{\text {lat}} )\) is consistent with leading-order perturbation theory (black dashed line) for \(\lambda _{\text {lat}} \le 2\), comparing U(N) gauge groups with \(2 \le N \le 4\) and \(L^3\!\times \! N_t\) lattice volumes with \(L \le 16\) and \(N_t \le 32\). Right: A scale-dependent effective anomalous dimension \(\gamma _{\text {eff}} (\Omega ^2)\) approaches the expected \(\gamma _*(\lambda ) = 0\) in the IR limit \(\Omega ^2 \rightarrow 0\). Considering \(16^4\) lattices for gauge group \(\text {U(}2\text {)}\), discretization artifacts become more significant as the coupling \(\lambda _{\text {lat}}\) increases

Full size image

The publicly available SUSY LATTICE parallel software package implements the improved action for \({\mathcal {N}} = 4\) SYM [77, 122], and is currently being used to study a wide range of interesting observables. For example, the left plot in Fig. 3 considers the static potential V(r), which is correctly seen to be coulombic at all accessible ’t Hooft couplings [76, 81, 154, 155]. By incorporating tree-level improvement into the lattice analyses of the static potential [155], and fitting the resulting data to the Coulomb potential \(V(r) = A - C / r\), we obtain the results for the Coulomb coefficient \(C(\lambda _{\text {lat}} )\) shown in the plot for several U(N) gauge groups and lattice volumes. With \(2 \le N \le 4\) and \(\lambda _{\text {lat}} \le 2\), the lattice results are consistent with the leading-order perturbative relation \(C(\lambda _{\text {lat}} ) \propto \lambda _{\text {lat}}\), which isn’t surprising given that higher-order perturbative corrections are suppressed by powers of \(\frac{\lambda }{2\pi ^2}\) [156,157,158]. In the strong-coupling planar regime \(\lambda \rightarrow \infty\) with \(\lambda \ll N\), there is a famous holographic prediction that \(C(\lambda ) \propto \sqrt{\lambda }\) up to \({\mathcal {O}} \left( \frac{1}{\sqrt{\lambda }}\right)\) corrections [159, 160], and more general analytic results have been obtained in the \(N = \infty\) planar limit [161]. Efforts are underway to search for this behavior, by building on Ref. [81] to access stronger couplings.

As a conformal field theory, \({\mathcal {N}} = 4\) SYM is characterized by its \(\lambda\)-dependent spectrum of scaling dimensions, which are a more challenging target for lattice calculations to predict. The right plot in Fig. 3 considers the analogue of the mass anomalous dimension, extracted from the eigenmode number of the lattice \({\mathcal {N}} = 4\) SYM fermion operator D [162]. Unlike the Dirac operator of the QCD-like theories where this approach was developed [163,164,165], the fermion operator is skew-symmetric, \(\Psi ^T D \Psi = \chi _{ab} {\mathcal {D}} ^{(+)}_{[a} \psi _{b]} + \eta {\mathcal {D}} ^{\dag (-)}_a \psi _a + \frac{1}{2}\epsilon _{abcde} \chi _{ab} {\mathcal {D}} ^{\dag (-)}_c \chi _{de}\), and the corresponding anomalous dimension \(\gamma _*(\lambda ) = 0\) for all couplings. The eigenmode number \(\nu (\Omega ^2)\) counts the number of \(D^{\dag } D\) eigenmodes with eigenvalues \(|\lambda _k|^2 \le \Omega ^2\), and scales \(\nu (\Omega ^2) \propto \left( \Omega ^2\right) ^{2 / (1 + \gamma _*)}\). Stochastically constructing and integrating a Chebyshev expansion of the spectral density [164, 165] produces numerical results for the eigenmode number. Fitting these results to a power-law within windows \(\left[ \Omega ^2, \Omega ^2 + \ell \right]\) of fixed length \(\ell\) provides a scale-dependent effective anomalous dimension \(\gamma _{\text {eff}} (\Omega ^2)\), which approaches \(\gamma _*(\lambda ) = 0\) in the IR limit \(\Omega ^2 \rightarrow 0\). Deviations from zero indicate the scale of lattice artifacts related to the breaking of conformality by the finite lattice volume and non-zero lattice spacing, which become more severe at the ’t Hooft coupling \(\lambda _{\text {lat}}\) increases.

This experience with a trivially vanishing anomalous dimension aids ongoing investigations of the non-trivial scaling dimension \(\Delta _K(\lambda ) = 2 + \gamma _K(\lambda )\) of the simplest conformal primary operator of \({\mathcal {N}} = 4\) SYM, the Konishi operator \({\mathcal {O}} _K = \sum _I \text{ Tr }\left[ \Phi ^I \Phi ^I \right]\), where \(\Phi ^I\) are the six real scalar fields of the theory. As for the static potential Coulomb coefficient, there are predictions for the Konishi scaling dimension from weak-coupling perturbation theory [166,167,168], from holography at strong couplings \(\lambda \rightarrow \infty\) with \(\lambda \ll N\) [169], and for all couplings in the \(N = \infty\) planar limit [170]. In addition, due to the conjectured S-duality of the theory, which predicts an invariant spectrum of anomalous dimensions under the interchange \(\frac{4\pi N}{\lambda } \longleftrightarrow \frac{\lambda }{4\pi N}\), the perturbative results are also relevant in the alternate strong-coupling regime \(\lambda \gg N\) [171]. Finally, the superconformal bootstrap program has been applied to analyze the Konishi anomalous dimension, with initial bounds on the maximum value \(\gamma _K\) could reach across all \(\lambda\) [172, 173] recently being generalized to \(\lambda\)-dependent constraints [174].Footnote 1 Preliminary lattice results with \(\lambda _{\text {lat}} \le 3\) in Ref. [1], obtained from Monte Carlo renormalization group (MCRG) stability matrix analyses, again appear consistent with perturbation theory.

Numerical lattice analyses of \({\mathcal {N}} = 4\) SYM clearly remain in their early stages, with many opportunities for both technical improvements as well as generalizations to other interesting targets. A great deal of effort is currently focused is on accessing stronger ’t Hooft couplings, both to make more direct contact with holographic predictions and also to investigate the behavior of the system around the S-dual point \(\lambda _{\text {sd}} = 4\pi N\). We will revisit this issue in Sect. 5.2. Another direction proposed by Ref. [177] is to adjust the scalar potential so as to study the theory on the Coulomb branch of the moduli space, where its U(N) gauge invariance is higgsed to \(\text {U(}1\text {)} ^N\). In this context S-duality relates the masses of the \(\text {U(}1\text {)}\)-charged elementary ‘W bosons’ and the magnetically charged topological ’t Hooft–Polyakov monopoles [151], each of which may be accessible from lattice calculations with either C-periodic or twisted BCs, for values of \(\lambda _{\text {lat}}\) that have already been studied successfully. The behavior of lattice \({\mathcal {N}} = 4\) SYM at non-zero temperatures will also be interesting to explore. In particular, there is motivation [14] to study the free energy, with the aim of using non-perturbative lattice calculations to connect the weak-coupling perturbative prediction [178] and the strong-coupling holographic calculation [169], which differ by a famous factor of \(\frac{3}{4}\).

5 Challenges for the future

While the recent progress in the three areas discussed above is substantial, many compelling directions remain for further research, ranging from improving large-N continuum extrapolations in lower dimensions, to revisiting \({\mathcal {N}} = 1\) SYM with Ginsparg–Wilson fermions, and pursuing stronger ’t Hooft couplings in \({\mathcal {N}} = 4\) SYM calculations. There are also many other aspects of supersymmetric gauge theories that have proven more challenging to tackle on the lattice. In this section we’ll conclude this brief review by commenting on two particular challenges — lattice analyses of superQCD and the possibility of sign problems in supersymmetric lattice systems.

5.1 Supersymmetric QCD

Considering first \({\mathcal {N}} = 1\) theories in four dimensions, the generalization from SYM to superQCD involves adding ‘matter’ multiplets — ‘quarks’ and ‘squarks’, which could in principle transform in any representation of the gauge group. This generalization would enable investigations of many important phenomena, including (metastable) dynamical supersymmetry breaking, conjectured electric–magnetic dualities and RG flows to known conformal IR fixed points. Of course, the presence of the scalar squarks implies many more relevant supersymmetry-violating operators, making the necessary fine-tuning far more challenging. Even exploiting the continuum-like symmetries preserved by overlap or domain-wall fermions, Ref. [11] counts \({\mathcal {O}} (10)\) operators to be fine-tuned, depending on the gauge group and matter content. The simplifications offered by the Ginsparg–Wilson relation appear particularly crucial in the context of superQCD. In particular, Refs. [11, 22] argue that this good control over the fermions might allow the scalar masses, Yukawas and quartic couplings to be fine-tuned “offline” through multicanonical reweighting, significantly reducing computational costs.

While the recent overlap investigation of \({\mathcal {N}} = 1\) SYM in Ref. [128] provides a starting point for generalization to superQCD, most recent lattice explorations use Wilson fermions and are confronted with the full fine-tuning challenge. One tactic is to proceed by using lattice perturbation theory to guide numerical calculations [147,148,149, 179,180,181]. There is also an initial investigation of a superQCD gradient flow that is consistent with supersymmetry in Wess–Zumino gauge [182]. Another approach is to omit the scalar fields at first, and initially study gauge–fermion theories including both adjoint gauginos and fundamental quarks [183, 184]. A similar strategy is also being applied to \({\mathcal {N}} = 2\) SYM, using overlap fermions [185]. These scalar-less studies also provide useful connections to investigations of near-conformal composite Higgs models, recently reviewed by Refs. [186, 187]. All of these studies in four dimensions remain in their early stages.

As an alternative, following the logic of Sect. 2, it can prove advantageous to investigate superQCD in the simpler setting of fewer than four space-time dimensions. Going all the way down to 0+1 dimensions, for example, Refs. [188,189,190] consider the Berkooz–Douglas matrix model [191], which adds \(N_f\) fundamental multiplets to \(Q = 16\) SYM QM in such a way as to preserve half of the supercharges in the continuum. In \(d = 2\) and \(d = 3\) dimensions, most effort so far has focused on constructing clever lattice formulations of superQCD [192,193,194,195,196,197,198].

Existing numerical calculations [199], and work in progress [120], employ a quiver construction of \(Q = 4\) superQCD in two dimensions [192, 193], based on the twisted formulation introduced in Sect. 2.2. Starting from \(Q = 8\) SYM in three dimensions, the system is reduced to have only two slices in the third direction. The twisted formulation is then generalized to have different gauge groups \(\text {U(}N\text {)}\) and \(\text {U(}F\text {)}\) on each of these two-dimensional slices. Decoupling the \(\text {U(}F\text {)}\) slice then produces a two-dimensional U(N) theory with F massless fundamental matter multiplets and \(Q = 4\) supersymmetries, one of which is preserved at non-zero lattice spacing. This same approach can be applied to construct \(Q = 8\) superQCD in two and three dimensions [196, 197], and may be generalizable to higher representations [198]. The main numerical result so far has been to compare \(\text {U(}2\text {)}\) superQCD with \(F = 3\) against \(\text {U(}3\text {)}\) superQCD with \(F = 2\), observing dynamical supersymmetry breaking when \(N > F\) and confirming that the resulting goldstino is consistent with masslessness in the infinite-volume limit [199].

5.2 Sign problems

Another challenge confronting some supersymmetric lattice field theories is the possibility of sign problems, at least in certain regimes. Ref. [200] provides a brief general introduction to sign problems. Here we will return to focusing on SYM, which involves only Majorana fermions, and hence the pfaffian of the fermion operator D, which can fluctuate in sign even when the determinant would be positive. Separating a generic complex pfaffian into its magnitude and phase, \(\text{ pf }\, D = |\text{ pf }\, D| e^{i\alpha }\), the results presented in Sect. 4 all come from ‘phase-quenched’ RHMC calculations that use only \(|\text{ pf }\, D|\) to carry out importance sampling. The resulting phase-quenched observables \(\left\langle {\mathcal {O}} \right\rangle _{\text {pq}}\) then need to be reweighted, \(\left\langle {\mathcal {O}} \right\rangle = \left\langle {\mathcal {O}} \right\rangle _{\text {pq}} / \left\langle e^{i\alpha } \right\rangle _{\text {pq}}\), with a sign problem appearing when \(\left\langle e^{i\alpha } \right\rangle _{\text {pq}} = Z / Z_{\text {pq}}\) vanishes within statistical uncertainties. In particular, in lattice calculations employing fully periodic BCs, the partition function Z is the Witten index and must vanish for any theory that can exhibit spontaneous supersymmetry breaking [201], implying a severe sign problem.

For lattice \({\mathcal {N}} = 1\) SYM, clover-improved Wilson fermions produce a real pfaffian whose sign can be computed efficiently [202]. Ref. [132] reports that \(\left\langle e^{i\alpha } \right\rangle _{\text {pq}} \approx 1\) in these studies, improving as the lattice spacing decreases. The twisted-mass approach to \({\mathcal {N}} = 1\) SYM allows the pfaffian to be complex. While Ref. [139] finds \(\left\langle e^{i\alpha } \right\rangle _{\text {pq}} \approx \cos (\alpha ) > 0.965\) for \(16^3\!\times \! 32\) lattices (extrapolating from smaller volumes), it also reports that \(\left\langle e^{i\alpha } \right\rangle _{\text {pq}}\) decreases exponentially in the lattice volume V, as expected [200]. This work has to extrapolate from smaller volumes because it directly evaluates the pfaffian, which can be extremely expensive, with computational costs scaling \(\propto\) \(N_{\Psi }^3\), where \(N_{\Psi } \propto N^2 V\) is the number of fermionic degrees of freedom in the lattice system. Such direct evaluations are more common in lower dimensions, where all works so far observe well-behaved \(\left\langle e^{i\alpha } \right\rangle _{\text {pq}} \rightarrow 1\) in the fixed-\(\widehat{T}\) continuum limit [34, 42, 52, 84, 98, 99, 107, 116, 120]. Lower-dimensional systems can also be used as testbeds for different algorithmic approaches, such as the complex Langevin method used to explore spontaneous supersymmetry breaking in Ref. [5].

Fig. 4
figure 4

Results for the pfaffian phase \(\left\langle \text{ Re } \left( e^{i\alpha }\right) \right\rangle _{\text {pq}} \approx \left\langle e^{i\alpha } \right\rangle _{\text {pq}}\) from four-dimensional lattice \({\mathcal {N}} = 4\) SYM calculations using the gauge-invariant improved action [153, 203]. Left: For a fixed ’t Hooft coupling \(\lambda _{\text {lat}} = 0.5\), only per-mille-level fluctuations are observed for U(N) gauge groups with \(N = 2\), 3 and 4, up to the largest accessible volumes. Right: For a fixed \(4^4\) lattice volume, the \(\text {U(}2\text {)}\) phase fluctuations increase significantly for stronger couplings, which obstructs studies of \(\lambda _{\text {lat}} \gtrsim 4\) with this lattice action

Full size image

Turning to lattice \({\mathcal {N}} = 4\) SYM in four dimensions, Fig. 4 presents results for the pfaffian phase using the gauge-invariant improved action [153, 203]. Despite implementing a parallelized pfaffian computation in SUSY LATTICE, only small N and small lattice volumes are computationally accessible, with each \(\text {U(}2\text {)}\) \(4^4\) pfaffian measurement requiring approximately 50 hours on 16 cores. In the left plot, only small per-mille-level phase fluctuations are observed on all accessible volumes with fixed ’t Hooft coupling \(\lambda _{\text {lat}} = 0.5\). In particular, the expected exponential suppression of \(\left\langle e^{i\alpha } \right\rangle _{\text {pq}}\) with the lattice volume is not visible, which is encouraging if not yet fully understood. However, the right plot shows that \(\left\langle e^{i\alpha } \right\rangle _{\text {pq}}\) decreases rapidly as the ’t Hooft coupling increases, obstructing studies of \(\lambda _{\text {lat}} \gtrsim 4\) with this lattice action. For gauge group \(\text {U(}2\text {)}\), Ref. [81] presents mixed-action evidence that sacrificing gauge invariance may improve control over pfaffian phase fluctuations at much stronger couplings \(\lambda _{\text {lat}} \sim {\mathcal {O}} (10)\), raising the possibility of more directly probing holography and S-duality.

5.3 Final remarks

Non-perturbative lattice investigations of supersymmetric QFTs are important and challenging, making this an area that should attract even more attention in the near future. It is encouraging that there has been so much recent progress in lattice studies of four-dimensional \({\mathcal {N}} = 1\) SYM and \({\mathcal {N}} = 4\) SYM, along with their dimensional reductions to \(d < 4\), where the consequences of supersymmetry breaking due to the discrete lattice space-time can be kept under control. There have also been advances in other areas of lattice supersymmetry that this brief review omits, in particular lattice studies of theories without gauge invariance, such as Wess–Zumino models and sigma models [4,5,6,7,8, 111,112,113,114,115, 204,205,206,207], as well as lattice investigations of the Green–Schwarz superstring worldsheet sigma model [208,209,210,211]. While supersymmetric QCD and sign problems present challenges that may be difficult to overcome, the overall prospects of lattice supersymmetry are bright, with many compelling opportunities for future work.