Sp(4) gauge theory on the lattice: towards SU(4)/Sp(4) composite Higgs (and beyond)

The Sp(4) gauge theory with two Dirac fundamental flavours provides a candidate for the microscopic origin of composite-Higgs models based on the SU(4)/Sp(4) coset. We employ a combination of two different, complementary strategies for the numerical lattice calculations, based on the Hybrid Monte Carlo and on the Heat Bath algorithms. We perform pure Yang-Mills, quenched computations and exploratory studies with dynamical Wilson fermions. We present the first results in the literature for the spectrum of glueballs of the pure Sp(4) Yang-Mills theory, an EFT framework for the interpretation of the masses and decay constants of the lightest pion, vector and axial-vector mesons, and a preliminary calculation of the latter in the quenched approximation. We show the first numerical evidence of a bulk phase transition in the lattice theory with dynamical Wilson fermions, and perform the technical steps necessary to set up future investigations of the mesonic spectrum of the full theory.

Abstract: The Sp(4) gauge theory with two Dirac fundamental flavours provides a candidate for the microscopic origin of composite-Higgs models based on the SU (4)/Sp(4) coset. We employ a combination of two different, complementary strategies for the numerical lattice calculations, based on the Hybrid Monte Carlo and on the Heat Bath algorithms. We perform pure Yang-Mills, quenched computations and exploratory studies with dynamical Wilson fermions.
We present the first results in the literature for the spectrum of glueballs of the pure Sp(4) Yang-Mills theory, an EFT framework for the interpretation of the masses and decay constants of the lightest pion, vector and axial-vector mesons, and a preliminary calculation of the latter in the quenched approximation. We show the first numerical evidence of a bulk phase transition in the lattice theory with dynamical Wilson fermions, and perform the technical steps necessary to set up future investigations of the mesonic spectrum of the full theory. 1 Introduction and motivation The study of Sp(2N ) gauge theories has the potential to unveil many new phenomena of general relevance in field theory and its phenomenological applications within high energy physics. The recent progress in lattice gauge theory makes it an ideally suited nonperturbative instrument for this type of investigation. The literature on the subject is somewhat limited (see for instance [1]). There is a large number of questions that we envision can be answered with dedicated lattice studies, and in this introduction we list them and discuss the long-term research programme that this work initiates, before specialising to the investigations and results we will report upon in this paper.

