Symplectic lattice gauge theories in the grid framework: Approaching the conformal window

gauge theories coupled

Symplectic gauge theories coupled to matter fields lead to symmetry enhancement phenomena that have potential applications in such diverse contexts as composite Higgs, top partial compositeness, strongly interacting dark matter, and dilaton-Higgs models.These theories are also interesting on theoretical grounds, for example in reference to the approach to the large-N limit.A particularly compelling research aim is the determination of the extent of the conformal window in gauge theories with symplectic groups coupled to matter, for different groups and for field content consisting of fermions transforming in different representations.Such determination would have far-reaching implications, but requires overcoming huge technical challenges.
Numerical studies based on lattice field theory can provide the quantitative information necessary to this endeavour.We developed new software to implement symplectic groups in the Monte Carlo algorithms within the Grid framework.In this paper, we focus most of our attention on the Sp(4) lattice gauge theory coupled to four (Wilson-Dirac) fermions transforming in the 2-index antisymmetric representation, as a case study.We discuss an extensive catalogue of technical tests of the algorithms and present preliminary measurements to set the stage for future large-scale numerical investigations.We also include the scan of parameter space of all asymptotically free Sp(4) lattice gauge theories coupled to varying number of fermions transforming in the antisymmetric representation.

I. INTRODUCTION
Gauge theories with symplectic group, Sp(2N ), in four space-time dimensions have been proposed as the microscopic origin of several new physics models that stand out in the literature for their simplicity and elegance.We list some compelling examples later in this introduction.Accordingly, lattice field theory methods have been deployed to obtain numerically a first quantitative characterisation of the strongly coupled dynamics of such gauge theories [1][2][3][4][5][6][7][8][9][10][11][12][13][14][15][16][17][18][19].Different regions of lattice parameter space have been explored; by varying the rank of the group, N , the number, N f,as , and mass, m f,as , of (Dirac) fermions transforming in the fundamental (f) and 2-index antisymmetric (as) representation, one can tabulate the properties of these theories.And, after taking infinite volume and continuum limits, the results can be used by model builders, phenomenologists, and field theorists working on potential applications.
A prominent role in the recent literature is played by the theory with N = 2, N f = 2, and N as = 3.It gives rise, at low energies, to the effective field theory (EFT) entering the minimal Composite Higgs model (CHM) that is amenable to lattice studies [20], 1 and also realises top (partial) compositeness [85] (see also Refs.[86,87]).It hence provides an economical way of explaining the microscopic origin of the two heaviest particles in the standard model, the Higgs boson and the top quark, singling them out as portals to new physics.
A special open problem is that of the highly non-trivial determination of the extent of the conformal window in strongly coupled gauge theories with matter field content.It has both theoretical and phenomenological implications, of general interest to model-builders, phenomenologists, and field theorists alike.Particular attention has been so far paid to SU (N c ) theories, more than Sp(2N ) (with N > 1) ones.Let us pause and explain what the problem is, on general grounds.Robust perturbation-theory arguments show that if the number of matter fields is large enough-but not so much as to spoil asymptotic freedom-gauge theories can be realised in a conformal phase.This is the case when long distance physics is governed by a fixed point of the renormalisation group (RG) evolution [157,158], and the fixed point is described by a conformal field theory (CFT).It is reasonable to believe that such fixed points may exist also outside the regime of validity of perturbation theory, when the number of matter fields is smaller.What is the smallest number of fermions for which the theory still admits a fixed point, rather than confining in the infrared (IR), is an open question.While gaining some control over non-perturbative physics is possible in supersymmetric theories (see Ref. [159] and references therein), the non-supersymmetric ones are the subject of a rich and fascinating literature [160][161][162][163][164][165][166][167][168], part of which uses perturbative instruments and high-loop expansions [169][170][171][172][173][174][175][176][177][178][179][180][181], but there is no firm agreement on the results-we include a brief overview of work in this direction, in the body of the paper.
Knowledge of the extent of the conformal window also has relevant phenomenological implications.Various arguments suggest that at the lower edge of the conformal window, the anomalous dimensions of the CFT operators might be so large as to invalidate naive dimensional analysis (NDA) expectations for the scaling of observable quantities [161,182].And it has been speculated that this might affect even confining theories that live outside the conformal window, with applications to technicolor, CHMs, top (partial) compositeness, SIMP dark matter (e.g., see Refs.[24][25][26][183][184][185][186][187][188] and references therein).
Lattice studies of the extent of the conformal window have mostly focused on SU (N c ) groups, with fermion matter in various representations of the gauge group. 2 Closely related to these studies is the emergence, in SU (3) gauge theories with eight (Dirac) fermions transforming in the fundamental representation [242][243][244][245][246][247][248][249][250], or (Dirac) fermions transforming in the 2-index symmetric representation [251][252][253][254][255][256], of numerical evidence pointing to the existence of a light isosinglet scalar state, that is tempting to identify with the dilaton, the PNGB associated with dilatations.
The aforementioned lattice studies of symplectic theories, motivated by CHMs and SIMPs, can be carried out with comparatively modest resources, and using lattices of modest sizes, because they require exploring the intermediate mass range for the mesons in the theory.By contrast, the study of the deep-IR properties of Sp(2N ) gauge theories requires investigating the low mass regime of the fermions, for which one needs lattices and ensembles big enough to overcome potentially large finite size effects and long autocorrelation times.The supercomputing demands (both on hardware and software) of these calculations are such that a new dedicated set of instruments, and a long-term research strategy, is needed to make these investigations feasible.With this paper, we make the first, propaedeutic, technical steps on the path towards determining on the lattice the extent of the conformal window in theories with Sp(2N ) group, for N > 1.
To this end, we elected to build, test, and make publicly available new software [291], that supplements previous releases of the Grid library [292][293][294][295], by adding to it new functionality specifically designed to handle Sp(2N ) theories with matter fields in multiple representations.The resulting software takes advantage of all the features offered by the modularity and flexibility of Grid, in particular its ability to work both on CPU-as well as GPU-based architectures.We present two types of preliminary results relevant to this broader endeavour: technical tests of the algorithm and of the physics outcomes are supplemented by preliminary analyses, conducted on coarse lattices, of the parameter space of the lattice theory.The latter set the stage for future large-scale numerical studies, by identifying the regions of parameter space connected to continuum physics.The former are intended to validate the software, and test its performance for symplectic theories on machines with GPU architecture.Unless otherwise specified, we use the Sp(4) theory, coupled to N as = 4 Wilson-Dirac fermions transforming in the 2-index antisymmetric representation, as a case study.The lessons we learn from the results we report have general validity and applicability.
This paper is organised as follows.We start by defining the Sp(2N ) gauge theories of interest in Sect.II, both in the continuum and on the lattice.We also summarise briefly the current understanding of the extent of the conformal window in these theories.Section III discusses the software implementation of Sp(2N ) on Grid, and the basic tests we performed on the algorithm.In Sect.IV we concentrate on lattice theories in which the fermions do not contribute to the dynamics, focusing both on the Yang-Mills theory and the quenched approximation.New results about the bulk structure of all the Sp(4) theories coupled to (Wilson-Dirac) fermions transforming in the 2-index antisymmetric representation can be found in Sect.V, while Sect.VI discusses scale setting (Wilson flow) and topology.A brief summary and outlook concludes the paper, in Sect.VII.Additional technical details are relegated to the appendix.

