Glueballs and Strings in $Sp(2N)$ Yang-Mills theories

Motivated in part by the pseudo-Nambu Goldstone Boson mechanism of electroweak symmetry breaking in Composite Higgs Models, in part by dark matter scenarios with strongly coupled origin, as well as by general theoretical considerations related to the large-N extrapolation, we perform lattice studies of the Yang-Mills theories with $Sp(2N)$ gauge groups. We measure the string tension and the mass spectrum of glueballs, extracted from appropriate 2-point correlation functions of operators organised as irreducible representations of the octahedral symmetry group. We perform the continuum extrapolation and study the magnitude of finite-size effects, showing that they are negligible in our calculation. We present new numerical results for $N=1$, $2$, $3$, $4$, combine them with data previously obtained for $N=2$, and extrapolate towards $N\rightarrow \infty$. We confirm explicitly the expectation that, as already known for $N=1,2$ also for $N=3,4$ a confining potential rising linearly with the distance binds a static quark to its antiquark. We compare our results to the existing literature on other gauge groups, with particular attention devoted to the large-$N$ limit. We find agreement with the known values of the mass of the $0^{++}$, $0^{++*}$ and $2^{++}$ glueballs obtained taking the large-$N$ limit in the $SU(N)$ groups. In addition, we determine for the first time the mass of some heavier glueball states at finite $N$ in $Sp(2N)$ and extrapolate the results towards $N \rightarrow +\infty$ taking the limit in the latter groups. Since the large-$N$ limit of $Sp(2N)$ is the same as in $SU(N)$, our results are relevant also for the study of QCD-like theories.

Recent years have seen a resurgence of interest in gauge theories based upon symplectic groups, driven by theoretical as well as phenomenological motivations, related to modelbuilding in the context of physics beyond the Standard Model (SM). In comparison with the SU (N ) and (limited to (2 + 1)-dimensional) SO(N ) cases [1][2][3][4][5], the literature on lattice studies of Yang-Mills theories with Sp(2N ) gauge groups is limited in its extent, scope and reach (see for instance Ref. [6]). In a recent publication [7] (see also Refs. [8][9][10][11]), some of us announced the intention to carry out a long-term, systematic lattice exploration of the strong-coupling dynamics of the theories based on Sp(2N ) gauge symmetry, and proposed a research programme that includes as one of its crucial steps the study of the dynamics of glueballs and strings in the pure gauge theory.
In the same paper [7], we presented our first comprehensive study of the Sp(4) pure gauge theory, and computed the spectrum of masses and decay constants of mesons consisting of two fundamental (Dirac) fermions, treated in the quenched approximation. In the context of Composite Higgs Models (CHMs) [12][13][14] (see also Refs. [15? -54]), that initial step provided an important source of quantitative information about the underlying dynamics. In particular, we started to explore and exploit the dynamical origin of low-energy Effective Field Theories (EFT) based upon the SU (4)/Sp(4) coset, which have a prominent role in the CHM context (see for instance Refs. ), as well as for related models of dark matter with strong-coupling origin [76][77][78][79]. More recently [80], some of us presented the first continuum results of the lattice study of the Sp(4) theory with dynamical Wilson fermions, hence making the treatment of the dynamics more realistic and useful in the CHM context. A first set of exploratory studies of the quenched theory with valence fermions in multiple representations has been published in Ref. [81].
In the present paper (see also Ref. [82,83]), we take major steps in a complementary direction, by focusing on the pure gauge theory without matter content, but extending the analysis to different Sp(2N ) gauge groups. Our specific objective is to obtain for the Sp(2N ) Yang-Mills theories in D = 3 + 1 dimensions a comparable level of control over the spectra of strings and glueballs as achieved for the previously studied SU (N ) and SO(N ) gauge theories [1][2][3][4][5]. On a theoretical side, this endeavour will allow us to study the approach towards the common large-N limit via an alternative sequence of groups in respect to SU (N ) and SO(N ). In turn, this will provide an alternative set of numerical tests for such conjectural behaviours as those put forwards for example in Refs. [82,84,85], as well as allowing comparison to calculations performed within the context of gauge-gravity dualities (see for instance Refs. [86][87][88][89][90][91][92][93][94][95][96]) or with alternative field theoretical methods [97][98][99][100][101]. In pragmatic terms, we will also set the stage for future studies in quenched theories realising the SU (4)/Sp(4) coset, based upon generic Sp(2N ) groups, extending the results obtained in [81] for the Sp(4) gauge theory.
In our investigation, we adopt a unified approach to the study of Sp(2N ) gauge theories, by applying the same heat bath (HB) algorithm exploited in Ref. [7] for the Sp(4) theory to the whole Sp(2N ) sequence. In addition to reconsidering Sp(2) ∼ SU (2), which allows to test our algorithm and procedures by comparing to existing results in the literature and to extending the N = 2 results discussed in Ref. [7] with new calculations, we consider the N > 2 cases. For the latter, with the exception of our study in Ref. [82] (focussing on a discussion of the two lowest-lying glueball states and on a remarkable universality property of their ratio) and Ref. [83] (presenting some preliminary numerical results, further discussed in the current work), no detailed calculation of the glueballs has been reported in the literature so far. From an operational perspective, we first compute the effective string tension and glueball masses in the large-volume limit for fixed lattice spacing and N = 1, 2, 3, 4. Then, after taking the continuum limit of the glueball spectrum at each investigated value of N , we perform a critical analysis of the large-N extrapolation and compare to other results in the literature, as appropriate.
The paper is organised as follows. In Section 2 we introduce the basic definitions and conventions adopted in the lattice calculations. In Section 3 we describe the spectral observables of interest. In Section 4 we present our numerical results, including also the extrapolations to continuum and large-N limits. Section 5 summarises our conclusions and suggestions for future further enquiries. We have relegated some important technical details to the Appendices.