The Sp(2N ) research programme
In the context of physics beyond the standard model (BSM), the discovery of the Higgs particle [2,3], combined with the absence of evidence for new physics at the TeV scale from LHC direct searches, exacerbates the little hierarchy problem. If the mass of the Higgs particle has a common dynamical origin with hypothetical new physics at multi-TeV scales, the low-energy effective field theory (EFT) description of the system is in general unnatural (fine-tuned). The framework of Higgs compositeness we refer to in this paper  addresses this problem by postulating the existence of a new underlying strongly-coupled theory, in which an internal global symmetry is broken spontaneously by the dynamicallygenerated condensates, resulting in a set of parametrically light pseudo-Nambu-Goldstone bosons (PNGBs). One writes their EFT description in terms of scalar fields, and weakly couples it to the standard-model gauge bosons and fermions. Four of the PNGBs of the resulting EFT are interpreted as the Higgs doublet, hence providing an elegant symmetry argument for the lightness of the associated particles. Particular attention has been devoted to models based on the SU (4)/Sp(4) coset [8, [26][27][28][29][30][31][32][33][34][35][36][37], as EFT arguments suggest that the resulting phenomenology is both realistic and rich enough to motivate a more systematic study of the underlying dynamics. This coset emerges naturally in gauge theories with pseudo-real representations, such as Sp(2N ) gauge theories with two massless Dirac fermions in the fundamental representation of the gauge group. Phenomenological arguments-ultimately related to the fact that if fundamental fermions carry SU (3) c (colour) quantum numbers, one can further implement partial top compositeness-select SU (2) ∼ Sp (2) and Sp(4) as most realistic viable candidates for BSM physics. A number of studies has considered the dynamics of SU (2) (see for instance [28][29][30][31][32][33][34]), while here we focus on Sp (4).
The primary objectives of our research programme include the study of the mass spectrum of mesons and glueballs, and the precise determination of decay constants and couplings of all these objects by means of lattice numerical techniques. 1 Eventually, we want to gain quantitative control over a large set of measurable quantities of relevance for phenomenological (model-building) considerations, which include also more ambitious determinations of the width of excited states, and of the values of the condensates in the underlying dynamics.
A separate set of objectives relates to the physics of baryons and composite fermions. In Sp(2N ) with fundamental matter, baryons are bosonic objects, and hence not suitable as model-building ingredients for top (partial) compositeness. Composite fermions could be realised for example by adding 2-index representations to the field content of the dynamics (see for example [26,27,35]). We envision to start soon a non-trivial programme of study of their dynamical implementation on the lattice.
A third set of dynamical questions that our programme wants to address in the long run pertains to the thermodynamic properties of the Sp(2N ) theories at finite temperature T and chemical potential µ. It is of general interest to study the symmetry-restoration pattern of these models at high temperature (see [34] for a step in this direction in the case of SU (2)). Furthermore, the pseudo-real nature of Sp(4) makes it possible to study the phase-space of the theory, while avoiding the sign problem.
Finally, there is a different field-theoretical reason for studying Sp(2N ) gauge theories. It is known that the Yang-Mills theories based on SU (N ), SO(N ) and Sp(2N ) all agree with one another on many fundamental physical quantities when the limit of large N is taken. Lattice results allow for non-trivial comparisons with field-theory and string-theory studies in approaching the large-N limit. While there is a substantial body of literature on SU (N ) theories on the lattice [39], for example for the calculation of the glueballs, and some literature on SO(N ) models (see for instance [40][41][42] and references therein), there is no systematic, dedicated study of the Sp(2N ) gauge theories. We aim at comparing with results in Yang-Mills theories based on other groups, and with conjectures such as those put forwards in [43] and [44].

Laying the foundations of the Sp(4) lattice studies
With this paper, we start the programme of systematic lattice studies of the dynamics of such gauge theories. We focus here on the Sp(4) gauge group, which is of relevance for the phenomenology of composite-Higgs models. We perform preliminary studies of the lattice theories of interest, a first exploratory computation of the meson spectrum in the quenched approximation and a first test of the same calculation with dynamical fermions.
We aim at gaining a quantitative understanding of the properties of the bound states, possibly describing them within the EFT framework. Starting from the leading-order chiral-Lagrangian description of the PNGBs, we extend it to include heavier mesons, aiming at providing dynamical information useful for collider searches. As is well known, the description of the spin-1 composite particles is weakly-coupled only in the large-N limit: we do not expect the EFT to fare particularly well for Sp (4), yet it is interesting to use it also in this case, in view of possible future extensions to Sp(2N ). We also begin the analysis of the next-to-leading-order corrections to the chiral Lagrangian of the model, as a first preliminary step towards understanding realistic model building in the BSM context. While still beyond the purposes of this paper, we find it useful also to briefly summarise the main goals of the exploration of the dynamics of composite fermions emerging from introducing on the lattice matter in different representations.
We devote a significant fraction of this paper to the study of the pure Yang-Mills theory. We perform our Sp(4) lattice calculations in such a way that the technology we use can be easily generalised to any Sp(2N ) theories, in view of implementing a systematic programme of exploration of the large-N behaviour. We present the spectrum of glueballs, and study the effective string-theory description, for Sp(4) pure Yang-Mills, with no matter fields. Our results have a level of accuracy that is comparable to the current state-of-the-art for SU (N ) gauge groups in four dimensions. This both serves as an interesting test of the algorithms we use, but also nicely complements existing literature on related subjects.
The paper is organised as follows. In Sec. 2 we define the field theory of interest, and introduce its low-energy EFT description in terms of PNGBs. We also extend the EFT to include the lightest spin-1 states in the theory (see also Appendix A and B). We define the framework of partial top compositeness for these models, and the lattice programme that we envision to carry out in the future along that line.
In Sec. 3 we describe in details the lattice actions, as well as the Heat Bath (HB) and Hybrid Monte Carlo (HMC) algorithms used in the numerical studies. In Sec. 4 we focus on scale setting and topology. These two technical Sections, together with Appendix A and C, set the stage not only for this paper, but also for future physics studies we will carry out. In Sec. 5 we present the spectrum of glueballs for Sp (4). We also explain in details the process leading to this measurement, that will be employed in the future also for the spectrum of Sp(2N ) with general N . Section 6 is devoted to the spectrum of mesons of Sp(4) in the quenched approximation, the extraction of the masses and decay constants and a first attempt at comparing to the low-energy EFT. Preliminary (exploratory) results for the full dynamical simulation are presented in Sec. 7. In particular, we exhibit the first (to the best of our knowledge) evidence that a bulk phase transition is present for Sp(4) with fundamental matter. We conclude the paper with Sec. 8, summarising the results and highlighting the future avenues of exploration that this work opens.
2 Elements of field theory, group theory and effective field theory The Sp(4) gauge theory of interest has matter content consisting of two (massive) Dirac fermions Q i a , where a = 1 , · · · , 4 is the colour index and i = 1, 2 the flavour index, or equivalently four 2-component spinors q j a with j = 1 , · · · , 4. The Lagrangian density is where the summations over flavour and colour indices are understood, and where the fieldstrength tensors are defined by In the m → 0 limit, the global symmetry is U (1) A × SU (4). The presence of a finite mass m = 0 is allowed within the context of composite-Higgs models, and may play an ≡ Ω ab q n a TC q m b , (2.3) so that in 2-component spinor language the mass matrix is M ≡ m Ω. We collect in Appendix A some useful elements of group theory. For the most part we ignore the anomalous U (1) A . We display the field content in Table 1, where we list also the transformation properties of the composite field Σ 0 , and the (symmetry-breaking) spurion M . The vacuum is characterised by the fact that 0 = Σ 0 ∝ Ω, hence realising the breaking SU (4) → Sp(4). In the absence of coupling to the SM fields, the vacuum structure aligns with the mass term, which hence contributes to the masses of the composite states, by breaking the global SU (4) while preserving its global Sp(4) subgroup. As a result, the lightest mesons of the theory arrange themselves into irreducible representations of Sp(4): the PNGBs π and axial-vectors a 1 transform on the 5 representation of the unbroken Sp(4), while the scalars a 0 and the vectors ρ on the 10 of Sp(4). There exist also the corresponding scalar, pseudo-scalar, vector and axial-vector Sp(4) singlets, but we will not discuss them in this paper.

EFT analysis
The EFT treatment of the lightest mesons depends on the coset, with only numerical values of the coefficients depending on the underlying gauge group. We summarise here some useful information about two different EFTs. Some of the material collected in this subsection can also be found in the literature [8, 11-18, 26-32, 34]. We construct the chiral Lagrangian for the SU (4)/Sp(4) coset, and its generalisation in the sense of Hidden Local Symmetry (HLS) [45][46][47][48][49] (see also [50][51][52][53]). The former assumes that only the pions are dynamical fields in the low-energy EFT, while the latter includes also ρ and a 1 as weakly-coupled fields. We will comment in due time on the regime of validity of the two.

EFT description of pions.
The low-energy EFT description of the pions π is constructed by introducing the real composite field Σ that obeys the non-linear constraint ΣΣ † = I 4 , and transforms as the antisymmetric representation Σ → U ΣU T , under the action of an element U of SU (4). The antisymmetric vacuum expectation value (VEV) Σ ∝ Ω breaks SU (4) to the Sp(4) subgroup, and as a result five generators T A , with A = 1, · · · , 5, are broken, while 10 other T A , with A = 6, · · · , 15, are not. We normalise them all as Tr T A T B = 1 2 δ AB . The field Σ contains the PNGB fields π = π A T A , conveniently parametrised as in terms of which, at leading-order, the EFT has the Lagrangian density (2.6) The pion fields are canonically normalised, hence f = f π is the pion decay constant. If it were promoted to a field, the spurion M would transform as M → U * M U † , so that the combination Tr M Σ would be manifestly invariant under the SU (4) global symmetry. The (symmetry-breaking) mass term is hence written as which confirms that the 5 pions are degenerate in the presence of the explicit breaking given by the Dirac mass for the fermions, because of the unbroken SO(5) ∼ Sp(4) symmetry. The GMOR relation is automatically satisfied, and takes the form: As is the case for the chiral-lagrangian description of low-energy QCD, we are making use of two expansions: the derivative expansion, that suppresses terms of higher dimension, and that is reliable provided we consider observables at energies E 4πf π , and the chiral expansion, reliable when the mass of the pion satisfies m π 4πf π . At the sub-leading order, we could for example add to the chiral Lagrangian the contribution (2.10) (2.11) The first term of the expansion comes from setting Σ = Ω, and amounts to a m-dependent rescaling of the interacting L 0 . No correction to the mass term appears, but just an overall rescaling of both f and π, so that the GMOR relation is respected. However, an additional quartic pion coupling is generated, that contributes to the ππ scattering amplitude. In the massless limit f π can be equivalently extracted from either 2-point functions or from the ππ scattering amplitude. For m = 0 there are subtleties: in the following we will always extract f π from the q 2 → 0 limit of the 2-point functions, and we will highlight this fact by denoting the pion decay constant (squared) as f 2 π (0) in the rest of the paper.

Hidden Local Symmetry description of ρ and a 1
Hidden Local Symmetry provides a way to include spin-1 excitations such as the ρ mesons into the EFT treatment, hence extending the validity of the chiral Lagrangian (see for instance [45][46][47][48][49] and also [50][51][52][53]). While very appealing on aesthetics grounds, when applied to QCD such idea shows severe limitations: the heavy mass and non negligible width of the ρ mesons imply that the weak-coupling treatment is not fully reliable. Yet, this description offers a nice way to classify operators and it is expected to become more reliable at large-N [49]. As we envision future studies with larger Sp(2N ) groups, it is useful to show the construction already in the programmatic part of this paper, and test it on Sp(4).  [45][46][47][48][49] (see also [50][51][52][53]). The two sites represent the SU (4) A × SU (4) B global symmetry. The scalar S transforms on the bifundamental representation, while Σ is antisymmetric. The SU (4) A group is gauged with coupling g ρ , while the SU (2) L ×U (1) Y SM subgroup of SU (4) B can be weakly gauged, with couplings g and g . In most of this paper we set g = 0 = g , and hence ignore the interactions of the strongly coupled dynamics with the standard-model fields.
The full set of ρ and a 1 mesons spans the adjoint representation of the global SU (4) symmetry. An EFT description of their long-distance dynamics can be built starting from the diagram in Fig. 1. The 15 spin-1 fields are introduced as gauge fields of SU (4) A . Two scalars, the antisymmetric Σ of SU (4) A , and the bi-fundamental S, transform as under the action of U A ∈ SU (4) A and U B ∈ SU (4) B . The VEV of S breaks the enlarged symmetry and provides mass for all the vectors. S is subject to the constraints S † S = I 4 , that are solved by parametrising S = e 2iσ F , with σ = σ A T A and F the decay constant. At the same time, we parametrise Σ = e 2iπ f Ω, in such a way that the two scalars together implement the breaking SU (4) A × SU (4) B → Sp(4), and describe the 15 exact Goldstone bosons that are higgsed away into the longitudinal components of the massive spin-1 states, as well as the remaining 5 (massive) PNGB, denoted asπ A in the following.
In composite Higgs models, the SM gauge group SU (2) L × U (1) Y is a subgroup of SU (4) B , and it is chosen to be a subgroup of the unbroken global Sp(4). The covariant derivative of S is given by where we have made use of the fact that Ω 2 = −I 4 . From this point onwards, we set g = 0 = g , and focus on the dynamics of the strongly-coupled new sector in isolation from the SM fields. We write the Lagrangian density describing the 15 gauge bosons A A µ , as well as the 20 pseudo-scalar fields π A and σ A , as In the expansion, we include two sets of operators. We call leading order (LO) the ones appearing in the first four lines, controlled by the parameters F , f , b, c, κ, g ρ and v. This is an exhaustive set of operators, at this order. We call next-to-leading order (NLO) those in the last four lines, controlled by the parameters v 1 , v 2 , y 3 , y 4 and v 5 . As we will discuss shortly, the list of NLO operators is incomplete. In total there are 12 parameters. This Lagrangian has to be used with caution. The appearance of ρ and a 1 fields in the EFT is fully justified only if the coupling g ρ is small, which must be discussed a posteriori, yet is expected to hold in the large-N limit, and as long as m is small.
The last four lines of Eq. (2.16) contain terms that are sub-leading in the powercounting. Because we are going to perform lattice simulations at finite mass m, a priori we do not know how important such terms are, and hence we include them. While nonvanishing values of m are allowed within the composite-Higgs framework, the EFT is useful only when m is small enough that truncating at this order is justified.
We do not include the full set of sub-leading four-derivative terms, because they are not important for our current purposes. These terms would become important when a complete analysis of 3-point and 4-point functions is performed, for example. We also omit topological terms. Furthermore, we do not include in L terms with the structure of Eq. (2.10), such as We will comment later in the paper on the implications of all these omissions. We conclude this subsection with a technical comment. Some of the terms in the Lagrangian density in Eq. (2.16) involve only nearest-neighbour interactions, in the sense of the diagram in Figure 1, while other couplings introduce non-nearest-neighbour interactions. Such additional interactions might for example emerge from the process of integrating out heavier degrees of freedom. One of the big limitations of the HLS language is that the number of independent, admissible such non-nearest-neighbour interactions grows rapidly with the number of fields in the theory, and hence by introducing more resonances the EFT Lagrangian density loses predictive power because of the proliferation of new free parameters. In the special Lagrangian we wrote, such interactions are controlled by the parameters κ, b, c, as well as v, v 1 , v 2 , y 3 , y 4 and v 5 . If only nearest-neighbour couplings were to be allowed, the set of parameters would be restricted to just f , F and g ρ , at this order in the expansion.

2-point functions
To compute masses and decay constants of the mesons, we use the language of the SU (2) t L × SU (2) t R symmetry that would be of direct relevance if we were to treat this as a Technicolor model. In particular, this symmetry is not a subgroup of the unbroken Sp(4) global symmetry, and the condensate breaks it. We treat this as a technical tool, that is convenient in order to extract physical quantities from the correlation functions. Yet our results hold also for finite m, and apply as well to the composite-Higgs scenario, as we never include in the calculations the effects of the couplings to the external (weakly-coupled) SM fields.
The Left-Left current-current correlator is (see Appendix B) from which one can read that the masses and decay constants are given by and that the pion decay constant is As anticipated, the notation explicitly specifies that f 2 π (0) is extracted from 2-point functions evaluated at q 2 = 0. The fact that f 2 0 = f 2 π (0) + f 2 ρ + f 2 a 1 is independent of m is the accidental consequence of the truncation we made, in particular of the omission of the operator in (2.17). Whether or not this is justified, depends on the range of m considered and on the size of the EFT coefficients, as emerging from lattice data.
One can compute the right-hand-side of the first and second Weinberg sum rules, within the EFT, to obtain hence showing explicitly that the non-nearest-neighbour couplings b, c, κ, y 4 , y 3 , v 1 and v 2 yield to direct violations of the Weinberg sum rules, within the EFT. As anticipated, this is not surprising: non-nearest-neighbour interactions are expected to emerge from integrating out heavy degrees of freedom, and result in the violation of the Weinberg sum rules because their rigorous derivation involves summing over all possible resonances. The additional couplings, in effect, parameterise the contribution to the sum rules of heavier resonances that have been omitted.

Pion mass and g ρππ coupling
To compute the physical mass and couplings of the pions, it is convenient to fix the unitary gauge, by setting σ A = 0 along the unbroken A = 6 , · · · , 15 generators, and for A = 1 , · · · , 5 along the broken generators, with the normalisation factor N chosen so that the physicalπ A are canonically normalised: The 5 degenerate pions have mass which modifies the GMOR relation to read The g ρππ coupling is conventionally defined by the Lagrangian density so that the width (at tree level) is . We find that (2.33) We conclude with a comment about unitarity. While the calculations performed here make use of the unitary gauge, we must check that the kinetic terms of all the Goldstone bosons be positive before setting to zero the linear combinations providing the longitudinal components of the vectors. We call the relevant normalisations k 10 , k 5 and k 5 , coming from the kinetic term of σ A with A > 5, as well as the trace and the determinant of the kinetic matrix mixing σ A and π A with A < 6. Such combinations are explicitly given by: We require that k 10 , k 5 , k 5 > 0. Furthermore, for the kinetic terms of the vectors to be positive definite one must impose κ + m y 4 < 1 and κ + m y 3 > −1.

On the regime of validity of the EFT
In the EFT we wrote to include the ρ and a 1 particles, we are making use of several expansions. Besides the derivative expansion and the expansion in the mass of the fermions m, appearing also in the chiral Lagrangian, there is a third expansion, that involves the coupling g ρ and deserves discussing in some detail. From lattice calculations of 2-point functions, one extracts the decay constants of π, ρ and a 1 , in addition to the masses. In the m = 0 limit, the expressions for the five quantities M ρ , M a 1 , f ρ , f a 1 and f π depend on the six free parameters f , F , b, c, κ and g ρ , that hence cannot all be determined. Let us choose to leave κ undetermined, for example, and solve the algebraic relations for the other five parameters in terms of the physical quantities. The g ρππ coupling can then be written as If one were restricted to the massless theory, only by gaining access to 3-point functions could one measure κ. Yet, detailed information about the m-dependence of 2-point functions can be used to predict g ρππ , and the width of the ρ meson. In principle, the width of the ρ meson Γ ρ could be compared with the physical width extracted from lattice calculations [55,56]. In this way we would be able to adjudicate explicitly whether the weak-coupling assumption that underpins this EFT treatment is justified. However, the direct extraction of Γ ρ from lattice data is highly non-trivial, and will require a future dedicated study.
There is no reason a priori to expect that g ρ , or g ρππ , be small, except in the large-N limit. The fact that from 2-point functions we can infer some of the properties of the EFT that enter the 3-point functions holds only provided the coupling is small, with g 2 ρππ /(48π) 1. Furthermore, if M ρ and M a 1 happen to be large in respect to f π , bringing them close to the natural cut-off set by the derivative expansion, it would again signal a break-down of the perturbative expansion within the EFT.
Nevertheless, even in the regime of large g ρ , we can still learn something from the expansion in small mass m. In particular, we should be able to use the EFT to reproduce the m-dependence of masses and decay constants, at least in the small-m regime. In the future, we envision repeating the study performed in the following sections for Sp(2N ) theory with dynamical fermionic matter, and with larger values of N , and hence track the N -dependence of the individual coefficients.

Spin−1/2 composite fermions and the top partner
In composite-Higgs models, the generation of the SM fermion masses is often supplemented by the mechanism of partial compositeness (PC). The SM fermions, in particular the top quark, mix with spin-1/2 bound states emerging from the novel strong-interaction sector (the composite sector), and phenomenologically this allows both to enhance the fermion mass (as in precursor top-color models) and to trigger electroweak symmetry breaking via vacuum (mis)alignement. As an example, we borrow some of the construction in [8] and [27]. So many other, equally compelling, examples exist in the literature, that we refer the reader to the review [16] and to the references therein. 3 Let us assume that the microscopic theory admits the existence of Sp(4)-colour singlet operatorsΨ i andΨ c i , that have spin-1/2, carry SU (3) c colour and, combined, span vectorial representations of the SM gauge group. The index i = 1, 2 refers to the SU (2) L singlets and doublets, respectively, and the notation refers to the fact that we write the operators as 2component fermions. Let us now consider the low-energy description of the lightest particles excited from the vacuum by such operators, and write it in terms of new 2-component spinorial fields Ψ i and Ψ c i with the same quantum numbers asΨ i andΨ c i . Coarse-graining over model-dependent details, Ψ i and Ψ c i have the correct quantum numbers to couple to the SM quarks, in particular to the SM top quark, represented by the 2-component Weyl fermions t and t c , provided Ψ i transforms on the fundamental of SU (3) c and Ψ c i on its conjugate.
Below the electroweak symmetry-breaking scale v W , the mass terms take the form where λ 1 , λ 2 , λ and y are dimensionless couplings, M * represents the typical scale of the masses of composite fermions in the Sp(4) gauge theory and Λ represents the underlying scale at which (third-generation) flavour physics arises (see also [8]). d Ψ = d Ψ c is the dimension of the operatorsΨ andΨ c in the underlying theory. Diagonalisation of the resulting mass matrix, under the assumption that yv W be small in respect to the other scales, yields two heavy Dirac masses approximately given by and finally the mass (squared) of the top is given approximately by (2.39) In order to assess the viability of these models, one needs to provide a microscopic origin for all of the parameters appearing in Eq. (2.36). To do so, one must specify the (modeldependent) microscopic details controlling the nature of the composite fermions. Spin-1/2 composite Sp(4)-neutral particles arise in the presence of fermions in higher-dimensional irreducible representations. As an example, ref. [27] proposes to extend the field content of the microscopic theory in Table 1 to include 2-component elementary fermions χ (and χ c ) in the antisymmetric representation of the gauge Sp(4), transforming as singlets of the global SU (4), and on the fundamental (and anti-fundamental) representation of the SU (3) c gauge symmetry of QCD.
The χ and χ c fermions carry QCD colour charge, which allows to construct coloured composite states in the antisymmetric, six-dimensional representation of the global SU (4) group, by coupling them to a pair of fundamental fermions q. For example, the operatorŝ Ψ andΨ c aforementioned can be obtained aŝ where summations over Sp(4) gauge indices are understood, while we show explicitly the (antisymmetrised) global SU (4) indices a and b, and the SU (3) c colour index α.
One of our long-term goals is to study the PC mechanism with lattice simulations, which requires generalising the lattice study we will report upon in the following sections to the case in which the field content contains at least two species of fermions transforming in different representations of the fundamental gauge group. The example we outlined here, though incomplete, immediately highlights how, from the phenomenological perspective, the determination of the masses of the top partners (the scale M * and couplings such as λ, as a function of the elementary-fermion mass m) in the PC mechanism are of direct interest, as they represent a way to test the theory. At the same time, they are accessible on the lattice, even without introducing (model-dependent) couplings to the SM fields.
The other additional, essential, input from non-perturbative dynamics of the microscopic theory is the anomalous dimension of the top-partner operators, such asΨ andΨ c in the example. For the PC mechanism to be valid, in principle one needs the operator dimensions to be small, for example d Ψ ≤ 5/2, which implies that the operatorΨ T 1C t c is relevant in the IR, and that the anomalous dimensions of the candidate operators have to be non-perturbatively large. In practice, since Λ/M * is not infinity, this requirement may be relaxed, at the price of admitting some degree of fine-tuning.
Finally, the (model-dependent) extension of the field content, required by the PC mechanism, also implies the enlargement of the global symmetry, and additional light PNGB's, some of which are neutral, some of which carry SU (3) c colour, and many of which may be lighter than the typical scale of the other composite particles. Lattice calculations of the masses of such particles would offer the opportunity to connect with the phenomenology derived from direct particle searches at the LHC.

Numerical lattice treatment
In this Section, we present the discretised Euclidean action and Monte Carlo techniques used in the numerical studies. We adapt the state-of-the-art lattice techniques established for QCD to the two-flavour Sp(4) theory. Pioneering lattice studies of Sp(4) gauge theory without matter can be found in [1]. Numerical calculations are carried out by modifying the HiRep code [57].

Lattice action
For the numerical study of Sp(2N ) gauge theory on the lattice, we consider the standard plaquette action where β = 4N/g 2 is the lattice bare gauge coupling, and N = 2 in the Sp(4) case of this paper. The plaquette P µν is given by where the link variables U µ (x) are Sp(4) group elements in the fundamental representation, whileμ andν are unit vectors along the µ and ν directions.
In the dynamical simulations with two Dirac fermions in the fundamental representation, we use the (unimproved) Wilson action where the massive Wilson-Dirac operator is given by where a is the lattice spacing and m 0 is the bare fermion mass.

Heat Bath
As a powerful way to perform calculations in the pure Sp(4) gauge theory, we implemented a heat bath (HB) algorithm with micro-canonical over-relaxation updates, to improve the decorrelation of successive configurations. As in the case of SU (N ) [58], the algorithm acts in turn on SU (2) subgroups, the choice of which can be shown to strongly relate to the ergodicity of the update pattern. A sufficient condition to ensure ergodicity is to update the minimal set of SU (2) subgroups to cover the whole Sp(2N ) group. This condition can be suitably translated to the algebra of the group and generalised to any Sp(2N ). In the Sp(4) case, of relevance to this paper, we choose to update a redundant set of subgroups, in order to improve the decorrelation of configurations. We provide below a possible partition of the generators used to cover all of the Sp(4) gauge group, written with the notation of [34].
• SU (2) R subgroup, with generators T i R in Eq. (B.7) of [34].   [1]. The symbols P and P H.P.W. denote the measurements from this work and those reported in Ref. [1], respectively. Our lattice volume is V = 8 4 , and both calculations use the HB algorithm with over-relaxation, as explained in Sec. 3.2. Each point has been obtained from 5000 measurements, and the errors are corrected for autocorrelations.
• SU (2) τ subgroup, with generators expressed in terms of B.4 in [34]: • SU (2)τ subgroup, with generators expressed in terms of B.4 in [34]: The set of 10 generators T i L , T i R , τ 1 , 2 andτ 1 , 2 spans the whole Sp(4). The minimal set of 5 elements that generate the whole group by closure consists for example of any two elements T i L , any two elements T j R and one additional element among τ 1 , 2 andτ 1 , 2 . As a check of correctness of the algorithm we employed, we compared the average of the elementary plaquette to the results obtained in [1], as shown in Fig. 2, confirming that they are compatible within the statistical errors.

Hybrid Monte Carlo
In the study of N f = 2 dynamical Dirac fermions, we make use of the hybrid Monte Carlo (HMC) algorithm. As Sp(4) is a subgroup of SU (4), most of the numerical techniques used for SU (N ) with an even number of fermions can straightforwardly be extended to our study. However, there are two distinguishing features.  First of all, in contrast to the HB algorithm, the explicit form of the group generators of Sp(4) is necessary for the molecular dynamics (MD) update. For instance, the MD forces for the gauge fields are given by where V µ (x) is the sum of forward and backward staples around the link U µ (x). The generators T A with A = 6 , · · · , 15 are given in Appendix B of [34]. The group invariant T F is defined as Tr (T A T B ) = T F δ AB , which in our case yields T F = 1/2, so that for Sp(4) the normalization is 2N T F = 2.
Secondly, due to machine precision, it is not guaranteed that the link variables stay in the Sp(4) group manifold. In analogy with the re-unitarization process implemented in SU (N ) studies, we perform a re-symplectisation at the end of each MD step. We describe in Appendix C the procedure, based on Sp(4) projection that makes use of the quaternion basis.
As a further test of this implementation of the HMC algorithm, we calculated the expectation value of the difference of the auxiliary Hamiltonian at the beginning and the end of a MD trajectory ∆H for various values of the integration step size , in the case with β = 6.9 and am 0 = −0.85 on a 24 × 12 3 lattice. We found that ∆H is proportional to 4 , as expected for a second order Omelyan integrator [59], and Creutz's equality exp(−∆H) = 1 [60] is satisfied. We also checked that the average plaquette values are consistent with each other for all values of the step size.
The HiRep code [57] is designed to implement SU (N ) gauge theories with a generic number of colours and flavours, with fermions in any two-index representation. One of its crucial features is that the gauge group and the representation can be fixed at compile time by using preprocessor macros. This provides us with great flexibility in implementing the aforementioned features of Sp(4).
As a nontrivial test of the HMC code, we first calculate the expectation value of the plaquette of the theory with two degenerate, very heavy fundamental fermions (a m 0 = 10.0) and compare the results with the pure Sp(4) results from [1]. In Fig. 3, we plot the differences of the average plaquette values between the two calculations for various values of β. The two series are compatible with each other, within the statistical errors.
It is known that the pure Sp(4) theory in 3 + 1 dimensions exhibits a first-order deconfinement phase transition [1]. Although a finite-size scaling analysis is needed to confirm the existence of the first-order phase transition, for the purpose of a consistency check of the code it is worth showing numerical evidence of the coexistence of the confined and deconfined phases. To this end, we calculate the expectation value of the Polyakov loop averaged over the space-like points, defined by The temperature T in lattice units is identified with the inverse of the extent of the temporal lattice, 1/N t . Near the critical temperature T c , the probability distribution of Φ indeed shows the coexistence of two phases, as in the second panel of Fig

Lattice calibration
This Section is devoted to discuss two lattice technicalities that are important in order to extract the correct continuum physics: we address the problem of setting the scale, using the gradient flow, and study the topology of the ensembles generated by our numerical process, to verify that there is no evidence of major problems in the lattice calculations.

Scale setting and gradient flow
Lattice computations are performed by specifying dimensionless bare parameters in the simulation, and all dimensionful results are extracted in units of the lattice spacing. These results have to be extrapolated to the continuum limit to make impact on phenomenology. It is also desirable to express them in natural units. Such demands make the scale setting an important task in lattice calculations. To carry out this task, the most straightforward approach is to compute a physical quantity on the lattice, and then compare with its experimental measurement. In the absence of experimental results for the Sp(4) gauge , roughly corresponding to temperatures T below, at, and above the critical temperature T c , respectively. theory, one can still accomplish reliable continuum extrapolations by employing a scale defined on theoretical grounds, such that one can determine the ratio a 1 /a 2 of the lattice spacings in two simulations performed at different choices of the bare parameters.
The gradient flow in quantum field theories, as revived in recent years by Martin Lüscher in the context of the trivialising map [61], is a popular method for scale setting [62,63]. To study the gradient flow in a field theory, one first adds an extra dimension, called flow time and denoted by t. An important point articulated by Lüscher is that a field theory defined initially with a cut-off can be renormalised at non-vanishing flow time. In addition, choosing carefully the bulk equation governing the gradient flow, the theory does not generate new operators along the flow time (counter-terms), keeping the renormalisation of the five-dimensional theory simple. 4 The Yang-Mills gradient flow of the gauge field B µ (t, x) is implemented via the equation the corresponding covariant derivative, and A µ (x) the initial gauge field in the four-dimensional theory. Noticing that Eq. (4.1) describes a diffusion process, the flow time t therefore has length-dimension two. It has been shown that, to all orders in perturbation theory, any gauge invariant composite observable constructed from B µ (t, x) is renormalised at t > 0 [63]. In particular, Lüscher demonstrated that the action density can be related to the renormalised coupling, α(µ), at the leading order in perturbation theory through The dimensionless constant k α is analytically computable [62]. Equation (4.2) can actually serve as the definition of a renormalisation scheme: the gradient-flow (GF) scheme. Furthermore, since t 2 E(t) ≡ E is proportional to the GF-scheme coupling, this quantity can be used to set the scale. In other words, if one imposes the condition where E 0 is a constant that one can choose, then √ t 0 should be a common length scale, assuming lattice artefacts are under control. In practice, one measures √ t 0 in lattice units.
That is, one computes √ t 0 /a ≡ t 0 . This allows the determination of the ratio of lattice spacings from simulations performed at different values of the bare parameters. It is worth mentioning that the diffusion radius in Eq. (4.1) is √ 8t, and it is convenient to define the ratio where L is the lattice size.
Given that the right-hand side of Eq. (4.1) is the gradient of the Yang-Mills action, the most straightforward way to latticise it is 5 is the Wilson plaquette action. The gradient flow serves as a smearing procedure for the gauge fields. This means the larger the flow time, the smoother the resultant gauge configurations will be. In other words, the larger the flow time is, the smaller the ultraviolet fluctuations of flown observables. On the other hand, it also means the gauge fields become more extended objects as the flow time grows. This results in longer autocorrelation time, and makes the statistics worse. Furthermore, having c τ > 0.5 can lead to significant finite-volume effects. These are issues one would have to consider carefully when choosing a value for the constant E 0 in Eq. (4.4).
The action density E(t) at non-vanishing flow time is obtained from the diffusion process in Eq. (4.6), starting from the bare gauge fields. To further reduce the cut-off effects in the scale-setting procedure, an alternative quantity was proposed in Ref. [65]. Define Then the scale can be set by where W 0 is again a dimensionless constant that one can choose. On the lattice, the calculation of E(t) depends on a definition of G µν , for which a variety of choices are available. The most obvious is to associate it with the plaquette P µν ; an alternative is to define a four-plaquette clover, which has a greater degree of symmetry [62].
In the continuum, all definitions should become equivalent, and at finite lattice spacing the relative difference between the two decreases at large t. The shape of E(t) at very small t is dominated by ultraviolet effects, and so differs strongly between the two methods; this introduces further constraints into the choice of E 0 . Figure 5 shows E(t) and W(t), calculated both via the plaquette and the clover. As anticipated from [65], the discretisation effects are smaller in W(t) than E(t); this is visible in the splitting between plaquette and clover curves being smaller in the W(t) case. 6 In the continuum theory, B µ (x) are elements of the Sp(4) gauge group; however, it is possible that the finite precision of the computer could introduce some numerical artefact that would cause the integrated B µ (x) to leave the group. Since the integration is an initial value problem, any such artefact introduced would compound throughout the flow, giving potentially significant distortions at large flow time. For this reason we have introduced the re-symplectisation procedure described in Appendix C after each integration step. We find no appreciable difference between the flow with and without re-symplectisation.
We now proceed to set the values of E 0 and W 0 , such that t 0 and w 0 avoid both the regions of finite lattice spacing and finite volume artefacts. In order to obtain a single value for the lattice spacing corresponding to a particular value of β, and allow comparisons to be drawn with pure gauge theory, we must also be in the vicinity of the chiral limit. Findings in [65] indicate that when the fermions are light enough any mass dependence in w 0 should be small. Figure 6 shows the fermion-mass dependence in √ 8t 0 and w 0 at β = 6.9, choosing The discretisation effects are significantly smaller in w 0 . We see a relatively strong dependence on the fermion mass in both t 0 and w 0 ; this is contrary to expectations from studies of QCD with light quarks [65]. Presently we are studying this fermion-mass dependence. Results of this study will be detailed in future publications. It should be noted that if the behaviour highlighted here persists also in proximity of the chiral limit, extra care will be needed in the process of taking the continuum limit.

Topological charge history
As the lattice extent is finite in all directions, a given configuration will fall into one of a number of topological sectors, labelled by an integer (or, at finite a, near-integer) topological charge Q, which is expected to have a Gaussian distribution about zero. Since it is probabilistically unfavourable to change a discrete global observable using a small local update, Q can show very long autocorrelations; as the continuum limit (i.e. the limit of integer Q) is approached, Q can "freeze", ceasing to change at all.
It is necessary to check that Q is not frozen, and instead moves hergodically, for two reasons. Firstly, the exponential autocorrelation time of the Monte Carlo simulation as a whole scales as one of the longest autocorrelation time in the system (see e.g. [69]). Secondly, the values of physical observables depend on which topological sector a configuration is in [70]; sampling a single Q or an unrepresentative distribution of Qs will introduce an uncontrolled systematic error. It is therefore necessary to verify that Q not only moves sufficiently rapidly, but also displays the expected Gaussian histogram.
The topological charge Q is computed on the lattice as and where x runs over all lattice sites. For gauge configurations generated by Monte Carlo studies, this observable is dominated by ultraviolet fluctuations; therefore it is necessary to perform some sort of smoothing to extract the true value. The gradient flow (described in the previous subsection) is used for the purposes of this work. We have examined the topological charge history for all our existing ensembles, including both pure gauge and those with matter. In most cases, Q is found to move with no noticeable autocorrelation, and shows the expected Gaussian distribution centred on Q = 0. Samples of these histories are shown in Fig. 7. Some marginal deviation is visible for example in the second of the three series.

The spectrum of the Yang-Mills theory
In this Section, we focus our attention on the Sp(4) Yang-Mills theory. We start by reminding the reader about several technical as well as conceptual points related to the physics of glueballs and to the description of confinement in terms of effective string theory. We then summarise the specific methodology we adopt in the process of extracting physical information from the lattice data. We conclude this section by presenting our numerical results on the glueballs, and commenting on their general implications.

Of glueballs and strings
At zero temperature, Sp(4) Yang-Mills theory is expected to confine. The particle states are colour-singlet gluon bound states, referred to in the literature as glueballs, and labelled by their (integer) spin J and (positive or negative) parity P quantum numbers as J P . 7 To distinguish between states with the same J P assignment but different mass, in the superscript of the n-th excitation we add n asterisks ( * ). For instance, 2 + * * denotes the second excitation with J = 2 and P = +, while 2 + is the ground state in the same channel.
The calculation of the mass spectrum of glueballs requires a fully non-perturbative treatment of the strongly interacting dynamics. We follow the established procedure that extracts glueball masses from the Monte Carlo evaluation of two-point functions of gaugeinvariant operators O(x). The operators O(x) transform according to irreducible representations of the rotational group and either commute or anti-commute with the parity operator, hence having well-defined J P . Given O(x) defined at any spacetime point x = (t, x), we separate the space-like and time-like components x and t, 8 and define the zero-momentum operator O(t) as where the sum runs over all spatial points x at fixed t. The lowest-lying glueball mass in the J P channel is then given by Assuming only contributions from poles (an hypothesis that certainly holds at large N ), we can insert a complete set of single-glueball states |g n (J, P ) carrying the same quantum Glueballs are not the only interesting observables in Yang-Mills theory. In the presence of infinitely massive, static quarks, the spectrum contains also confining flux tubes. While flux tubes are exposed by the static probes, their physics is fully determined by the Yang-Mills dynamics and plays a crucial role in the study of confinement. Consider a static quark Q S and the corresponding antiquarkQ S , a distance ∆x apart. In a confining theory, the static quark-antiquark pair is bound by a linearly rising potential where the quantity σ (having dimension of a mass squared) is the (confining) string tension, and provides a measurement of the dynamically generated confinement scale. In Yang-Mills theory there is only one dynamically generated dimensionful quantity, hence the square root of the string tension also sets the scale of the glueball masses, besides providing a fundamental test of confinement. The semiclassical cartoon associated with linear confinement explains the latter as arising from the formation of a thin (flux) tube in which the conserved colour flux is being channeled. Over distances much bigger than the transverse size of the confining flux tube, the latter can be represented by a string of tension σ binding quark and antiquark together.
At zero temperature, a signature of confinement is the area law: where the Wilson loop W (∆x, ∆t) is defined as The contour integral of the gauge field A µ extends over a rectangular path R of sides ∆x along one spatial direction and ∆t in the temporal direction. In Eq. (5.6), g is the coupling, Tr indicates the trace and the exponential is path-ordered along R. The potential is then obtained as At finite temperature, the temporal direction of size τ is compactified on a circle, and the resulting thermal field theory has temperature T = 1/τ . The order parameter for confinement can be identified with the expectation value of the Polyakov loops: with C being the circle (of circumference τ ) at fixed spatial point x. 9 The expectation value of this quantity vanishes in the confined phase. This observable has the advantage that it makes transparent the fact that the transition is associated with the breaking of the centre symmetry of the gauge group. In this respect, the Sp(2N ) theories play a useful 9 To avoid confusion with the average over spatial directions Φ in Eq. (3.8), when referring to Polyakov loops we explicitly indicate the x-dependence, being it understood that the average over the other space-like coordinates is taken. For example, we will later indicate as Φ(z) the average of Φ( x) over two space-like directions x and y.
complementary rôle with respect to SU (N ), the centre of the former being Z 2 for every N , as opposed to the Z N centre of the latter. In this set-up, the propagation of a pair of static conjugated quarks is represented by two oppositely-oriented Polyakov loops and their correlator Φ † ( 0)Φ( x) probes strings attached to two static lines at 0 and x. In the language of string theory, the confining string stretching between static sources is an open string.
Yet, in Euclidean space we can reinterpret the zero-th direction as a compact spatial dimension and (for instance) the third direction as Euclidean time. From this point of view, the string is not attached to any static source but closes upon itself. For this reason, it can be also interpreted as a closed string. Choosing x = (0, 0, z) and inserting a complete set of eigenstates |l n of the transfer matrix (the time-translation operator) in the third direction yields with c l n the overlap between the state Φ( 0)|0 and the n-th eigenstate of the Hamiltonian along z and E n the corresponding energy eigenvalue. In this case the Polyakov loop correlator probes (closed) string states wrapping along the compact direction, created at 0 and annihilated at x.
The fact that the same correlator can be interpreted in terms of either propagating closed or open strings expresses the open-closed string duality, a key observation that has profound physical implications. Among them, the most direct and practically relevant for our study is the fact that the string tension can be extracted in the closed string channel from correlators of Polyakov loops. This is related to the fact that the topology of the world-sheet swept by the string is cylindrical. 10 If we instead consider the operator obtained by averaging Φ( x) along two dimensions where the sum runs over the two spatial coordinates in the directions orthogonal to z, for the correlator we obtain Φ † (0)Φ(z) = n c l n e −m l n z , and open-closed string duality implies that The state corresponding to m l 0 (where the subscript l stands for loop) can be interpreted as the ground state mass of a torelon, a stringy (flux tube) state that wraps around the compact direction. In general, torelon states can be labelled by their length τ , the absolute value of their angular and longitudinal momenta J and q, their (transverse) parity P t in a plane transverse to their symmetry axis, and their longitudinal parity P t along the wrapping direction. As for glueballs, the gauge group being pseudo-real, charge conjugation is always positive, and furthermore we are interested only in torelons with both transverse momenta equal to zero and both positive parities.
The quantum fluctuations around the classical world-sheet solution corresponding to the area law in Eq. (5.5) generate a spectrum of modes for the flux tube that can be computed using an effective string theory description. The relevant degrees of freedom are identified as the D − 2 Goldstone bosons living in the 2-dimensional world-sheet of the flux tube that breaks the D-dimensional Poincaré group ISO(D) according to: If the theory has a mass gap, as is the case for Yang-Mills theory, and no other degrees of freedom are present on the world-sheet, the most general effective action S eff [X] describing the dynamics can be expressed as an expansion in derivatives of X µ = {ξ a , X i } with respect to the world-sheet parameters (ξ 0 , ξ 1 ), where a, b, c = 0, 1, while i, j = 2, . . . , D − 1 and summation over repeated indices is understood. This action can be naturally recast as an expansion in powers of 1/(στ l), as a low momentum expansion around an infinitely long string. This expansion is meaningful as long as lτ 1/σ. In turn, a flux tube is string-like provided the long-string expansion is valid (l τ ), and hence provided l 1/ √ σ.
In lattice calculations, spacetime is a box of finite extent. When taking limits such as the r.h.s. of Eq. (5.12), the extraction of m l 0 is contaminated by short-distance contributions that can be non-negligible. Let us specialise to the closed string channel, for which the world sheet is a cylinder with time direction collinear to its axis. The mass can be systematically approximated as The dimensionless quantities d l k -related to the C k in Eq. (5.14)-enclose all the shortdistance effects that have been integrated out in this effective description. They are not completely independent: the request that S eff inherits the same space-time symmetry as the underlying Yang-Mills field theory imposes constraints (open/closed string duality being an example), and only few of the d l k -equivalently, the C k -are left as free parameters. This striking effective string-theory feature determines the universal nature of the quantum corrections to the area law, to which we devote the remainder of this subsection.
In one of the earliest works on the subject [71], it was shown that the first non-null coefficient appears at O 1 στ 2 , and has a universal nature: 11 The correction is known as the Lüscher term [73], and can be interpreted as the Casimir energy of the string massless modes. The coefficient c captures the central charge of the effective string theory, hence encoding information on its nature. For example, the purely bosonic string theory in non-critical dimension has c = 1, and there is evidence that it can be used to describe SU (N ) Yang-Mills gauge theory (see e.g. [74][75][76]). The next-to-leading-order correction coming from the effective field theory has been studied in [77,78] (see also [79][80][81][82] for a related derivation in the Polchinski-Strominger approach). It arises at O(1/(στ 2 ) 2 ). Duality arguments (in particular the annulus-annulus duality [78]) can be used to prove the universality of this correction. At this order, the effective mass appearing in Polyakov loop correlators is given by Interestingly, the same result can be obtained by expanding in powers of 1 στ 2 the light-cone spectrum obtained from the Nambu-Goto action in non-critical dimensions (D = 26): Going for a moment beyond the specific purposes of our paper, we recall that establishing the nature of the string that forms between a static quark-antiquark pair in a confining gauge theory is a very interesting programme in itself, as it can shed some crucial light on the nature and on the mechanism of colour confinement. The recent revival in its interest has resulted in new fundamental advances, among which the key observation that the constraints mentioned above are in fact particular cases of a more general viewpoint allowing to extend universality considerations to higher orders. Because the effective action is still Poincaré invariant (despite spontaneous symmetry breaking), the difference between the number of derivatives and the number of fields (called weight) is an invariant. The expansion of S eff can be organised according to the weight, and (non-linearly realised) Poincaré invariance imposed upon each of the terms. The result is the emergence of recurrence relations among d l k terms of the same weight. The unique weight-0 invariant action is precisely the Nambu-Goto action, with the leading correction appearing at weight-4. This explains, for example, why up to order 1/τ 5 and in D = 4, the predictions of the light-cone spectrum for the ground state energy are fully universal. For a more detailed analysis, we refer the reader to [83,84] and references therein.
Coming back to our lattice calculation, we build on the results available for SU (N ) and assume that the nature of the confining string is reproduced also in Sp(4). We use (rather than derive) the expressions for the bosonic string in order to extract numerical values for the string tension σ. Subject to the validity of such working assumption this enables us to remove large finite-size corrections from the extraction of the string tension.

A variational approach to the mass spectrum
In this subsection, we outline the lattice methods used to extract both σ and the spectrum of glueball masses [85][86][87]. It is important to note from the outset that because of centre symmetry, glueball and flux-tube states do not mix in the confined phase. Glueballs and flux tubes are sourced by products of link matrices around contractible and non-contractible loops C, respectively, with their geometrical symmetry properties determining J P for the former and size L, momentum q, angular momentum J and parities P l and P t for the latter. We limit our study of flux tubes to the state with J = q = 0 and P t = P l = +. We refer to [76] for a detailed analysis of creation of excited flux tubes in various channels.
The isotropic lattice breaks the continuum rotational group to the octahedral group (the 24-elements group of the symmetries of the cube). At finite lattice spacing, glueball states are classified by the conventional names R = A 1 , A 2 , T 1 , T 2 , E of the irreducible representations of the octahedral group, so that we label the glueballs as R P , instead of using the continuum J P .
The irreducible representations of the octahedral group can be decomposed into irreducible representations of the continuum rotational group. Since the octahedral group is finite, different continuum spins are associated with a given octahedral irreducible representation. For instance, the continuum J = 0 spectrum is found in the A 1 representation, which also contains (among others) J = 4 states. While on physical grounds one can assume that the lightest A 1 state corresponds to a J = 0 glueball (at least when P = +), distinguishing different continuum channels in the excited spectrum measured on the lattice is not an easy task. Guidance is provided by the degeneracies that are expected in different octahedral representations, where different polarisations of the same state can appear. This is for instance the case for the continuum J = 2 states, two polarisations of which are to be found in the octahedral E representation, and the other three in T 1 . Hence, close to the continuum limit, states that are degenerate in the E and T 1 channel can be interpreted as would-be continuum states with J = 2.
Given a lattice path C with given shape and size, located at reference coordinates (t, x) on the lattice, that transforms in an irreducible representation R of the octahedral group and is an eigenstate of parity, the lowest-lying mass in the R P channel can be extracted from the asymptotic behaviour of the correlator where m i is the mass of state |i and O C (t) is the zero momentum operator: 12 The decay of each exponential appearing in the spectral decomposition is controlled by the squared normalised amplitudes In practice, since the statistical noise is expected to be constant with t, the signal-tonoise ratio decays exponentially, eventually defying attempts to isolate the ground state. Hence, for a generic choice of C, the mass that is extracted at large but finite t suffers from contamination from excited states. We notice that, as a consequence of unitarity, this results in an overestimate of the mass.
To improve accuracy in the extraction of m i , in principle one should choose the operators O to maximise the overlap of O C (t)|l with the desired state |l . While such an operator is not known a priori, we can operationally construct a good approximation by performing a variational calculation involving loops C of various shapes and sizes. The size λ of the loop C must be chosen appropriately: in order for it to capture the infrared physics, it should have a size of the order of the confinement scale. This means that in practice the size of C in lattice units must grow as the lattice spacing goes to zero. Over the years, various methods to circumvent these potential problems have been suggested. In this work, we shall use a variational calculation involving an operator basis constructed with a combination of smearing and blocking operations. We briefly review the method used, and we refer to [40] for more details, before presenting our results in Sec. 5.3.
Given a set of N operators O i , defined as in Eq. (5.20) for paths C i of different shape and sizes labelled by i, we compute the N × N , normalised correlation matrix Assuming maximal rank, C ij can be diagonalised, and we callC ii the N resulting functions of t. The special linear combination i α i O i (t), corresponding to the maximal eigenvalue, has the maximal overlap with the ground state in the given symmetry channel. Assuming its mass is the only one present in the given channel, we obtain it by fitting the data with the functionC where |c i | 2 and m i are the fit parameters, and where the appearance of the cosh is due to the inclusion of backward propagating particles through the periodic boundary. In general the data is contaminated by contributions from states with higher mass. Hence we must restrict the fit to the range for which we see the appearance of a plateau in the quantity In order to include operators which extend beyond the ultraviolet scale, following [40], we subjected the lattice links to a combination of smearing and blocking transformations.
These are iterative procedures similar to block transformations in statistical mechanics, except that we restrict them to space-like links. The operators O i defined as in Eq. (5.20) for paths C i are evaluated using the output links from each iterative smearing and blocking step. After S iterations, one has a collection of S ×N such operators, where N is the number of basis paths in the given channel. We chose to build the operators O i by starting from a large set of basic lattice paths. In this set, we include all the closed paths with length λ up to ten in units of the lattice spacing, appropriately symmetrised to transform according to an irreducible representation of the octahedral group and to have definite parity.
We implement the process of smearing along the directions orthogonal to the direction of propagation, by starting from the link U s=0 where andî refer to the unit-length displacements along the lattice directions j and i, respectively, while the positive parameter p a controls how much smearing is taking place at each step. The smeared matricesŨ s>0 i (x) are not elements of the gauge group. We project those matrices to the target group by finding the Sp(4) matrix U s>0 . This is done in two steps: a crude projection is operated by using one of the re-symplectization algorithms presented in Appendix C, and afterwards a number of cooling steps [88] (15 in our case) is performed on the link.
Similarly, blocking is implemented by replacing simple links U b=0 i (x) ≡ U i (x) with super-linksŨ b>0 i (x) that join lattice sites 2 b spacings apart, where b is the number of blocking iterations, as described bỹ Again, each such step yields a matrixŨ b+1 i (x) that does not belong to the Sp(4) group, and hence must be projected onto U b+1 i (x) within the group in same way as for the smearing. In practical terms, when performing numerical lattice studies blocking allows to reach the physical size of the glueball in fewer steps, while at the physical scale smearing provides a better resolution. Due to the fact that the identification of the physical scale is a dynamic problem, an iterative combination of n smearing steps (generally, n = 1, 2) with a blocking step generally proves to be an efficient strategy [89].

Lattice results
In this Subsection, we report the results of our numerical analyses of the glueball spectrum and the string tension of the pure Sp (4)   performed on fully isotropic lattices of various sizes and lattice spacings. To investigate the finite size effects, we first consider the coarsest lattice with β = 7.7. Based on the estimate of the critical coupling of the bulk phase transition [1], the choice of this value should provide a prudent yet reasonable compromise between the practical necessity of performing a detailed calculation at a lattice coupling at which the physically large volumes can be reached on a moderately coarse grid and the physically paramount request that the lattice gauge theory be in the confining phase connected to the continuum theory as a → 0. Indeed, we have found evidence in our calculations that at β = 7.7 the lattice theory is in the physically relevant confined phase. We started with this β = 7.7 value, and increased the lattice size, starting from L = 10a, until we found the best economical choice at which the (exponentially suppressed) finitesize effects became much smaller than the statistical errors. Assuming scaling towards the continuum limit, this analysis provides a lower bound for the physical volume of the system such that finite-size effects are negligible with respect to the statistical errors, and hence ensures that the calculations cannot be distinguished from infinite volume ones.
We repeated the same set of measurements on progressively finer lattices (larger β), always making sure that the physical volume were large enough for the calculation to be considered at infinite volume for all practical purposes, and extrapolated the glueball spectrum to the continuum limit.
The whole procedure is illustrated more quantitatively in the following two sub-subsections. Our parameter choices for the continuum extrapolation are reported in the first two columns of Table 2. For each lattice setup 10000 configurations were generated, and a binned and bootstrapped analysis of errors was carried out to take care of temporal autocorrelations. Operators blocked to the level N b ≤ L and with 15 cooling steps were used, resulting in a variational basis of ∼ 200 operators.

The string tension
To extract the string tension from measurements of masses of closed flux tubes, we turn to effective string theory. A finite overlap with flux-tube states can be obtained with lattice operators defined on non-contractible loops, as described earlier. We produce two measurements of the string tension, that we denote as σ t and σ s . The former is obtained as follows. We first consider loops that wrap the time-like direction as in Eq. (5.8) and average them along two space-like dimensions as in Eq. (5.10). We then compute the correlators as in Eq. (5.11), with an additional statistics-enhancing average over interchanges of (x, y, z), to  Table 3. Masses obtained from Polyakov loop correlators winding in the time direction (m t ) and in a spatial direction (m s ), together with the corresponding string tensions σ t and σ s extracted from the Nambu-Goto (NG) prediction for the ground state energy of the flux tube of length L at β = 7.7. In the last column, we report the result of the weighted average in Eq. (5.29).
extract the lowest mass m l 0 (which in this case we refer to as m t ). Finally, we determine the string tension as in Eq. A similar procedure is performed for obtaining σ s , which is extracted from the mass m s associated with correlators in time of Polyakov loops winding one of the spatial directionsexcept that in this case there is no average over interchange of equivalent directions. Because the lattice used is isotropic, we expect σ t to be compatible with σ s , since the latter could be obtained from the former by interchanging the roles of the time direction and one of the space directions used for defining the correlation functions associated with m t . Notice that because of the averaging over interchanges of spatial directions, the statistical error on σ t is reduced in respect to σ s .
In order to carefully assess finite size effects, we show the results of the analysis for β = 7.7 in Table 3. We perform a best fit analysis of the data for m t and m s by using Eqs. (5.16)-(5.18). We start from the largest flux length L = 24a and gradually add in the fit lower-length values, until the value of the χ 2 /d.o.f becomes larger than a fixed threshold that we conventionally set at 1.2. We find the following best fit results for σ t , obtained with the largest possible range for which χ 2 /d.o.f < 1.2:   For both m t and m s , our requirement for the acceptability of the fit is verified down to L = 12a for all the three proposed functional forms, except for the leading-order ansatz in the case of σ t , which requires L = 14a. We regard this last case as a warning that at L = 12a the effective string description at leading order might break, although the results for the NLO and NG descriptions give us confidence that higher orders cure this problem. The fits provide very good indications that the description we adopted is robust for L/a ≥ 14. Indeed, all the reported fitted values are compatible, regardless of the ansatz used and of the direction of correlation of the Polyakov loops from which we extract the relevant mass. Conversely, when we try to extend the fit down to L = 10a, we typically find a significantly larger χ 2 /d.o.f, of order 3-10, indicating that the effective string description cannot be trusted between L = 10a and L = 12a. The only exception is the NG description of σ s , for which we get χ 2 /d.o.f 1.75. While it would be tempting to interpret this result as evidence that the NG ansatz provides a better description of the data, in the absence of confirmation of this hypothesis in the σ t case (for which an extension down to L = 10a leads to χ 2 /d.o.f 8. 22) and given also the scope of our calculation, we prefer to take a cautious attitude towards our results and assume that a safe lower bound for all effective string models analysed to work (and to produce compatible results) is L = 14a.
Given all of these considerations, and taking into account all our estimates of a 2 σ, a safe infinite-volume value for the latter quantity that encompasses the spread of the fits is σa 2 = 0.5179 (50), which translates into √ σa = 0.2276 (11). Using this result, in physical units one finds that L √ σ = 2.731(13) for L = 12a and L √ σ = 3.186(15) for L = 14a. The fact that effective string description works remarkably well for Polyakov loop correlator masses down to at least L = 14a is consistent with the picture of confinement through the formation of thin flux tubes. It is of practical relevance for numerical studies to assess how well the infinite-volume value of the string tension is represented by the result extracted inverting Eqns. (5.16-5.18) at a single finite size L, and how this would be affected by varying L. To provide information about this, we report the results of our procedure in Table 3 and in Fig. 8. As we see from the figure, the value of σa 2 is constant for a wide range of L. This holds for the LO, NLO and the NG extractions, with the corresponding values being always compatible within errors. Based on these observations, we use the NG approximation to extract our best estimate. We detect finite-size effects for the smallest lattice volumes L = 10a. Though we do not present the detailed results, we also detect a discrepancy at L = 24a between the space-like and time-like string tensions. This discrepancy may be a consequence of the systematic error coming from the difficulty in extracting the asymptotical behaviour of the correlator for very large masses.
Our final estimate for the value of √ σ as a function of L is obtained from the first two columns of Table 3 by computing the weighted average where the error is given by The resulting values are reported in the last column. From the behaviour of √ σ we conclude that finite size effects are certainly smaller than the statistical error for L ≥ 14a, and we take as a final estimate of √ σ at this coupling the value at L = 16a. We also note that compatible results are obtained for L = 12a. Assuming scaling towards the continuum, from our finite-size study we obtain firm evidence that all lattices for which L √ σ ≥ 3 provide an estimate of the string tension that is compatible within the statistical errors with the infinite-volume value. Hence, we conclude that finite-volume effects are negligible once L √ σ > ∼ 3. In particular, we have verified that the condition L √ σ > ∼ 3 is safely fulfilled when carrying out calculations on lattice ensembles with larger β, starting from the finite-size analysis at β = 7.7. Table 2 reports the lattice parameters of the calculations we have used to extrapolate to the continuum limit and the corresponding results for √ σ.

The mass spectrum of glueballs
As for the string tension, we began our analysis of glueball masses with a study of finitesize effects for lattice coupling β = 7.7. We aimed at estimating finite-volume effects as a function of the lattice size L, and bound L such that the systematic error due to the finite size be negligible with respect to statistical error on the measurement of the masses.  (15)   Our results for the mass spectrum at β = 7.7 for various volumes are reported in the rows of Table 4. While this particular choice of β is suitable for finite-size studies, as it allows us to reach large lattices in physical units with a relatively small computational effort, the coarseness of the lattice spacing pushes most of the masses above the lattice cut-off, making their extraction numerically challenging. For this reason, we observe a systematic effect related to the isolation of the ground state on all irreducible representations other than the lowest-lying one. While in our tables we quote only the statistical error, for higher excitations the systematic error coming from the ground state isolation in a given channel is expected to have a comparable size.
Another systematic effect that affects our calculation of the glueball masses is contamination of the spectrum by multi-glueball scattering states and torelon states (the latter being lattice artefacts associated with state propagation of pairs of oppositely directed Polyakov loops). Separating the physical spectrum from those unwanted states requires a more demanding calculation 13 that goes beyond the scope of this first exploration of the glueball masses in Sp (4). What results is a hard-to-control error related to mixing with spurious states, which is strongly dependent on the volume and manifests itself in occasional sudden jumps and irregularities of the extracted masses. Indeed this is visible for some of the most massive states we report in Table 4. In view of all these considerations, and to focus the discussion to the main purposes of this paper, we limit the analysis of finite-size effects at β = 7.7 to the ground state (found in the A1 + channel) and to the would-be continuum 2 + glueball (expected to appear in the E + and in the T 1 + channels).
Starting with the A1 + state, we fit the behaviour of the extracted mass to the size dependency parameterised as follows:  and we find its infinite volume limit to be am = 0.746 (6) , This value for am is compatible within errors with all the measured masses in the A1 + channel for L/a ≥ 12. Hence, as in the case of the string tension discussed earlier, the systematic error due to neglecting the finite size of the lattice is found to be comfortably less than the statistical error as long as L √ σ ≥ 3. Thus, having taken care that this bound be satisfied at all β values simulated for the continuum extrapolation, in the extraction of the continuum limit we shall use the value measured at one lattice size. It is also interesting to look at the masses in the E + and in the T 1 + channels, as the lowest-lying states in these two channels are expected to be degenerate in the continuum limit, since they correspond to different polarisations of the continuum 2 + state. From Table 4 we infer that this degeneracy is satisfied for lattice sizes other than the smallest L/a = 10 (for which large finite-size effects are expected to arise differently in the two channels) and the largest L/a = 24, where some uncontrolled systematics (possibly due to mixing with spurious states) creates visible anomalies in this measurement. We take the agreement in the intermediate region as a good indication that β = 7.7 is sufficiently close to the continuum limit to justify its inclusion in the continuum extrapolation. Going towards the continuum, we have measured glueball masses in the various lattice channels for the couplings and volumes reported in the first two columns of Table 2. Continuum extrapolations for the ratio m G / √ σ are obtained from the expectation that the corrections due to discretisation are linear in σa 2 . The results are reported in Table 5 and represented in Figs. 9 and 10. Some excited states are shown in Fig. 11. As expected, in the continuum limit the states T ± 1 and E ± are degenerate in pairs. The channels T ± 2 show strong fluctuations and discretisation effects that are caused by the difficulty in extracting their masses in lattice units. The values obtained for the 0 + and 2 + states are at the same order of magnitude as those obtained in SU (N ) theories [90].

Epilogue
In this Section, we have reported on what is (to the best of our knowledge) the first controlled calculation in the continuum limit of the glueball spectrum and of the string tension for the Sp(4) pure Yang-Mills theory. The main purpose of our study is to gain some understanding of the glue dynamics in this theory, and progressively aim at providing an interpretation of the results emerging from the theory with dynamical quarks. Nevertheless, the outcomes of our investigation in the pure gauge case are interesting in their own right: for instance, they provide a first step towards a systematic calculation of the pure Yang-Mills spectrum in the large N limit of Sp(2N ) gauge theories and give us an opportunity to contrast the non-perturbative dynamics of Sp(2N ) with that of SU (N ) gauge theories. Albeit expected, the first remarkable outcome of our calculation is that the pure gauge quantities behave not dissimilarly from those of SU (N ). The mass associated to the decay of Polyakov-loop correlators follows the properties expected from a confining flux tube, hence supporting the confining nature of the pure gauge dynamics in Sp(4). The lowestlying state in the glueball spectrum is the 0 + glueball. The ratio m 0 + / √ σ 3.55 is not far from the large-N value from SU (N ) groups, m 0 ++ / √ σ 3.3 [39], and is indeed compatible within errors with the SU (3) value of m 0 ++ / √ σ = 3.55(7) [40], with the value for SU (4) of m 0 ++ / √ σ = 3.36(6) being appreciably different, albeit close.
This behaviour is consistent with the above ratio being a mildly decreasing function of the number of generators-while SU (4) has 15 generators, SU (3) has 8 and Sp(4) has 10and deserves further enquiry by performing numerical studies of Sp(2N ) gauge theories at larger N . An interesting observation has been put forward in [44], by suggesting that where C 2 (A) and C 2 (F ) are the quadratic Casimir of the adjoint and of the fundamental representation, respectively, and η is a universal constant, in the sense that it depends on the dimensionality of the spacetime, but not on the gauge group. Noting that for Sp(2N ) we find that our determination of the proportionality constant η = 5.27 (15) , is compatible with the value η = 5.41(12) extracted from SU (N ) in 3+1 dimensions [44]. The rest of the glueball spectrum also follows a pattern that is broadly similar to that of SU (N ). Another interesting quantity in the glueball sector is the ratio m 2 ++ /m 0 ++ . Using universality arguments, it has been argued in [43] that for confining theories where the dynamics does not yield large anomalous dimensions, as in pure Yang-Mills, one should find m 2 ++ /m 0 ++ = √ 2. Our numerical results give m 2 ++ /m 0 ++ = 1.425 (32), a value that is fully compatible with the conjecture of [43].
Besides being relevant for models of electroweak symmetry breaking based on a Pseudo-Nambu-Goldstone interpretation of the Higgs field, the investigation of which is the central Table 6. Interpolating operators used for the measurement of the properties of flavoured mesons (i = j) in the pseudo-scalar (PS), vector (V) and axial-vector (AV) cases, associated with the π, ρ and a 1 mesons, respectively. We also report the Lorentz-group quantum numbers J P . The summations over colour and spinor indices are understood.
leitmotif of this paper, studies of Sp(2N ) pure gauge theories provide new relevant information on universal aspects of Yang-Mills dynamics. We shall develop this latter line of research in future numerical investigations.

Of quenched mesons: masses and decay constants
In this Section, we present our results for the masses and decay constants of the lightest flavoured mesons in the quenched approximation. Our main purpose is to illustrate the process that we envision we will carry out once simulations with dynamical quarks are available. Although we are aware of the fact that the quenched results may not capture in full the features of the theory, we still expect it to provide some useful information about its qualitative features. Experience on QCD with light quarks suggests that several quantities are well captured by the quenched approximation, although we already cautioned the reader about the fact that such considerations may not extend to these dynamical theories. The EFT in Sec. 2.1, within the limitations discussed therein, describes the continuum limit of the dynamical simulations, not the quenched ones. In principle, one could make more sense of the comparison by adopting the approach of quenched chiral perturbation theory [91,92] or of partially-quenched chiral perturbation theory [93][94][95], but for present purposes our strategy will suffice, though we invite the reader to use caution, in particular for quantities such as the g ρππ coupling, that are certainly affected by the quenching procedure.

Observables
As discussed in Sec. 2.1, the observables of most direct phenomenological relevance, and at the same time the most directly accessible to lattice calculations, are the masses, m M , and the decay constants, f M , of pseudo-scalar, vector and axial-vector mesons. This numerical study focuses on flavoured particles. The interpolating operators and quantum numbers are summarised in Table 6.
We define the ensemble average of two-point meson correlators in the Euclidean space as where O M denotes any of the meson interpolating operators in Table 6. In the limit in which the Euclidean time t is large, and for vanishing three-momentum p, the correlation function is dominated by the ground state. Its exponential decay is controlled by the meson mass m M , and can be approximated as where T and L are the temporal and the spatial extent of the lattice, respectively. In the meson states |M , we define M = M i T i , with T i the generators of the group. Using this convention, the mesonic matrix elements are parameterised in terms of decay constants f M and masses m M as 14 where µ is the polarisation vector obeying µ p µ = 0 and µ µ = +1. Using Eq. (6.2) and Eq. (6.3), for vector and axial-vector mesons we can rewrite the correlation functions To calculate the decay constant of the pseudo-scalar meson, we additionally consider the following two-point correlation function The pion mass m π and matrix element 0|O P S |P S are obtained from the pion correlator, Meson decay constants computed on the lattice are matched to the continuum. In this work, we perform one-loop matching in lattice perturbation theory. Because we are using Wilson fermions, the axial and vector currents are not conserved in the lattice theory, and hence receive (finite) renormalization, that we write as −7.75 −3.00 Table 7. Results of one-loop integrals for the matching coefficients in Eq. (6.8) at the choice of the Wilson parameter r = 1, and taken from [96] .
The pion decay constant f π is renormalised by Z A , as the axial current is used to define f π in Eq. (6.3). In the continuum limit, the renormalization constants are expected to be unity. The one-loop matching coefficients taken from [96] are given by where g is the coupling constant and the eigenvalue of the quadratic Casimir operator is C(F ) = 5/4 for Sp(4). The coefficients ∆ I , relating the lattice computation with the continuum MS regularisation scheme, result from one-loop integrals performed numerically, and are summarised in Table 7. We notice that the coefficients reported here have been obtained by restricting the integrals within the first Brillouin zone. We verified explicitly that the errors in numerical evaluation for these integrals are 2% or less, and that there are no discernible finite-volume effects. Therefore we neglect the uncertainty on the ∆ i coefficients in the rest of this paper. The coefficient ∆ Σ 1 is taken from the wavefunction renormalization of the external fermion lines, without taking into account the power-divergence contribution, while the coefficients ∆ γµ and ∆ γ 5 γµ are extracted from the one-loop computations of the vertex functions.
Wilson fermions receive quite large renormalisation and thus the perturbative expansion with the bare coupling is reliable only at very large values of the lattice coupling β = 4N/g 2 . As in the continuum case, in lattice perturbation theory, the appropriate expansion parameter is the renormalised couplingḡ rather than the bare counterpart g. Following [97], we use the simple tadpole improved coupling defined bȳ where P is the plaquette operator.

Numerical results and EFT
As a first exploratory step towards understanding the qualitative features of mesons of the Sp(4) theory, as well as to test the low-energy EFT and illustrate its use, we calculate the mesonic observables described in the previous section in the quenched approximation at  Table 9. Masses and decay constants (squared) of pseudo-scalar (π), vector (ρ), and axial-vector (a 1 ) mesons, as obtained from the quenched calculations described in Sec. 6.2, for β = 7.62 on lattice of size 48 × 24 3 . g), and the gradient flow scales t 0 and w 0 . The numerical results are summarised in Table 8. The gradient flow scales are determined by setting E 0 = W 0 = 0.35. As discussed in Section 4, we use two definitions of G µν denoted by the superscript p for a plaquette and c for a (four-plaquette) clover, with the difference of the two measuring the size of finite lattice spacing artefacts. We find that, compared to t 0 , the difference between the clover and the plaquette regularisations for w 0 is significantly smaller, and thus we choose to use w c 0 to convert lattice units to physical ones.
The mesonic two-point correlation functions in Eq. (6.1) are measured using stochastic wall sources [98] at various quark masses. While at sufficiently large time we use the asymptotic expressions of two-point correlation functions in Eqs. (6.2), (6.4) and (6.5) to extract the meson masses and decay constants, we perform multi-exponential fits when the time extent is not large enough to reach the asymptotic region, according to where E 0 = m M < E 1 < E 2 < · · · . We find that the two-exponential fits are good enough  fit range (m 0 ) to describe the numerical data in most cases. The numerical results (not renormalised) are summarised in Tables 9 and 10. The statistical errors are estimated by using a standard bootstrapping technique (with about 250 bootstrap samples) and a correlated fit with χ 2 minimisation.
To analyse the numerical results using the continuum EFT developed in Sec. 2.1, we reinstate the correct dimensionality by expressing all mass and decay constants in units of w 0 = w c 0 . The decay constants are further renormalised using perturbative one-loop matching.
Let us first restrict our attention to the pseudo-scalar mesons and check the GMOR  relation. In the upper panel of Fig. 12 we plot m 2 P S f 2 P S , against the quark mass, where the latter is defined by the difference between the bare quark mass and the critical mass at which the mass of pseudo-scalar mesons vanishes. Even for the lightest masses m 0 available we do not find the expected leading-order linear behaviour. Therefore, we fit the data to the NLO results including the v 5 term in Eq. (2.16). The solid and dashed lines in Fig. 12 are the fit results for β = 8.0 and 7.62, respectively, where the fitting ranges and the resulting fit parameters are found in Table 11. As the quark mass is scheme-dependent, to compare the lattice results obtained with different coupling we would have to properly renormalise the quark mass. Instead of doing so, we consider the GMOR relation with respect to the pseudo-scalar meson mass, m P S , which is a scheme independent (physical) quantity, and show the result in the bottom panel of Fig. 12. As shown in the fit results in Tab. 12, we find that the parameters of the NLO GMOR relation between two lattices are statistically consistent, implying that finite lattice artefacts affect this observable only in a negligible way.
One of the interesting quantities in the EFT is f 2 0 , the sum of the squared decay constants of pseudo-scalar, vector, and axial-vector mesons, which is independent of the mass m 0 if one adopts Eq. (2.16). As shown in Fig. 13, the numerical results are in good agreement with such expectation. In the figure, blue circles and orange crosses represent the choices β = 7.62 and for β = 8.0, respectively. Notice that we derived this relation from our NLO EFT at the tree level. In principle, one has to consider the contribution from the one-loop corrections, including chiral logarithms. We performed a constant fit to the results and obtained f 2 0 (β = 8.0) = 0.0926 (17) and f 2 0 (β = 7.62) = 0.0944(18) denoted by blue and orange bands, respectively. As for the case of the GMOR relation, we find that lattice spacing artefacts in f 2 0 are negligible, in the sense that no appreciable difference is measured.
We finally use the NLO EFT relations to construct a global fit to the meson masses and decay constants. The results are illustrated in Figs. 14 and 15. We perform an uncorrelated fit to the data, restricted to the eight and six lightest masses m 0 for β = 7.62 and 8.0, respectively. The fitting range is shown as the shaded region in the figures. There are two main technical difficulties in this fit procedure. First of all, the parameter space is too large to determine the actual global minimum. In the NLO EFT we have thirteen fit parameters (including the critical bare mass m * 0 ). The standard χ 2 minimisation is not stable, and it typically yields two qualitatively very different results, one fairly linear and one exhibiting highly nonlinear behaviours with respect to the quark mass. We take the former as our best fit as the stability of the fits is better than the latter when we vary the fitting range. Secondly, statistical uncertainties vary widely for different types of mesons, and as a result the pseudo-scalar mesons tightly constrain the fit results. Furthermore, the numerical data suffer from several systematics such as quenching and discretising effects.
Undeterred by all these limitations and difficulties, since the purpose of this explorative study with quenched calculations is to show how our EFT works, we attempted to perform We explicitly checked that all of these fits satisfy the unitary constraints. Notice that the values of χ 2 /d.o.f are very reasonable, as also shown by the figures. Yet, the fitting procedure is very rough, the comparison between quenched data and continuum EFT is not rigorous, and hence we do not include uncertainties on these values of the couplings as a way to stress the fact that they should be used just for illustrative purposes. With this illustrative values of the parameters, one finds that g 2 ρππ /(48 π) ∼ 0.76 (for β = 7.62), and g 2 ρππ /(48 π) ∼ 1.0 (for β = 8.0). Such large values of g ρππ are affected by the uncontrolled systematics originating from quenching effects, and hence should not be used beyond the illustrative purposes of this exercise.
In future studies with dynamical calculations, a dedicated examination of the statistics and the systematics of the EFT fits will be required to determine the corresponding lowenergy constants in a meaningful way. Furthermore, we anticipate that it will be more involved to apply the continuum EFT result to the dynamical simulation, because the scale-setting procedure becomes more subtle. For instance, we observed in Sec. 4.1 that the scale w 0 used in our analysis changes visibly as the quark mass is varied in the dynamical case. Finally, it would be interesting to compare the results for masses and decay constants to the calculation presented in [38], and a possible extension that includes the dependence on the quark mass m 0 .
In order to investigate finite lattice spacing artefacts, as we observed earlier it is more effective to plot the mesonic observables with respect to the pseudo-scalar meson mass, rather than the bare quark mass. The results are shown in Fig. 16. In this case, only the masses of vector and axial-vector mesons are plotted. As we already learned from the GMOR relation, discretisation effects for the pseudo-scalar mesons are negligible in their spectroscopy. Similar conclusion can be drawn for the axial-vector mesons, given the current uncertainties. On the other hands, the masses and decay constants of vector mesons are affected significantly by lattice artefacts.