II. GAUGE THEORIES WITH SYMPLECTIC GROUP
The Sp(2N ) continuum field theories of interest (with N > 1), written in Minkowski space with signature mostly '−', have the following Lagrangian density (we borrow notation and conventions from Ref. [4]): The fields They can be written in terms of the gauge field A µ ≡ A a µ t a , where t a are the generators of Sp(2N ), normalised so that Tr t a t b = 1 2 δ ab , to read as follows: where g is the gauge coupling.The field-strength tensor is given by where [• , •] is the commutator.The form of Eq. ( 1) makes it easy to show that the SU (N f ) L × SU (N f ) R and SU (N as ) L × SU (N as ) R global symmetries acting on the flavor indexes of Q i and Ψ k , respectively, are enhanced to SU (2N f ) and SU (2N as ).By rewriting Eq. ( 1) in terms of 2-component fermions [4,296] , ( C = −iτ 2 , τ 2 is the second Pauli matrix) we get the Lagrangian where the global symmetries are manifest where we defined σµ ≡ 1 2×2 , τ The mass terms break the symmetries to the maximal Sp(2N f ) and SO(2N as ) subgroups.Bilinear fermion condensates arise non-perturbatively, breaking the symmetries according to the same pattern, and hence one expects the presence of N f (2N f − 1) − 1 PNGBs in the (f) sector (for N f > 1), and N as (2N as + 1) − 1 in the (as) sector.
The main parameters governing the system are hence N , N f , and N as , and in most of the paper we refer to the theory with N = 2, N f = 0, and N as = 4 as a case study.The running coupling, g, obeys a renormalisation group equation (RGE) in which the beta function at the 1-loop order is scheme-independent, and is governed by the coefficient b 1 , which for a non-Abelian theory coupled to Dirac fermions can be written as and, specifically for Sp(2N ) groups, becomes The coefficients C 2 (G), C 2 (f), C 2 (as) are quadratic Casimir operators in the adjoint, fundamental and antisymmetric representations, while d G , d f , d as are the dimensions of these representations, respectively.We restrict attention to asymptotically free theories, for which b 1 is positive.For Sp(2N ) theories with N f = 0, this requirement sets the upper bound N as < 11(N +1) 4(N −1) , which for N = 2 yields N as < 33/4-perturbatively, as-type fermions make double the contribution of f-type ones, in Sp(4).The spectrum of mesons depends on the mass, m f,as , of the fermions, by varying which we can test which of the following three possible classes the theory falls into.
(i) The theory confines, similarly to Yang-Mills theories.One expects to find a gapped spectrum, and a set of PNGBs that become parametrically light in respect to other states, when m f,as → 0. The small mass and momentum regime is described by chiral perturbation theory (χPT) [297][298][299][300].
(ii) The theory is IR conformal.In this case, a gap arises only because of the presence of the mass terms, and would disappear into a continuum for m f,as → 0. The spectrum and spectral density exhibit scaling, in the form described for example in Refs.[195,[301][302][303][304][305][306]-see also Ref. [307].

A. The conformal window
The three possible classes of gauge theories described above are determined by whether the theory is, respectively, far outside, inside or just outside the boundary of the conformal window.The determination of the conformal window is tantamount to showing the existence of the IR fixed point at non-zero coupling so that the theory is interacting and IR conformal.We provide here some more detail and information about this challenging endeavour and what is known to date, starting from perturbative arguments.The coefficient of the (scheme-independent) 2-loop RG beta function, b 2 , which is found to be, for generic non-abelian gauge theories, and for Sp(2N ) groups reduces to When b 2 is negative, one finds that for a positive and sufficiently small value of b 1 , a perturbative IR fixed point at coupling α IR ≃ α BZ = −4πb 1 /b 2 ≪ 1 arises.This is referred to as a Banks-Zaks (BZ) fixed point [157,158].The upper bound of the conformal window therefore coincides with that of asymptotically free theories, given by b 1 = 0.The determination of the lower bound of the conformal window is hindered by the vicinity of the strong coupling regime.To see this, one can fix the value of N and decrease the number of flavors N f,as .The coefficient b 2 then becomes less negative and eventually approaches zero, while b 1 remains finite and positive.Accordingly, the coupling at the (perturbative) BZ fixed point, α BZ , becomes larger and larger and the perturbative analysis of the β function is no longer reliable.Despite such inherent limitations, several (approximate) analytical methods have been proposed to estimate the critical value N cr f,as corresponding to the lower edge of the conformal window.We now briefly summarise known results, for the theories of interests, that can be used to guide dedicated studies using non-perturbative numerical techniques, such as those based on lattice field theory.
FIG. 1: Estimates of the extent of the conformal window in Sp(4) theories coupled to N f Dirac fermions transforming in the fundamental and Nas in the 2-index antisymmetric representation.The black solid line denotes the upper bound of the conformal window, while different colored and shaped lines denote alternative analytical estimates of the lower bound, obtained with different approximations.The dashed line is obtained by imposing the constraint b2(N f , Nas) = 0.The dot-dashed line is the result of the all-order beta function with the assumption that the anomalous dimensions of the fermion bilinears are γ ΨΨ = γ QQ = 1.The dotted line is the result of the SD analysis.The BZ expansion leads to the lower (blue) solid line.Details about these approximations can be found in the main text and in the reference list.
Let us start by setting N f = 0 and varying N as .A naíve estimate can be derived by taking the perturbative 2-loop beta function to hold beyond perturbation theory, using it to compute N BZ, cr as , and assuming that the fixed point disappears when α BZ → ∞, or equivalently by looking for solutions of the condition b 2 → 0. Doing so yields N BZ, cr as ≃ 3.7 for Sp(4).This approach can be systematically improved by including higher-order loops, up to ℓ max > 2, in the expansion of the beta function β(α).One then seeks values of N as for which α IR → ∞, with α IR determined by solving β(α) ≡ −2α ℓmax ℓ=1 b ℓ α 4π ℓ = 0.In particular, one finds N4−loop, cr as ≃ 4.1 from the perturbative beta function at four loops in the MS-scheme [311].It should be noted, however, that the results are affected by uncontrolled systematics, since the coefficients, b ℓ , of the beta function, β(α), depend on the renormalisation scheme at three or higher loops, when ℓ ≥ 3.
An alternative approach makes use of the Schwinger-Dyson (SD) equation in the ladder approximation, in which case conformality is assumed to be lost when α IR ≡ α cr , with α cr = π/3C 2 (R), which yields N SD as ≃ 6 for Sp(4).Going beyond the perturbative coupling expansion, a conjectured all-orders beta function β all−orders (α) [164], which involves the first two universal coefficients of β(α) and the anomalous dimension of fermion bilinear operator, γ ψψ (α), has been proposed. 3In this case, the conformal window is determined by solving the condition β all-orders = 0 with the physical input for the value of γ ψψ at the IR fixed point.For γ ψψ = 1, one finds N all−orders, cr as ≃ 5.5 for Sp(4). 4 More recently, the scheme-independent BZ expansion in the small parameter ∆ Nas = N AF as −N IR as has been extensively applied to the determination of physical quantities such as the anomalous dimension, γ ψψ , at the IR fixed point-see Ref. [176] and refs.therein.In Ref. [167], the authors determined the lower edge of the conformal window by imposing the critical condition of γ ψψ (2 − γ ψψ ) = 1.This condition is identical to γ ψψ = 1 at infinite order, but displays better convergence at finite order in the ∆ Nas expansion.The 4th order calculation yields N γcc, cr as ≃ 5.5 for Sp(4) [168].These analytical approaches can be extended to determine the conformal window for the theory containing fermions in the multiple representations, {R 1 , R 2 , • • • , R k }, in which case the upper and lower bounds of the conformal window are described by (k − 1)-dimensional hyper-surfaces.For the Sp(4) theories of interest with N f Dirac fermions transforming in the fundamental and N as in the 2-index antisymmetric representation, the results are summarised in Fig. 1. 5 The upper bound is determined by the condition b 1 (N f , N as ) = 0.The various alternative determinations of the lower bound are estimated as follows.The dashed line is obtained by setting b 2 (N f , N as ) = 0.The dot-dashed line corresponds to the result of the all-order beta function with the input of γ ΨΨ = γ QQ = 1.The dotted and solid lines are the results of the SD analysis and the BZ expansion of γ ΨΨ at the 3rd order in ∆ N f (nas) [179] with the critical conditions applied to the antisymmetric fermions, α BZ = α cr as = π/3C 2 (AS) and γ ΨΨ (2 − γ ΨΨ ) = 1, respectively, as fermions in the higher representation are expected to condense first, resulting in the larger values of α cr and γ IR [313].It might be possible to make use of the five-loop computations in Refs.[169,170], to further improve these estimations of the conformal window, but this goes beyond the purposes of this discussion.For the purpose of phenomenological applications, the most interesting physical quantities one would like to determine within the conformal window are the anomalous dimensions of fermion bilinear operators (mesons) and chimera baryon operators.Perturbative calculations of the former are available in the literature, up to the 4th order of the coupling expansion [314] and at the 3rd order of the BZ expansion [179], while that of the latter is only available at the leading order in α [62].All of these considerations, summarised in Fig. 1, offer some intuitive guidance for what can be expected, but non-perturbative instruments are needed to test these predictions and put Fig 1 on firmer grounds.

B. The lattice theory
In presenting the lattice theory, we borrow again notation and conventions from Ref. [9].The theory is defined on a Euclidean, hypercubic, four-dimensional lattice with spacing a, with L/a sites in the space directions and T /a in the time direction.The generic lattice site is denoted as x, and the link in direction µ as (x, µ).The total number of sites is thus Ṽ /a 4 = T × L 3 /a 4 .Unless stated otherwise, in the following we set L = T .The action is the sum of two terms where S g and S f are the gauge and fermion action, respectively.Among the several choices for the former-the Iwasaki, Symanzik, DBW2 and Wilson gauge actions-for simplicity, we show our results using the Wilson action, defined as where is known as the elementary plaquette operator, U µ (x) ∈ Sp(2N ) is the link variable defined on link (x, µ), and β ≡ 4N/g 2 0 , where g 0 is the bare gauge coupling.For the fermions, we adopt the massive Wilson-Dirac action, where Q j and Ψ j are the fermion fields transforming, respectively, in the fundamental and 2-index antisymmetric representation and j is a flavor index, while color and spinor indices are omitted for simplicity.The massive Wilson-Dirac operators in Eq. ( 15) are defined as where m f 0 and m as 0 are the bare fermion masses in the fundamental and 2-index antisymmetric representation, and The link variables U (as) µ (x) are defined as in Ref. [9], as follows: where e (ab) are the elements of an orthonormal basis in the (N (2N −1)−1)-dimensional space of 2N ×2N antisymmetric and Ω-traceless matrices, and the multi-indices (ab) run over the values 1 ≤ a < b ≤ 2N .The entry ij of each element of the basis is defined as follows.
while for b = N + a and 2 ≤ a ≤ N , e i,i+N = −e (ab) , for i < a , , It is easy to verify that each element of this basis satisfies the Ω-traceless condition Tr(e (ab) Ω) = 0, where the symplectic matrix Ω is defined in Eq. (A1).Finally, we impose periodic boundary conditions on the lattice for the link variables, while for the fermions we impose periodic boundary conditions along the space-like directions, and anti-periodic boundary conditions along the time-like direction.

III. NUMERICAL IMPLEMENTATION: GRID
Our numerical studies are performed using Grid [292][293][294]: a high level, architecture-independent, C++ software library for lattice gauge theories.The portability of its single source-code across the many architectures that characterise the exascale platform landscape makes it an ideal tool for a long-term computational strategy.Grid has already been used to study theories based on SU (N c ) gauge groups with N c ≥ 3, and fermions in multiple representations [315,316].In this section, we describe the changes that have been implemented in Grid in order to enable the sampling of Sp(2N ) gauge field configurations.With the aim of including dynamical fermions in future explorations of Sp(2N ) gauge theories, we focused our efforts6 on the Hybrid Monte Carlo (HMC) algorithm and on its variation, the Rational HMC (RHMC), used whenever the number of fermion species is odd.
The (R)HMC algorithms generate a Markov chain of gauge configurations distributed as required by the lattice action described in Sect.II B. The ideas underpinning these two algorithms can be summarized as follows.Firstly, bosonic degrees of freedom ϕ and ϕ † , known as pseudofermions, are introduced replacing a generic number n f of fermions.Powers of the determinant of the hermitian Dirac operator, Q R m = γ 5 D R m , in representation R can then be expressed as where flavor and color indices of ϕ and ϕ † have been suppressed for simplicity.For odd values of n f , the rational approximation is used to compute odd powers of the determinant above, resulting in the RHMC.Second, a fictitious classical system is defined, with canonical coordinates given by the elementary links and Liealgebra-valued conjugate momenta π(x, µ) = π a (x, µ) t a , where t a are the generators of the sp(2N ) algebra in the fundamental representation.The fictitious hamiltonian is where H g = S g and H f = S f .The molecular dynamics (MD) evolution in fictitious time τ is dictated by where F (x, µ), known as the HMC force, is defined on the Lie algebra sp(2N ), and can be expressed as The detailed form for F g (x, µ) and F f (x, µ) [191], the gauge and fermion force can be found in Section IIIA of Ref. [191].Numerical integration of the MD equations thus leads to a new configuration of the gauge field, which is then accepted or rejected according to a Metropolis test.The update process can hence be described as follows: • pseudofermions distributed according to the integrand in Eq. ( 21) are generated with the Heat Bath algorithm, • starting with Gaussian random conjugate momenta, the MD equations in Eqs. ( 23) are integrated numerically, • the resulting gauge configuration is accepted or rejected by a Metropolis test.
In this section we provide details on the implementation of the operations listed above, focusing on the alterations made to the pre-existing structure of the code designed for SU (N c ) gauge theories.We then describe and carry out three types of technical checks, following Ref.[191].We test the behaviour of the HMC and RHMC algorithms.We produce illustrative examples of the behaviour of the molecular dynamics (MD).Finally, we carry out a comparison between HMC and RHMC algorithms.The purpose of these tests is to verify that the dynamics is implemented correctly.

A. Software development
As in the case for the pre-existing routines handling the theories with gauge group SU (N c ), our implementation of Sp(2N ) allows for a generic number of colors.The starting point of the MD is the generation of random Lie-algebravalued conjugate momenta.The generators of the sp(2N ) Lie Algebra in the fundamental representation, as they appear in Grid, are provided by the relations described in Appendix B, where conventions for their normalisation are also established.Generators in higher representations of the gauge group can be derived from the fundamental ones [191,315].In particular, the generators of the algebra of Sp(2N ) in the antisymmetric representation can be obtained from the definition in Eq. (18), by Taylor expanding to first order around the unit transformation, (t a as ) (ab)(cd) = Tr e (ab) T t a f e (cd) + e (ab) T e (cd) t a T f .
In the numerical integration of Eq. ( 23), it is required to project the HMC force on the Lie algebra of the gauge group.In Grid, the embedding of the force-projection within the integrator requires the forces to be anti-hermitian.Hence, a projection operation to the matrices of the algebra sp(2N ) must be defined.This can be done in analogy with the projection to su(N c ), defined for a generic matrix M as where P tr M ≡ M −1 Nc Tr(M )/N c and P aH M ≡ (M −M † )/2 are the projectors to its traceless and to its anti-hermitian parts, respectively.For sp(2N ), the projection is instead defined as, where Notice that P − Sp returns an anti-hermitian matrix, while P + Sp projects on a space of hermitian matrices.
The resympleticisation of gauge links to the Sp(2N ) group manifold has also been implemented in Grid.The algorithm [1] is a modification of the Gram-Schmidt process designed to take into account the condition in Eq. (A3).After normalising the first column of the matrix U , the (N + 1)-th column is set to The second column is then obtained by orthonormalisation with respect to both the first and the N + 1-th column.An iteration of this process leads to a Sp(2N ) matrix.This procedure, performed after every update, prevents the gauge fields from drifting away from the group manifold due to the finite precision of the simulation.

B. Basic tests of the algorithm
In this subsection, we follow closely Sects.III and IV of Ref. [191].Our MD evolution is implemented using a second-order Omelyan integrator [319].However, in this work, the inversion of the fermion matrix is treated without preconditioning [9,320].
We now restrict attention to the theory with N = 2, N f = 0, and N as = 4, and perform a set of preliminary checks on the algorithms we use.We present the results in Figs. 2, 3, 4, 5, and 6, obtained, for convenience, setting the lattice parameters to β = 6.8, and am 0 = −0.6, on an isotropic lattice with volume Ṽ = (8a) 4 .
The first test pertains to Creutz equality [321]: by measuring the difference in Hamiltonian, ∆H, evaluated before and after the MD evolution, one should find that This is supported by our numerical results: Fig. 2 shows the value of exp − ∆H for five different choices of the time-step used in the MD integration, with ∆τ = τ /n steps , and the choice τ = 1.The numerical results are obtained by considering a thermalised ensemble consisting of 3400 trajectories, that we find has integrated auto-correlation time τ c = 6.1(2), measured using the Madras-Sokal windowing process [322].A closely related test is shown in Fig. 3: the value of the ensemble average of the plaquette is independent of ∆τ .
A third test pertains to the dependence of ⟨∆H⟩ on ∆τ , which for a second-order integrator is supposed to scale as ⟨∆H⟩ ∝ (∆τ ) 4 [323].In Fig. 4 we show our measurements, together with the result of a best-fit to the curve log⟨∆H⟩ = K 1 log(∆τ ) + K 2 , with K 1 = 3.6(4) determined by minimising a simple χ 2 .We find good agreement, as quantified by the value of the reduced χ 2 /N d.o.f.= 0.6, and K 1 is compatible to 4. A closely related test is displayed in Fig. 5, confirming the prediction that the acceptance probability of the algorithm, P acc , obeys the relation [324]: The final test of this subsection is displayed in Fig. 6.Following Refs.[191,325], we also want to ensure that the reversibility of our updates is respected.Reversibility is one of the fundamental properties required in order to pursue a correct HMC update.Our update algorithm, based on leapfrog, is reversible analytically.Yet, when using this algorithm numerically on computers, because of the finite precision, exact reversibility is lost.It is therefore  important to verify that implementation of the fundamental steps of the algorithm can be considered as reversible to good approximation, in order to avoid that rounding errors introduce a significant bias in our calculations.One can show that the quantity |δH|-the average difference of the Hamiltonian evaluated by evolving the MD forward and backward and flipping the momenta at τ = 1-doesn't change significantly in our simulations.Since the Hamiltonian in these tests is of order ∼ 10 6 and the typical δH ∼ 10 −11 , the results show that the violation of reversibility is consistent with having |δH|/H of the order of the numerical accuracy.This is the expected relative precision for double-precision floating-point numbers.Moreover, the violation |δH| is independent of ∆τ .As a "microscopic" and related effect, reversibility violations may occur while updating the gauge link variables and momenta updates during the MD evolution.To ensure this doesn't occur we update the gauge links through the exponentiation of the momenta, so that U (π) U (−π) = 1.Moreover, thanks to the double-precision nature of the variables we use, the entity of relative violations for momenta results to be within machine-accuracy in our simulations, as in the global reversibility violation case.

C. More about the Molecular Dynamics
For illustration purposes, we find it useful to monitor the contribution to the MD of the fields, and how this changes as we dial the lattice parameters.We focus on the theory with N = 2, N f = 0, and N as = 4, and consider a few ensembles with isotropic lattice with Ṽ = (8a) 4 , and lattice coupling β = 6.8, but vary the mass am as 0 .We show in Fig. 7 the force, F , as defined in Eq. ( 23), split in its contribution from the gauge and fermion dynamics, the latter computed using the HMC for all fermions.The results are normalised so that the gauge contribution is held constant.As can be clearly appreciated, for large and positive values of am as 0 the fermions can be neglected, as for these choices of the mass, one expects to be in the quenched regime.When decreasing the mass, the fermion contribution increases.For large, negative values of the Wilson bare mass (close to the chiral limit), the fermion contribution is even larger than the contribution of the gauge part of the action.

D. Comparing HMC and RHMC
While in this paper we are mostly interested in the theory with N = 2, N f = 0, and N as = 4, and hence we can use the HMC algorithm, for the general purpose of identifying the extent of the conformal window in this class of lattice gauge theories it may be necessary to consider also odd numbers of fermions, for which we resort to the RHMC algorithm.The latter relies on a rational approximation in the computation of the fermion force, but the presence of a Metropolis accept-reject step ensures that the algorithm is exact.Thus, a preliminary test must be made to check the consistency of the implementation-as was done for SU (3) theories, see for instance Ref. [318]. 7To gauge whether the numerical implementation is working at the desired level of accuracy and precision, we performed the exercise leading to Fig. 8.We computed the average plaquette, ⟨P ⟩, where P is defined as for ensembles having lattice volume Ṽ = (8a) 4 and coupling β = 6.8, for a few representative choices of the bare mass −1.4 ≤ am as 0 ≤ 0.0.We repeated this exercise three times: at first, we treated all fermions with the HMC, then we treated them all with the RHMC, and finally we used a mixed strategy, treating two fermions with the HMC, and two with the RHMC.We display, in the two plots in the figure, the differences of the second and third approaches to the first one, respectively.For most of the data, the differences are compatible with zero within the statistical uncertainties.More generally, given the number of observables, the probability to find a deviation larger than 3σ from the null value results to be ∼ 12%.

IV. THE N = 2 LATTICE YANG-MILLS THEORY
In this section, we start to analyse the physics of the Sp(4) theory of interest.We begin from the pure Yang-Mills dynamics, with N f = 0 = N as .We verify that centre symmetry, (Z 2 ) 4 , is broken at small volumes, but restored at large volumes, by looking at the (real) Polyakov loop, in a way that is reminiscent of Ref. [218].Following Ref. [315], we then consider the spectrum of the Dirac operator in the quenched approximation, both for fundamental and 2-index antisymmetric fermions, to verify the symmetry-breaking pattern expected from random matrix theory.
The results for the first of these tests are shown in Fig. 9.At a coupling β = 9.0, we generate four ensembles in the pure Sp(4) theory, at different values of the space-time volume, Ṽ = (2a) 4 , (4a) 4 , (12a) 4 , (20a) 4 .For each configuration, we compute the spatial averaged (real) Polyakov loop, defined as where U 0 (t, ⃗ x) is the time-like link variable.For our current purposes, we choose the lattice to be isotropic in all four directions, N t = N s = L/a.For each ensemble, we display the frequency histogram of the values of Φ.The expectation is that the zero-temperature Sp(4) lattice theory should preserve the (Z 2 ) 4 symmetry of the centre of the group in four Euclidean space-time dimensions.This is indeed the case for sufficiently large volumes, as shown by the bottom-right panel of Fig. 9, for which N t = N s = 20, that exhibits a Gaussian distribution centred at the origin.But for small enough lattice volumes, this expectation is violated.This is visible in the other three panels in Fig. 9, in which the distribution is non-Gaussian, and two other peaks emerge.In the extreme case of N s = N t = 2, the two peaks at finite value of Φ dominate the distribution, which is otherwise symmetrical around zero.Interestingly, Polyakov loops can be used also to perform the more physical study of the finite-temperature confinement/deconfinement phase transition.
In this case, one would consider N t ̸ = N s , vary the coupling β to identify the transition temperature, and then perform continuum and infinite volume extrapolations.The characterisation of the deconfinement phase transition is of interest for both theoretical as well as phenomenological reasons-see Ref. [156] and the review in Ref. [19]-but requires a dedicated, extensive programme of numerical work, and possibly cutting-edge new technology to address some of the difficulties faced by conventional Monte Carlo sampling methods (see, e.g., the discussions in Refs.[326] and [327]), while the simpler analysis performed here suffices for the more modest purposes of this paper.
Ensembles of gauge configurations without dynamical fermions can also be used to verify that our implementation of the Dirac operators is correct.To this purpose, following Ref.[315] (and Ref. [9]), we consider the theory with quenched fermions in either the fundamental or 2-index antisymmetric representation, and compute the spectrum of eigenvalues of the hermitian Wilson-Dirac operator Q m = γ 5 D m .The numbers of configurations are N conf,f = 88 and N conf,as = 47, while the number of eigenvalues in each configuration used is 3696 for fundamental fermions and 5120 for antisymmetric fermions.Then, we compute the distribution of the folded density of spacing, P (s).Finally, we compare the results to the exact predictions of chiral Random Matrix Theory (chRMT) [328,329].Because the spectrum captures the properties of the theory, in particular the pattern of chiral symmetry breaking [330], the distribution P (s) differs, depending on the symmetry-breaking pattern predicted.The folded density of spacing is where β is the Dyson index.This index can take three different values: β = 2 corresponds to the symmetry breaking pattern The latter two are the cases we are interested in, corresponding to fundamental and 2-index antisymmetric fermions for the symplectic theory.
In order to make a comparison with the chRMT prediction in Eq. ( 33), we compute the eigenvalues of Q m for N conf configurations.This process yields a set of eigenvalues λ The constant N is defined so that the density of spacing has unit average over the whole ensemble, ⟨s⟩ = 1.Finally, the (discretised) unfolded density of spacings, P (s), is obtained by binning numerical results for s and normalising it.In Fig. 10, we show an example of the folded distribution of eigenvalues of the Wilson-Dirac operator, computed numerically.As it can be seen, in the case of fermions in the fundamental representation, one finds a distribution that is compatible with the symmetry breaking pattern leading to the coset SU (2N f )/Sp(2N f ).Conversely, for fermions in the 2-index antisymmetric representation, our numerical results reproduce the prediction associated with the coset SU (2N as )/SO(2N as ).The spectacular agreement with chRMT confirms that there are no inconsistencies in our way of treating fermions.The size of the lattices we have considered has been chosen in order to reduce finite-size effects.These effects, as shown in Ref. [9], can become evident in smaller lattices and they lead to discrepancies due to some abnormally large spacings for the smallest and largest eigenvalues.This was interpreted to be an artefact due to the finiteness of the lattice size.We observe that, as in previous studies, the antisymmetric representation already matches the predictions in a 4 4 volume, while for the fundamental to reproduce the predictions chRMT, we had to remove the 200 lowest and highest eigenvalues (reducing the number of eigenvalues from 4096 to 3696).In this fashion, the differences with chRMT are no longer visible to the naked eye even for lattices with modest volume, Ṽ = (4a) 4 .

V. THE N = 2 THEORIES COUPLED TO FERMIONS: BULK PHASE STRUCTURE
In this section, we present our main results for the theory with N = 2, N f = 0, and varying number of fermions transforming in the antisymmetric representation, starting from N as = 4-for which we apply the HMC algorithm.We performed a coarse scan of the lattice parameter space, to identify phase transitions in the (β, m 0 ) plane, by studying the average plaquette, ⟨P ⟩, its hysteresis, and its susceptibility.We provide an approximate estimate of the upper bound coupling for the bulk phase, β * , above which there is no bulk phase transition, and hence one can safely  4) theory with Nas = 4 fermions transforming in the 2-index antisymmetric representation, with ensembles generated from a cold start, using the HMC.We show the value of the average plaquette, ⟨P ⟩, as a function of the bare mass, for a few representative values of the coupling.The lattice size is Ṽ = (8a) 4 , and each point is obtained by varying the lattice coupling β = 7.0, 6.8, 6.6, 6.5, 6.4, 6.3, 6.2, 6.0, 5.8, 5.6 and the bare mass −1.4 ≤ am as 0 ≤ 0.0.
perform lattice numerical calculations at finite lattice spacing, yet confident that the results can be extrapolated to the appropriate continuum limit.Figure 11 displays the average plaquette, ⟨P ⟩, obtained in ensembles generated using a cold start.The lattice size is V = (8a) 4 , and each point is obtained by varying the lattice coupling β = 7.0, 6.8, 6.6, 6.5, 6.4, 6.3, 6.2, 6.0, 5.8, 5.6 and the bare mass −1.4 ≤ am as 0 ≤ 0.0.The figure shows that, for small values of β and large, negative values of the bare mass, the average plaquette displays an abrupt change at a particular value am as * 0 , while being a smooth, continuous function elsewhere.This is a first indication of the existence of a first-order bulk phase transition.
To better understand whether a first-order phase transition is taking place, we study the effect of adopting two different strategies in the generation of the ensembles, repeating it using of thermalised (hot) starts, and redoing the measurements.Figure 12 shows the comparison of the average plaquette, ⟨P ⟩, computed for several fixed choices of the coupling β, while varying the bare mass −1.4 ≤ am as 0 ≤ 0.0.The two curves in the plots represent the behaviour measured in ensembles obtained from a cold and hot start configuration.The effects of hysteresis are clearly visible for β < 6.4 and are an indication of the presence of a first-order phase transition taking place at a critical value of the bare mass am as * 0 .The final test of the nature of the phase transition is shown in Fig. 13.For illustration purposes, we choose two values of the coupling for which we have evidence of a phase transition (β = 6.2), or of smooth behaviour of ⟨P ⟩ for all value of am as 0 (β = 6.5), respectively.We compute the plaquette susceptibility, defined as and compare the numerical results obtained with ensembles having two different volumes, Ṽ = (8a) 4 and Ṽ = (16a) 4 .
The results indicate that the peak height scales as the 4-volume when β is small, in which case the position of the peak also moves to a different value of am as 0 .These are indeed the expected signature of a first order phase transition.For large β, we observe a shift in values of am as 0 and no clear and large change in entity for the peak heights between the two volumes.This is a clear indication of a smooth crossover.We hence conclude that, in the theory with N = 2, N f = 0, and N as = 4, there is numerical evidence of a line of first-order phase transitions turning into a crossover at β > β * = 6.4.