Numerical simulations of the Lattice model
In four Euclidean dimensions, the Sp(2N ) gauge theory is defined by the following action where g 0 is the gauge coupling, the trace is over colour indices, the field-strength tensor F µν ≡ A F A µν τ A is defined by and the gauge fields are A µ = A A A µ τ A , with the indices taking the values A, B, C = 1 , · · · , N (2N + 1), for Sp(2N ). The 2N × 2N matrices τ A are the generators of the algebra associated with the Sp(2N ) group, written in the fundamental representation, and normalised according to Tr τ A τ B = 1 2 δ AB . The structure constants of the algebra are defined as the commutation relations We regularise the theory on a lattice, in which the continuum coordinates are discretised with lattice spacing a. The four dimensional Euclidean hypercubic lattice consists of sites that are denoted by their position x in the lattice. The sites are connected by links that are characterised by the position x and direction µ, where µ, ν = 0, .., 3 label the four space-time coordinates. The elementary variables of the lattice regularised Sp(2N ) gauge theory are the link variables, defined as withμ the unit vector in direction µ. The 2N × 2N matrices U µ (x) transform according to the fundamental representation of the Sp(2N ) group. Gauge transformations take the form U µ (x) → g(x)U µ (x)g † (x +μ), with g(x) a group element. The simplest gauge invariant operator is the trace of the product of link variables around an elementary square of the lattice, The matrices P µν (x) are called the elementary plaquette variables or just plaquettes for short.
The Sp(2N ) lattice gauge theory (LGT) we adopt in this paper is defined by the Wilson action, Tr P µν (x) . (2.6) In this expression, Tr P µν (x) is the real part of the trace of P µν (x). The inverse coupling β is related to g 0 by the request that, when the lattice spacing a → 0, Eq. (2.6) tends to the continuum Yang-Mills action in Eq. (2.1), at leading order in a. From this requirement, one finds β = 4N g 2 0 . (2.7) Monte Carlo numerical evaluations of the integrals appearing in the definitions allow us to explore the long-distance regime of the Sp(2N ) (pure) Yang-Mills theories, capturing nonperturbative phenomena that are not accessible to perturbation theory. For any quantity O(A µ ) that depends on the gauge fields, the physical observables are estimated as ensemble averages, which are schematically given by , (2.8) where the denominator is Z(β) ≡ DU µ e −S W . (2.9) These expressions can be computed numerically by sampling the space of configurations of U µ (x), according to the probability distribution e −S W . This can be achieved by defining a Markovian process that evolves a particular configuration according to an update algorithm.
The algorithm must respect detailed balance and reproduce the correct equilibrium distribution. Then, if i labels the M configurations produced sequentially, the ensemble average can be obtained as the simple average where O i is the value that the observable O(U µ ) takes on configuration i. The algorithm adopted in this work to produce successive configurations is a combination of local heat bath (HB) and overrelaxation (OR) updates, adapted to Sp(2N ) from the SU (2N ) implementation provided in Ref. [102] (see Appendix A for further details). Configurations are updated sequentially, one link at a time, with one HB update followed by four OR updates. An update of all the links on the lattice is called a lattice sweep. Successive configurations produced in this manner are correlated; to reduce the effects of autocorrelation, the ensemble averages used for physical calculations are restricted by sampling the history in steps that are separated by 10 lattice sweeps. Our implementation of the algorithms above is based on the HiRep code [103], originally designed for the treatment of SU (N ) theories with matter fields in general representations. 1 The lattice size being finite, we impose periodic boundary conditions in all directions. In the continuum, it is known that resulting configurations of gauge fields are characterised by an integer topological number [104], defined as Q ≡ 1 32π 2 µνρσ d 4 xTr F µν F ρσ . (2.11) The associated susceptibility can be related to the large mass of the η particle [105]. The configuration space is thus divided into sectors, each characterised by an integer value of the topological number Q, and separated from each other by potential barriers. Because of the lattice discretisation, the topological charge Q takes nearly integer values [106][107][108]. There are many microscopic lattice definitions of the topological charge that reproduce the same, correct long-distance results in the a → 0 limit. In this work we adopt the definition Q ≡ x q(x) , (2.12) with q(x) ≡ 1 32π 2 µνρσ Tr {U µν (x)U ρσ (x)} , (2.13) and where x runs upon all lattice sites. Since these definitions make use of the shortdistance degrees of freedom, calculations are affected by short-range fluctuations. These effects can be reduced by the use of smoothing operations such as the Gradient (or Wilson) Flow [109], which we will introduce below. As in the continuum, also on the lattice the different topological sectors are separated by potential barriers. If these barriers are not too steep, in simulations a sufficient number of tunnelling events between sectors will occur, and the resulting measured topological charge will be Gaussian distributed around zero. However, superselection of topological sectors can be shown to emerge close to the continuum limit [106,107]. As a consequence, Monte Carlo update algorithms tend to become trapped inside one of the topological sectors. Hence, close to the continuum limit, the topological charge has a long autocorrelation time. This phenomenon is referred to in the literature as topological freezing. Due to large-N suppression of small-size instantons, which are crucial for changing the topological charge in numerical simulations [110], topological freezing becomes more severe as N increases. We shall discuss implications of this algorithmical trapping more extensively later in the paper, focusing on the effects of topological freezing on the observables that are of interest to us.
To remove ultraviolet fluctuations that would otherwise dominate the signal in the extraction of the topological charge, we employ the Gradient Flow [109,111] of the Wilson action ( i.e. the Wilson flow). The Gradient Flow provides a first-principles approach to the smoothening of configurations with efficiency comparable to that of the more empirical and time-honoured cooling methods (see for instance Ref. [112]). Moreover, the evolution of observables under the Gradient Flow can be determined with numerical procedures that can be easily implemented. For this reason, this method has gained a prominent role in lattice studies in recent years.
With t the coordinate in an additional fifth dimension (referred to as flow time) and x a point in four-dimensional space, the Gradient flow B µ (t, x) is defined by the following differential equations and boundary conditions: (2.14) Here A µ (x) is the continuum gauge field, while the covariant derivative is which yields the field-strength tensor: On the lattice, the Gradient Flow for the action in Eq. (2.6) is defined by Here, S flow is the Wilson plaquette action for V µ . The Gradient Flow describes a diffusion process with time t. At the leading order in the coupling g 0 , the flow to time t acts on the gauge fields as a Gaussian spherical smoothing operation, with root-mean-square radius √ 8t, the flow time t having the dimension of a length squared. Furthermore, to all orders in perturbations in g 0 , any gauge invariant composite operator constructed from B µ (t, x) is renormalised at t > 0, and thus directly encodes physically observable properties. Using a value of the flow time τ such that a √ 8τ R, where R is a typical hadronic scale, provides four-dimensional smoothed configurations V (τ, x) that are not affected by ultraviolet fluctuations and still encode the correct infrared behaviour of the theory.

The spectrum
Non-Abelian Yang-Mills theories confine, and their spectra consist of massive colour-neutral states called glueballs. If a non-Abelian gauge theory is formulated on a space with one or more compact directions, wrapping torelon states arise. The validity of the confinement picture for the specific case of Sp(4) has been confirmed explicitly in the numerical calculations reported in Refs. [6,7]. The main objectives of this work are to show through lattice calculations that, as one would expect, confinement arises also in Sp(6) and Sp (8), to measure the resulting glueball mass spectrum, and to determine the large-N limit of the latter.
Before discussing our numerical results, in this section we review the methodology we shall adopt. The methodological material presented in this section is based upon notions that have been tested and are well established in the literature. Details beyond our exposition can be found, e.g., in Refs. [113][114][115][116][117][118], from which we draw heavily in what follows.

The variational method
Let H be a Hamiltonian of the 3-dimensional system of volume L 3 defined on a lattice with L t time slices. 2 Let |n and E n be the eigenstates and eigenvalues of H, i.e., H|n = E n |n . (3.1) The transfer matrix, is the operator that evolves one time slice of the system into the next. Note that in this section, for simplicity, we reabsorbe β in the definition of H. In terms of T, the partition function in Eq. (2.9) can be expressed as Masses of particle states can be obtained from the large time decay rate of (normalised) two-point correlators of interpolating operators, where |Ω is the vacuum state, normalised so that |Ω = T|Ω , and O(t) is an interpolating operator that produces the single-particle state |Ψ by acting on the vacuum, with Ω|Ψ = 0. Inserting a complete set of eigenstates of H in Eq. (3.4), we obtain where the coefficients c n , given by are called overlaps. If E 0 < E 1 < · · · , then, This equation implies that, in principle, E 0 can be obtained by fitting an exponential to the large t values of C(t) as measured from the lattice. When O(t) creates a zero-momentum state, the energies E i are identified with particle masses m i . In our calculation we will restrict to this case.
Following from Eq. (3.9), we define the effective mass m ef f (t) as (3.10) If a one-particle eigenstate of the Hamiltonian were propagating, m eff (t) would be constant with respect to t with a value equal to the mass of that state. In the presence of other states contributing to the correlation function, we expect this effective mass to be an upper bound for the true asymptotic mass at any finite t. In numerical studies, a t min can be identified such that, for t ≥ t min , only the ground state (or, more precisely, the smallest mass eigenstate with non-zero overlap) contributes to C(t) within the statistical precision, and hence am eff (t) becomes constant. The plateau value of m eff (t) provides an estimate of the ground state mass m 0 , which can be extracted by fitting a single exponential to the data for C(t) for t ≥ t min .
While this programme is at the basis of standard techniques for extracting masses from correlators, its direct implementation is not straightforward, and requires a careful treatment of numerical data. The first difficulty one encounters stems from the statistical noise affecting the measurements. In fact, while the statistical fluctuations of C(t) are roughly independent of t, the magnitude of correlation functions decays exponentially. This gives an exponentially-suppressed signal-to-noise ratio which is prohibitively hard to improve upon with an increase in the measurement sample size alone. In addition, the value t min of the onset of the single-exponential asymptotic regime is not known a priori; it is a model-dependent feature, sensitive to the mass spectrum in the given channel and to the choice of the operator O, as well as to the precision of the numerical calculation. The time t min is extracted from the simulations. Moreover, simple arguments based on asymptotic freedom show that for a given operator and in a given channel, t min grows exponentially as the continuum limit is approached.
The discussion above highlights the necessity to go to large times to isolate the ground state, but then the signal-to-noise ratio degrades and this makes it difficult to estimate m 0 in a reliable way. If one could find operators with correlators that provide single exponential behaviours, one could perform fits at small times, when the signal is still well visible above the noise. Although this ideal situation can not be reproduced in numerical investigations, since the knowledge of operators giving rise to single exponential correlators would only arise from an at least partial solution of the theory, one can try to engineer the calculation in such a way that in each relevant channel t min is as small as possible. For this purpose, at each value of a we construct interpolating operators that maximize the overlaps with the spectral states of interest. The main idea is to approximate the (unknown) exact eigenfunctions of H with an appropriate linear combination of a set of states {|Ψ i }, chosen on the basis of symmetry considerations, as trial wave functions. Then, in the given channel, the mass of the lowest lying state above the vacuum can be bound as where Ψ denotes any linear combination of the variational basis Ψ i , subject to the constraint Ψ|Ω = 0, and τ is a time chosen for minimisation, which is performed across the linear combinations of our basis operators. This bound is saturated by the lowest-lying eigenstate of the Hamiltonian in the chosen channel, which can be obtained using a complete set of variational states {|Ψ i }. Since variational bases used in calculations are necessarily finite, the bound is in general not saturated when the variational method is used in practice. Nevertheless, with a suitably large variational basis, the extracted variational mass m var will eventually be compatible within the statistical errors with m 0 . This framework, referred henceforth as the variational technique, can be implemented algorithmically in order to extract both the glueball and the torelon spectrum in various channels [117]. The success of this approach and the quality of the results obtained with this technique crucially depend on the nature of the operators that we include in the variational basis. For this reason, particular attention needs to be paid to its construction. We will review in the following two subsections the approach we followed to construct trial states to be used in the variational calculation, by discussing how gauge invariant states are created on the lattice in Sect. 3.2, and how one obtains the irreducible representations of the symmetry group of the lattice in which these states must transform in Sect. 3.3. In Sect. 3.4 we will show how to perform the extremisation provided in Eq. (3.11) in an effective way, in order to obtain robust estimates of m 0 . The effective description of torelon states as closed fluxtubes will be summarised in Sect. 3.5. The estimates of m 0 will be affected by systematic errors of different origins, which will be discussed in Sect. 3.6.

