SU(2) Gauge Theory with Two Fundamental Flavours: a Minimal Template for Model Building

We investigate the continuum spectrum of the SU(2) gauge theory with $N_f=2$ flavours of fermions in the fundamental representation. This model provides a minimal template which is ideal for a wide class of Standard Model extensions featuring novel strong dynamics that range from composite (Goldstone) Higgs theories to several intriguing types of dark matter candidates, such as the SIMPs. We improve our previous lattice analysis [1] by adding more data at light quark masses, at two additional lattice spacings, by determining the lattice cutoff via a Wilson flow measure of the $w_0$ parameter, and by measuring the relevant renormalisation constants non-perturbatively in the RI'-MOM scheme. Our results for the lightest isovector states in the vector and axial channels, in units of the pseudoscalar decay constant, are $m_V/F_{\rm{PS}}\sim 13.1(2.2)$ and $m_A/F_{\rm{PS}}\sim 14.5(3.6)$ (combining statistical and systematic errors). In the context of the composite (Goldstone) Higgs models, our result for the spin-one resonances are $m_V>3.2(5)$ TeV and $m_A>3.6(9)$ TeV, which are above the current LHC constraints. In the context of dark matter models, for the SIMP case our results indicate the occurrence of a compressed spectrum at the required large dark pion mass, which implies the need to include the effects of spin-one resonances in phenomenological estimates.


