Sp (4) gauge theories on the lattice: Nf = 2 dynamical fundamental fermions

We perform lattice studies of the gauge theory with Sp(4) gauge group and two flavours of (Dirac) fundamental matter. The global SU(4) symmetry is spontaneously broken by the fermion condensate. The dynamical Wilson fermions in the lattice action introduce a mass that breaks the global symmetry also explicitly. The resulting pseudo-Nambu-Goldstone bosons describe the SU(4)/Sp(4) coset, and are relevant, in the context of physics beyond the Standard Model, for composite Higgs models. We discuss scale setting, continuum extrapolation and finite volume effects in the lattice theory. We study mesonic composite states, which span representations of the unbroken Sp(4) global symmetry, and we measure masses and decay constants of the (flavoured) spin-0 and spin-1 states accessible to the numerical treatment, as a function of the fermion mass. With help from the effective field theory treatment of such mesons, we perform a first extrapolation towards the massless limit. We assess our results by critically comparing to the literature on other models and to the quenched results, and we conclude by outlining future avenues for further exploration. The results of our spectroscopic analysis provide new input data for future phenomenological studies in the contexts of composite Higgs models, and of dark matter models with a strongly coupled dynamical origin.


Introduction
The Large Hadron Collider (LHC) recently discovered a new scalar particle [1,2], which has experimental properties compatible with those of the Higgs boson and mass ∼ 126 GeV. This discovery stands out against the current absence of clear evidence of new physics beyond the Standard Model (SM), both in direct and indirect experimental searches, up to and beyond the TeV scale-evidence of a little desert in high energy physics. Composite Higgs Models (CHMs) implement the symmetry-based mechanism first proposed in Refs. [3][4][5]. They provide a compelling framework that allows to soften the level of fine-tuning required to accommodate the little desert within an Effective Field Theory (EFT) description of Electro-Weak Symmetry Breaking (EWSB).
In the SM, EWSB is induced at the scale v W ∼ 246 GeV by the dynamics of a complex doublet of weakly-coupled Higgs fields. Composite Higgs models reinterpret its components as some of the interpolating fields appearing in the EFT that describes the low energy, weakly coupled dynamics of a set of composite pseudo-Nambu-Goldstone bosons (pNGBs). They emerge from a more fundamental-possibly strongly coupled and UV completetheory, that dynamically drives the spontaneous breaking, at scale f > v W , of an extended approximate continuous global symmetry. After coupling this EFT to the SM gauge bosons and fermions, what results is a natural and stable dynamical origin for the little hierarchy between the masses of the SM particles (the Higgs boson in particular) and the higher scale beyond which new phenomena arise. The ratio v W /f < 1 is determined by the interplay of strong-coupling dynamics and weakly coupled, symmetry-breaking perturbations. Electroweak symmetry breaking is triggered via what, in the jargon of the field, is often referred to as vacuum misalignment-a phrase that highlights the intrinsic differences with other models of EWSB with a strongly coupled dynamical origin [6].
A wide variety of implementations of these ideas has been proposed in recent years, and their model-building features, phenomenological implications, and dynamical properties are the subject of a rich, diverse, and rapidly evolving literature (see for instance Refs. ). Particular attention has been devoted to models and dynamical theories characterised by the SU (4)/Sp(4) coset (see for instance Refs. [39][40][41][42][43][44][45][46][47][48][49][50][51][52][53][54][55][56][57][58]). The low-energy EFT Lagrangian contains, in this case, five scalar fields, that describe the long-distance dynamics of the five pNGBs. With an appropriate choice of embedding for the SU (2) × U (1) gauge group of the SM, four such fields reconstruct the complex Higgs doublet familiar to the Reader from the Standard Model, with the fifth scalar a new real, neutral singlet, extending the SM Higgs sector.
The investigations summarised in these pages contribute to the study of the strong dynamics underlying the SU (4)/Sp(4) theory, under the working assumption that it originates from a Sp(2N ) gauge theory with two fundamental (Dirac) fermions, which is amenable to lattice numerical treatment. We ultimately aim at computing the many free parameters of the low-energy EFT, by starting from fundamental principles. This paper summarises the findings of the second stage of development of the programme outlined in Ref. [59] (see also [60][61][62][63]), by focusing attention on the case in which the matter field content consists of two dynamical Dirac fermions transforming in the fundamental representation of the Sp(4) gauge group.
The main step forwards we make here is that we move beyond the quenched approximation adopted in earlier explorative work [59], as we implement in its stead fully dynamical (Wilson) fermions on the lattice. In the presence of two fundamental lattice parameters, the scale-setting process has to be reconsidered-especially in the dynamical regime away from the chiral limit. We compute the spectrum of the lightest mesons, which are organised in irreducible representations of the unbroken global Sp(4) group. We mostly discuss scalar S, pseudoscalar PS, vector V and axial-vector AV composite flavoured particles. For completeness we also analyse the mass of states sourced by both the antisymmetric tensor T and its axial counterpart AT, although only the latter sources genuinely independent states-the tensor T and vector V operators source the same states. We extract masses and (appropriately renormalised) decay constants from 2-point functions of the relevant interpolating operators. We assess the size of finite volume effects and perform continuum limit extrapolations.
An important limitation of this work is that the physical masses of the pNGBs are large enough that the vector mesons are effectively stable. While we attempt an extrapolation towards light masses of relevance to phenomenologically viable CHMs, the large-mass regime is interesting in itself, being relevant for models of dark matter with a strongly coupled dynamical origin, along the lines discussed in Refs. [64][65][66]. A crucial piece of dynamical information in this context turns out to be the strength of the coupling g VPP between V and two PS mesons, which plays an important role in controlling the relic abundance.
The g VPP coupling can in principle be extracted by careful analysis of 4-pNGB amplitudes [67][68][69] (see also Ref. [70] for an application in the context of new physics), but the amount of data generated for this paper is not sufficient to even approach this gargantuan task, which we leave for future dedicated studies. We instead perform a first, preliminary extrapolation of our results towards the massless limit, with help from low-energy EFT instruments. The EFT treatment we proposed in Ref. [59] is based on the ideas of hidden local symmetry (HLS), adapted from Refs. [71][72][73][74][75] (see also Refs. [76][77][78][79]), and supplemented by some additional, simplifying working assumptions. This process allows us not only to estimate the masses and decay constants of the spin-1 states in the regime relevant to electroweak models, but also to extract an estimate of g VPP , hence providing a first, possibly rough measurement of its size based on a numerical, dynamical calculation. We also discuss several non-trivial features of the spectra, and compare them to previously published results obtained in other related gauge theories as well as the results of quenched calculations, that are reported elsewhere [80].
The paper is organised as follows. In Section 2 we introduce the model by defining its lattice action. We recall some useful notions about the Hybrid Monte Carlo (HMC) algorithm we employ, and present our choices of lattice parameters. In Section 3 we employ the gradient flow method to set the physical scales. The mass dependence of the flow scale can be understood in EFT terms [81]. We also discuss the size of finite-volume artefacts. In Section 4 we present our numerical lattice results for the spectra of mesons and renormalised decay constants. We define the corresponding interpolating operators and analyse their 2point correlation functions. Details on the HMC algorithm, diquark operators and the fits of correlation functions are presented in Appendix A. We present our strategy and perform the continuum limit extrapolation by employing a mass-dependent prescription [82] introduced through the Wilson chiral perturbation theory (WχPT) [83,84] (see also Ref. [85], and the literature on improvement [86,87]), and report the results in Section 4.3. Given the extensive amount of information we are communicating, we find it useful to conclude this part of the paper with a short summary of the lattice numerical results in Section 5.
From Section 6 onwards, we restrict the discussion to the study of the continuum limit extrapolations obtained by using only a set of eleven ensembles selected in Sec. 4.3. We deploy our EFT tools and perform extrapolations towards the massless limit to determine the Low Energy Constants (LECs). We also critically discuss implications, applications, and limitations of the resulting numerical fits. Details of the numerical results are presented in the histograms of Appendix B. Section 7 is devoted to comparing our results to the analogous observables in other theories, by borrowing published data available in the literature, as well as to the results obtained within the quenched approximation [80]. As we shall see, besides providing an important sanity check, the latter also allows us to assess the impact of quenching on 2-point correlators, information that might be of general value, as it provides guidance towards future studies of Sp(2N ) theories with N > 2. We provide a summary of the most important results in the continuum in Section 8, based upon the detailed information provided in Sections 6 and 7. We conclude the paper with a short list of open avenues for future exploration in Section 9.