States on the lattice
In this section we explain how to create gauge invariant states out of the vacuum and their interpretation in terms of glueball and torelon states. Consider traced path ordered products of links, defined in Eq. (2.4), around closed spacelike loops C, where C can be defined as a set of successive displacements, where each f j is one of the elementary vectors of the lattice { e i }. The sequence f 1 , f 2 , · · · , f L is defined up to cyclic permutations. The closeness of the path C implies that i f i = 0. A generic gauge invariant operator O such that Ω|O|Ω = 0 can be obtained as a sum of products of operators O C , each defined as for specified choices of C.
Single trace operators create states called glueballs when C is contractible and torelons states when C wraps around a spatial direction of the space-time hypertorus, and is thus non-contractible. These two classes of states transform in different representations of the center of the group and hence do not mix. We will start our analysis from the contractible loops. Most of our arguments are applicable also to non-contractible loops, which will be analysed more specifically in Sect. 3.5.
Multitrace operators are monomials involving products of at least two of the operators in Eq. (3.12). Operators in this class can be used to generate multi-particle states. Some of these states have the same quantum numbers as single particle states we are interested in, and can thus mix with them. This mixing can result in a systematic error in the extraction of masses of single-particle states. In our calculation, we will neglect mixing of genuine glueball states with multi-particle states. The justification for neglecting multitrace operators resides in the fact that matrix elements involving them go to zero in the large-N limit.

Symmetries
In the lattice theory, the Poincaré symmetry of the continuum is explicitly broken to the discrete subgroup of symmetries of the hypercubic lattice and discrete translations by an integer number of lattice spacings. In particular, on an infinite lattice, for a time slice, this is the semi-direct product of discretised translations T d and of the point groups of invariance of the elementary (cubic) cell of the lattice: the Octahedral group O h (see e.g. Ref. [119]). The study of the representations of this symmetry group of the lattice is simplified by the fact that T d is an invariant abelian subgroup. The one-dimensional representations of O h (related to momentum) can thus be studied separately from those of T d .
In a finite box of size L, and with periodic boundary conditions, the momentum is quantised in every direction as p n = 2πn/L. On a lattice its value must also lie in the Brillouin zone containing p = 0. Operators at fixed momentum can be obtained as Fourier sums of their counterpart in coordinate space, (3.15) Zero momentum combinations (to which we restrict ourselves in this study) can be simply obtained as sums over fixed time slices of operators of the type given in Eq. (3.14), We now briefly describe the irreducible representations of O h and their relation with the representations of the Poincaré group. The Octahedral group is the symmetry group of a cube. This group has 24 elements divided in 5 conjugacy classes. Accordingly, it has 5 inequivalent irreducible representations, labelled by R = A 1 , A 2 , E, T 1 , T 2 , of dimensions 1, 1, 2, 3, 3, respectively. The spatial parity P has two eigenstates, which we label by an additional ± sign, depending on whether they remain invariant (+) or are reflected (-) under a parity transformation. We will label the states of the lattice theory with R = R P and their mass with m R P . Asterisks will denote excitations of the ground state: A + * 1 will denote the first excited state of A + 1 , A + * * 1 the second, etc. The states generated from the vacuum by gauge invariant operators U (C) will transform in the same representation as the paths on which they were defined according to Eq. (3.12). In general, single trace operators belong to reducible representations of the octahedral group. Under the action of an element r of the group, the operators O C transform in representation U (r) in the following way where the law of transformation of C can be inferred from its definition in Eq. (3.13), The decomposition of U (r) in terms of its irreducible components can be obtained from the orthonormality property of characters, supplemented by a choice of orthonormal bases for each of the irreducible representations R P of O h . For this, the projector method borrowed from Ref. [119] has been used.
In the continuum limit, we expect the Poincaré symmetry to be recovered. The relationship between the representations of the octahedral group defined above and those of the Poincaré group enables us to decompose the former in their continuous spin components. The representations of the Poincaré group are labelled by the mass m and the quantum numbers J P C , where J is associated to irreducible representations of the rotation group, P to spatial parity and C to charge conjugation. Owing to the pseudo-reality of the representations of Sp(2N ), C is always positive. Hence, we will drop this quantum number from now on.
If we restrict the elements of the rotation group in a representation J to the discrete rotations that lie in O h , we obtain the subduced representation J ↓ O. We report in Tab. 1 the subduced representations for the lowest values of J, adapted from Ref. [117]. In O h , these representations are reducible in terms of A 1 , A 2 , E, T 1 and T 2 . Thus, degenerate states with the same spin but different polarisation of the continuum spectrum might have a different mass on the lattice. In the continuum limit, nevertheless, the restoration of continuum rotational invariance implies that these states become degenerate. For instance, the E and T 2 representations of the octahedral group contain respectively two and three of the five polarisations of spin-2 particles. Hence, corresponding states extracted in the E ± and T ± 2 channels must become degenerate in the continuum limit. The degree of degeneracy of these states at finite lattice spacing will thus provide an important measure of the effect of lattice artefacts.

Extraction of masses
Let us now consider a specific irreducible representation R P and build a generic linear combination Φ of basis elements O R P at time t, which we denote as Table 1. Subduced representations R of the continuum rotation group and their components labelled with the spin J, up to J = 4. 20) where, in general,

The 2-point correlation function is
with c a i = a|O R P i (0)|Ω . As a result, Eq. (3.10) can be rewritten as .
The matrix C ij (t) is positive-definite (see Eq. (3.21)), and its eigenvalues are given by λ a (t) = e −mat . Hence, extracting the spectrum is equivalent to the diagonalisation of C ij (t).
Unfortunately, due to statistical fluctuations, eigenvectors and masses of the measured C ij (t) do depend on t. In order to resolve this dependency, we seek a solution to the generalised eigenvalue problem by diagonalising C −1 (0)C(τ ) for some τ > 0. The eigenvectors of C −1 (0)C(τ ) provide us with a practical choiceΦ i of the optimal operators. The corresponding masses m i can be obtained from fits of the correlators of theΦ i (which we refer to asC i ) at t > t 0 , using the ansatzC over ranges of t for which reaches a plateau value. We still denote as am eff the effective mass, although we adopt from now on a definition that differs from the one in Eq. (3.10). The reason for the discrepancy, which is visible only away from the large volume limit, is a consequence of adopting periodic boundary conditions in time, which allows for both forward and backward propagating states. The mass of the ground state, m 0 , is obtained from a fit of the largest eigenvalue λ 0 . The masses of higher energy states can be obtained in the same manner from the diagonal correlators of eigenvectors associated to the other eigenvalues computed in the generalised eigenvalue problem. As discussed earlier, a crucial ingredient for an efficient variational calculation is the preparation of trial states that have an extension comparable to the target glueball state. A priori, we have no information about the physical size characterising glueball states. To determine an efficient linear combination, we shall insert in our variational set operators obtained from prototypical paths of different sizes and shapes, and also operators obtained from the original basis at each of S iterations of smearing and blocking operations, with the combination obtained in Ref. [2] (to which we refer the Reader also for the definition of the operations of blocking and smearing and for specific details on the particular paths used to define the basis operators). In this way, we obtain a variational basis that finely scans the propagating states from length scales corresponding to the lattice spacing all the way up to the lattice extent L.
From the technical perspective, the only procedural change to the methodology employed for Sp(4) in Ref. [7] (to which we refer for further details) lies in the projection and cooling routines, that had to be adapted to the case of Sp(2N ). With M elementary paths and S smearing steps used for constructing the basis, our variational basis is formed by S × M operators in total, andC ij (t) will accordingly have (S × M )(S × M + 1)/2 elements. In our calculations, we perform the maximum number of blocking steps N b allowed by the finite size, provided by the maximisation of the l.h.s. in the inequality 2 N b ≤ L/a. At each blocking step we perform 2 smearing steps and 15 cooling steps to reproject on the group. In general, our variational basis contains approximately 200 elements.