I. INTRODUCTION
New composite dynamics is often invoked to construct extensions of the Standard Model (SM) physics that can address one or several of the SM shortcomings.
In TC models the Higgs boson is the lightest scalar excitation of the fermion condensate responsible for electroweak (EW) symmetry breaking [14][15][16][17][18].The physical Technicolor Higgs mass can be light due to near conformal dynamics [14,19] and the interplay between the TC sector and the SM fermions and electroweak gauge bosons [20].
In composite Goldstone Higgs models [12,13], the new sector has an underlying fundamental dynamics with larger global symmetry group than the one strictly needed to break the EW symmetry.In this case the Higgs state can be identified with one of the additional Goldstone Bosons (GB), and it is therefore naturally light.However, to break the EW symmetry, typically radiative corrections are not enough and yet another sector is required to induce the correct vacuum alignment for the EW gauge bosons and for the Higgs to acquire the observed mass.
The underlying fundamental theory studied here constitutes the ultra minimal composite template for any natural UV completion that simultaneously embodies both the TC and composite Goldstone Higgs models [21][22][23][24].It is also well known that fermion mass generation constitutes a challenge for any composite dynamic extension.For the present theory an extension that makes use of chiral gauge theories [25][26][27] has been put forward recently in [28].The constructions yield distinctive experimental signatures and can be used universally for both types of model building.
Novel composite dynamics has also been advocated to construct natural candidates for DM stemming from a composite EW sector.In fact, several asymmetric DM candidates were put forward which are stable baryons in TC models [29,30] or Goldstone bosons of a new strong sector [22,[31][32][33].
Another interesting class of DM models, unrelated to the composite EW scenario, was recently revived in [34].Here an alternative mechanism [35,36] is employed for achieving the observed DM relic density.It uses 3 → 2 number-changing processes that should occur in the dark sector involving strongly interacting massive particles (SIMPs).Compared to the WIMP paradigm, where the dark matter particles typically are expected to be around the TeV scale, this model can yield dark matter particles with masses of a few 100 MeVs.In [37,38] a realisation of the SIMP paradigm was introduced in terms of composite theories for which the model investigated here again provides the minimal template.In this realisation, the pions constitute the dark matter particles and the topological Wess-Zumino-Witten (WZW) term [39][40][41] introduces a 5point pion interaction, making it an ideal candidate for the 3 → 2 annihilation process.
The most minimal realisation of this breaking pattern comes indeed from the underlying Sp(2)=SU(2) gauge group (but can be generalised to any Sp(N c ) gauge group).The first consistent investigation of the phenomenological viability of this construction, that properly takes into account important next-to-next-leading-order corrections via chiral perturbation theory, appeared in [38].Here it was shown that higher order corrections substantially increase the tension with phenomenological constraints.Because the energy scale of the SIMP is very light, it is especially relevant to know at which energy scale dark spin-one resonances will appear, or more generally to understand its spectroscopy [42].
Furthermore the new states will modify the scattering at higher energies introducing possible interesting resonant behaviours [43] and, as it is the case for ordinary QCD, will impact on a number of dark-sector induced physical observables.
In this work we investigate the SU(2) gauge theory with N f = 2 flavours of Dirac fermions in the fundamental representation.One important feature of this minimal SU (2) template model is that, due to the pseudo-reality of the fundamental representation, the flavour symmetry is upgraded to an SU(4) (locally isomorphic to SO( 6)) symmetry which is expected to break spontaneously to Sp(4) (locally isomorphic to SO(5)), thus leading to 5 Goldstone bosons.
The theory has previously been studied on the lattice, and in particular, it has been shown that the expected pattern of spontaneous chiral symmetry breaking is realised [44].
A first estimate, affected by large systematic errors, of the masses of the vector and axialvector mesons, in units of the pseudoscalar meson decay constant, have been obtained in [1].The scattering properties of the Goldstone bosons of the theory have also been considered [45], and the model has furthermore been investigated in the context of possible DM candidates related to the EW scale [46,47].Other groups have also investigated the spectrum of this model on the lattice [48,49] concluding that chiral symmetry is broken, although no continuum extrapolation was attempted as the focus of both works was on the comparison with the six flavours theory to understand the approach to the conformal window in SU(2) gauge theories.
Here we extend our previous analyses by improving our control on the systematics.
Our simulations achieve smaller fermion masses, include two additional lattice spacings, and we also perform a precise determination of the lattice spacings used.Finally we determine the relevant renormalisation constants non-perturbatively.
The paper is organised as follows.We first describe the lattice setup in section II and the procedure to set the lattice spacing through the Wilson Flow observable w 0 in section III.In section IV we discuss the calculation of the renormalisation constants using the RI'-MOM scheme.Finally we provide in section V an improved estimation of the spectrum of the theory.

II. LATTICE SET-UP
We simulate the SU(2) gauge theory with two Dirac fermions in the fundamental representation discretised using the (unimproved) Wilson action for two mass-degenerate fermions u, d and the Wilson plaquette action for the gauge field.The numerical simulations have been performed using an improved version of the HiRep code first described in [50].The fermionic part of the action reads: where U µ is the gauge field, ψ is the doublet of u and d fermions, and am 0 is the 2 × 2 diagonal mass matrix proportional to the identity.
Our simulation are performed at four values of the inverse lattice gauge coupling β, for various fermion masses and on several volumes.This is needed in order to perform the  All the others runs are referred to in the text as "large volume runs".required extrapolations to the chiral limit and infinite volume and to give an estimate of the systematic errors stemming from such extrapolations.We detail the procedure used in the following sections.
The bare parameters of our simulations are summarised in Table I.We have extended our previously published dataset considerably, in particular towards the chiral regime and by adding two additional lattice spacings at β = 1.8, 2.3.As we will discuss in more detail below, note that the lightest quark masses now reach, in some cases, the decay threshold for the vector meson resonance.The simulations in Table I denoted with an asterisk, are only used to study the systematic errors due to finite size effects.The remaining runs will be referred to as "large volume runs" in this paper.This is justified as for all these lattices we have m PS L ≥ 5 which implies a systematic error of about 5% for the quantities studied here [1].
For convenience, we define the following operators: where Γ denotes any product of Dirac matrices.
We extract the meson masses from zero-momentum two-point correlation functions The quantities of interest in this study are the pseudoscalar Γ = γ 5 , vector Γ = γ k (k = 1, 2, 3), and axial vector Γ = γ 0 γ 5 γ k mesons.We use Z 2 × Z 2 single time slice stochastic sources [51] to estimate the meson 2-point correlators.From those, we define an effective mass m eff (Γ) (t) as in [52,53] by the solution of the implicit equation: where T is lattice time extent.At large Euclidean time, m eff (Γ) (t) approaches the value of the mass of the lightest state with the same quantum numbers as the operator O (Γ)  ud .In the following, we will denote the pseudoscalar meson mass m PS , and the isovector vector and axial-vector meson mass m V and m A respectively.
In addition to the meson masses above, we will use in the present analysis two other quantities: the current quark mass m PCAC and the Goldstone boson decay constant F PS .
We define the quark mass through the Partially Conserved Axial Current (PCAC) relation where The Goldstone boson decay constant can be calculated as: where G PS is obtained from the asymptotic form of f γ 5 (t) at large t: On a lattice of finite temporal extent, we use the same definitions as in [52,53].
The (bare) values in lattice units for m PCAC , m PS , F PS , m V and m A corresponding to the large volume lattices considered in this paper are reported in table VII in appendix A.
To convert the lattice quantities to physical units, we determine the lattice spacing for our simulations and the appropriate non-perturbative renormalisation constants.
It is well known that for Wilson fermions, the pseudoscalar decay constant renormalises multiplicatively with the scale independent renormalisation constant Z A and that the bare PCAC mass renormalises with the ratio Z A /Z P (µ 2 ).
The lattice spacing, in a generic composite model, is fixed by the requirement that the renormalised Goldstone boson decay constant has a given value specified for the physical model considered.For example in the case of composite dynamics at the electroweak scale, a value of 246 GeV yields the correct mass for the electroweak gauge bosons.For the more general fundamental composite Goldstone Higgs scenario described in [24] the scale is still set by the same requirement, but the constraint on the renormalised Goldstone boson decay constant now reads F PS sin(θ) = 246 GeV.The actual value of the parameter θ in this model depends on the electroweak gauge bosons corrections, the top corrections as well as the effects of other possible sources of explicit breaking of the initial SU(4) symmetry.The Technicolor limit is recovered for θ = π/2 while the composite pGB Higgs case corresponds to small, but non-vanishing θ.Any other value of θ is also allowed and the resulting model thus interpolates between these two limits.For the details we refer to [24].
Another case of immediate interest is the SIMPlest composite model [37] for DM where, as shown in [38], it is important to control the underlying dynamics.By stretching chiral perturbation theory to its limit of validity, the interesting phenomenological values for the pion decay constant would be as low as 10 MeV with pion masses of the order of 100 MeV.Besides the rescaling the pion decay constant, another major difference, when compared to composite dynamics at the electroweak scale, resides in the fact that the SIMP requires quite massive pions.
For definiteness, below we present our results in units of the EW scale with sin(θ) = 1 but the dependence on θ can be reinstated when needed.At the end we will also comment on the results for the SIMPlest case.

III. SCALE SETTING
Following [54], we consider the following "Wilson flow" equation for the gauge fields: where t denotes the fictitious flow "time", U(x, µ) are the gauge links, and S G is the plaquette gauge action.One important property is that correlation functions at flow time t > 0 are finite, when the four-dimensional theory is renormalised as usual, and the flow thus maps gauge fields into smooth, renormalised gauge fields [55].Observables at non-zero flow time can, in particular, be used to define a scale, as shown in [54].
Two different scale-setting observables have been introduced in the literature, known as t 0 [54] and w 0 [56].In terms of E(t), the action density at flow time t, they are defined through the following equations: where In order to extrapolate to the chiral limit, we use the NNLO expansion in terms of m 2 PS which reads [57]: where F is the pseudoscalar decay constant and k 1 , k 2 are dimensionless low energy constants.Note that the chiral logarithm enters only at NNLO.In practice we fitted our data at each β with the following ansatz : where A, B and w χ 0 are free parameters with the choice w χ 0 µ = 1.The fit is performed for each of the four β-values independently and the gray bands indicate the 1σ error regions.The best fit parameters and their statistical errors are reported in Table II.In the left panel of Fig. 1, the red dotted vertical line indicates the upper limit of the y 2 region used in the NNLO fit.
For three of our data sets we have also performed a fit to the NLO expression.The black vertical dotted line indicates the upper limit of the y 2 region included in the NLO fit.
Due to lack of data, we cannot perform this fit for β = 2.3.For the three remaining lattice spacings available, the results of the NLO and NNLO fits agree well within uncertainties.
In the right panel of Fig. 1 we show w 0 /w χ 0 for all four lattice spacings.The deviation from a universal curve of such a quantity is a measure of lattice discretisation errors.As can be seen, these are small in the w 0 observable for our three finest lattice spacings.The same conclusion can also be reached by looking at the dimensionless coefficient A and B as determined from the fits, given in Table II.

A. RI'-MOM scheme
In this section we describe the method used to determine the non-perturbative renormalisation constants of the isovector vector (V), axial (A), and pseudoscalar (P) bilinear operators.They are needed for the renormalisation of the pseudoscalar decay constant F PS and of the quark mass m PCAC .
We use the RI'-MOM scheme (regularisation invariant momentum scheme) as in [58].
We define the following bilinear operators : and the fermion propagator : S(x, y) = ψ(x)ψ(y) , and S(p) = p e ip(x−y) S(x, y).(15) Note that we have omitted to write explicitly spin and color indices.We also define the following Green's function: and we will denote the corresponding vertex function by: where S −1 (p) is the inverse propagator in spin and color space.The RI'-MOM scheme [58] is then defined by imposing the conditions that in the chiral limit and at a given scale p 2 = µ 2 , the inverse propagator and amputated Green's function Π Γ (p) satisfy the following equations: where the trace is over spin and color indices and the projectors P Γ are defined as follows: For convenience we define: such that in the chiral limit:

B. Evaluation of the correlators
Following the approach introduced in [59], we use momentum sources.This approach has the advantage to be computationally inexpensive and to have a high statistical accuracy.We will shortly summarise the procedure.
The vertex functions defined in Eq. ( 17) are not gauge invariant, and must be computed in a fixed gauge.We chose the Landau gauge by minimising a functional proposed in [60].
We introduce S(y, p) defined to be the solution of the following linear equation where ½ stands for the identity matrix in spinor and color indices.It is straightforward to obtain that and

C. Twisted boundary conditions
In order to interpolate easily between the lattice momenta we use twisted boundary conditions [61,62] by imposing : where L µ=1,2,3 = L and L 4 = T and θ is the twist angle.The boundary conditions are imposed by modifying the Dirac operator in the valence only.The accessible momenta are then p µ = 2π L µ n µ + π L µ θ µ .Eq. ( 25) and ( 26) of sect.IV B can be generalised in the case of twisted boundary conditions.
In practice the propagator S(p) and the Green's function G Γ (p) are evaluated for for every pair (l, ].Note that we also use negative values for l ′ in order to obtain the same values of p 2 from twisting with different initial momentum.This is useful in order to estimate cut-off effects.From Fig. 2 it is clear that they are small.Finally, note that we choose "non-democratic" momenta in Eq. ( 28).

D. Results & Analysis
The vertex functions Λ X for X ∈ {P, V, A, P/S} at a fixed quark mass, as a function of momentum (ap) 2 , are shown for β = 2.0 and β = 2.2 in Fig. 2. The filled symbols are obtained with twist angle θ = 0, while the empty symbol denotes the results obtained for θ 0.
In order to determine the value of the renormalisation constants, the first step is to extrapolate the result in the chiral limit.At fixed p 2 the behaviour of the vertex functions, which do not involve the pseudoscalar density, is expected to be polynomial in (am bare PCAC ) 2 .Concerning the pseudoscalar vertex functions, it is well known that special care must taken due to the presence of the Goldstone bosons pole [58,63,64].In that case, we use the following ansatz to perform pion-pole subtraction : where A, B and C are functions of p2 .The subtraction is performed for each p 2 by fitting the data at a given β, and we will denote m PCAC the subtracted vertex function at a given fermion mass.
We illustrate the chiral extrapolation at fixed p 2 in Fig. 3, where we show Λ X (p 2 = 1/a 2 ) as a function of (am bare PCAC ) 2 .In the plot we also included the Goldstone boson subtracted vertex function for X = P and P/S.The chiral extrapolation is obtained by fitting a second order polynomial in (am bare PCAC ) 2 to the data.The vertical dashed-dotted line indicates the extent of the region included in the fit.The best fit curve and its statistical error are included in the figure.The typical χ 2 /ndof for these fits are larger than one, because of the small statistical error bars on Λ X .Given our target accuracy of a few percent, those effects are negligible, however.
We show in Fig. 4 the dependence of the chirally extrapolated vertex function Λ X as a function of (w 0 p) 2 .In the continuum, Z V , Z A and Z P /Z S are renormalisation scale independent, the observed scale dependence is a manifestation of discretisation effects. 1n order to have meaningful estimates of Z X (p 2 ), one relies on the existence of a renormalisation window: Λ < p < O(a −1 ).The lower bound guarantees that the Goldstone pole contamination is small and that the Wilson coefficient entering in the operator product  III.Renormalisation constant obtained using (w 0 p) 2 = 7 as a reference scale expansion, which relates the physical process and the matrix element, can be computed in perturbation theory.The upper bound guarantees small lattice artefacts.In our case, reformulating the inequality in unit of w 0 , and setting w χ 0 Λ ∼ w χ 0 m V ∼ 1 we have: Since the smallest value w χ 0 /a obtained at β = 1.8 is w χ 0 /a ∼ 2, we would have 1 < w χ 0 p 2 < O (4).This is the famous window problem occurring at coarser lattice spacing.We thus have to relax the upper bound of the inequality and introduce larger cut-off effects for our coarser lattices.In practice we chose (w χ 0 p) 2 = 7, which corresponds to the lattice cutoff at β = 2.0.
In the following we will check that this particular choice of the reference scale does not affect scale-independent quantities, by using a second reference momentum, at the higher end of the sensible momenta range, namely: (w χ 0 p) 2 = 17.As shown below, our final results are very stable and do not depend, within errors, on the particular choice of reference momentum.
We summarise the values of the renormalisation constants, defined at our reference scale (w χ 0 p) 2 = 7, for the four β values, in Table III.

V. SPECTROSCOPY A. Effective Masses
We compute the mass of the lightest (isovector) pseudoscalar, vector and axial-vector meson resonances using two-point correlators.As explained in Section II, the mass can be extracted using the large time behaviour of the effective mass as decribed by Eq. (4).This approach is justified if the state is stable.We illustrate effective masses for various ensembles in Fig. 5, 6, 7 and 8.
The effective masses are fitted on a given plateau range, which is determined for each state by individual inspection.Systematic errors introduced by the choice of the plateaux are small for the pseudoscalar and vector resonances, and for this reason we will neglect them in the following.The best fit value for the effective mass is plotted for each state in the figures together with its statistical error.The masses of the vector and pseudoscalar mesons are clearly determined for all ensembles.For the axial vector correlator we do not observe long plateaux, due to the much worse signal-to-noise ratio as a function of Euclidean time separation.This results in significantly larger systematic errors, which are not yet fully under control.
In each plot, we also show the two-and three-pion thresholds.This shows that the vector meson resonance, whose main decay channel is expected to be the decay in two pions, is stable for almost all of our simulations.In a few cases, our most chiral points at β = 1.8 and β = 2.0 are at kinematical threshold.A similar conclusion can be drawn for the isovector axial-vector meson, whose main decay channel is expected to be three pions.

B. m PS and F PS
The continuum expressions for m PS and F PS have been worked out in [65] at next-toleading order in chiral perturbation theory: where x = 2Bm f (4πF) 2 and m f is the renormalised fermion mass at a given scale.In the conventions of [65], the condensate is given by Σ ≡ −2BF 2 .Note that F and B appear in both expressions.The range of applicability of the effective theory is not known a priori.In order to make the fits more stable, we will rewrite the expansion in a new parameter, x = m 2 PS (4πF) 2 .At this order Eq. ( 31) and Eq. ( 32) remain unchanged (this is, however, not true  at NNLO) and read: From this result we observe that the expansion of F PS now is independent of B, which will allow us to perform the fit in two steps: first a fit to F PS to obtain F and then using it as an  must then be taken into account.In order to obtain a reliable estimate, we will use two different strategies.
The first strategy (strategy I) is based on fitting the pseudoscalar mass and decay constant using several lattice spacings simultaneously together with a given model for the lattice discretisation effects: Here the new fitting parameters δ M,F and γ M,F control the discretisation effects.Note that the two coefficients a F,M are fixed in the continuum, but here we consider them as free parameters.
To control the stability of the fit, we consider four subsets of our data S 1 = {β = 2.0, 2.2}, 2} and S 4 = {β = 1.8, 2.0, 2.2, 2.3} and perform the fit on each of these subsets.The result of the fit for the S 2 subset is shown in Fig. 9 and 10.
The second strategy (strategy II) consists of fitting each of the lattice spacings independently, to obtain the coefficients B, F, a F,M and b F,M , while setting to zero the coefficients δ M,F , γ M,F in Eq. ( 35) and (36).In a second step, lattice discretisation effects can be assessed by studying the dependence of the coefficients as a function of the lattice spacing.
In all fits we use w χ 0 µ = 1 as a scale.The results of the fits, including their χ 2 per degrees of freedom, are summarised in Table IV for strategy I and Table V for strategy II.The fits are performed on a given range of values for (w χ 0 m PS ) 2 below the "cut" given in the tables.Strategy II allows us to extract an estimate of w χ 0 F and w χ 0 B for each lattice spacing.This is shown in Fig. 11, where the value of B has been re-scaled by a factor of 20 for convenience.The scaling towards the continuum limit is compatible with a linear behaviour and no O(a 2 ) effects are visible.On the plot we also show the results obtained directly in the continuum using the first strategy for the subset of gauge ensembles S 1 and S 2 .The results obtained with strategy I for the subsets S 3 and S 4 have a χ 2 /ndof ∼ 10 and thus do not describe the data well.
Our final estimates for the chiral parameters are w χ 0 B = 2.88(15) (17) and w χ 0 F = 0.078(4) (12).The central value and statistical error comes from the linear extrapola- We repeated a similar analysis using p = √ 17/w χ 0 as reference scale, which is shown in appendix B. As claimed in the previous section, we do not observe any statistically significant change in the continuum values of F and B.