Towards dynamical fermions
The study of strongly-coupled gauge theories on discretised Euclidean space-time assumes the existence of a proper continuum limit as the lattice spacing a decreases, so that the field-theory dynamics is recovered. In order to avoid uncontrollable systematic effects in the continuum extrapolation, one must explore the lattice parameter space to identify any singularities or bulk transitions before carrying out detailed numerical studies of physical observables. In particular, a bulk phase where lattice discretisation effects dominate the behaviour of the system is expected to be present at strong coupling, with the interesting physical region separated by this bulk region by a first-order phase transition, or a very sharp cross over.
The identification of the associated (pseudo-)critical coupling is also strongly desired for practical purposes, because with finite numerical resources one cannot reduce the lattice spacing to arbitrarily small values, while at the same time using large enough volumes. For the pure Sp(4) lattice theory with the standard plaquette action the numerical study in [1] shows the absence of bona-fide bulk phase transitions. To the best of our knowledge, no such a study for dynamical simulation with N f = 2 Wilson fermions exists in the literature.
A possible choice of order parameter associated with the lattice bulk transition is the expectation value of the plaquette. To have a rough mapping of the transition, we first scan the parameter space over the range of β = [6.0, 7.0] using a 4 4 lattice, and show the results in Fig. 17. For each lattice coupling, we calculate the average plaquette values and vary the bare mass in steps of 0.1, over the range 0.0 ≤ −am 0 ≤ 1.4. Near the region in which the change of the plaquette value is large as a function of am 0 , we add twice more data points, to increase the resolution. The abrupt change of the plaquette expectation value, visible at smaller values of β, strongly suggests the presence of a bulk phase transition.
To find more concrete evidence of the bulk phase and determine the phase boundary, we increase the volume to 6 4 and 8 4 for two lattice couplings, β = 6.6 and 6.8. In Fig. 18 we plot the trajectories of the average plaquette measured on a 8 4 lattice with various values of bare quark mass close to the transition. All configurations are generated from a cold start-the individual link is the unit matrix. The top panel for β = 6.6 shows evidence of metastability at the critical mass which is expected for a first-order bulk phase transition. By comparison, in the bottom panel of Fig. 18, obtained for β = 6.8, the plaquette values varies smoothly. These results are further supported by measuring the plaquette susceptibility χ = ( P 2 − P 2 ) V , and investigating the dependence of its maximum on the lattice fourvolume V , as shown in Fig. 19. The peaks of χ from the largest two volumes roughly scale with V for β = 6.6, while they are almost constant for β = 6.8, indicating cross-over behaviour.
From the combination of all these numerical results, we find that a conservative estimate of the minimum value of β that ensures the existence of the continuum limit is β ≥ 6.8. Based on this finding we perform fully dynamical simulations at β = 6.9 as a very preliminary study of the meson spectrum. We generate six ensembles, am 0 = −0.  algorithm.
In Fig. 20, we show the trajectories of the plaquette values. The asymptotic value of plaquette gradually increases as the bare fermion mass decreases. The typical thermalisation time appears to be n traj. ∼ 300, while the typical autocorrelation time ranges from 12 to 32, depending on the ensembles. The MD time steps for gauge and fermion actions are optimized such that the acceptance rate in the Metropolis test is in the range of 75 − 85%. Fig. 21 shows the effective masses for pseudo-scalar, vector, and axial-vector mesons at am 0 = −0.91, as an illustrative example. The pseudo-scalar and vector mesons clearly exhibit plateaux at large time, starting around t = 10 and persisting over six time slices. These two mesons are much lighter than the UV cutoff and we find that mπ mρ ∼ 0.8. The best estimation of the axial-vector meson mass is carried out by fitting the effective mass for t = [7 , 9], but the resulting mass is already at the scale of the UV cutoff ∼ 1/a, as visible from Fig. 21. More interesting numerical studies including the EFT fits are beyond the aims of this paper: a dedicated investigation of systematic effects, such as those relevant to the scale-setting and the continuum extrapolation, has to be carried out, before we perform a detailed analysis of confronting the lattice data with the EFT.