Effective String Theory
Torelon states are generated from the vacuum by path ordered products of link variables along non-contractible paths, i.e. paths that wrap the periodic lattice along a given direction. These states have the same quantum numbers as physical states in which a wrapping closed loop of glue with fixed length is propagating in the system. We refer to this configuration as a fluxtube. When the fluxtube is long enough, it can be described by an effective string theory. This classical effective theory is written in terms of a single physical parameter, the string tension σ, that governs the energy of the fluctuations. In order to extract it from the data with the highest precision, we will make use of effective string theory, as briefly summarised in this section.
Effective string theory is based on approximating the fluxtube as a one dimensional fluctuating object -a string -with constant energy per unit length. Classically, the mass m and the length L of the fluxtube are proportional, This classical string description becomes exact in the long string limit L 2 σ → ∞.
At finite length, quantum corrections become relevant. The energy of the fluxtube is obtained as a power expansion in 1/(σL 2 ) around the long string limit. In general where the dimensionless coefficients d k , which are in principle calculable, can be determined by matching the power series to the results of numerical measurements. Universality theorems allow to fix some of these coefficients on the basis of symmetry arguments. The formation of the fluxtube can be described as a process of spontaneous breaking of some of the generators of Poincaré symmetry. We omit details, for which we refer the Reader to the literature [120].
The ground state mass m of a torelon wrapping along one direction of extent L is given, in a spacetime of dimension D, by where we included the leading order correction in an expansion in 1/σL 2 , and which describes m(L) up to the next-to-leading-order correction. At this order, one can show that these predictions are universal, i.e. the coefficients are fixed by Poincaré invariance and certain geometric dualities. The only physical parameter to consider is thus σ.
In general, for the ground state, no deviations with respect to the Nambu-Goto formula are allowed up the term 1/(σ 2 L 5 ). These results will allow us to compute σ from the mass of torelons, keeping under control the effects of working at finite L.

Sources of systematic errors
There are several sources of systematic errors that affect the computation of glueball and torelon masses. In this section, we discuss the most relevant ones for our study.
As explained in Sect. 3.4, the variational technique depends on our choice of basis of operators. A potential source of error is the choice to only include single trace operators in our variational set. By doing so, we are neglecting scattering states and multitorelon states that share the quantum numbers of single glueball states of interest. In the case of scattering, we deal with states with two or more glueballs. Neglecting the interactions (an approximation that holds at large N ), these states have masses that are about twice as large as the smallest glueball mass. Thus, below this threshold, we can safely neglect the effect of scattering states. Even above that threshold, scattering states decouple at large N .
Multitorelon states have a mass that is in general an increasing function of L. Therefore, at large enough values of L, they decay quickly in correlators and can hence be neglected as well. The effect of these states can in principle be controlled by including the corresponding operators in the variational basis and evaluating their overlaps, as done in Ref. [117].
As a consequence of the fact that we are simulating a finite lattice, all our physical estimates will be affected by finite size effects. These effects have been reported in Ref. [121], where it is shown that they obey the relationship with m(L) and m the masses in volume V = L 3 and at infinite volume, respectively. b is a coefficient that, a priori, depends on the symmetry channel. Under the assumption that these corrections are independent of the lattice spacing a, we will be able to compute them at one value of a and use the same prediction for all others. More so, we will be able to neglect them altogether once we find that, at a certain combination of a and L, these effects are much smaller than the statistical error. Discretisation errors come from the dependence of the masses on the lattice spacing. A trivial dependence can be inferred from dimensional analysis. The lattice combination ma is dimensionless. Since all masses obtained on the lattice depend on the lattice spacing a in this way, we consider ratios of dimensionful objects where the trivial dependence simplifies in the ratio. As a reference scale for the ratio, we use the square root of the string tension √ σ. This choice is motivated by the fact that, thanks to the results discussed in Sect. 3.5, we can measure the string tension more accurately than any other quantity of interest. Hence, the use of the string tension reduces significantly the systematic error due to the scale setting process, which is a necessary step to provide quantities in physical units.
Beyond the overall dependence of the mass on the lattice spacing a, we know, by computing the naïve continuum limit of the theory described by the lattice action in Eq. (2.6), that the leading corrections to mass ratios start at order a 2 . Therefore, close to the continuum limit, for a glueball state R P , we approximate To conclude this overview of systematic effects, we return to mentioning topological freezing. Near the continuum limit, the Monte Carlo updates tend to get trapped in a sector at fixed topology. This topological trapping becomes more pronounced at larger N [110]. Restricting the gauge theory to a sector at fixed topology generates power-law corrections in the inverse volume that delay the onset of the large volume regime [122,123]. Both large-N reduction arguments [124] and the large-N scaling prescription of the θ angle [125] suggest that finite volume corrections due to topological freezing are suppressed at large N . The decreased severity of topological freezing as N increases has been verified explicitly in Ref. [126]. In Sect. 4, we will show that topological freezing affects only a small subset of our calculations. When discussing the relevant ensembles, we shall describe how topological freezing has been accounted for in those specific cases.

Numerical Results
In this section, we present and discuss our main numerical results. In Sect. 4.1 we perform calibration and validation studies of the underlying algorithm. We also select the values of the coupling β at which to compute the masses of torelons and glueballs for the Sp(2N ) Yang-Mills theories with N = 3, 4. In Sect. 4.2, we compute the ground state mass of torelons of various lengths at fixed lattice spacing. We compare the results to the predictions discussed in Sect. 3.5. We also evaluate finite size effects, alongside exposing our strategy for extracting the string tension using one lattice size, in the asymptotic regime. In Sect. 4.3 we report the results of the continuum limit extrapolations of the glueball spectrum for N = 1, 2, 3, 4, while we cover in detail in Appendix C all pertinent technicalities. The continuum limit values of the masses are then used to extrapolate towards the large-N limit, in Sect. 4.4.

Preliminary tests and calibration studies
We compute the expectation value of the plaquette P , which is defined as We consider several values of β, and focus attention on N = 3 and N = 4. Independent ensembles are generated for each chosen value of β, with either unit (cold), or random (hot) starting configurations in the Monte Carlo update algorithm. We calculate P each 5 sweeps, record 10 4 individual measurements of this quantity out of the 5 · 10 4 sweeps performed for each β value. By comparing the history of P starting separately with unit and random configurations, we are able to identify and discard the initial transient due to thermalisation. We have verified explicitly that the integrated autocorrelation times are less than 1.5 for all values of β. We finally bin and bootstrap the measurements of P . For Sp(2N ) Yang-Mills theories with N = 3, 4, the results are shown in Fig. 1 for lattices with size (L/a) 4 = 16 4 . Our algorithm is based on a heat bath update of Sp(2) subgroups that when combined provide a covering of the whole Sp(2N ) group (see Appendix A for a detailed explanation). For validation purposes, we obtained alternative, independent estimates for the average plaquette using the simpler (and slower) Metropolis-Hastings update algorithm. For both N = 3 and N = 4, and at every relevant value of β, the estimates obtained with the two different algorithms are compatible with each other, within one standard deviation. For N = 3, independent numerical results are also available through Ref. [6], and our results are compatible with theirs within one standard deviation, when comparisons are possible. Finally, the limits of weak (β → ∞) and strong (β → 0) coupling can be controlled analytically [6]. It is expected that and at weak and strong coupling, respectively. In Fig. 1, we compare the leading terms of these analytical predictions to our numerical data in the relevant regime. The combination of all these tests supports the robustness of the algorithm we are employing.
The behaviour of P as a function of β is also used to detect the potential presence of a bulk phase transition that separates the weak and the strong coupling regimes of the theory. While the latter is dominated by strong lattice artefacts, the former is relevant to continuum physics. The pseudo inflection point visible in Fig. 1 (for both N = 3 and N = 4) is a potential signature of such a phase transition. We study the nature of this change of regime in Appendix B, where we conclude that our numerical data are compatible with a crossover, confirming the findings of Ref. [6] for N = 3 and extending this conclusion to N = 4.
In principle, the absence of a genuine phase transition may allow extraction of physical observables by performing an extrapolation to the continuum limit that makes use of generic values of β. Nevertheless, by restricting our choices of β to the weak coupling regime we maintain better control over the approach to the continuum. Our choice of the values of β for the simulations results from a pragmatic compromise aimed at reducing discretisation errors while deploying the finite amount of available computational resources.