Lattice model
A distinctive feature of Sp(2N ) gauge theories with N f massless Dirac fermions in the fundamental representation is the enhancement of the global symmetry to SU (2N f ), which originates from the pseudo-real character of the representation.
The lattice formulation in terms of Wilson fermions introduces an operator that breaks the global symmetry, and introduces a (degenerate) mass term for the fermions. The global symmetry is expected to be spontaneously broken by the formation of a non-zero fermion bilinear condensate. With N f = 2, both explicit and spontaneous breaking follow the (aligned) SU (4) → Sp(4) pattern. The resulting low-energy dynamics is governed by five pNGBs. They describe the SU (4)/Sp(4) coset, and have degenerate masses.

Lattice action
The four-dimensional Euclidean-space lattice action contains the gauge-field term S g , together with the fermion matter-field term S f : We use the standard Wilson plaquette action for the discretised gauge fields, with the gauge links U µ group elements of Sp(4) in the fundamental representation: The plaquette P µν is defined by The trace Tr is over the fundamental of Sp(4), and the lattice coupling is given by β = 8/g 2 . We define the fermion sector by using the (unimproved) Wilson action for two massdegenerate Dirac fermions Q in the fundamental representation where a is the lattice spacing and am 0 is the bare mass in lattice units.

Numerical Monte Carlo treatment
We use the lattice action in Eq. (2.1) to study the Sp(4) theory with N f = 2 Dirac fundamental fermions, with the gauge configurations generated by the standard Hybrid Monte Carlo (HMC) algorithm. In Ref. [59] we extensively discussed the relevant numerical techniques adopted, including the need to project back onto the symplectic group after each HMC update of the configurations, and the associated modifications to the HiRep code [88]. During this study, we have further improved the code to implement arbitrary values of N ≥ 2 and to reduce the storage size of an individual gauge configuration by a factor of two-details are presented in Appendix A.1.
Pioneering lattice studies of Sp(2N ) Yang-Mills showed that a bulk phase transition is absent in the Sp(4) theory, implying that one can in principle take the continuum limit by choosing any values of β [89]. By contrast, in the case of dynamical simulations with two Wilson-Dirac fermions, the preliminary study of the average plaquette value P ≡ (24N t N 3 s ) −1 Re x µ<ν Tr P µν (x) detected evidence of a first-order bulk phase transition [59]-N t and N s are the temporal and spatial extents of the lattice, respectively. Hybrid Monte Carlos trajectories of P started from cold (unit) and hot (random) configurations at small lattice volume show signs of hysteresis. Careful study of the volume dependence of the plaquette susceptibilities indicates that the continuum extrapolation can be carried out safely when β > ∼ 6.8. In this regime, the desired continuum limit can be reached safely, and subsequently the fermion mass can be lowered smoothly, avoiding the unphysical Aoki phase near the massless limit [90]. For all ensembles considered in this work, we operate far enough from the massless limit and no sign of the Aoki phase is visible.
The parameters characterising the ensembles generated by the dynamical simulations are summarised in Table 1. In the table we present the values of lattice coupling β and bare fermion mass am 0 . The former is chosen to be in the range 6.9 ≤ β ≤ 7.5. The choices of the latter will be discussed later in the paper, where we will see that some of the larger choices of am 0 will not be used in the analysis.
The four-dimensional Euclidean lattice has size N t ×N 3 s , and we impose periodic boundary conditions in all directions for the gauge fields. The physical volume  Table 1: List of ensembles generated for this study. The lattice parameters β and am 0 are, respectively, the bare coupling and bare fermion mass. The lattice sizes are denoted by N t × N 3 s , separately highlighting the time-like and space-like dimensions. The number of configurations used to estimate the average plaquette P and the gradient flow scale w 0 /a is denoted by N configs . The separation of trajectories between adjacent configurations is denoted by δ traj . In the results for P and w 0 /a, in parenthesis we indicate the statistical error.
obtained by setting T = N t a and L = N s a. For the Dirac fields we implement periodic boundary conditions in the spatial directions, but anti-periodic boundary conditions in the temporal one.
We anticipate that all lattice volumes satisfy the condition m PS L > ∼ 7.5, where am PS denotes the mass of the lightest pseudo-scalar meson (expressed in lattice units). The latter is extracted from the two-point correlation functions, as will be discussed in details in Section 4 (see also Table 5). As we shall see later in the paper, this choice guarantees that the volumes are large enough that the finite-size effects are under control. In the table we also present the results for the average plaquette P and the gradient flow scale w 0 defined in Section 3.1, which are measured from N configs configurations separated by δ traj trajectories. Throughout this work we estimate the statistical uncertainties by using a standard bootstrapping method for resampling [91].

Scale setting, topology and finite volume effects
Lattice calculations yield dimensionless numbers. The inverse of the lattice spacing a acts as a hard momentum cut-off Λ cut , and all lattice measurements in lattice units can be written as a function of a. For example, a dimensionless mass m lat can be expressed as m lat = ma, with m having canonical units. But this is not sufficient to take the non-trivial limit a → 0, as the lattice spacing does not have a precise counterpart in the continuum theory. A physical quantity that can be measured both in the continuum as well as in the discretised theory must be used to set a common reference scale, and yield a scale setting procedure within the continuum extrapolation.
We adopt Lüscher's Gradient Flow (GF) technique [92]. Besides achieving good accuracy with modest numerical effort, this procedure has two other major advantages. The reference scale is defined on fully theoretical grounds, which is convenient for theories that have not been tested experimentally. Furthermore, it preserves the topological charge Q, while strongly suppressing ultraviolet (UV) fluctuations. In this section, after carrying out the scale-setting programme, we also discuss the topology of the ensembles. We conclude by studying finite volume effects and arguing that they are smaller than the statistical uncertainties.

Gradient flow and scale setting
We denote by A µ (x) the four-dimensional non-abelian gauge field evaluated at the spacetime coordinates x. The gradient flow is defined in the continuum theory by a diffusion equation (in Euclidean five-dimensional space) for a new gauge field B µ (t, x) at the fictitious flow time t. The equation reads: where D ν is the covariant derivative in terms of B ν , while G µν = [D µ , D ν ] is the fieldstrength tensor. Repeated indices are summed over. Along the flow time the gauge fields evolve into renormalised gauge fields, smoothed over a radius of √ 8t, the characteristic scale of the diffusion process. As shown in Ref. [93], at t > 0 the correlation functions of the renormalised fields are finite to all orders in perturbation theory. In particular, the following gauge-invariant observable does not require any additional renormalisation other than that at zero flow time (t = 0): The expectation value of E(t) is proportional to the inverse of the flow time squared. We consider two different proposals for defining the gradient flow scale, and denote them by t 0 [92] and w 0 [94]. We first define the dimensionless observables at positive flow    Table 1. and Then the scales are set by imposing the conditions and Here E 0 and W 0 are common, dimensionless reference values. In numerical studies, we measure the dimensionless quantities t 0 /a 2 and w 0 /a, which determine the relative size of the lattice spacing between ensembles obtained by using different (bare) lattice parameters. In this project, consistently with our previous work [59], we employ the Wilson-flow method [92] to proceed with the lattice implementation of Eq. (3.1). In our previous publication [59], we performed detailed numerical studies of the GF scheme for the quenched theory, as well as full dynamical calculations for β = 6.9. We found that w 0 shows smaller cut-off-dependent effects, compared to t 0 . In particular, no significant deviation was found between the values of w 0 obtained by using the action density at non-zero flow time E(t) constructed from the average plaquette and from the symmetric four-plaquette clover, as defined in [92].
In this study, we consider a finer lattice with β = 7.2. The results are presented in Fig. 1. We find that while the values of t 0 show significant discrepancies, the measured values of w 0 from the two definitions of E(t) are in good agreement over the wide range of W 0 and m 0 we considered, in particular for W 0 = 0.3 ∼ 0.4. The agreement in the flow scales has improved with respect to the results from coarser lattices in [59]. In Table 1, and in subsequent calculations, we elect to use the gradient flow scale w 0 , which we compute with the reference value of W 0 = 0.35, on the four-plaquette clover action density-for which smaller lattice artefacts are observed. For convenience, we introduce the following notation:m ≡ m lat w lat 0 = mw 0 denotes the dimensionless quantity corresponding to a mass. We useâ ≡ a/w 0 when we discuss lattice-spacing artefacts in Section 4.2. Figure 1 shows that the scales √ 8t 0 /a and w 0 /a depend on the fermion mass am 0 . The title of this subsection is borrowed from Ref. [81], to reflect the fact that we employ the EFT treatment suggested in this reference and we apply it to our numerical results. The EFT treatment assumes that the square root of the flow scale t 0 is much smaller than the Compton wavelength of the pseudoscalar meson.