Summary and Outlook
With this paper (see also [99][100][101]), we have started a programme of systematic lattice studies of the dynamics of Sp(2N ) gauge theories with N f = 2 fundamental Dirac fermions and N > 1. As explained in the Introduction, we envisage developing this programme along several distinct lines, of relevance in the contexts of composite-Higgs phenomenology, of composite fermions at strong coupling, and of thermodynamics at finite T and µ. We are also interested in the pure Yang-Mills theories (N f = 0), and in studying how the properties of these theories evolve as a function of N . We outlined here the whole programme, and took some important steps along these lines. We focused for the time being on Sp(4), but constructed our numerical algorithms, data-analysis procedures and EFT treatment in such a way that they generalise straightforwardly to larger N . We conclude the paper by summarising our findings, and how they allow us to proceed to the next steps in the near future.
We performed preliminary, technical analyses of the Sp(4) lattice gauge theory, and the results are shown in Sections 3 and 4. We used two different lattice algorithms; we successfully checked that the Hybrid Monte Carlo and Heat Bath both yield results that are compatible with each other as well as with those reported previously in the literature, having introduced an adequate re-symplectization procedure. We tested that the topological charge moves across sectors with dynamics suggesting good ergodicity properties and that the distributions of the ensembles for various choices of the lattice parameters do not show appreciable indications of severe autocorrelation. We also addressed the question of scale setting, by studying the gradient flow associated with the quantities t 0 and w 0 defined in the main body of the paper, and we found visible signals of quark-mass dependence. From the field-theory perspective, this may not come as a surprise, given that the RG flow is twodimensional and non-trivial. Yet, it shows that in calculations with dynamical fermions one has to use extra caution in the process of extrapolating to the continuum limit. We expect that at least when the physical mass is small compared to the confinement scale, and for lattice calculations performed close enough to the continuum limit, the RG flow be driven mostly by the gauge coupling, with small dependence on the mass, as is the case of QCD [65]. We would like to collect evidence of this with values of lattice parameters beyond those employed in this paper (but see also [30,102]).
The main body of this paper mostly focused on the dynamics of the glue. In Sec. 5 we studied in detail the pure (Yang-Mills) Sp(4) theory, showed that it confines in a way that is compatible with the effective string description, and performed the first detailed study of the spectrum of glueballs. We were able to perform the continuum-limit extrapolation of the latter, hence providing a set of determinations for the physical masses that is of quality comparable to the current state-of-the-art for other gauge theories. We could hence compare the spectrum of Sp(4) to that of SU (N ) gauge theories, and in particular we found novel numerical support for two general expectations from the literature: the ratio R between the masses of the lightest spin-2 and spin-0 glueballs is independent of the gauge group [43], and the ratio of the mass of the lightest spin-0 glueball to the string tension obeys Casimir scaling [44].
The long-term objective of this programme is the investigation of whether composite-Higgs models of new physics based upon the SU (4)/Sp(4) coset are realistic and predictive. Having discussed the main features (and limitations) of the low-energy effective field theory description of pions, ρ, and a 1 mesons in Section 2, and postponing to the future the study of dynamical fermions, we performed a first, exploratory calculation of the masses and decay constants of the mesons in the quenched approximation, and reported the results in Sec. 6. The main purpose of this study is to show that the whole technology works effectively. We took particular care of precisely defining the operators of interest on the lattice, and of renormalizing the decay constants with one-loop matching coefficients. We performed the calculations by varying the value of the bare mass over a large range, while considering only two values of the lattice coupling β and one choice of lattice volume.
We found several potentially interesting results for the mesons, although the preliminary nature of this quenched study implies that much caution has to be adopted. While we found that the spectra and decay constants of pions, ρ and a 1 mesons can be fitted satisfactorily with the EFT description we provided, the spin-1 mesons are heavy with respect to the pion decay constant, and their coupling to the pions is large, hence bringing into serious question the reliability of the EFT description itself. We do not know whether this feature persists also with dynamical fermions, yet we expect an improvement with larger values of N , and hence find it encouraging that the fits within quenched Sp(4) work well. We also found several other interesting features. For example, the special combination f 0 (defined in Sec. 2) of the decay constants of the mesons appears to be independent of the quark mass. It will be interesting to test such features beyond the quenched approximation, and possibly explaining them within field theory.
Finally, we uncovered evidence of a first-order (bulk) phase transition in the lattice theory with dynamical fermions, and presented in Sec. 7 the first results of the coarse scanning of the parameter space, hence identifying regions that are safely connected to the field theory in the continuum limit. We exemplified the calculations of the meson spectrum in the full dynamical theory for one choice of such parameters. A much more extensive study of the spectra would be needed to match to the expectations from field theory, particularly because of the subtleties involved in taking the continuum limit for generic, non-trivial values of the fermion mass. Having shown the feasibility of such a study with the instruments we put in place, we postpone to the future this extensive task.
The next steps of our programme will involve the following studies.
• We will compute the spectrum of glueballs in pure Yang-Mills for generic Sp (2N ).
The HB algorithm has been already generalised to any N and tested [101], and the process leading to the extraction of masses has been shown here to be robust. This will allow us to put the level of understanding of the spectra of Sp(2N ) Yang-Mills theories on the same level as the SU (N ) ones.
• The mass spectrum and decay constants of mesons will be studied with dynamical fermions, hence providing quantitative information of direct relevance to modelbuilding and phenomenology in the context of composite-Higgs models of new physics.
• We want to extend the present study to be of relevance to the context of composite fermions, by generalising the underlying action to include fermionic matter in different representations. This is a novel direction for lattice studies, the very first such attempts having appeared only recently [102,103]. We envision to perform a preliminary study, possibly quenching part of the fermions, before attacking the non-trivial (and model-dependent) problem of analysing the properties of fermionic composite states.
Further in the future, we intend to extend the study of the mesons to other non-trivial dynamical properties of relevance to composite-Higgs models, such as the width of the excited mesons, and the value of the condensates. We are also interested in extending to Sp(2N ) the study of the high temperature behaviour of the theory, along the lines followed for SU (2) ∼ Sp(2) in [34], and to introduce non-trivial chemical potential. Combinations of all these studies will provide a coherent framework within which to gain new insight of relevance for field theory, model building, and thermodynamics in extreme conditions.