C. Heavier states
In this section we report our results for the mass of two heavier isotriplet meson resonances, namely the vector in Fig. 12 and the axial-vector in Fig. 13.All the masses are presented in units of w χ 0 as functions of (w χ 0 m PS ) 2 .In each figure we present a global fit, including all the available data at four lattice spacings, to the following fit ansatz:   VI.
For the vector meson the fit describes our data well and the observed cutoff effects are small.We find w χ 0 m χ V = 1.01(3) with a χ 2 /ndof = 23/16.Note that for our data m V is always less than 2m PS , except maybe for the most chiral point used in the fit, so that the vector meson is expected to be stable and its mass can be reliably extracted from the large (Euclidean) time behaviour of the appropriate two-point function.
For the mass of the axial-vector meson, our data is more noisy already at the level of the effective masses and we therefore have larger systematic uncertainties.The ansatz Eq. ( 37) fits the data well, within large errors, and the resulting value for the mass is: 1) with χ 2 /ndof = 20/16.In units of F PS we have m V /F PS ∼ 13.1(2.2) and m A /F PS ∼ 14.5 (3.6).

VI. CONCLUSION
We analysed the SU(2) gauge theory with N f = 2 flavours of fermions in the fundamental representation using lattice techniques.Dynamical simulations have been performed at four lattice spacings and a number of volumes and masses to asses systematic effects and to carry out the necessary extrapolations.We determined non-perturbatively, in the RI'-MOM scheme, the relevant renormalisation constants and performed a detailed analysis of the mass and decay constant of the pseudoscalar Goldstone bosons, including an extrapolation to the chiral and continuum limits to take into account the lattice cutoff effects present in our computation.We use a conservative estimate of all systematic uncertainties to obtain a reliable estimate of F PS .Finally we analysed the mass of the spin-1 bound states and determined the ratios m V /F PS = 13.1(2.2) and m A /F PS = 14.5 (3.6) for the continuum theory in the chiral limit, using similar extrapolation methods.Our final results are consistent with, and improve upon, previous results for this model, which were performed with only two lattice spacings, at much larger quark masses and using a perturbative estimate of the renormalisation constants.
In the context of the fundamental composite (Goldstone) Higgs dynamics [24] which are beyond the present LHC constraints, even in the Technicolor limit [6] where In the context of dark matter models, in paticular for the SIMPlest case, because the dark pion is estimated to be around ten times its decay constant [38], we cannot use the estimate above.Nonetheless, a preliminary result can be obtained from our simulations reported in the first line of Table VII, at β = 2, which yield m PS /F PS ≈ 7.5, m V /F PS ≈ 8.3 and m A /F PS ≈ 13.7.Although these results need crucial refinement they immediately show that for such large values of the dark pion mass one cannot neglect the effects of higher mass states since the overall spectrum is much more compressed than in the case of the chiral limit.