A. Varying Nas
We repeat the parameter scan for other choices of N as , using the RHMC for all fermions when N as is odd, and the HMC algorithm otherwise.The purpose of the exercise is to study the dependence of the upper bound coupling for the bulk phase β * on the number of fermions, N as .Indeed, it is expected that for small N as we expect the theory to confine, while for larger values of N as ∼ N c as the theory should approach the lower end of the conformal window, and eventually lose asymptotic freedom-we recall that the latter requires to impose the bound N as < 33/4 in Sp(4), while setting the stage for a first truly non-perturbative determination of the former is the main motivation for this For pure gauge, we just vary the value of β.All the fermions are treated with the HMC/RHMC algorithms.The lattice size is Ṽ = (8a) 4 and the base mass is chosen in the range −1.4 ≤ am as 0 ≤ 0.0 for Nas ≥ 2, and −1.5 ≤ am as 0 ≤ 0.0 for Nas = 1.For the pure gauge theory, the coupling is chosen to be 1.0 ≤ β ≤ 16.0.For Nas = 1, we have chosen β = 7.1, 7.0, 6.9, 6.8, 6.7, 6.6, while for Nas = 2 we have β = 6.8, 6.7, 6.6, 6.5, 6.4, 6.2.For Nas = 3, the coupling is β = 6.8, 6.7, 6.6, 6.5, 6.4, 6.2, 6.0, 5.8, while for Nas = 5 we've chosen β = 6.6, 6.5, 6.4, 6.3, 6.2, 6.1, 6.0, 5.8.For Nas = 6, β = 6.4,6.3, 6.2, 6.1, 6.0, 5.8.For Nas = 7, β = 6.4,6.2, 6.1, 6.0, 5.9, 5.8 and for Nas = 8, β = 6.3, 6.1, 6.0, 5.9, 5.8, 5.7.