Chiral perturbation theory for gradient flow observables
Following [81], we use the leading order (LO) relation in the chiral expansion m 2 PS = 2Bm f (where m f is the fermion mass), to write the next-to-leading-order (NLO) result for the GF scale w NLO 0 as where w χ 0 is the GF scale and f PS is the pseudoscalar decay constant, both defined in the chiral limit. It is convenient to rescale this expression by writing wherek 1 , w χ 0 /a andm PS = w 0 m PS are dimensionless parameters 1 . We also find it convenient to report here on the extraction of the constantsk 1 and w χ 0 from the dynamical ensembles in Table 1. This requires that we anticipate the use of numerical data for the measurement of m PS /a that will be discussed extensively in Section 4 and Appendix A.3, and will be reported in Table 5. Figure 2 shows data from Table 1 combined with Table 5, together with the result of the two separate fits (for β = 6.9 and 7.2) to Eq. (3.8) of the five ensembles with smallest m PS . The resulting values of the fit parameters are reported in Table 2. The values of χ 2 /N d.o.f. at the minima indicate that chiral perturbation theory for w 0 well describes the data. Deviations from linear mass dependence of w 0 /a appear at aroundm 2 PS ∼ 0.4. We 1 The difference between w0 mPS and w χ 0 mPS is a sub-leading effect, which would appear at next-tonext-to-leading order (NNLO).  Figure 2: Gradient flow scale dependence on the pseudoscalar meson mass for β = 6.9 (left panel) and β = 7.2 (right panel). Numerical data from Table 1 and from  (10) 0.5 Table 2: Results of the NLO fits for w 0 /a from Table 1 andm PS from the combination with Table 5. The fit uses the five ensembles with smallest mass to extract the parameters k 1 and w χ 0 /a in Eq. (3.8).
anticipate that this scale is in broad agreement with the upper bound inferred by studying the pseudoscalar decay constant, and will be discussed in Section 4. We observe that the value ofk 1 is smaller for the finer lattice. We have too few ensembles for other lattice couplings to extract the values ofk 1 , yet the generic trend is consistent with expectations, as visible in Fig. 3, with the mass-dependence becoming milder for larger choices of β. 2 Later in the paper, we will perform simultaneous continuum and massless extrapolations via a global fit of all measurements of physical quantities-masses and decay constants of mesons-expressed in units of w 0 .

Topology
By analogy with the continuum definition ( 1 , the lattice topological charge of a gauge configuration is defined by summing over lattice sites i as   Figure 3: Inverse of the gradient flow scales w 0 with respect to the pseudoscalar meson mass squared, for data taken from Table 1 and from Table 5. Blue, purple, green, red and brown colours (top to bottom, approximately linear series) have β = 6.9, 7.05, 7.2, 7.4 and 7.5, respectively.
The HMC algorithm yields gauge configurations in which ultraviolet fluctuations have typical sizes that are orders of magnitude larger than the desired signal. The resulting large cancellations prevent a reliable extraction of Q. A smoothing procedure must be applied, that preserves the topological charge of a configuration while removing UV fluctuations, in order to measure Q. The gradient flow described in the earlier subsections can be used for this purpose. For each ensemble, we measure Q from uncorrelated, thermalised configurations, after flowing to a flow time t/a 2 = N 2 t /32 (equivalent to √ 8t = T /2). With our choice of boundary conditions, at finite lattice volume, Q is quantised, and assumes only integer values. Infinitesimal changes in the field configuration cannot alter the topological charge. The change at each Monte Carlo time-step required to yield a good acceptance rate in the Metropolis step of the HMC becomes smaller for fine lattice spacings and small fermion masses. For this reason, at large volumes and small masses, there is a risk that the topological charge freezes at a single value, for the finite HMC trajectories that one implements in practice. Since values of observables may depend on which topological sector is being observed, this effect may introduce a systematic error. To ensure that our measurements are not (heavily) affected by this type of systematic effect, we plot and inspect histories of Q for each ensemble used.
We expect the measurements of the charge Q to obey a Gaussian distribution about zero in the limit of infinite simulation time. We use the fitting form in order to fit the histogram data. All ensembles used to generate the results presented in this paper were found to be free of significant topological freezing, with Q 0 always  compatible with zero within errors (at the 1σ level). A sample of topological charge histories  Table 3: Ensembles and numerical results used to estimate the size of finite volume effects. The number of configurations and the separation of trajectories between adjacent configurations are given by N configs = 200 and δ traj = 12, respectively, for all these ensembles. The pseudoscalar masses in the infinite volume limit m inf PS are estimated from the exponential fits as discussed in the text.
is shown in Fig. 4.

Finite size effects
The finiteness of the lattice volume represents an inherent source of systematic uncertainties in lattice calculations. In a confining gauge theory, finite-volume (FV) contaminations of the masses and decay constants are expected to be exponentially suppressed as a function of the length of the spatial lattice L, at least when the volume is much larger than the inverse of the Compton wavelength of the lightest state, i.e. as long as m PS L 1. The size of FV effects can be estimated in a systematic way within chiral perturbation theory, as long as the volume is larger than the hadronic scale. For f PS L 1, the dominant contribution is expected to arise from one-loop tadpole integrals of the pseudoscalar mesons [95,96]. The leading-order FV corrections to the meson masses can be written as where m inf M(PS) is the meson (pseudoscalar) mass in the infinite volume limit. In principle, the coefficient A M can be determined within chiral perturbation theory without introducing any new parameters in the infinite-volume chiral Lagrangian. Nevertheless, since our data is far from the massless limit, we treat A M as free.
To quantify the size of FV effects, we calculate masses of pseudoscalar PS and vector V mesons am PS(V) on three spatial lattice volumes of N 3 s = 16 3 , 20 3 , 24 3 at two sets of lattice parameters, (β, a m 0 ) = (7.2, −0.77) and (β, a m 0 ) = (7.2, −0.79). Numerical results are reported in Table 3. At the smallest volume m inf PS L ∼ 5, the masses deviate from those at the largest volume by 4 ∼ 5%, more than expected from statistically uncertainties. At m inf PS L ∼ 6.3 the deviations decrease to the level of 1 ∼ 2%, compatible with the statistical uncertainties for am V , but not for am PS . We performed a fit of the pseudoscalar masses to Eq. (3.11). As shown in Fig. 5, the data are well described by the exponential fit function: Figure 5: Volume dependence of the pseudoscalar masses, as in Table 3 the blue dashed lines correspond to the best fits, while the black solid lines indicate the masses in the infinite volume limit. From this analysis we conclude that the size of FV effects is no larger than ∼ 0.3%-smaller than the typical size of statistical uncertainties in the spectroscopic measurements reported in the following sections-if we require that m PS L 7.5. In the rest of this work we restrict attention to ensembles that satisfy this requirement-see Table 5.
The gradient flow scale also receives a correction from the finite size of the lattice volume. The flow along the fictitious time t can be understood as a smearing procedure with scale √ 8t, hence FV effects are controlled by the dimensionless ratio c τ = √ 8t/L. At the reference value of E 0 = 0.35 we find that c τ does not exceed 0.2 in most ensembles. Using the results in Ref. [97] we estimate the size of FV corrections to Eq. (3.4) to be at most at a few per-mille level. The only exception are the ensembles with β = 7.2 on a 36 × 16 3 lattice and the one with β = 7.5. These ensembles do not play an important role in the analyses that follow, because of the large physical masses associated with them.

Meson spectroscopy and decay constants
In this section we summarise the main results of our lattice study, by focusing on the properties (masses and decay constants) of the mesons in the dynamical theory. We start in Section 4.1 by defining the operators we are interested in, the correlation functions we measure, and the renormalisation procedure we apply. We then present all the main spectroscopy results in Section 4.2, and discuss the continuum limit extrapolation in Section 4.3. Useful supplementary details are relegated to Appendix A.2 and A.3.

Two-point correlation functions
Following established procedure, we extract the masses and the decay constants of flavoured mesons by studying the behaviour of the relevant two-point correlation functions at large To avoid mixing with the flavour singlets, we restrict to i = j the flavour indices of the Dirac fermions, while colour and spinor indices are summed and omitted. For completeness, we also show the J P quantum numbers and the corresponding particle in the QCD classification of mesons. Notice that two of the operators source the same particles (ρ meson) because of the breaking of chiral symmetry. In the last column we report the irreducible representation of the unbroken global Sp(4) spanned by the meson (see also [42]). In brackets are irreducible representations of Sp(4) that are sourced by operators with the same Lorentz structure, but that we do not discuss in this context.
Euclidean time t. The interpolating operators which carry the same quantum numbers with the desired meson states take the generic form where i, j = 1, 2 are the flavour indices, and Γ M refer to the Dirac structures summarised in Table 4. Summations over spinor and colour indices are understood. The lightest spin-0 and spin-1 mesons, denoted by PS, V and AV for pseudoscalar, vector and axial-vector mesons, respectively, appear in the low-energy EFT described in [59]. In Section 6 we will use their masses and decay constants to test the EFT and to extrapolate towards the chiral limit. We also consider additional interpolating operators with the Dirac structures 1, γ 0 γ µ , γ 5 γ 0 γ µ , referred to as scalar S, (antisymmetric) tensor T, and axial tensor AT, though from the related correlation functions we only extract the meson masses. We restrict our attention to the flavoured mesons (by choosing i = j in flavour space)-the analogous mesons in QCD are π, ρ, a 1 , a 0 , and b 1 . As the global symmetry is broken, the states created by the interpolating operators denoted by V and T mix, with the low-lying state corresponding to the ρ meson in QCD (see [98] and references therein, as well as Fig. 1 of both Refs. [99] and [100]).
For all the meson interpolating operators O M listed in Table 4, we define the zeromomentum Euclidean two-point correlation functions at positive Euclidean time t as In the numerical study, the resulting mesonic two-point correlation functions are studied by replacing the point-like sources in Eq. (4.1) with Z 2 × Z 2 single time slice stochastic sources [101], with number of hits 3.
Because of the pseudoreal nature of the representation of the symplectic gauge group, diquark operators are indistinguishable from the mesonic operators. We report in Table 4 also the multiplicity of each state (the size of the irreducible representation of Sp (4)). For instance, five pNGBs form a multiplet of the unbroken Sp(4), in the SU (4) → Sp(4) enhanced symmetry pattern of the gauge theory considered here. Compared with what happens with gauge group SU (N ), these five pNGBs include the three associated with the breaking SU (2) L × SU (2) R → SU (2) V , together with two diquarks 3 . In Appendix A.2 we explicitly show the equivalence of meson and diquark correlators by using the lattice action in Eq. (2.4).
At large Euclidean time t the correlation functions in Eq. (4.2) are dominated by the lowest excitation at zero spatial momentum so that the mass m M appears in the asymptotic expression: where T is the temporal extent of the lattice. The decay constants f M are determined from the matrix elements, which are parameterised as The polarisation vector µ is transverse to the momentum p µ and normalised by * µ µ = 1.
The meson states |M are conventionally defined by the self-adjoint isospin fields, as in M = M A T A , where T A are the generators of the group. We adopt conventions such that in QCD the analogous experimental value of the pion (pseudoscalar) decay constant is f π 93 MeV. In Eq. (4.4), the pseudoscalar decay constant f PS is defined via the local axial current. To calculate the decay constant f P S , we introduce an additional two-point correlation function . In practice, we calculate m PS and f PS by performing a simultaneous fit to the numerical data for C O PS (t) and C Π (t).
The details of the fit of the meson correlators, including the effective masses and best-fit ranges, are provided in Appendix A.3.
The matrix elements in Eq. (4.4), calculated from the lattice at finite lattice spacing a, must be converted to those renormalised in the continuum. For Wilson fermions the decay constants in the continuum are determined from lattice ones via where Z V and Z A are the renormalisation factors for vector and axial-vector currents which are expected to approach unity in the continuum. Since the pseudoscalar decay constant f PS is defined using the axial current as in Eq. (4.4), it receives renormalisation with the factor of Z A . The renormalisation factors are determined by the one-loop renormalisation procedure in lattice perturbation theory for Wilson fermions, and the expressions for the matching factors are the following [103]: (4.7) The eigenvalue of the quadratic Casimir operator with fundamental fermions is C(F ) = 5/4 for the Sp (4) .7) is defined via the mean field approach to the link variable, which effectively removes the tadpole diagrams, asg 2 = g 2 / P [104], where P is the average plaquette value and g the bare gauge coupling.