Torelons and strings
In this subsection, we discuss the methodology we use to extract the string tension from ground state mass of torelons of length L for N = 3, 4, while also testing the predictions of Sect. 3.5. We first perform an analysis of the L dependence of the mass at a fixed value of the lattice coupling, in order to identify the regime in which the string effective description Sp(6), β = 16.5 Sp (8) is applicable. We then extract the string tension from torelon masses measured at one asymptotic value of L for each choice of β. This procedure allows us to obtain accurately the string tension as a function of the finite lattice size, using the known functional form of the torelon mass. We retained 10 4 thermalised configurations for post-production analysis. The variational basis we adopt includes the elementary Wilson line winding around a compact spatial dimension, and averaged over all three spatial directions, alongside its blocked and smeared improved versions, up to blocking level N b such that in the inequaliy 2 N b ≤ L/a the l.h.s. is maximised. Following Ref. [7], to which we refer for further details, we performed either one (for the coarsest lattices) or two (for the finest lattices) smearing steps in-between one blocking step. For the study of the finite size dependence of the torelon mass, we generated configurations with β = 16.5 for Sp (6), and β = 26.7 for Sp(8), on the lattice sizes listed in Tab. 2. These values of β are chosen to be small enough that large physical volumes are reached with moderate computing cost, while still remaining within the weak coupling regime. The values of the masses thus obtained, denoted by am s , are reported in Tab. 2 and plotted in Fig. 2. In order to extract am s , we performed a maximum likelihood analysis based upon Eq. (3.24). The value of the resulting χ 2 /N d.o.f. is usually below or around one; exceptions to this are mostly restricted to the largest lattice studies in Sp (8), where am s becomes of order one and as a consequence the signal decays quickly.
We now test the predictions of Sect. 3.5. From Fig. 2, we see that, at the largest values of L/a, m s a is an approximately linear, increasing function of the length, both in N = 3 and N = 4. This behaviour supports the intuitive description of a torelon state as a closed fluxtube with constant energy per unit length. In order to extract the string tension σ, as a first approximation we use Eq. (3.26) applied to the largest value of L/a, treating the fluxtube as a classical string. We call σ cl the resulting string tension. For Sp(6), we find σ cl a 2 = 0.0212(4) at L/a = 28 and for Sp(8), we obtain σ cl a 2 = 0.0520 (14) at L/a = 20. The large-L expansion is expected to be well-behaved when σL 2 1. At a given value  Table 3. Measurements of the string tensions σ, based upon applying Eqs. for N = 4. From this analysis, we observe that at the chosen values of β the string picture provides a good description of the torelon mass down to lattice size L = 16a for Sp(6) and L = 12a for Sp (8). These values correspond to L √ σ 2.4 for Sp(6) and L √ σ 2.8 for Sp (8), somewhat more generous than the generic, conservative estimates we anticipated. We hence impose the bound L √ σ ≥ 3, in order to control the extraction of the string tension through an asymptotic large-L expansion including at the least the LO correction. This observation will be used in the following to extract the string tension at other β values for both N = 3 and N = 4, when we will apply the NG expression to torelon masses obtained at a single size L and test a posteriori that L/a fulfils the condition L √ σ ≥ 3.

The Glueball spectrum
In this Section, we report on the the spectrum of glueballs in Sp(2N ) gauge theories for N = 1, 2, 3, 4, for each fixed value of N , focusing on the results obtained in the continuum. Calculations of the masses in units of the lattice spacing, finite size effects studies and more technical details on the continuum extrapolation can be found in Appendix C.
We report the glueball masses in Tab. 4, in units of √ σ (top half of the table) as well as in units of the mass of the E + state (bottom half). The spectra at the various values of N are also presented in Fig. 3, where, together with the lattice quantum numbers, we display the continuum quantum numbers of glueball states. The latter have been obtained from the decomposition presented in Tab. 1 under the assumption that lighter states correspond to lower continuum spin. For N = 1, since Sp (2) SU (2), results for the spectrum are already present in the literature, see for example Ref. [2]. In this case, the results obtained in our study are useful for comparison and as a test of our procedure. For SU (2), Ref. [2] finds the values m A + 1 / √ σ = 3.78 (7) and m E + / √ σ = 5.45 (11). These values are compatible (84) 3.577 (49) 3.430 (75) 3.308 (98) 3.241(88) A + * 1 5.22 (33) 6.049 (40) (33) 0.711 (21) (14) for the E ++ channel (data taken from Ref. [2]). As expected, at least for these three channels, which are the only ones for which we can compare, the masses of Sp(N → ∞) and SU (N → ∞) theories are compatible.
within one standard deviation with the values obtained in this work (see Tab. 4). For N = 2, some of us already obtained first results for the spectrum in Ref. [7]. We combine our new measurements for Sp(4) with our earlier results, and in Tab. 4 we report the weighted averages of the two. The available data sets for Sp(4) are discussed more in detail in Appendix D.
A look at Fig. 3 shows that, while specific details depend on N , there are common patterns across the investigated values of N . As expected, the A + 1 channel is consistently the lightest, followed by the (T + 2 , E + ) (degenerate) pair. At a slightly larger mass we find the A − 1 channel and the T − 2 and E − (degenerate) pair. As explained in Sect. 3.3, the degeneracy of these pairs provides evidence that the rotation invariance of the continuum theory is recovered as a → 0. The remaining channels, A ± 2 and T ± 1 , are also almost degenerate in pairs and their masses are larger than those of all other states. Since the smallest masses in the A ± 2 and T ± 1 channels are comparable with twice the ground state mass of the lowestlying A + 1 state, numerical results for these masses may be affected by systematic errors due to mixing with scattering states, as discussed in Sect. 3.6. An indication of this is the fact that the error bars for the masses of those heavier states are visibly larger. Large error bars are also the result of the higher level of noise affecting the extraction of masses of heavier states.
We were able to extract masses of excited states for the A ± 1 channels at all values of N . These masses are reported in Tab. 4 and displayed above the corresponding ground states in Fig. 3. The error bars of the A + 1 states are comparable to those of the ground state in the A − 1 channel, while for the A − 1 states they are similar to those found in heavier channels.
Finally, we note that, where determined in both calculations, corresponding states obtained from a recent SU (3) study [127] are in broad agreement with the spectrum resulting from our investigation.