Appendix A: Numerical results
We report in this section our numerical results for the main spectroscopy quantities studied in this article.The column "stat" reports the number of thermalised configurations used in the analysis, while the column "N rep " is the number of "replicas" runs used, i.e. number of independent runs with the same bare parameters.

Appendix B: Systematic error due to the choice of renormalisation scale
In this appendix we report the dependence of our continuum extrapolation results for the low energy constants F and B on the choice of renormalisation scale (w χ 0 p) 2 .The main result in the text are obtained using (w χ 0 p) 2 = 7 as the reference momentum scale.Here we present the same analysis for another value of the reference momentum scale: (w χ 0 p) 2 = 17.This corresponds to a much higher scale where lattice cutoff effects are expected to become more relevant.We show below in ref and W ref are two dimensionless reference values.In this work we will use w 0 to set the scale.The value of w 0 obtained for each quark mass needs to be extrapolated to the chiral limit to obtain a scale w χ 0 for each lattice spacing.We investigated finite volume errors in w 0 at the chosen reference value W ref by comparing two simulations performed on spacial sizes L = 16 (m PS L ∼ 5.1) and L = 32 (m PS L ∼ 8.4) at bare parameters m 0 = −0.75 and β = 2.2 .These values of the bare parameters were chosen to correspond to one of the lightest points in our dataset, at a fine lattice spacing.The values of w 0 obtained are w 0 (L = 16) = 3.39(6)a and w 0 (L = 32) = 3.36(10)a which agree well within statistical errors, indicating that finite volume effect for w 0 can be safely neglected for m PS L > 5 within our numerical precision.A. Determination of w χ 0 In Fig. 1 (left panel) we show our results for w 0 /a for the four lattice spacings considered in this study as a function of y 2 , where y = w 0 (m PCAC ) m PS .Here the reference value chosen is W ref = 1.For all the points in Fig. 1 we have m PS L > 5.5 and are thus safe from finite volume effects.