study.
The results of these studies are shown in Fig. 14, which displays our measurements of the average plaquette, ⟨P ⟩, as a function of the bare parameters of the theories.For the pure gauge Sp(4) theory, we get plaquette values that are in agreement with the ones shown in Ref. [156].The corresponding upper bound value of the coupling is roughly estimated to be β * ≃ 7.2.
For theories with dynamical fermions, we vary both the masses and the coupling of the theories.As can be seen from Fig 14, for N as = 1 the upper bound is β * ≃ 6.7.For N as = 2 the upper bound is β * ≃ 6.7, and for N as = 3 it is β * ≃ 6.5, in agreement with the values found in Ref. [2].At a larger number of fermions species, we obtain progressively smaller values of β for the upper bound of the bulk phase β: for N as = 5, we get β * ≃ 6.3.For N as = 6, the upper bound coupling is β * ≃ 6.2.For N as = 7, we get β * ≃ 6.1 ÷ 6.2 and for N as = 8, β * ≃ 6.1.Overall, we notice a trend according to which the more fermion flavors are present in the Sp(4), the smaller the upper bound value of the coupling we find and the bigger is the corresponding critical bare mass am (as) * 0 .

VI. SCALE SETTING AND TOPOLOGY
We return now to the theory with N = 2, N f = 0, and N as = 4.We discuss a scale setting procedure that uses the Wilson flow.We also monitor the evolution of the topological charge, to show that topological freezing was avoided.We focus the discussion on a few representative examples, although we checked that our conclusions have general validity for all choices of parameter relevant to this study.
The gradient flow [331], and its discretised counterpart, the Wilson flow [332], are useful for two complementary purposes.On the one hand, the Wilson flow provides a universal, well defined way to set the scale in a lattice theory, that is unambiguously defined irrespectively of the properties of the theory and of model-dependent considerations.On the other hand, the process we will describe momentarily consists of taking gauge configurations and evolving them with a flow equation, which results in the smoothening of such configurations, and the softening of short-distance fluctuations.The former property is beneficial because it allows to compare to one another different theories for which no experimental information is available (yet), and that might have different matter content.The latter characteristic allows, in practical terms, to reduce the short-distance numerical noise and the effects of discretisation in the lattice We follow Refs.[1,11] (and references therein).One introduces the flow time, t, as an additional, fifth component of the space-time variables, and solves the defining differential equation subject to the boundary conditions B µ (x, 0) = A µ (x).Here A µ (x) are the gauge fields, and the covariant derivatives are As anticipated, the main action of the flow is to introduce a Gaussian smoothening of the configurations, with mean-square radius √ 8t.
In order to use this object to introduce a scale, one defines the quantities and introduces a prescription that defines the scale on the basis of a reference value for either of the two.Two common choices in the literature are the scale, t 0 , defined by setting or the scale, w 0 , defined implicitly by the condition Both E 0 and W 0 are set on the basis of theoretical considerations.For example, Ref. [11] advocates to set W 0 = c w C 2 (F ), where C 2 (F ) = (1 + 2N )/4 is the quadratic Casimir operator of the fundamental representation in Sp(2N ) theories, and one sets c w = 0.5, though other choices are possible.On the discretised lattice, one replaces the gauge field, A µ (x), with the link variable, U µ (x), and the flow equation is rewritten by replacing B µ (x, t) with the new V µ (x, t) (with V µ (x, 0) = U µ (x)).There are then at least two ways to replace G µν with a discretised variable.We introduced the elementary plaquette P µν when defining the lattice action in Eq. ( 14).The clover-leaf plaquette operator, C µν , provides an alternative to the elementary plaquette, and can be seen as a simple form of improvement.We borrow the definition from Refs.[339,340], that for generic link variables U µ (x) reads: In principle, one would like to set the scale in a way that does not depend crucially on microscopic details.The scale setting using Wilson flow depends on the way the flow equation in Eq. ( 36) is latticised and how the observable Tr[G µν G µν ] in Eq. ( 37) is discretised.Therefore, different choices lead to different values for the scale at a given cutoff, but choosing a suitable flow time t allows us to set the scale while reducing drastically these effects.To this purpose, in Fig. 15 we consider the Sp(4) theory with N f = 0 and N as = 4, for two representative choices of β, and a representative choice of volume, Ṽ , and bare mass, am as 0 , and we show E(t) and W(t) as functions of the flow time, t, by comparing explicitly the results obtained by adopting either the elementary or the clover-leaf plaquette as defining the lattice regularisation of the action.The plots illustrate the general trend evidenced elsewhere in the literature, according to which the function W(t) displays a milder dependence on the short distance regulator.In the following, we set the scale w 0 by conventionally setting W 0 = 1 2 C 2 (F ).Recently, several studies for the usage of the Wilson flow observables were performed to define a non-perturbative running coupling-see, e.g., Refs.[333][334][335][336][337] and [338] for a review.This is a very intriguing application of these tools, but it requires a dedicated study with non-negligible effort and large lattices, which exceed the lattice sizes we explored in this introductory paper.
The topological charge density is defined as and the topological charge is Q L (t) ≡ x q L (x, t), where, again, t is the flow time.In general, the topological charge on the lattice is not quantised, and in cases where it is the physical quantity of interest-for example because one is working towards a determination of the topological susceptibility, as in Ref. [11] and references therein-one needs to evolve to large t, and introduce a rounding process.
For the current purposes, we do not need a discretisation algorithm: what we want to verify is that there is no evidence of topological freezing, and to this purpose we perform three simple tests.In Fig. 16 we display the value of Q L (t = w 2 0 ) in the Sp(4) theory coupled to N f = 0 and N as = 4 fermion species, for two values of the coupling, β, and a common value of the bare mass.We show how the topological charge evolves along the trajectories, and supplement it with a histogram displaying its distribution.Both visual tests confirm that there is no evidence of topological freezing.We can make these tests more quantitative by applying the standard Madras-Sokal windowing algorithm [322], and provide estimates of the integrated autocorrelation time τ Q of the topological charge, which in both examples, as shown in Fig. 15, turns out to be many orders of magnitude smaller than the number of trajectories.Furthermore, fits of the histograms are compatible with a Gaussian distribution centered at ⟨Q L (t = w 2 0 )⟩ = 0.The main message from this section is that the behaviour of the Wilson flow and of the topological charge, computed using the new software based on Grid, and tested on GPU architecture machines, to examine the properties of the lattice Sp(2N ) gauge theory with N = 2, N f = 0, and N as = 4, provide results that are broadly comparable to those in the literature for related, though different, field theories.This suggests that the implementation of the simulation routines and of the observables are both free from unwanted effects.