Masses and decay constants
Using the techniques described in the previous subsection, we calculate meson masses and decay constants for the ensembles in Table 1. The resulting values in lattice units are summarised in Table 5 for mesons sourced by PS and S operators, and in Table 6 for those sourced by V, AV, T, and AT operators. The decay constants in the tables are renormalised as in Eq. (4.6). As expected, the pseudoscalar mesons are the lightest states for all ensembles. In Table 5 we also present the numerical values of m PS L and f PS L. The lattice volumes in all ensembles are large enough that m PS L 7.5, and, as we discussed in Section 3.4, finite volume effects are negligible. Furthermore, all the values of f PS satisfy the condition f PS L > 1, to ensure that the lattice volume is large enough to capture the chiral symmetry breaking scale. The ensembles to be used for the massless and continuum extrapolations are further restricted to f PS L 1.5.
Mesons with higher mass (and spin) can in principle decay into 2 and/or 3 pseudoscalars [53], but for all the ensembles we considered they cannot decay due to large pseudoscalar masses. Table 6 shows that the masses of mesons sourced by V and T operators are in good agreement, within current statistical uncertainties. As is the case within QCD, low-lying states identified by a given set of quantum numbers (I, J P C ) can result from admixture 0.8344 (11) 0.1431 (7) 1.52 (4) 13.351 (17)    of more than one possible operator when the global symmetry is strongly broken [98]. We also find that the masses of three heavier states, sourced by S, AV, and AT operators, and that have different quantum numbers, are close to one another. They are affected by large statistical and systematical uncertainties when their masses are extracted from the two-point correlation functions (see Appendix A.3 for technical details). In order to compare to one another results obtained from ensembles at different bare parameters, we must relate them to the corresponding renormalised quantities in the continuum limit. We do so by adopting the GF scheme as explained in Section 3.1: we define the meson masses and decay constants in units of w 0 using the notation We can compare measurements at different lattice couplings β, and use the EFT to remove residual artefacts due to the discretisation. This procedure will be explained in detail in section 4.3.
Ensemble (17)    On the lattice, the bare fermion mass m f is a free parameter. It can vary from being very small (yielding very light pseudoscalars) to assuming somewhat large values (yielding stable vector mesons). A detailed study of the renormalisation of the fermion mass m f requires a dedicated study that goes beyond the purposes of this work. In its stead, we replaced m f by the pseudoscalar mass squared m 2 PS , which is a physical quantity. In the low-energy EFT, the two are related through the LO relation m 2 PS = 2Bm f , with B one of the coefficients in the chiral Lagrangian.
Before proceeding, we must check what is the regime of parameters for which the LO mass relation between fermion mass and pseudoscalar mass is a good approximation to the data. To illustrate this point, in the top panel of Fig. 6 we present the pseudoscalar masses various bare masses. 4 The critical massm c 0 is determined numerically, extrapolating from the linear fit to the lightest five data points to the value for whichm 2 PS = 0. As shown in the figure, deviation from linearity appear form 2 PS > 0.4. The decay constant squared also shows a linear behaviour, and its slope is not negligible. The bottom panel of Fig. 6 shows deviations from linearity of the dependence of the combination m 2 PS f 2 PS on the fermion mass. The Gell-Mann-Oakes-Renner (GMOR) relation [105], m 2 PS f 2 PS = m f ψ ψ , would imply a mass dependence of the fermion condensate over the range of mass considered. We alert the reader that a rigorous discussion of the GMOR relation would require first to determine the values ofm bare mass and coupling, while keeping the lattice spacing in units of w 0 fixed), while in this simplified discussion we kept the lattice coupling β = 7.2 fixed. This is adequate for the purposes of this subsection, but we refer the reader to Section 6.1 for an assessment of the validity of the GMOR relation in the continuum limit.
In Fig. 7 we show the numerical results of the decay constants squared of PS, V and AV mesons, while in Figs. 8 and 9 we present the masses squared of mesons sourced by the operators V and T, and by S, AV, and AT, respectively, as functions ofm 2 PS . Our first observation is that discretisation effects inf 2 PS andm 2 V (orm 2 T ) are significant, given the visible difference between data collected at different lattice couplings (and denoted by different colours). For other quantities the deviations are no larger than the statistical uncertainties. Also, the masses and the decay constants decrease as we approach the massless limit, with the exception off 2 AV . Overall, the masses and decay constants show linear