The spectrum towards the large-N limit
As shown for instance in Ref. [128], while corresponding quantities in SU (N ) and Sp(2N ) Yang-Mills theories converge to a common large-N limit, the 1/N expansions around this limit are different: in the case of SU (N ), only even powers of 1/N are present, while for Sp(2N ) the power expansion is genuinely in 1/N . Following the strategy that has been implemented in the large-N extrapolation of the SU (N ) glueball masses, we shall investigate whether the lowest order correction to the large-N limit is sufficient to describe the large-N glueball spectrum in Sp(2N ) for all the simulated values of N . Therefore, we fit the finite-N spectrum with the ansatz where c R P is a constant (expected to be of order 1 in a well-behaved expansion) that depends on the glueball channel. If the ansatz provides a sufficiently accurate description of the data, is a reliable infinite-N extrapolation of the ratio of the mass in the channel R P normalised to the square root of the string tension.  (81) 1.49 Table 6. Resulting values of the universal constant η for the Casimir scaling described in Eq. (4.7) for Sp(2N ) and SU (N ) groups. The combined fit to both groups is also reported.
for a range N > 1. Although generally the χ 2 /N d.o.f. are smaller, in this latter case the parameters c R P and m R P / √ σ(∞) are estimated from three data points only and thus only one degree of freedom remains to assess the goodness of the fits. For this reason, we opt to present the extrapolations including the (generally still acceptable) N = 1 data points, postponing to future studies that investigate larger N the question of whether N = 1 is captured by a simple leading correction with the current precision of the data. For the time being, in the absence of any evidence that would suggest otherwise, we assume that indeed this is the case. For some of the lightest states for which the continuum mass in the large-N limit is available in the literature, we can verify that the large-N extrapolation of the Sp(2N ) and of the SU (N ) values are compatible. In Fig. 3 the spectrum in the large-N limit is represented together with the finite-N one, to allow for such a comparison. Recalling that charge conjugation is always positive in Sp(2N ), for the sake of comparing corresponding states, we temporarily reintroduce the corresponding index in the notation for glueball states for the rest of the current subsection. With the second + superscript identifying positive charge conjugation, we borrow the values of the A ++ 1 , A ++ * 1 and E ++ channel masses for SU (∞) from Ref. [2]. Fig. 3 shows that the large-N extrapolations of the A ++ 1 , the A ++ * 1 and the E ++ in Sp(2N ) and SU (N ) are compatible with each other, as predicted by general large-N arguments.
Armed with the results of the mass calculation of the A ++ 1 , we can provide further support to the conjecture put forward in Ref. [85] that the quantity η in the relationship is a universal constant depending only on the dimension of space-time. In this equation, C 2 (F ) and C 2 (A) are the quadratic Casimir operators in the fundamental and adjoint representations respectively, whose ratio in Sp(2N ) is given by After performing the standard identification of the A ++ 1 with the lowest-lying scalar glueball, we tested this conjecture by performing a fit of Eq. (4.7) to the data, using η as a fitting parameter. The result can be found in Tab. 6 and is represented in the top panel of Fig. 5. The behaviour of η as a function of N is compatible with a constant for both the Sp(2N ) and the SU (N ) sequence. Moreover, the values of η extracted in each sequence are compatible with each other within one standard deviation, as reported in Tab. 6. As an additional test of Eq. (4.7), the behaviour of m 2 0 ++ /( √ ση) is represented in Fig. 5 for both Sp(2N ) and SU (N ) groups, along with the ratio of the quadratic Casimir operators. The weighted mean of the values of m 2 0 ++ /( √ ση) obtained in each series is also reported in Tab. 6 and represented in Fig. 5. This analysis provides further indications of the validity of the conjectured Casimir scaling, at least within current accuracy and precision. Another remarkable universal property is the independence on the gauge group of the ratio between the mass of the tensor glueball and the mass of the scalar glueball. This has been the subject of the investigation reported in Ref. [82], that makes use of the measurements reported here. We do not repeat the details of that analysis, but refer the interested Reader to Ref. [82].    2N ). Note that the naming convention for the symplectic group has been altered, using the variable N c = 2N , to better accomodate the data into the plots; fits are also shown for the Sp(N c ) family, the SU (N c ) family and the combination of Sp(N c ) and SU (N c ) results. Bottom: measured ratios m 2 0 ++ /σ further divided by the fitted universality constant η plotted as a function of 1/N c ; lines are the ratios of the quadratic Casimir operators of the adjoint representation over the corresponding ones of the fundamental representation as N c varies. (we note that, for the sake of the visualisation, in this figure we have represented N c as a real number).

Conclusions
We have performed a numerical study of the low-lying spectrum in Sp(2N ) Yang-Mills gauge theories. We have considered the lattice theory formulated with N = 1, 2, 3, 4 and we have measured numerically its torelon and low-lying glueball spectrum as a function of the lattice coupling β. After estimating finite size effects on the target observables by using effective-string-theory motivated predictions applied to torelon masses at N = 3, 4, we have extracted the string tension as a function of β and N from the latter quantities. As a byproduct, through this calculation we have explicitly verified the realisation of the confinement scenario in Sp(6) and Sp(8) Yang-Mills theories, by exposing one of its most typical signatures: the presence of stringy states wrapping compact directions. While this is hardly surprising, direct validation of the expected behaviour in these two gauge theories had never been done before in the literature. We have then extrapolated to the continuum limit the measurements of the adimensional ratios between the glueball masses and the square root of the string tension. Finally, we have obtained the large-N limit of the glueball spectrum in the Sp(2N ) sequence of groups through an extrapolation in a power series in 1/N . For the lowest-lying masses, the leading corrections O(1/N ) to the large-N limit appear to be sufficient to describe the N -dependence down to the smallest value N = 1. We have assessed the size of systematic errors connected with the continuum and the large-N extrapolations and showed that they are negligible at the level of precision of our data.
We have found that, for the states for which the large-N extrapolation in the SU (N ) sequence has been measured, their masses in the large-N limit agree with the ones we have extracted taking the same limit in the Sp(2N ) sequence, as expected. The other states we have determined in this calculation extend our knowledge of the continuum large-N spectrum, therefore providing a more complete set of masses to compare to analytic methods that naturally work at N = ∞, such as gauge-gravity duality techniques. Through an analysis of the ratio of the lowest-lying glueball mass squared and the string tension as a function of N , we have provided further support to the conjecture put forward in Ref. [85], that this ratio is proportional to the ratio of the quadratic Casimir of the adjoint over that of the fundamental representation of the gauge group. Indeed, we have verified that the m 2 0 ++ /σ ratio normalised with the appropriate ratio of quadratic Casimir operators is constant within the Sp(2N ) and the SU (N ) family, and takes compatible values in the two. Our calculation bounds possible N -dependent corrections to this constant to be less than 10%, the latter being the minimum precision with which we have measured the ratio as a function of N .
In addition to the glueball spectrum at finite N ≤ 4, our study has also provided a preliminary investigation of the topological charge in Sp(2N ) gauge theories, in relation to systematics effects in the generation of configurations and in the extraction of spectral masses. An extended analysis of topological observables and a more thorough analysis of topological freezing effects at large N is currently in progress, and will be reported elsewhere.
We envision a number of possible future avenues for exploration, in order to improve and extend this study. Beside the obvious increase in precision that can be obtained by simulating at larger N and smaller lattice spacing (both of which, however, are affected by increased autocorrelation times near the continuum limit and as N grows), one could investigate the effect of including double-torelon and scattering states in the operator basis, in order to have a better resolution of genuine glueball masses. A study of glueball scattering would also provide an extension to the physical reach of our current investigation. Indeed, a scenario in which Sp(2N ) glueballs may play a central role is gluonic dark matter [129,130]. In order to assess the viability of a dark matter scenario based on Sp(2N ) glueballs, one would have to compute the cross section for the decay of the higher-spin glueballs into two scalar glueballs. This (very challenging) calculation would require a dedicated study of multi-point glueball functions. A study of correlators describing glueball scattering would be a natural starting point for such an investigation.
Finally, it is worth reminding the reader that the main motivation for our work has been provided by our ongoing investigation of the pseudo-Nambu-Goldstone-Boson mechanism of electroweak symmetry breaking based on the SU (4) → Sp(4) global symmetry breaking pattern in Sp(2N ) gauge theories with two fundamental Dirac flavours, following the programme outlined in Ref. [7]. In this context, the obtained glueball masses in Yang-Mills provide a reference energy scale to compare with the fermionic matter spectrum, both directly, in the quenched approximation [81], or indirectly, as a proxy for masses of equivalent states in the full theory, for which calculations performed in QCD suggest that the presence of dynamical fermions does not shift significantly the mass of gluonic bound states (see, e.g., Ref. [131]). The mild N dependence in the gluonic observables provides a first indication that no large variations will emerge across corresponding relevant physical observables evaluated in different Sp(2N ) gauge theories, as long as the theory is dominated by gluon dynamics, with small numbers of matter fields.