3 FIG. 1 .
FIG. 1. Chiral behaviour of w 0 as a function of y 2 in unit of the lattice spacing (left panel) and in unit of w χ 0 (right panel) for W ref = 1.The data at four lattice spacings are displayed.

5 FIG. 2 .
FIG. 2. Λ X=P,V,A,P/S as a function of (ap) 2 for the most chiral point at β = 2.0 (left panel) and β = 2.2 (right panel).The filled data points are obtained without twisted boundary conditions while the empty symbols denote the use of a non vanishing θ.The vertical lines indicate where (ap) 2 = 1.

3 FIG. 9 .
FIG. 9. F PS versus m 2 PS for the four lattice spacings.The curves correspond to the best fit parameters obtained fitting only β = 2.0, β = 2.2 and β = 2.3 (subset S 2 ) and drawn for the corresponding lattice spacing.The black curve indicate the continuum results.

3 FIG. 10
FIG. 10. m 2 PS /m f versus m 2 PS for the four lattice spacings.The curves correspond to the best fit parameters obtained fitting only β = 2.0, β = 2.2 and β = 2.3 (subset S 2 ) and drawn for the corresponding lattice spacing.The black curve indicate the continuum results.

input for a second fit to m 2 PS
/m f to obtain B. The renormalised values for F PS and m 2 PS /m f at four values of the lattice spacing are shown as function of m 2 PS in Fig. 9 and 10.All the lattices included in the fit satisfy m PS L ≥ 5.6.The fermion mass is given by m f (p 2 ) = m PCAC Z A /Z P (p 2 ) and the renormalised pseudoscalar decay constant is F PS = F (bare) PS Z A .As a reference scale for the renormalisation constants we use p = √ 7 /w χ