Continuum extrapolation
We are now in a position to perform the continuum extrapolation, and to eliminate discretisation artefacts in the meson masses and decay constants. In order to do so, we introduce the important tool of WχPT [83,84] (see also Ref. [85], and [86,87]). It extends the continuum effective field theory by a double expansion, both in small fermion mass m, as well as lattice spacing a, as both of them break chiral symmetry and can be introduced in the EFT as spurions. We denote the lattice spacing in units of w 0 bŷ a ≡ a/w 0 = 1/w lat 0 . wheref χ = f χ w χ 0 is the pseudoscalar decay constant in the massless and continuum limit. The fermion masses used in this study are comparatively large, and hence it is legitimate to omit from Eq. (4.10) the chiral logs, which are important for small values ofm PS . The coefficientsb χ f andŴ χ f control the size of corrections due to finite mass and finite lattice spacing, respectively. In principle one should measure all observables in units of w χ 0 , while we instead use the mass-dependent w 0 , as measured at the finite mass of the individual ensembles, hence avoiding the need to extrapolate to the massless limit [82]-we collected enough data to attempt such extrapolation only for two values of β. The replacement of w χ 0 by w 0 does not affect the NLO EFT, the difference appearing at higher orders in m 2 PS . Compared to the continuum NLO expression in Ref. [106], it results in a shiftb χ f byk 1 in Eq. (3.8), due to fitting the measurements of f PS w 0 .
Among the underlying assumptions of WχPT is the requirement that the measurements it describes satisfy where Λ χ is the symmetry breaking scale, for which we take the estimate Λ χ = 4πf PS . In Section 3.2, we found the NLO EFT to describe the measurements ofŵ 0 up tom 2 PS < ∼ 0.4. The numerical results for the pseudoscalar decay constant squared in Fig. 7  The first estimate is roughly consistent with the scale separation in Eq. (4.11). But we find that aΛ χ evaluated on some of the ensembles is greater than unity. We further constrain the  analysis to the ensembles that satisfy the conditionâ 1, which is satisfied by DB1M5−7, DB2M1−3, DB3M5−8, and DB4M2. Only these ensembles will be used in the continuum and massless extrapolations that employ the tree-level NLO WχPT. We will report the values of the χ 2 /N d.o.f. in the analysis, and we anticipate here that the quality of the fits of the data supports the results of this exclusion process, otherwise based upon a set of estimates.
As    (4.14), and by using the fit parameters summarised in Table 7.
In the table, the numbers in the two parentheses denote the statistical and systematic uncertainties associated with the numerical fits. The latter is estimated by varying the fitting range to include or exclude the coarsest or the heaviest ensemble. We find acceptable values of χ 2 /d.o.f for all the fits, in support of the applicability of the tree-level NLO EFT to describe our data.

Lattice results: summary
We provide here a list of lattice highlights (from Sections 2, 3, and 4). We will be very schematic: precise definitions and detailed discussions can be found in the main text.
L1. The lattice ensembles used in the dynamical study, and the values of the lattice parameters that characterise them, are discussed in Section 2, and listed in Table 1, including the average plaquette and the gradient flow scale w 0 .
L2. The gradient flow scale w 0 is defined and studied in Section 3.1. Its mass dependence is illustrated by Fig. 1, 2, and 3 and Table 2. Throughout this work we employ the gradient flow scale to set the physical scale of the lattice calculations. L4. Lattice finite size effects are studied in Section 3.4. They are negligibly small provided the condition m PS L > ∼ 7.5 is satisfied-see also Table 3 and Fig. 5. In the analysis presented in this paper, all ensembles satisfy this condition.
L5. In Section 4, we define the operators and the two-point functions used for the measurement of decay constants and masses of the mesons-see Table 4. The decay constants are perturbatively renormalised according to Eqs. (4.6) and (4.7).
L6. The masses and (renormalised) decay constants (in units of the lattice spacing a) of the spin-0 mesons are reported in Table 5, and those of the spin-1 mesons in Table 6.
L7. In Sec. 4.2 we introduce the notationm ≡ m M w 0 , and equivalent for all dimensional quantities expressed in terms of the gradient flow scale. We discuss the GMOR relation and illustrate it in Figure 6. Linearity with the fermion mass holds for m 2 PS < ∼ 0.4. From this point on, we eliminate the dependence on the fermion mass, and study our observables only as a function of the mass (squared) of the PS state.
L8. Masses and decay constants, expressed in units of the gradient flow scale, are plotted for all ensembles in Figs. 7, 8, and 9. We notice that the V and T masses are compatible with each other, and so are the AV and AT-although the latter two states are physically distinct. We do not further discuss the T and AT states in the paper.
L9. We apply tree-level WχPT via Eqs. (4.13) and (4.14), that we use to perform the continuum and massless extrapolations. We impose the conditionsâ < ∼ 1, m PS L > ∼ 7.5, and m 2 PS < ∼ 0.4, which identify a subset of eleven ensembles to be used for the continuum extrapolation off 2 PS . For the other observables we relax the condition on the mass to m 2 PS < ∼ 0.6. The results are shown in Table 7 and in Figs. 10 and 11.
L10. As reflected in the fit results for W 0 presented in Table 7,f 2 PS andm 2 V (and similarlŷ m 2 T ) are affected by discretisation effects, the sizes of which are about 4 ∼ 13% and 10 ∼ 35%, respectively, over the considered ranges of bare parameters. In all other observables, discretisation effects are no larger than the statistical errors.

Low-energy phenomenology: EFT and sum rules
In this section, we study some of the long distance properties of the Sp(4) gauge theory that may be useful inputs for phenomenological studies performed in the context of composite Higgs models based upon the SU (4)/Sp(4) coset. We restrict our attention to the PS, V and AV states, and consider only ensembles withâ < ∼ 1 andm 2 PS < ∼ 0.4. We perform a simplified continuum extrapolation, by subtracting from their masses and decay constants the last term of Eqs. (4.13) and (4.14), with the coefficients W 0 f,M and W 0 m,M obtained by a fit similar to the one presented in Sec. 4.3, but restricted to include only the eleven ensembles that satisfy all the constraints. We then use the resulting eleven independent measurements of the three masses and three decay constants, that we report in Table 8, to perform a simplified global fit based upon the NLO EFT in Ref. [59]. The results of the fit allow us to provide an extrapolation to the massless limit, and to assess the validity  of the EFT itself by providing a first estimate of the size of the g V P P coupling between a vector state V and two pseudoscalar states PS. We also discuss several phenomenological quantities of general interest, connected with classical current algebra results, such as the aforementioned GMOR relation [105], and the saturation of the Weinberg sum rules [107] by the lightest spin-1 states.