A Some useful elements of group theory
We choose the generators of SU (4) and of its Sp(4) maximal subgroup as in Appendix B of [34]. We summarise some useful properties of the symplectic groups of interest, which we conventionally refer to as Sp(2N ), and the real algebra of which is denoted C N in [104].
The group Sp(2N ) is defined as the set of 2N × 2N unitary matrices U with complex elements that satisfy the relation U ΩU T = Ω, where Ω is the symplectic form, written-consistently with Eq. (2.2)-in N × N blocks as As suggested by Eq. (A.1), U may be written in block form as where A and B satisfy A † A + B † B = I and A T B = B T A. From these relations, we can deduce many properties of Sp(2N ) matrices. Having unit determinant, the matrices of Sp(2N ) can be shown to form a compact and simply connected subgroup of SU (2N ). Moreover, the structure in Eq. (A.3) implies that the centre of the group is isomorphic to Z 2 for any N . Lastly, since U * = ΩU Ω T and Ω ∈ Sp(2N ), every representation of the group is equivalent to its complex conjugate. Thus Sp(2N ) has only pseudo-real representations, and charge conjugation is trivial. In model-building as well as numerical applications, a prominent role is played by the subgroups of Sp(2N ), especially those isomorphic to some SU (N ). In particular, one notices that Sp(2N ) ⊂ SU (2N ), and that Sp(2(N − 1)) ⊂ Sp(2N ). Starting with Sp(2) ∼ SU (2), this allows us to use the machinery already developed for the Monte Carlo simulation of SU (2N ) groups to the case of Sp(2N ). Particular attention has to be given to the choice of subgroups, as we further discuss in Appendix C in the HB context.
The subgroup structure of Sp(2N ) can be understood in terms of its algebra, to the study of which we now turn. Locally, one can represent a generic group element with the exponential map U = exp(iH) and impose the constraints of Sp(2N ). This is equivalent to taking only the generators of SU (2N ) that satisfy Eq. (A.1), i.e. the hermitian traceless matrices with H * = ΩHΩ, from which a block structure for H follows, The properties A = A † and B = B T are a consequence of H † = H. These conditions leave a total of 2N (N + 1) degrees of freedom for H, which is also the dimension of the Sp(2N ) group. The choice of generators that we use in this work is explicitly stated in [34]. The rank of the group is N , thus in Sp(2N ) we can find N independent SU (2) subgroups. Once the elements of the algebra have been chosen, the SU (2) subgroups of Sp(2N ) follow from their matrix structure (see, once more, Appendix C).