Acknowledgments
The work of EB has been funded by the Supercomputing Wales project, which is partfunded by the European Regional Development  As briefly mentioned in Sect. 2, ensembles of Sp(2N ) configurations distributed according to Eq. (2.6) are obtained from lattice sweeps of single link heat bath (HB) and overrelaxation (OR) updates. In our implementation of these algorithms, we have used an adaptation of the Cabibbo-Marinari method [102] to the case of Sp(2N ). The Cabibbo-Marinari algorithm updates a group matrix via subsequent updates of SU (2) subgroups covering the whole target gauge group. The choice of the set of SU (2) subgroups to update is crucial. For SU (N ), an efficient implementation can be obtained starting from all the Cartan generators (i, j) having 1 on the i-th diagonal element, −1 on the j-th diagonal element (with 1 ≤ i < j ≤ N ) and 0 everywhere else, along with their eigenvectors under conjugate action. Each generates an SU (2) subgroup of SU (N ). Since Sp(2N ) is a subgroup of SU (2N ), the desired set of subgroups can be obtained from the set found for SU (2N ) by excluding the SU (2) subgroups that alter the block structure in Eq. (A.3) of the Sp(2N ) matrices. Chosing a larger set improves the decorrelation of the algorithm. In this work, we used N 2 subgroups.
To better understand how these subgroups are embedded in Sp(2N ) matrices, we reformulate the considerations above in the language of group representations. Each choice of Cartan generators, along with its eigenvectors, exponentiates to a SU (2) subgroup of SU (N ). The elements of the matrices of this subgroup are different from a unit matrix only at the positions {(i, j), (j, i), (i, i), (j, j)}. A SU (2) matrix is thus embedded into a SU (N ) matrix. We denote this embedding with (i, j). The different embeddings (i, j) can be seen as completely reducible representations of SU (2) that are unitarily equivalent to R ⊕ 1 N −2,N −2 , i.e. to the (1, 2) embedding, where R is the 2 × 2 irreducible representation of SU (2). The unitary transformation that maps one embedding into another is the exchange of rows and columns i and j with 1 and 2 respectively. If [N ] SU is the fundamental representation of SU (N ), {2} the fundamental representation of SU (2), all the embeddings above can be decomposed as For the Sp(2N ) case, the allowed SU (2) embeddings must respect the block structure Eq. (A.3). These embeddings can be split in three classes that are not unitarily equivalent. The embedding (1, 2) is unitarily equivalent to the embeddings (i < N, j < N ). Embeddings in this class can be denoted by with a, b ∈ C such that |a| 2 + |b| 2 = 1 and a * b − b * a = 0. There are N (N − 1)/2 of those embeddings. The embedding (1, 2) is unitarily equivalent to the embeddings (i < N, j < N ) . Embeddings in this class can be denoted by with a, b ∈ C such that |a| 2 + |b| 2 = 1 and a * b − b * a = 0. There are N (N − 1)/2 of those embeddings. The embedding (1, 1 + N ) is unitarily equivalent to the embeddings (i, i + N ). These can be denoted by Examples are, for N = 3, with a, b ∈ C such that |a| 2 + |b| 2 = 1 and a * b − b * a = 0. There are N of those embeddings. With our construction, we have identified N 2 embeddings that cover the whole of For this work, we performed one HB iteration followed by 4 OR iterations for each link variable. Repeating these iterations for all the links of the lattice is a lattice sweep. To prevent the desymplectisation and deunitarisation of the configuration caused by the accumulation of numerical error, we reprojected each link of the configuration on the group after each 10 lattice sweeps with a modified Gram-Schmidt algorithm that preserves the Sp(2N ) structure [7].

B Searching for the bulk phase transition
A phase transition taking place in the (d+1)-dimensional classical canonical system defined by Eq. (2.9) is called a bulk phase transition. This transition separates the strong and weak coupling regimes of the theory, limiting the range of β that is analytically connected to the continuum limit. The identification of bulk phase transition points is hence a crucial step for extrapolating numerical data to the a → 0 limit in a controlled way.
In general terms, a bulk phase transition takes place at values of β for which one of the derivatives of Z(β) with respect to β is singular (in the L → ∞ limit). For a system defined by Eq. (2.9), with the action in Eq. (2.6), the first derivative of ln Z(β) corresponds to the expectation value of the average plaquette, and its second derivative to the plaquette susceptibility, As we observed in Section 4.1, P (β) shows a pseudo-inflection point at some value β c of the lattice coupling. This pseudo-inflection corresponds to a maximum χ c (L) of the plaquette susceptibility. If the latter is associated to the smoothing of a proper phase transition, we expect χ c (L) → ∞, and β c (L) → β c , as L → ∞. Conversely, if χ c (L) stays finite when L → ∞, the change from the strong coupling regime to the weak coupling one happens not through a phase transition, but via a crossover. 3 To study the scaling of the maximum of the plaquette susceptibility with the size of the system, we focused our attention on the neighbourhood of the identified pseudocritical coupling β c , computing at various values of β near this coupling for L/a = 8, 12, 16, and collecting measurements of P (β). A total of 3000 data points at each investigated value of β and L were collected, one every five sweeps. The corresponding results for P (β) are reported in Tab 7. The Monte Carlo histories of our simulations were searched for any sign of metastability, which would have signalled a first order phase transition, with negative results. This allows us to exclude the presence of a discontinuous phase transition for both N = 3 and N = 4. The plaquette susceptibility χ p (β) was computed at each β. The results are presented in Fig. 6. At each volume, β c and χ c were estimated using the multihistogram method. The obtained values are reported in Tab. 8. No appreciable scaling of the peak values can be detected as L is increased. Thus, from our data we can conclude that no phase transition is present for N = 3, 4, with the change of behaviour from strong to weak coupling being described by a crossover. These conclusions are in agreement with the findings in Ref. [6] for the case N = 3.
Even if a phase transition is excluded, the presence of a crossover can still affect physical observables near the change of regime. An example of a similar effect in SU (4) Yang-Mills with a fundamental Wilson action is described in Ref. [118], where the effect of the presence of a crossover reflects in a dip of the measured scalar glueball masses at the corresponding values of β. Similar results emerge also in SU (2) with a mixed fundamentaladjoint action [134,135]. Therefore, to extrapolate lattice observables to the continuum limit with a simple and controlled dependence in a √ σ, it is still necessary to be in the weak coupling regime. We achieved this by ensuring that our data points were far enough from the inflection points and then verifying that there was no visible signal of bulk phase contamination in our observables.

C Continuum and infinite volume extrapolations
As mentioned in Sect. 3.6, estimates of glueball and torelon masses obtained from the lattice are affected by systematic errors. We focus on the systematic error caused by working on a lattice of finite size in Sect. C.1 and on the error caused by the discretisation in Sect. C.2.

C.1 Finite-size effects (FSEs)
The spectrum of a theory defined in a finite box of linear size L with periodic boundary conditions depends on L. The problem was studied for instance in Ref. [136], and this dependence was found to be described by Eq. (3.31). The magnitude of the leading finitesize effects decays exponentially as a function of mL, where m is the lightest excitation in the spectrum.
If ma is estimated to a given finite precision, a value L min /a exists such that for L > L min the FSEs on the spectrum are negligible-by which we mean that their size is much smaller than the statistical error. For L > L min , the measured spectral masses can thus be considered as an estimate of the infinite size spectrum at fixed lattice spacing. In the scaling regime, mL is also a constant as a → 0, and once L min /a is found for a value of a, the FSEs will remain negligible as a is taken to 0, provided the physical volume is kept approximately constant in the process. The precise value of L min /a depends on the precision of our estimates and on the theory under scrutiny, and must be determined empirically.
To determine L min and obtain an estimate for the spectrum at infinite size for N = 3, 4, we used the ensembles described in Sect. 4.2. For each ensemble, we determined the glueball spectrum and the string tension. The results are reported in Tab. 9 for N = 3 and Tab. 10 for N = 4. In this Appendix we focus on the lightest channel, which is consistently found to be A + 1 and suffers from the largest FSEs. (Exceptions to this rule can be found, but they can only occur in the small L/a regime, in which we are not interested.) The estimates of am A + 1 are presented in Fig. 7, for N = 3 in the top panel, and N = 4 in the bottom panel. From these figures we see that am A + 1 rapidly settles on a plateau as L/a is increased. This means that, as expected, FSEs become negligible as L is increased. A rough estimate yields L min /a = 20 for Sp(6) at β = 16.5 and 10 for Sp(8) at β = 26.7. As an a posteriori check, note that m A + 1 L min ∼ 9.76 for Sp(6) and m A + 1 L min ∼ 6.94 for Sp (8). The infinite size spectrum can then be estimated by any one of the results at L > L min . We fitted Eq. (3.31) to the data using b and m(∞) as fitting parameters. The fitted curves and the related χ 2 /N d.o.f. are displayed in Fig. 7.
From the analysis above, we conclude that FSEs are negligible when L > 20a, for N = 3, and L > 10a, for N = 4. On these lattices, the condition L √ σ ≥ 3, which identifies the large volume regime of torelons, is also fulfilled. Hence, we choose this condition throughout as an indicator that finite volume effects can be neglected.  (14) Table 10. Estimates of the masses of glueballs ma in all symmetry channels R P in units of the lattice spacing at β = 26.7 for various L/a and for N = 4.