Global fit and low-energy constants
The EFT presented in Ref. [59] describes the low-energy behaviour of the lightest PS, V, and AV mesons, by adopting the ideas of HLS [71][72][73][74][75]. The extension of the chiral Lagrangian is achieved by enhancing the SU (4) global symmetry to SU (4) A × SU (4) B , weakly gauging the SU (4) A group factor, and implementing the spontaneous breaking to the global Sp(4) by the VEVs of two spin-0 fields subject to non-linear constraints, and transforming one on the bi-fundamental of SU (4) A × SU (4) B and the other on the antisymmetric of SU (4) A , respectively. In this way, besides the 5 pNGBs identified with the PS states of the theory, the spectrum also contains 15 massive spin-1 particles, identified with the 10 V and 5 AV states. Explicit breaking of the global symmetry due to the fermion mass m f is implemented by the familiar spurion analysis. The results for physical quantities in terms of the 12 free parameters (f , F , b, c, g V , κ, v, v 1 , v 2 , v 5 , y 3 , y 4 ) of the EFT are summarised in Eqs. (2.19)-(2.24), together with Eqs. (2.29), (2.30) and (2.33) of Ref. [59]. 5 The simplified global fit in Ref. [59] suffered from uncontrolled systematics such as quenching effects and discretisation artefacts. In this paper, we overcome these limitations  Figure 12: Decay constants squared in the continuum limit-subtracting lattice artefacts due to a finite lattice spacing (see Table 8)-as a function of the pseudoscalar massm 2 PS . The bars represent the statistical uncertainties. The global fit results are denoted by shaded bands whose widths indicate the statistical errors.
with the continuum extrapolations of the data obtained by dynamical simulations, as discussed in previous sections. A further technical difficulty of the global fit is due to the large parameter space and the limited number of observables available. Because we employ Wilson fermions, we have one more unknown fit parameter, the critical bare fermion mass m c 0 , if we use the bare mass m 0 for the fermion mass in the EFT. A better way to determine the parameters in the EFT might be to consider a more sophisticated definition of the fermion mass, such as the one calculated via partially conserved axial current (PCAC). However, this would require to carry out a more involved computation of the renomalisation factors associated with the pseudoscalar operators, which is beyond the scope of this work.
Along the lines discussed in Section 4.2, here we take a different approach. We eliminate from the analysis any direct reference to the fermion mass, effectively replacing its role with the mass of the PS state. Because the linear relation m 2 PS = 2Bm f holds only at leading-order, we restrict our attention to ensembles for whichm 2 PS 0.4 (see Table 8). Furthermore, we expand the dependence of the other observable masses and decay constants [59] Figure 13: Masses squared in the continuum limit-subtracting lattice artefacts due to a finite lattice spacing (see Table 8)-as a function of the pseudoscalar massm 2 PS . The bars represent the statistical uncertainties. The global fit results are denoted by shaded bands whose widths indicate the statistical errors.
Up to O(m 4 PS ), we find having make use of the redefinitionŝ This linearised ansatz involves 10 unknown parameters (f ,F , b, c, g V , κ,v 1 ,v 2 ,ŷ 3 ,ŷ 4 ) to be determined from 5 measurements. The remaining two parameters (v and v 5 ) are associated with the GMOR relation, as they enter at NLO the definition ofm 2 PS . We will return to this point later on.  Table 9: Coefficients appearing in Eqs. (4.13) and (4.14), determined by using the results of the global fit described in the main text. The results are compatible with those obtained with the alternative fit process in Table 7, except for the coefficient L 0 f,AV .
We perform a global fit using the functions in Eqs. (6.1)-(6.5). The χ 2 function is defined in a simplified manner, by summing the individual χ 2 obtained from the five independent fit equations, and hence ignoring correlations among the fit equations. Besides subtracting the discretisation effects W 0 m,Mâ and W 0 f,Mâ from the original data, and restricting consideration tom 2 PS 0.4 andâ 1, as anticipated, we also constrain the fit by implementing the following conditions, which are the consequences of unitarity [59]: In Figs. 12 and 13 we present the continuum values of the masses and decay constants from Table 8, along with the results of the global fits, represented by blue bands, obtained through a constrained bootstrapped χ 2 minimisation with 10 parameters, with the constraints guided by an initial minimisation of the full data set. The fit functions describe the data well, with χ 2 /N d.o.f ∼ 0.4. But we also find that some of the EFT parameters are not well determined (see discussion in Appendix B, in particular the histograms of 200 resampled data obtained by bootstrapping, for the LO parameters in Fig. 21 and for the NLO ones in Fig. 22), indicating the presence of flat directions in the χ 2 .
The fact that the EFT yields a good quality fit of the masses and decay constants justifies using its results to calculate other physical quantities. The coupling g VPP is responsible for the decay of the vector meson into two pseudoscalar mesons. The NLO EFT expression can be found in [59] (where g ρππ denotes g VPP ) and in the massless limit we have that the g χ VPP coupling is We find that g χ VPP = 6.0(4)(2)-see also Appendix B, in particular the histogram of the coupling in Fig. 23, which shows a nice gaussian distribution.
As a consistency check, we also calculate the relevant coefficients in Eqs. (4.13) and (4.14), by using the results of the global fit. We present the results in Table 9: they are consistent with the ones in Table 7, with the only exception of L 0 f,AV which is over-constrained by Eqs. (6.1)-(6.5), but also affected by larger systematic as well as statistical uncertainties.
The results of this exercise have to be taken with caution, for a number of reasons. The fermion masses considered are not small enough to fully justify the linearisation of the fit equations, nor the truncation of the EFT itself to exclude higher-order terms. This is also clear from the fact that the masses considered are such that the V mesons are stable, because the 2-PS channel is closed. These objections can be addressed in the future, by lowering the fermion mass studied on the lattice. Also, the coupling g V (or alternatively the coupling g VPP ) turns out not to be small, bring into question the truncation of the EFT. This coupling should be suppressed if we consider future lattice studies of Sp(2N ) theories with N ≥ 3.

GMOR relation and Weinberg sum rules
There are several consequences of the HLS EFT (and of current algebra) that can be tested on our lattice data. If we first restrict our attention to the PS sector, we have the GMOR relation, which at NLO takes the corrected form where v and v 5 are the two low-energy constants defined in Ref. [59] that we removed from our global fit. As noted in Section 4.2, we cannot extract v and v 5 without renormalisation of the fermion mass. In Fig. 14 we plot the numerical results ofm 2 PSf 2 PS with respect tô m 2 PS , as an illustration of Eq. (6.10). Going beyond the pseudoscalar sector, the first non-trivial prediction of the NLO EFT-related to some reasonable assumptions for the truncation of the series of operatorsis that the sum of the decay constants (at zero momentum) is expected to be independent on the fermion mass. In the top-left panel of Fig. 15 we plot the numerical results off 2 0 with respect tom 2 PS , with the continuum extrapolation shown in the top-right panel. As seen in the figure, at the level of the present precision the numerical results support the mass independence of f 2 0 over the range of mass considered. While the Weinberg spectral theorems, which involve integration over all momenta, are exact, because they reflect the properties of the underlying condensates in the theory [107], the saturation of the resulting sum rules over a finite number among the lightest spin-1 states is an approximation. The NLO EFT analysis captures such violations in the non-nearestneighbour interactions appearing in the HLS Lagrangian. These terms are generated by integrating out heavier resonant spin-1 states that contribute to the dispersion relations. Saturated over the lightest resonances, the first and second sum rules are  Table 8). and respectively. The aforementioned corrections give finite contributions to the right-hand-side of these two relations. The numerical results for these quantities, normalised by f 2 V and f 2 V m 2 V , respectively, are shown in the middle and bottom panels of Fig. 15. As the figures show, after continuum extrapolation the results are closer to saturating the Weinberg sum rules. The conservative estimate of the uncertainties includes both statistical and systematic error, and with its current size this comparison is somewhat inconclusive.