B EFT and Technicolor
The 2-flavour Sp(4) theory can be used as a technicolor (TC) model, provided the embedding of the SM symmetries is such that the condensate breaks them, and so can the EFT treatment we apply to the spin-1 states, provided we identify the natural SU (2) t L × SU (2) t R symmetries acting on the left-handed and right-handed components of Q i a as the SM global symmetry (following the conventions in Appendix B in [34]). We added the superscript t to the weakly-coupled gauge groups to distinguish them from the embedding used in the body of this paper, in the context of composite-Higgs models The embedding is chosen so that (SU (2) t L × SU (2) t R ) ∩ Sp(4) = SU (2), hence realising spontaneous symmetry breaking. One can then match this model to the familiar electroweak chiral Lagrangian and its extensions, based on the (SU (2) × SU (2))/SU (2) coset, which is practically advantageous as it allows to re-use well known results.
Gauge invariance of the electroweak theory in this case requires setting M = 0, so that the pions are massless, 15 and to gauge SU (2) t L × U (1) t Y , where the latter factor is the subgroup of SU (2) t R generated by the diagonal t 3 R . The gauge couplings areg andg , respectively. By doing so, the exact symmetry is reduced from SU (4) to SU (2) t L × U (1) t Y × U (1) B , with the last abelian factor denoting baryon number. The condensate then breaks SU (2) t L × U (1) t Y → U (1) e.m. , and the W and Z bosons acquire a mass via the usual Higgs mechanism.
We discuss this TC model mostly for completeness, and for technical reasons related to the calculation of the 2-point functions. Among the phenomenological reasons why this is not a realistically viable TC model are the following.
• Precision parameters such as S might exceed experimental bounds.
• The spectrum does not contain a light scalar (Higgs) particle.
• There is no high scale to suppress higher-dimensional operators and FCNC transitions.
• There is no natural way to enhance the mass of the top quark.
We address here only the first two points, within the low-energy EFT, mostly for technical reasons that are of interest also in the composite-Higgs framework. We introduce the adjoint spurion G, formally transforming as We add to the Lagrangian density the additional (symmetry-breaking) term in such a way as to induce the explicit breaking SU (4) B → SU (2) t L × SU (2) t R × U (1) B . By expanding explicitly we find that
We can now focus on the 2-point functions, by matching the theory defined by Eq. (2.16) onto the leading-order part of the effective Lagrangian for the transverse components of the SM gauge bosons V i µ = (L 1 µ , L 2 µ , L 3 µ , R 3 µ ), which reads where P µν = −q µ q ν /q 2 + η µν , q µ is the four-momentum, and all the dynamics is contained in the non-trivial functions Π ij (q 2 ). The functions Π ij are defined by matching the (gaussian) path integrals (at the treelevel). In practice, one takes the second derivatives in respect to the 15 + 3 + 1 = 19 gauge bosons V µ i of the original theory P µν Π ij T ≡ P µρ ∂ 2 ∂V ρi ∂V νj L = P µν q 2 δ ij − (M 2 ) ij (where M 2 is the complete mass matrix of the vectors), then one inverts it and retains only the 4 × 4 sub-matrix along the directions of SU (2) t L × U (1) t R , and finally inverts it again to obtain Π ij = Π −1 Focusing on the 1 and 2 components of Π ij , one can write Π 11 (q 2 ) = Π 22 (q 2 ) = q 2 − By expanding Π ij (q 2 ) in powers of q 2 one also obtains the precision electroweak parameters. They are defined by first normalising the fields V µi so that Π ii (0) ≡ d dq 2 Π ii (q 2 ) q 2 =0 = 1 for i = 1 , · · · , 4, and then defining [105] S ≡g g In the current case, one finds thatT = 0, because of the custodial SU (2), while the experimental bounds obtained from precision parameters areŜ < 0.003 (at 3σ c.l.) [105].
A subtlety should be noted here: the structure of the effective Lagrangian describing all the vectors in the case of relevance for TC should include additional terms in respect to Eq. (2.16), as in this case the condensate breaks the SU (2) t L symmetry. For example, this yields kinetic mixing between the W and B gauge bosons, which contributes directly to theŜ parameter, and can be thought of as the contribution toŜ coming from heavier composite states that have been integrated out. We do not include such terms, because they go beyond the purposes of this paper.
If one applies the formalism described here to the Lagrangian density in Eq. (2.16), what results is not the S parameter, 16 but rather the contribution to S coming from ρ and a 1 mesons only, which we called S 0 . We notice that the following sum rule (occasionally referred to as zeroth order in the literature) is verified exactly within the framework defined by Eq. (2.16): We can express this result explicitly in terms of the coefficients in the EFT Lagrangian density: We stress once more that this result holds only within the EFT in Eq. (2.16), but not in the fundamental theory, 17 for which one would have to replace the right-hand side of Eq. (B.8) with a summation over all possible spin-1 states, not just the ground state. To include such contributions in the EFT would require introducing explicit symmetry-breaking terms in L of the form of kinetic mixing terms that are forbidden within the composite-Higgs scenario, and are hence omitted here.
We conclude with an exercise. By making use of the quenched calculations reported in subsection 6.2, and of the fits of the low-energy couplings in Eq. (2.16), we find that S 0 /4π ∼ 0.049 (for β = 7.62), and S 0 /4π ∼ 0.041 (for β = 8.0). The large values of S 0 , in respect to QCD (for which the same treatment would yield S 0 /4π ∼ 0.025, obtained by just replacing the masses and decay constants from experiments [109]), might be due, at least in part, to the dimension of the gauge group Sp(4) being larger than that of QCD. Yet these numerical results show a significant decrease in approaching the continuum limit, and are affected by quenching, hence they must be taken as just for illustration purposes.