C.2 Continuum limit extrapolations
The behaviour of the discretisation error was studied in Sect. 3.6. The dimensionless ratio m R P / √ σ behaves, at leading order in a, as where c R P is a constant that depends on the symmetry channel and on the excitation number. The multiplicative term on the right hand side is the continuum limit of the ratio, while the term in σa 2 describes the deviation with respect to this limit for sufficiently small a.
The continuum limit of the spectrum of glueballs can be obtained from sets of estimates thereof obtained at finite lattice spacing with a fit of Eq. (C.1) to the data. This is the For each N = 1, 2, 3, 4, ensembles of 10, 000 thermalised configurations were obtained at several values of a and L/a and stored for later analysis. The values of L/a were always chosen so that FSEs could be safely neglected. This has been verified a posteriori from the measured values of m R P L. The glueball and torelon masses were estimated in units of the lattice spacing, as explained in Sect. 3, for each ensemble.
While not strictly related to the continuum extrapolation, a comment is in order regarding the estimation of the uncertainty on m R P a and to guide the Reader in navigating the tables of results. To determine t min , defined in Sec. 3.1, the quantity m eff (t) is computed on all the available range of t/a. If a plateau can be found, the fits of Eq. (3.24) over the range t > t min provide an estimate of m R P a together with the statistical error of the fit and the corresponding χ 2 /N d.o.f. . It is often the case, however, that the plateau is only 2a − 3a long, that an accurate determination of its preimage in t/a is hindered by the contamination from larger mass states, or that the mass itself is large. These difficulties in determining t min lead to a systematic error on m R P a that can be larger than the statistical error of the fitting procedure. In such cases, the statistical error of the fit cannot be trusted in describing the fluctuations of m R P a. Hence, we use a safe estimate of the mass and its error from the envelope of the points at plateau. A value for χ 2 /N d.o.f. cannot be defined and the corresponding entry in the table is left empty. Finally, there are extreme cases for which an estimate for the mass simply cannot be found, i.e. a plateau is absent. In that case, the corresponding entry in the table is left empty.
All our estimates are reported in Tabs. 11-22 for the ensembles with N = 1, 2, 3, 4. The values of β and L/a are found in the first row of each table; the subsequent rows correspond to the 10 symmetry channels, until the last row, that corresponds to the string tension. For each value of N , these estimates are plotted as a function of σa 2 in Figs. 8-11. In general, we found that the statistical errors and the confidence intervals are of the same order of magnitude, the latter being slightly larger in the majority of cases. This can be taken as an indication of the correctness of the method detailed above. In the following, we refer to the uncertainty in the determination of m R P a generically as its "error".
Let us now comment on the features of these finite-a estimates that are common across all the values of N . The fact that in every case m R P L > 3 allows us to safely neglect FSEs for all the ensembles, as anticipated. Moreover, all the estimates satisfy m R P a ≤ 2, except for the two roughest lattices when N = 1, corresponding to β = 2.2986 and β = 2.3726. Therefore, we can assert that our choices of β are well calibrated to study the flow to the continuum limit of the spectrum of these systems. In the glueball sector, the A + 1 is consistently the lightest channel, followed by the (E ± , T ± 2 ) degenerate pairs. The error of the estimates is larger for larger m R P a, as is to be expected on the basis of the discussion in Sect. 3.6. The E ± and T ± 2 are degenerate over the whole range of a probed at least at the 2σ level, with the mass difference being below one standard deviation in most of the cases. This is a non-trivial a posteriori test of the restoration of continuum rotational invariance and can be taken as a signal that we are in the regime for which Eq. (3.32) is valid.
An additional source of systematic error, the effects of which are difficult to account for, is the autocorrelation time of the system, which grows as the continuum limit is approached. Since the topological charge is one of the quantities with the longest autocorrelation time, studying the evolution of this observable yields a conservative estimator of these effects. Particular attention should be paid to cases in which the topology is (nearly-)frozen. These can be detected by analysing the time series of the topological charge. To this end, a subset of configurations were obtained from the N = 3 and N = 4 ensembles, by picking one configuration each 100. The Gradient flow was then used to smooth each of the configurations, and the regularised topological charge was computed using Eq. (2.11) at each smoothing step. The continuum topological charge is obtained when the regularised topological charge reaches a plateau under further smoothing operations. The results of this analysis are visible in Fig. 12 and Fig. 13.
For both N = 3 and N = 4, there is a value β min above which the topological charge barely changes with the Monte Carlo steps that we are able to perform. These ensembles β = 2.2986  Given that topological freezing affects only our two largest values of N , where the systematic effects it induces on measurements of masses are expected to become less severe (as discussed in Sect. 3.6), we included the estimates obtained from these frozen ensembles in the extrapolation to the continuum limit.
A related potential source of systematic error lies in the length of the initial thermalisation. Our earlier estimates of the continuum spectrum, especially for N = 3 and N = 4, presented a visible dip in the calculated masses for the smallest values of σa 2 . This urged us to perform an overall check of the invariance of the final result under the prolongation of the simulation trajectory. In Fig. 14, we show results of m R P / √ σ at finite lattice spacing as a function of the initial thermalisation time. The fact that these estimates are largely independent of this initial thermalisation time suggests that the Markov chains from which the final averages are obtained are long enough for the system to be at statistical equilibrium.
Let us now discuss the continuum extrapolations of the ratios m R P / √ σ for given values of a. These ratios can be easily formed for each ensemble from the estimates in Tabs Table 14. Estimates of glueball masses and string tensions for N = 2, in units of the lattice spacing a, on lattices of linear size L and lattice spacing determined by the inverse coupling β. The error in brackets are discussed in the main text.    For the new ensembles, the continuum and large-N extrapolated estimates can be found  (14) 1.02 A + *   Fig. 9, and the large-N extrapolation is shown in Fig. 16.
The old ensembles have been reanalysed following the approach used in this work-see  E On the inclusion of N = 1 in the large-N extrapolation In Sect. 3, the large-N extrapolation of the spectrum has been provided including the value N = 1 for all channels. For a handful of channels, this gives a value of χ 2 /N d.o.f. above 2, indicative of a lower statistical significance of the extrapolation. This may suggest that, for these specific channels, the N = 1 value is not captured by the large-N expansion. Indeed, in previous studies of SU (N ) gauge theories (e.g., Ref. [2]), the value of the χ 2 /N d.o.f. has been used as an indication of the reliability of the truncation of the large-N series at a given order for capturing results at some finite value of N . In our work, excluding the data for N = 1 generally improves the value of the χ 2 /N d.o.f. . However, this leaves only three points for the extrapolation, and hence creates a larger systematic bias on the latter. Likewise, adding a higher order correction will decrease the number of degrees of freedom and hence introduce more noise. Being faced with the necessity to make a choice, we have opted to systematically include N = 1 in all large-N extrapolations. This means that we interpret larger value of the χ 2 /N d.o.f. as results of fluctuations in the data or of some  These are obtained from the data generated for this work (new ensembles). In the right column, the extrapolation to N = ∞ obtained from fits of Eq. (4.6) to the data using the data in the left column for N = 2 and the same data as before for N = 1, 3, 4.
unknown systematics, rather than as stemming from the fact that N = 1 is not described by the expansion. The question is left open by this study. For completeness, we compare in Tab. 25 our results for the extrapolations with N = 1 systematically included and excluded.
Most of the results are compatible at the two sigma level.   Table 24. Calculations of the masses in the continuum limit for N = 2 and each channel, in units of √ σa and m E + , using only a reanalysis of the N = 2 data from Ref. [7] (old ensembles), as explained in the text.        and E ++ channels for SU (∞) (borrowed from [2]).   Only the data produced for Ref. [7] were used (old ensembles).
For each symmetry channel R P , the value at σa 2 = 0 is the continuum limit, obtained from a best fit of Eq. (C.1) to the data. The best fits lines are represented as solid lines. The extrapolated results are reported in Tab. 24  Figure 18. Spectrum of the Sp(2N ) theory in the continuum limit for N = 1, 2, 3, 4 and N = ∞, in units of √ σ from the data collected for Ref. [7]. For ease of comparison, we have reported also the masses of the A ++ 1 and E ++ channels for SU (∞) (borrowed from [2]). For N = 2 only the data created for Ref. [7] were used. The point corresponding to 1/2N = 0 is the value of m R P / √ σ(∞) obtained from the best fit of Eq. (4.6) to the data.