FIG. 11 .
FIG.11.Values of the chiral parameters B and F in units of the reference lattice scale w χ 0 as extracted using strategy II described in the text.In this plot B has been rescaled by a factor of 20 for graphical convenience.

FIG. 12 .
FIG.12.Combined chiral and continuum extrapolation of the vector meson mass m V .Our data for four lattice spacings is presented together with the best fit at each lattice spacing.The grey band is our result for the continuum extrapolation and its 1-σ confidence region.

FIG. 13 .
FIG.13.Combined chiral and continuum extrapolation of the axial vector meson masss m A .Our data for four lattice spacings is presented together with the best fit at each lattice spacing.The grey band is our result for the continuum extrapolation and its 1-σ confidence region.

3 FIG. 14 .
FIG. 14. F PS versus m 2 PS for four lattice spacing, using (w χ 0 p) 2 = 17 for the renormalisation scale.The curves correspond to the best fit parameters obtained fitting only β = 2.0, β = 2.2 and β = 2.3 (subset S 2 ) and drawn for the corresponding lattice spacing.The black curve indicate the continuum results.

3 FIG. 15
FIG. 15. m 2 PS /m f versus m 2 PS for four lattice spacing, using (w χ 0 p) 2 = 17 for the renormalisation scale.The curves correspond to the best fit parameters obtained fitting only β = 2.0, β = 2.2 and β = 2.3 (subset S 2 ) and drawn for the corresponding lattice spacing.The black curve indicate the continuum results.

TABLE I .
Parameters used in the simulations.Runs with * are used only to study finite size effects.

TABLE IV .
0 .As can be seen, significant cut-off effects are observed.In order to estimate the low energy constants F and B in the continuum, discretisation effects Results of the global fits of m 2 PS /m f and F PS on subset S 1,2,3,4 using (w 0 p) 2 = 7 as a reference renormalisation scale.

TABLE V .
Results of the fixed lattice spacing fits for each β value using (w 0 p) 2 = 7 as a reference renormalisation scale.

TABLE VII :
Numerical results for large volume runs used in the analysis presented in this paper.