VII. SUMMARY AND OUTLOOK
A number of new physics models based upon Sp(2N ) gauge theories has been proposed in the literature, in such diverse contexts to include Composite Higgs Models, top partial compositeness, dilaton-Higgs models, strongly interacting dark matter models, among others.It is essential to the development of all these new physics ideas to provide model builders and phenomenologists with non-trivial information about the non-perturbative dynamics.
The programme of systematic characterisation of Sp(2N ) theories is still in its early stages, though.Prominently, the challenging question of identifying the lower end of the conformal window in these theories coupled to matter fields in various representations of the group requires the non-perturbative instruments of lattice field theory.As a necessary step in this direction, we developed and tested new software, embedded into the Grid environment to take full advantage of its flexibility.In this paper we reported the (positive) results of our tests of the algorithms, that set the stage for future large-scale dedicated studies.We focused particularly on the Sp(4) theory coupled to N as = 4 (Dirac) fermions transforming in the antisymmetric representation, that might be close to the onset of conformality.
We performed a long list of non-trivial exercises.We both tested the effectiveness of the algorithm and software implementation, but also provided a first characterisation of lattice theories that had never been studied beforealthough for present purposes we used comparatively small and coarse lattices.We reported in this paper illustrative examples demonstrating that there are no obvious problems in the software implementation.We computed effectively such observables as the averages of the plaquette and (real) Polyakov loop, the plaquette susceptibility, the Wilson flow, and the topological charge.We catalogued the first measurements of the critical couplings in Sp(4) lattice theories with N as < 33/4-below the bound imposed by asymptotic freedom-hence identifying the portion of lattice parameter space connected with the continuum theories of interest.
are Dirac fermions that transform in the fundamental representation of Sp(2N ), as indicated by the index a = 1, • • • , 2N , while the Ψ k ab ones, with k = 1, • • • , N as , transform in the 2-index antisymmetric representation of the gauge group.The covariant derivatives are defined by making use of the transformation properties under the action of an element U of the Sp(2N ) gauge group, according to which

 0 1 1 0 