Comparison to other gauge theories
This section is devoted to comparing our results with those obtained in a few other lattice gauge theories with two Dirac fermions. We start by presenting our results in a different way, by normalising all masses and decay constants to the decay constant f PS of the lightest PS state. By doing so, we remove direct reference to the gradient flow scale w 0 , and make the comparison to other theories more transparent. We report the results of this exercise in Table 10 and Fig. 16. We restrict attention to the eleven ensembles for which we are able to take the continuum limit for all the observables, as described in Section 4.3. We then compare our results to those obtained in gauge theories with SU (2) [44], SU (3) [108], and SU (4) [82] with matter field content consisting of N f = 2 Dirac fermions in the fundamental. We do so by testing the validity of the two phenomenological Kawarabayashi-Suzuki-Riazuddin-Fayyazuddin (KSRF) relations [109,110]. We conclude this section with a comparison to the continuum extrapolation of the quenched Sp(4) theory; details of additional quenched computations performed beyond those presented in Ref. [59] to facilitate this analysis are presented elsewhere [80].  Figure 15: Top to bottom:f 2 0 , the first, and the second Weinberg sum rule, saturated over the lightest spin-1 states, as a function ofm 2 PS . Left panels: we do not perform the continuum extrapolation, colour coding represents the bare lattice coupling of the various ensembles. Blue, purple, green, red and brown colours have β = 6.9, 7.05, 7.2, 7.4 and 7.5, respectively. Right panels: continuum-extrapolation results as in Table 8 (in blue), and the results of the continuum and massless extrapolation (in black). Statistical and systematic errors are summed in quadrature.
The fermion masses in the selected eleven ensembles are moderately heavy, with 1.39 < ∼ m V /m PS < ∼ 1.87, so that both V and AV mesons are stable, as the decay channels to two and three PS mesons, respectively, are kinematically forbidden. Bearing this in mind, we nevertheless apply Eqs. (4.13) and (4.14) to the eleven ensembles, and we perform the massless (and continuum) extrapolations. In the last row of the Table 10 we report the resulting ratios of the masses and decay constants tof PS , which are rendered in black colour  : Ratios of meson masses and decay constants, extrapolated to the continuum limit by the subtraction method explained in Section 4.3. We restrict the data to the eleven ensembles for which the subtraction of effects due to the finite lattice spacing can be done for all measurable quantities. In the last row, we report the results of the continuum and massless extrapolation, obtained by applying Eqs. (4.13) and (4.14) to the eleven ensembles. Statistical and systematic uncertainties have been added in quadrature.
in Figure 16 . Notice that the ratio are independent of the gradient flow scale w 0 , so that f V /f PS =f V /f PS (and analogous for all dimensionless quantities in Table 10). The spectroscopy of the lightest meson states is captured to some approximation by the two KSRF phenomenological relations [109,110]. The first such relation states that which is rather close to real-world QCD, as taking the experimental value of the rho meson decay constant to be f ρ 148 MeV, yields f ρ /f π ∼ 1.6 [111]. We can compare our results for Sp(4) with the first KSRF relation by looking atf V /f PS in Figure 16. Form 2 PS < 0.4 the ratio monotonically increases fromf V /f PS ∼ 1.5 asm 2 PS decreases. A simple linear extrapolation yields the ratio in the massless and continuum limit to bef V /f PS ∼ 2.1. Therefore, our numerical results do not support the first KSRF relation: the resulting values not only depend on the fermion mass, but also become larger in the massless limit.
The second KSRF relation involves m V , f PS , and the on-shell coupling constant g VPP associated with the decay of a vector meson: In real-world QCD, the mass of ρ meson m ρ 775 MeV, expressed in units of the pion decay constant f π , yields roughly m ρ / √ 2f π ∼ 5.9 [111]. For comparison, we adopt the treelevel definition for the decay rate of ρ meson, Γ ρ =  Figure 16: Ratios of masses and decay constants, in the continuum limit, expressed in units of the decay constant of the PS meson (as in Table 10), shown in blue. We also show the extrapolation to the continuum and massless limit in black, except form PS , which vanishes in the limit of massless fermions.
experimental values Γ ρ 150 MeV and m π 140 MeV. We find g ρππ 6.0, which is in quite good agreement. By evaluating the right-hand side of the second KSRF relation for the Sp(4) gauge theory, computed with the lightest ensemble and in the massless limit, we findm V / √ 2f PS = 5.47 (11) andm V / √ 2f PS = 5.72(18)(13), respectively. By comparing with the independent measurement of g VPP = 6.0(4)(2) from Section 6, obtained from the global fit of the EFT, we conclude that the second KSRF relation holds.
It is interesting to compare the right-hand side of the second KSRF relation with other lattice results available in the literature on SU (N ) gauge theories with two fundamental Dirac fermions. We show the comparison in Fig. 17. For the lightest ensembles available, it is found that m V / √ 2f PS ∼ 11.4 (13) for SU (2) [44] and m V / √ 2f PS ∼ 5.2(3) for SU (4) [82], respectively. The general trend in SU (N ) theories is that the value of m V / √ 2f PS decreases as N increases, which complies with the expectation that g VPP decreases in the large-N limit. 6 Near the threshold of m PS /m V ∼ 0.5, the vector meson mass we find for Sp(4) in the continuum limit lies between the values for SU (3) and SU (4).
We close this section by comparing the dynamical results to the quenched ones. Experience from lattice QCD suggests that the quenched results, which require only a small amount of computing resources, capture well the qualitative features of the dynamical results for masses and decay constants, in spite of the fact that the associated systematic uncertainties are not under analytical control. Quenching effects are expected to be smaller as the number of colours increases with a fixed number of fundamental flavours. 6 It is also interesting to investigate the flavour dependence of the ratio mV/ √ 2fPS. A recent lattice study for SU (3) gauge theory coupled to N f fundamental fermions finds that the ratio is statistically independent on N f up to N f = 6, for which all theories considered are expected to behave in a way resembling QCD [112]. On the other hand, the ratio could depend on the group representation of the fermion matter fields. For instance, large-Nc arguments suggest that gVPP ∝ 1/ √ Nc and gVPP ∝ 1/Nc for single index and two-index fermion representations, respectively. Pioneering lattice results in SU (4) are consistent with this scaling [82]. In order to make the comparison possible, our preliminary exploration of quenched simulations [59] had to be extended by considering larger lattice volumes and various lattice couplings [80]. As in the dynamical case, the quenched ensembles satisfy the condition m PS L 7.5, that allows one to neglect finite volume effects. Continuum extrapolations use ensembles constructed for β = 7. 62, 7.7, 7.85, 8. the differences become more substantial as the fermion mass decreases. We estimate the discrepancies to be δf 2 PS /f 2 PS ∼ 20% and δm2 S /m 2 S ∼ 25%, in the massless limit. The mass of the V meson shows a somewhat milder discrepancy, at the level of ∼ 10%. For other quantities, quenching effects are not visible: the corresponding discrepancies are smaller than the uncertainties associated with the fits. Interestingly, the resulting values ofm V / √ 2f PS for the dynamical and quenched simulations, which may be used to estimate the coupling g VPP via the second KSRF relation, are found to be consistent with each other in the massless limit [80]. The general conclusion of the comparison with the quenched results is quite encouraging, although at present we do not know whether this conclusion is an indication that the quenched approximation adequately captures the information encoded in the two-point functions-possibly because of the proximity to large-N -or whether it is just a trivial consequence of the large fermion masses we studied.

Continuum results: summary
In this section, we briefly summarise the continuum extrapolation results for the dynamical theory, presented in Section 6 and 7. C1. Our continuum results for the decay constants and masses of PS, V, and AV states, for the ensembles satisfyingm 2 PS ≤ 0.4, are reported, in units of the gradient flow scale, in Table 8. C2. The global fit based on the EFT describing PS, V and AV states using hidden local symmetry yields the results in Table 9, illustrated in Figs. 12 and 13.
C3. Section 6.1 discusses the GMOR relation and three sum rules in the continuum limit.
The main results are shown in Figs. 14 and 15.
C4. In Section 7 the continuum limit results are discussed in units of the decay constant f PS of the PS states, that are summarised in Table 10, and illustrated in Fig. 16.
We include also the mass of the scalar flavoured state S. All our measurements are in the range 1.39 < ∼m V /m PS < ∼ 1.87, in which the V states cannot decay to states containing two PS particles. (Analogous considerations apply to the 3-body decay of the AV mesons.) This range may be of direct relevance in the context of dark matter models.
C5. We perform the extrapolation to the massless limit. The masses of V, AV, and S states are m V /f PS = 8.08 (32), m AV /f PS = 13.4(1.5), and m S /f PS = 14.2(1.7).
The decay constants of V and AV states in the continuum and massless limit are f V /f PS = 2.15 (8) and f AV /f PS = 2.3(4), respectively-see Table 10.
C6. We find that f V /f PS = 2.15 (8) is larger than expected on the basis of the first KSRF relation, which would yield f V = √ 2 f PS -see Section 7.
C8. We compare the continuum and massless limit results to the literature on theories with two Dirac fermions on the fundamental. Fig. 17 shows that the VPP coupling is smaller than in the SU (2) theory, but comparable to SU (3) and SU (4).
C9. We close Section 7 by comparing our results, obtained with dynamical fermions, to quenched calculations [80]. (See Figures 18 and 19.) We find that, in the massless limit, the decay constant squared of the PS statef 2 PS is ∼ 20% lower than the quenched result. The mass squared of the (flavoured) scalar is ∼ 25% lower than in the quenched result, while that of the vector is lower by ∼ 10%. In the other observables, the dynamical and quenched results are compatible with one another, given the current uncertainties.

Conclusions and outlook
Following along the programme outlined in Ref. [59], with this paper we have made a substantial step forwards in the study of the gauge theory with Sp(4) group and N f = 2 (Dirac) fundamental fermions, that at long distance yields pNGBs describing the SU (4)/Sp (4) coset, of relevance in the CHM context. We have adopted the Wilson-Dirac lattice action for gauge and fermion degrees of freedom, and performed numerical studies via the HMC method, with dynamical fermions. We have repeated the calculations at several values of the lattice bare parameters (fermion mass and gauge coupling). Our main result is the first continuum extrapolation of the Sp(4) measurements for the masses and decay constants of the lightest, flavoured, spin-0 and spin-1 mesons. While numerical studies are performed in a regime of fermion masses large enough to preclude the decay of heavier mesons to the pNGBs, we have presented also a preliminary massless extrapolation, and compared the results to those of the quenched approximation. Summaries of the lattice and continuum limit results can be found in the brief Sections 5 and 8, respectively. Other results of relevance to the broad research programme of lattice gauge theories with Sp(2N ) group will be reported in Refs. [80] and [113] (see also Ref. [63]).
In the CHM context, while it is not strictly necessary to have massless fermions, it would be important in the future to extend this study to reach closer to the massless regime, for pNGB masses that allow the heavier mesons to decay. This would allow to test whether the moderate departure of our results from the quenched ones is due to the large fermion mass in the action, to the approach to the large-N limit, to dynamical properties of the theory or to a combination of the above. Of great importance in the CHM context is also the problem of vacuum alignment. A first step towards addressing this point, without specifying the full set of couplings to SM fermions, would be to compute the coefficients of higher-order operators in the EFT obtained by retaining only the pNGB degrees of freedom. This information could be used as input in model building, to compute the (perturbative) Higgs potential responsible for EWSB. In general, it would also be interesting to track the dependence on N of the spectra of mesons in Sp(2N ) theories with the same fermion field content, and understand how closely it follows the expectations from the large-N analysis.
The mass regime we studied is suitable for direct application in the context of models with strongly-coupled origin of dark matter (along the lines of Refs. [64][65][66], for example). The information we provide here is a necessary first step towards adopting this theory as a candidate for the origin of dark matter. Further information about the dynamics feeds into the relevant cross-sections, the calculation of which requires knowing the interactions. The first coupling one would like to measure in the future is the g VPP coupling between one vector state V and two pseudoscalar states PS. We have performed a preliminary determination of this coupling based upon the EFT constructed to include also V and AV mesons, but the intrinsic limitations of this process imply that the results should be used with some caution.
Finally, CHM constructions often involve advocating for (partial) compositeness of some of the SM fermions, the top quark in particular. To this purpose, model building requires the existence of new composite fermion states with special properties, that can be realised in the Sp(4) gauge theory by introducing n f = 3 (Dirac) fermions transforming in the 2-index antisymmetric representation of the group [12,41]. Some results about the (quenched) spectrum of mesons in the theory with this field content will be presented in Ref. [80], but the implementation of fully dynamical calculations with multiple fermion representations will require a significant amount of development, with only few examples of calculations of this type existing in the current literature [27,31,36,82,114].