C Projecting on Sp(4)
The update of local gauge links must be supplemented by a projection onto the Sp(4) group manifold, to remove machine-precision effects that would bring the algorithm outside of the 16 The normalization of S by Peskin and Takeuchi [106] is such that S 4π ≡ 4 g 2Ŝ , which yields a result independent of the SM gauge couplingg, contrary to the phenomenologically more convenientŜ [105]. 17 We refer the reader to Ref. [107] for an example of such calculations on the lattice, specialised to the SU (3) gauge theories with n f = 2 and 6 fermions in the fundamental representation, and to [108] for an earlier calculation of S within QCD.
group manifold. In this Appendix we explain how we implement this re-symplectisation procedure.
To gauge the potential size of the violation of symplectic conditions and the importance of the Sp(4) projection, we plot in Fig. 22 the trajectories of the plaquette values with (top) and without (bottom) re-symplectisation. The result shows that the re-symplectisation is successfully working and the plaquette values are stable.
In the HB calculations, we use instead a variant of the (modified) Gram-Schmidt algorithm. A re-symplectisation algorithm that inherits the numerical stability of the Gram-Schmidt process and that can be applied to any N is obtained by noting that, owing to the general form of Sp(2N ) matrices, once the first N columns of the matrix are known, the remaining ones can be computed from After normalising the first column, one can obtain the (N +1)-th. The second column is then obtained by orthonormalisation with respect to the first and the (N +1)-th. Repeating the process for every column, one obtains an Sp(2N ) matrix. [7] K. Agashe, R. Contino