i and ω jk = ω jk ≡  Nas Nas .The antisymmetric matrix Ω has the same form, as defined in Eq. (A1) of Appendix A, for both the gauge and fundamental flavour symmetries, but the indices run with a = 1, • • • , 2N for the former and with j = 1, • • • , 2N f for the latter.

FIG. 8 :
FIG.8: Compatibility between plaquette averages ⟨P ⟩ obtained with HMC and RHMC algorithms for the theory with N = 2, N f = 0, and Nas = 4. ⟨P ⟩HMC is obtained running two couples of fermions with HMC.For ⟨P ⟩RHMC (top panel), RHMC was applied individually to each of the fermions.⟨P ⟩2HMC+2RHMC (bottom panel) is obtained running two fermions with HMC, while the other two were run with RHMC.The lattice coupling is β = 6.8, with the bare mass in the range −1.4 ≤ am as 0 ≤ 0.0.The lattice is isotropic and has volume Ṽ = (8a)4 .

5 SUFIG. 10 :
FIG.10: Distribution of the folded density of spacing between subsequent eigenvalues of the hermitian Dirac-Wilson operator Qm = γ5Dm, and comparison with predictions from chRMT, computed in the quenched approximation, with ensembles having β = 8.0, am0 = −0.2, and lattice volume Ṽ = (4a)4 , in the Sp(4) theory.The left panel shows the case of fermions transforming in the fundamental representation, and the right is for fermions in the 2-index antisymmetric one.
(c) i with c = 1, • • • , N conf .The eigenvalues are arranged in one long list, in which λ (c) i are ordered in ascending order.Any degeneracy that is present in the 2-antisymmetric case is discarded.Then, for each c = 1, • • • , N conf , a new list of values is produced, that contains n (c) i , the positive integer position of the eigenvalue λ (c) i in the long list ordered in ascending order, instead of λ (c) i .The density of spacing, s, is replaced by the expression

FIG. 11 :
FIG.11: Parameter scan of the Sp(4) theory with Nas = 4 fermions transforming in the 2-index antisymmetric representation, with ensembles generated from a cold start, using the HMC.We show the value of the average plaquette, ⟨P ⟩, as a function of the bare mass, for a few representative values of the coupling.The lattice size is Ṽ = (8a)4 , and each point is obtained by varying the lattice coupling β = 7.0, 6.8, 6.6, 6.5, 6.4, 6.3, 6.2, 6.0, 5.8, 5.6 and the bare mass −1.4 ≤ am as 0 ≤ 0.0.