Acknowledgments
We acknowledge useful discussions with Axel Maas and Roman Zwicky. We would like to thank Michele Mesiti and Jarno Rantaharju for their assistance on the modification and improvement of the HiRep code for this project.
The work of EB has been funded in part by the Supercomputing Wales project, which is part-funded by the European Regional Development Fund (ERDF) via Welsh Government. The

A Lattice action and numerical calculations: additional details
We collect in this Appendix some supplementary details, pertaining the lattice numerical study that ultimately leads to the results summarised in Tables 5 and 6, which provide the data upon which we perform the continuum and massless extrapolations. In Appendix A.1 we discuss the treatment of the Sp(2N ) matrices within the HMC algorithm. We briefly clarify the role of diquark operators in Appendix A.2. Finally, in Appendix A.3 we present details of the fitting procedure of the correlation functions used to extract masses and decay constants.

A.1 Hybrid Monte Carlo
We perform numerical simulations using a variant of HiRep code [88], which is designed to simulate SU (N ) and SO(N ) lattice gauge theories with fermions in higher representations. While a detailed description of the implementations of Sp(2N ) theories in the HiRep code is described in Ref. [59] (see also Refs. [60] and [62]), here we briefly summarise its main features and then report on some improvements we implemented for the purposes of this and future projects. This paper focuses on the spectroscopy of Sp(4) with two fundamental fermions, higher-dimensional representations will be discussed elsewhere. 7 To study the N f = 2 theory we use the standard hybrid Monte Carlo algorithm, a well established technique for lattice QCD.
In our preliminary study of the two-flavor Sp(4) theory [59], we used the specific form of group generators T A given in Ref. [48] with the normalisation of Tr(T A T B ) = δ AB /2 for the molecular dynamics (MD) evolution. We also implemented a resymplectisation process, to ensure that updated configurations stay inside the Sp(4) group manifold, that requires performing a (normalised) group projection onto the quaternion basis. This process works well for Sp(4) theories, but in order to enhance the capability of the software in view of future studies of dynamical Sp(2N ) theories with arbitrary N , in this study we further improve and generalise the code in the following way.
First, we remind the reader that the group elements U of Sp(2N ) satisfy the condition U is a unitary matrix, and it can be written as U = exp(ia A T A ), where T A are Hermitian traceless 2N × 2N matrices, and a N real numbers. Combining Eq. (A.1) and unitarity, one can also write the matrix U in the block-diagonal form: where V and W are complex N × N matrices. Because of its adaptability to any N , we perform the resymplectisation via a variant of the modified Gram-Schmidt algorithm, that has been tested for the pure gauge model with Heath Bath algorithm [59,62]. The basic idea is that the (N +j)-th column of U can be obtained from j-th via col j+N = −Ωcol j , while the Gram-Schimt procedure is used to find the (j+1)-th column through the orthonormalisation with respect to the j-th and (N + j)-th. We also modified the code to save two N × N matrices (V and W ), instead of the full 2N × 2N unitary matrix U , reducing by half the size of each configuration.
The MD update in HMC makes explicit use of the generators of the group, not just of the group elements. We write the group generators T A as follows. Eq. (A.1) can be rewritten in terms of the group generator T A , as the condition which allows to write T A in block-diagonal matrix form as where of the two N × N matrices, X is Hermitian and Y is complex symmetric. We use the definition in Eq. (A.4), supplemented by the normalisation Tr (T A T B ) = δ AB /2, for the generators implemented into the HiRep code. We conclude by summarising some technical details about the algorithm used in the Sp(4) gauge theory with two Dirac fermions in the fundamental representation. Gauge configurations are generated using the HMC with a second order Omelyan integrator for the MD update. We use different lengths of MD time steps δτ g and δτ f for gauge and fermions actions, respectively, which are optimised to keep the acceptance rate of the Metropolis test, performed at the end of each HMC update, in the range of 75 − 85%. Thermalisation and autocorrelation lengths are determined by monitoring the average value of the plaquette.

A.2 Of diquarks
In the theory studied in this paper, mesons and diquarks combine together in the lowenergy spectrum, to form irreducible representations after the spontaneous breaking of the enlarged global SU (4) symmetry breaking to Sp(4). Hence, we do not calculate the diquark correlators, as they are identical to the corresponding meson correlators. A general discussion of both real and pseudoreal representations can be found in Ref. [102], yet we think it is useful to explicitly write the diquark operators and show the identity at the level of correlators by using our lattice action in Eq. (2.4). Such an analysis for SU(2) theory with two fundamental fermions can be found in Ref. [42].
The diquark operators are defined by where C is the charge conjugation operator satisfying γ µT C = −Cγ µ , Γ is a generic matrix with spinor indices, spinor indices are understood and summed over, and the anti-symmetric matrix Ω defined in Eq. (A.1) acts on the (understood) colour indices. Then, the diquark correlation function is

A.3 Effective mass and fitting procedure
In order to extract the lattice measurements of masses and decay constants, we fit the numerical data for the relevant two-point correlation functions. In all ensembles we produced,  Table 11: Fitting intervals of the Euclidean time I fit = [t i , t f ] used in the single-exponential fit of the measured correlators for PS, V and AV mesons. The number of configurations and the separation between the adjacent configurations are denoted by N configs and δ traj , respectively. We carry out a correlated fit using standard χ 2 minimisation. We show the resulting values of χ 2 /N d.o.f. , as a way to assess the quality of the fit itself.
the Euclidean time is large enough to reach a plateau region in the plot of m eff versus time t, in which the correlation functions are well approximated by a single exponential function, with its decay rate identified by the mass of mesons-the relevant 2-point function asymptotically behaves as C O M (t) ∝ e −m M t . We use standard χ 2 minimisation to this functional form to determine the best fit parameters (mass and decay constant, in units of the lattice spacing) and their statistical uncertainties.
This process requires choosing, for each ensemble and each meson, an interval I fit = [t i , t f ] of Euclidean time, and restricting the fit to data collected between its minimum value t i and its maximum time t f . We choose the fitting intervals as follows: we first look at the effective mass plots and identify the emergence of the plateau. Over the plateau region, we perform a single-exponential fit. We vary the choices of t i and t f , and ultimately identify  Table 12: Fitting intervals of the Euclidean time I fit = [t i , t f ] used in the single-exponential fit of the measured correlators for T, AT and S mesons, where the same N configs and δ traj in Table 11 are considered. We carry out a correlated fit using standard χ 2 minimisation. We show the resulting values of χ 2 /N d.o.f. , as a way to assess the quality of the fit itself.
the best fit range with the choice that yields the smallest value for χ 2 /N d.o.f , with the largest possible range t f − t i . We provide the relevant details of this process in Tables 11  and 12 We notice that while the effective masses retained in the fit extend to the maximum length of the temporal directions T max for PS and V mesons, those for AV and S mesons typically cease at t < T max due to severe numerical noise problems, which in practical terms reduce the fitting ranges. As a result, we expect a comparatively large systematic error  Figure 21: Histograms of the distribution of the six low-energy constants appearing in the EFT of Section 6.1 at the leading order, as determined from the global fit. associated with the choice of the fitting range for AV and S states (analogous arguments apply to the AT states).

B Low-energy constants and global fit
In this Appendix, we present the numerical results for the LECs in Eqs. (6.1)-(6.5) obtained from the simplified global fit to the data discussed in Section 6.1. As anticipated, we find it instructive to explicitly show the histograms associated with the LEC distributions. Figs. 21 and 22 report the histograms for the LECs appearing in the EFT at the leading and the next-to-the-leading order, respectively. As seen in the figures, some fit parameters do  Figure 22: Histograms of the distribution of the four low-energy constants appearing in the EFT of Section 6.1 at the next-to-the-leading order, as determined from the global fit.  Figure 23 shows the histogram for the g VPP coupling as defined and discussed in Eq. (6.9) in Section 6.1.