The Yang-Mills gradient flow and SU(3) gauge theory with 12 massless fundamental fermions in a colour-twisted box

We perform the step-scaling investigation of the running coupling constant, using the gradient-flow scheme, in SU(3) gauge theory with twelve massless fermions in the fundamental representation. The Wilson plaquette gauge action and massless unimproved staggered fermions are used in the simulations. Our lattice data are prepared at high accuracy, such that the statistical error for the renormalised coupling, g_GF, is at the subpercentage level. To investigate the reliability of the continuum extrapolation, we employ two different lattice discretisations to obtain g_GF. For our simulation setting, the corresponding gauge-field averaging radius in the gradient flow has to be almost half of the lattice size, in order to have this extrapolation under control. We can determine the renormalisation group evolution of the coupling up to g^2_GF ~ 6, before the onset of the bulk phase structure. In this infrared regime, the running of the coupling is significantly slower than the two-loop perturbative prediction, although we cannot draw definite conclusion regarding possible infrared conformality of this theory. Furthermore, we comment on the issue regarding the continuum extrapolation near an infrared fixed point. In addition to adopting the fit ansatz a'la Symanzik for performing this task, we discuss a possible alternative procedure inspired by properties derived from low-energy scale invariance at strong coupling. Based on this procedure, we propose a finite-size scaling method for the renormalised coupling as a means to search for infrared fixed point. Using this method, it can be shown that the behaviour of the theory around g^2_GF ~ 6 is still not governed by possible infrared conformality.


I. INTRODUCTION
The discovery of a Higgs-like light scalar state at the Large Hadron Collider (LHC) has stimulated a significant amount of studies in various electroweak symmetry breaking (EWSB) models which can accommodate such a state. The Standard Model (SM) contains a fundamental scalar Higgs field, and is successful in explaining all the experimental results hitherto. Nevertheless, the scalar sector of the SM is widely believed to be trivial [1][2][3], therefore the cut-off is indispensable. This makes its predictions of low-energy quantities, such as the Higgs boson mass, sensitive to new physics effects which can appear as higher-dimensional irrelevant operators [4][5][6][7]. In view of this, it is desirable to find an alternative EWSB model in which the interaction is described by a relevant operator. There are various approaches to achieve this while having a light scalar state in the spectrum. Amongst these, one possibility is that this scalar particle is a dilaton [8][9][10][11] in the composite Higgs, or walking technicolour (WTC), scenario [12][13][14]. In this scenario, it is necessary to construct a gauge theory which exhibits asymptotic freedom and quasi scale invariance in the infrared (IR) regime. The Goldstone boson (the dilaton) resulting from the breaking of the IR scale invariance can be parametrically light compared to all the other states [15,16]. In addition to the possible existence of a light scalar state, WTC models also contain other appealing features. Any such model can incorporate the generation of fermion masses and the origin of flavours, withal dynamical suppression of flavour-changing neutral current (FCNC) processes.
The search for gauge theories that are viable for WTC model building has been a popular subject in the lattice community [17,18] 1 . This subject has led to studies of phase structure in lattice gauge theories with many flavours of fermions [20][21][22][23][24], broadening our understanding of the lattice regularisation. In addition to numerical calculations, it is also resulting in progress in analytic work related to lattice gauge theory [25][26][27][28][29][30][31]. One important task in this research avenue is to determine the conformal window, for a gauge theory coupled to N f fermions in a particular representation 2 , where N AF f is the number of fermions above which asymptotic freedom is lost, and N cr f is the "critical" number of fermions below which the theory is confining in the IR. A candidate theory for the WTC scenario is believed to have N f just below N cr f . This makes the determination of N cr f an endeavour with phenomenological importance.
Amongst the intensive investigations of candidate WTC theories using the lattice technique, the study of SU(3) gauge theory with twelve flavours of fermion in the fundamental representation has a long history. Several groups of authors found that this theory is conformal in the IR [22,[35][36][37][38][39][40][41][42][43][44][45][46]. Nevertheless, in Refs. [47,48] it was argued that chiral symmetry is broken in this theory. Amidst all these works, one popular approach is the step-scaling method for computing the running coupling constant, which was originally formulated in the Schrödinger-functional (SF) scheme [49,50]. This method was first used to determine the low-energy behaviour of SU(3) gauge theory with N f = 12 in Ref. [37], where the authors claimed the existence of an infrared fixed point (IRFP) by studying the coupling constant in the SF scheme. Recently, we adopted the same procedure in a different renormalisation scheme, namely the Twisted Polyakov Loop (TPL) scheme [51][52][53], in our investigation of the same theory, and found evidence for IR conformality as well [54] 34 .
With the lattice computation for the search of IR-conformal field theories becoming mature, the importance of controlling errors in such calculations is now receiving growing attention. To illustrate this point as relevant to our work presented here, in Fig. 1 we plot the perturbative predictions for the change of a generic renormalised coupling constant with respect to the doubling of the length scale in 12-flavours SU(3) gauge theory 5 . This plot indicates that the β−function is very small in this theory, and there can be an IRFP at strong coupling where perturbation theory may not be applicable. To obtain concrete evidence for confirming or ruling out the existence of this IRFP, one can resort to the lattice numerical implementation of the step-scaling method. However, Figure 1 also shows that in order to make any statistically-meaningful statement in such lattice studies, one would have to control the error to the sub-percentage level. Without this accuracy, it is difficult to demonstrate that the theory flows out of the vicinity of the asymptotically-free ultraviolet fixed point (UVFP), and then flows towards the IRFP when increasing the length scale. In previous step-scaling calculations for SU(3) gauge theory with twelve flavours, such statistical precision was not achieved. In addition, it was difficult to estimate the systematic effects in the continuum extrapolation in those calculations. This makes it challenging to draw reliable conclusions from these computations. For instance, in our previous step-scaling work employing the TPL scheme [54], the lattice simulations were performed to give statistical error around 2% in the extracted coupling constant 6 . The evidence for IR scale invariance discovered there was subsequently shown to be unreliable, after the addition of data at a larger volume that enables better estimation of the continuum-extrapolation error [55].
In this project, we perform lattice simulations with massless unimproved staggered fermions and the Wilson plaquette gauge action from which the Yang-Mills gradient flow [56] is implemented. As pointed out in Ref. [57], this method allows for the extraction of the renormalised coupling in the "gradient-flow (GF) scheme", via computing the energy density of the Yang-Mills field. The step-scaling approach is carried out with the step-size s = 2, on the lattice sizesL = L/a = (8, 10, 12) −→ 2L = (16,20,24), where L is dimensionfull, and a is the lattice spacing. Since our procedure only involves the calculation of the smallest Wilson loops on the lattice, we can determine the renormalised coupling with high statistical accuracy. Using about one-hundred thousand Hybrid Monte-Carlo (HMC) trajectories, we are able to control the statistical errors to be within 0.5% for the renormalised couplings computed on our lattice data. Furthermore, we adopt two different discretisations, namely the plaquette and the clover operators, in obtaining the energy density. This enables us to investigate the reliability of our results in the continuum limit. It has been well known that the continuum extrapolation is the main source of systematic error in the step-scaling study, therefore this feature of our analysis is welcome.
The current article is a report for our analysis of the GF-scheme coupling constant in SU(3) gauge theory with twelve fermions in the fundamental representation. In this work, we are able to probe the low-energy regime of this theory, with the bare coupling g 2 0 ∼ 1.45 and the corresponding renormalised coupling g 2 GF ∼ 6, using lattices as large asL = 24. In this regime, the coupling runs significantly slower than the two-loop perturbative prediction. However, our work does not allow us to draw definite conclusion regarding the existence of an IRFP in this theory. At bare couplings larger than 1.45, we begin to observe the onset of the bulk phase structure of the lattice theory. This means that the investigation of the theory at stronger renormalised coupling can only be achieved with simulations at larger lattice volume, such asL = 32. This is beyond the scope of this project.
In this article, we discuss in detail the application of the continuum extrapolation ansätz a'la Symanzik in lattice studies for the conformal windows in gauge theories. We point out that one has to be cautious in using this conventional method for confirming IR scale invariance. When the theory is tuned to be close to the possible IRFP, the continuum extrapolation may have to be conducted using a formula containing an unknown power of the lattice spacing. Since this kind of extrapolation is very challenging to carry out in practice, we propose a finite-size scaling method based on the same IR scaling property. We perform this finite-size scaling test in this work. The result of this test indicates that at g 2 GF ∼ 6, the behaviour of theory is not governed by IR conformality. This paper is organised in the following way. In Sec. II we discuss our strategy and lattice simulations. We then present the details of our analysis and results in Sec. III. Section IV contains our comment on the continuum extrapolation, and the proposal of a finite-size scaling test of the renormalised coupling in the strong-coupling regime. We compare our result with previous lattice computations in Sec. V, and conclude in Sec. VI.

II. STRATEGY AND THE LATTICE SIMULATION
A. The colour-twisted boundary condition In this work, we make use of twisted boundary condition (TBC) [58], where the gauge field is periodic up to a gauge transformation (µ, ν = 1, 2, 3, 4 are the Lorentz indices) where L ν is the linear size in the ν direction (withν being the unit vector). The SU (3) matrices Ω µ (x) are called twist matrices and must obey the consistency relation Where n µν is an anti-symmetric tensor of integers modulo 3 called twist tensor. The concrete choice of twist matrices is largely irrelevant, since, they change under gauge transformations. All the physical information about the twisted boundary conditions is contained in the twist tensor n µν . Our particular choice consists in using n µν = −n νµ = 1 µ = 1 and ν = 2 0 otherwise (5) On the lattice this particular choice can be realised by choosing the twist matrices constant in space and obeying the relations A concrete representation is Note that on the lattice the link variables, U µ (n) (µ = 1, 2, 3, 4 is the Lorentz index andn is the position of a lattice site), obey the boundary condition Fermions must also obey some specific boundary conditions in order to maintain gauge invariance and singlevaluedness of the action. In particular it is well known that the number of fermion flavours have to be an integer multiple of the rank of the gauge group [59]. The different flavours (usually called "smells") transform one in the other under translations of a full period of the torus according to where ν = 1, 2 and Latin (a, b) and Greek (α, β) indices are colour and smell indices, respectively. The factor e iπ/3 is introduced to lift the zero-momentum modes in these directions. In the directions ν = 3, 4, we impose PBC, ψ(n +νL ν /a) = ψ(n).
Using TBC in lattice calculations leads to various advantages. It removes the toron configurations. Therefore in the weak-coupling regime, the power laws in the coupling in finite-volume perturbation theory are the same as those in the infinite-volume, continuum case [60][61][62]. This boundary condition also lifts the zero-momentum modes, making it possible to perform simulations at vanishing fermion mass.

B. The gradient flow and the renormalised coupling constant
In recent years the Yang-Mills gradient flow [56,63] has become a standard technique to define renormalised couplings (see [64] for a recent review).
The basic idea is to add an extra coordinate to our gauge field, that we will call "flow time" and denote it by t. The evolution of the gauge field B µ (t, x) with respect to the flow time is given by the Yang-Mills gradient flow equation, where G µν is the field strength associated with B µ (t, x), and D µ = ∂ µ + [B µ , ·] the corresponding covariant derivative. The initial condition for Eq. (10) is given by is the fundamental gauge field. It has been proven to all orders in perturbation theory that gauge invariant composite observables made of the flow field B µ (t, x) are automatically renormalised for t > 0 [65]. In particular Ref. [56] suggests using the action density to define a renormalised coupling at a scale µ = 1/ √ 8t, with In the context of finite size scaling, the renormalisation scale is identified with the linear size of the box µ = 8t/L is a constant that defines our renormalisation scheme. Several finite volume renormalisation schemes have been defined by using different boundary conditions: periodic [57], Schrödinger functional [66], open-SF [67], and more directly related with this work, twisted boundary conditions [68].
There is quite some freedom when translating the flow equation Eq. (10) to the lattice. Since the r.h.s. of Eq. (10) is just the gradient of the Yang-Mills action, a straightforward option consists in evolving the lattice gauge links V µ (t, x) according to the gradient of a lattice action 7 Moreover one has also the freedom to choose amongst different discretisations to define the observable E(t), the most popular one being the clover and the plaquette. In general when evaluating the coupling these two choices of discretisations (i.e. the lattice flow equation and the observable) are, together with the action chosen to produce the configurations, the three sources of cutoff effects [69,70]. Although recently a discretisation of the flow equation and observables free of O(a 2 ) effects has been proposed [71], in this work we have used the Wilson flow (i.e. Eq. (13) with the Wilson plaquette action for S latt ), and two different discretisations of the observable (clover and plaquette).
In full glory our coupling reads,ḡ where is the discretised version of E(t) defined in Eq. (12), and the constantN (c τ , a/L) has been computed to tree-level with our choice of discretisations (see [68] for the concrete expressions).

C. The step-scaling method
In identifying gauge theories that are viable for the walking technicolour scenario, one may have to follow the evolution of the running coupling constant for a large range of scales. Although the lattice regularisation introduces both IR (the volume) and ultraviolet (the lattice spacing) cut-off scales which normally differ by only a factor of twenty to one hundred, such a task is made possible by employing the step-scaling method. In the following, we briefly review this method, in the context of our calculation of the coupling constant in the GF scheme.
To apply the step-scaling technique in our work, we first compute the coupling constant,ḡ latt , defined in Eq. (14). This technique relies on the use of the lattice size as the renormalisation scale. To ensure that this is feasible, one has to make certain that in addition to any intrinsic scale in the theory under investigation, L is the only length scale that is being introduced in the analysis. This means that we have to work with fixed c τ 8 , and remove the effects of the lattice spacing.
By performing this computation at various values ofL and β [Eq. (15)], we can tune these input parameters to determine the renormalised coupling that does not depend on the lattice spacing, In practice, one simulates at fixed values ofL. Since the lattice spacing, a, is determined by the bare coupling, β (or g 2 0 ), the above equation describes a procedure of adjusting a to achieve a constant physical length scale, L at various chosenL in the simulations. This tuned coupling, g GF (L) =ḡ cont (L), is the "input coupling" in the 7 For the precise definition of the Lie-algebra valued derivative ∂x,µ see [56]. 8 Different choices of cτ correspond to various renormalised schemes defined using the Yang-Mills gradient flow. determination of the step-scaling function. Following the conventional notation, we define the "input variable", u, as In practice, the tuning in Eq. (16) is carried out by choosing a few (as denoted by n 0 ) values ofL, and adjusting β accordingly to achieve the above condition. After this tuning,ḡ cont (L) does not depend the lattice spacing, therefore must be renormalised at the length scale L. In this work, our simulations for this "tuning procedure" have been performed at L/a = 8, 10, 12.
The details of the tuning for β will be given in Sec. III A.
Using the n 0 values of (β, L/a), as tuned via Eq. (16), we then calculate the step-scaling function, where s is the step size. With these n 0 results for Σ at the same physical volume, L, but different lattice spacings, we can perform the continuum extrapolation for the step-scaling function, We set the step size s = 2 in this work. Given the choices ofL in Eq. (18), we then have to compute the step-scaling function at sL/a = 16, 20, 24.
With these three values ofL, two different discretisations for the energy density in Eq. (12) are implemented. This allows us to investigate the systematic effects arising from the continuum extrapolation. We will give details of this issue in Sec. III B.
To make the presentation clear, we define The step-scaling function is simply the renormalised coupling, therefore its value depends on the choice of renormalisation scheme. A more suitable approach in demonstrating the existence of the IRFP is through the calculation of the ratio, This ratio becomes one when the β-function vanishes. The existence of zeros of the β−function is schemeindependent, although the values of the coupling at these zeros are not. In order to confirm that an asymptoticallyfree gauge theory contains an IRFP, we have to show that in this theory r σ (u) is one at both ultraviolet (UV) and IR regimes, while being positive in between. This demonstrates that when increasing the length scale, the theory flows out of the vicinity of the UV Gaussian fixed point, and then flows towards the IRFP at strong coupling.

D. Simulation parameters and the raw data
In this work, we perform simulations using the Wilson plaquette gauge action, and unimproved massless staggered fermions. As discussed in Sec. II A, the number of flavours in our calculation has to be a multiple of N c = 3, as a result of the introduction of the smell degrees of freedom. Staggered fermions contain four tastes, making it suitable for the lattice computation of SU(3) gauge theory with twelve flavours.
To implement the step-scaling method for computing the running coupling constant, as discussed in Sec. II C, we carry out simulations at the lattice volumes, L/a =L = 8, 10, 12, 16, 20, 24.
For each volume, we perform lattice calculations at several values of the bare coupling constant. As mentioned in the Introduction, the slow-running behaviour of the coupling in SU(3) gauge theory with twelve flavours results in the demand for controlling the statistical error to the subpercentage level. Furthermore, this behaviour makes it necessary to perform computations at a large range of bare coupling (lattice spacing), in order to trace the renormalised coupling from the UV to the IR regimes. We have the input bare coupling, g 2 0 , mostly in the range, in the simulation. For eachL, we perform computations at 20 to 34 choices of β, as summarised in Table I. In our analysis procedure, we have to interpolate the renormalised couplings extracted from simulations at different bare couplings, as detailed in Sec. III A. Having a large number of data points for this interpolation is essential for obtaining statistically independent results in the step-scaling study, as discussed at the end of Sec. III C. In this table, we also list the corresponding results ofḡ 2 latt at the largest g 2 0 for a few selections of c τ . The values of g 2 latt and g 2 0 differ significantly in the strong-coupling regime. On the other hand, as a consequence of asymptotic freedom, at g 2 0 ∼ 0.3, we observe that the bare and the renormalised coupling strengths are similar. To avoid strong cut-off effects, one should work with large enough c τ , such thatḡ 2 latt decreases with increasing c τ at fixed g 2 0 andL. For this purpose, it is clear from Table I that our analysis has to be carried out with data obtained at c τ ≥ 0.375. We implement the Yang-Mills gradient flow for many values of c τ . However, in order to illustrate our study in this work, it suffices to present results at in the remaining of this paper.
As discussed in Refs. [54,74], in lattice calculations of SU(3) gauge theory with twelve massless flavours, the Markov chains can exhibit tunnelling behaviour amongst local minima of the effective potential. Such behaviour may lead to artificially long autocorrelation time, and should be monitored carefully. To reduce the probability of having this tunnelling in our simulations, we thermalise the Markov chains by starting from a configuration with zero fermion mass, and U µ (n 1 ,n 2 ,n 3 ,n 4 ) = 1 elsewhere. This forces the system to be at the true vacuum, in which Polyakov loops in the untwisted directions contain nonvanishing imaginary parts [74]. The above prescription also produces the largest gap in the vicinity of zero modes in the fermion matrix. We observe that the autocorrelation time for the renormalised coupling,ḡ 2 latt , remains at about 20 to 100 HMC trajectories, and it increases with c τ . We have been able to obtain data with statistical uncertainties forḡ 2 latt below 0.5%, even for the largest flow time that corresponds to c τ = 0.5 in this work. This is achieved with about 100,000 HMC trajectories.
It is also important that we implement our simulations away from artificial bulk phases in the lattice theory [21]. For this purpose, we have checked the plaquette expectations values and confirmed that all our computations were performed in the weak-coupling phase.
In Appendix A, we tabulate all our raw data for the renormalised couplings, extracted using both clover and plaquette discretisations, at c τ = 0.5. Lattice data at other values of c τ can be obtained from the authors upon request.

III. ANALYSIS DETAILS FOR THE STEP-SCALING STUDY
In this section, we give the details of our analysis procedure. As shown in Table I in the last section, we perform simulations at 172 combinations of the bare coupling, g 0 , and the lattice size,L. It can also be seen in this table that typical statistical errors for our raw data at the largest couplings are between ∼ 0.2% (c τ = 0.375) to ∼ 0.45% (c τ = 0.500). This feature is common in all values of (g 0 ,L) in this work. Such error budget is achieved by having at least ∼ 100, 000 HMC trajectories, followed by carrying out measurement every 50 to 200 trajectories, and creating 100 to 200 bins for data at each simulation. One thousand bootstrap samples are created for our analysis.

A. Interpolation in the bare coupling
In principle, the step-scaling method described in Sec. II C has to be implemented by tuning the bare coupling, on the L/a = 8, 10, 12 lattices, to achieve the condition in Eq. (16). However, in this work we aim at tracing the coupling constant over a large range of scale, while computing the ratio, r σ (u), at many values of the input renormalised coupling, u, in the process. This renders the time-consuming tuning procedure impractical. In view of this, we make use of an interpolation method as the substitution for the tuning of the bare coupling. Namely, for each L/a we perform simulations at many values of β in a wide range, and then obtainḡ latt (β, L/a) at other bare couplings in this range through interpolation. In the rest of this section, we discuss this procedure in detail.
The feature of the input β−values in our simulations is given in Table I of Sec. II D. Since the chosen range of bare coupling constant straddles between the perturbative and non-perturbative domains, it is challenging to perform the interpolation in β using a well-inspired function. We begin by noting that in the perturbative (large−β) regime, one-loop approximation has to be applicable. Therefore, at fixed L/a, u latt (β, L/a) ≡ḡ 2 latt (β, L/a) ≈ with g 0 being the bare gauge coupling. Equation (28) also leads to the motivation for using polynomials in 1/β to carry out the bare-coupling interpolation. In this work, we perform simulations for many values of β at each L/a. This makes it possible to use high-degree polynomials which can in principle result in good fits. Nevertheless, having such high-degree polynomials in the interpolation procedure can introduce artificial oscillatory behaviour of the fit function, known as the Runge phenomenon. In order to avoid this artefact, we first note that the renormalised coupling, u latt , should be non-decreasing in 1/β at fixed L/a in this work, since our simulations are all performed in the weak-coupling phase of the lattice theory. In view of this, we can impose the non-decreasing constraint on the bare-coupling interpolation and use the function, where It can be seen in Fig. 2, where we plot u latt against g 2 0 , that this constraint is well-justified. Combining Eqs. (28) and (29), we further impose h 0 = 0, h 1 = 6 (then c 0 = √ 6).
This condition results in the number of fit parameters, with N deg and N h defined in Eq. (29).
In addition to reducing the artificial oscillatory behaviour in the bare-coupling interpolation, the non-decreasing polynomial fit function has the advantage that its inverse is singled-valued. This single-valuedness is an essential requirement in the implementation of the step-scaling method. It is also a necessary consequence resulting from the fact that all our simulations are performed in the weak-coupling phase of the lattice theory. Figure 2 shows the bare-coupling interpolation at c τ = 0.5 using Eq. (29). In this figure, N param is fixed to result in the best χ 2 /d.o.f. volume by volume (see Table II). The interpolation is performed with uncorrelated fits. It is clear from this plot that the resulting curves are smooth and they explain the data well. We find the same behaviour for all our other choices of c τ . In Table II, we give the values of N param corresponding to the best χ 2 /d.o.f. in these fits. In this interpolation procedure, we observe that χ 2 /d.o.f. is close to unity for all the optimal fits.
The results of this β−interpolation will be used in the step-scaling study of the renormalised coupling. It is obvious that this can lead to correlation amongst results of the step-scaling function and r σ determined at different input couplings. This issue will be addressed in Sec. III C.

B. Continuum extrapolation
Through the bare-coupling interpolation procedure, we can achieve the tuning of the input renormalised coupling in Eq. (16). This allows us to compute the corresponding lattice step-scaling function, Σ(β, L/a, u, s = 2), in Eq. (19). The next step in the analysis is the extrapolation of these lattice step-scaling functions to the continuum limit, as indicated in Eq. (20). Since our simulations are performed using unimproved staggered fermions and the standard Wilson plaquette gauge action, the lattice artefacts are polynomials 9 of (a/L) 2 . In this work, we carry out computations atL = 8, 10,12,16,20,24. This allows us to determine Σ(β, L/a, u, s = 2) at L = L/a (8, 10, 12) −→ 2L = (16,20,24) , then extrapolate to the continuum limit with the linear function Notice that this fit function is valid only when the effects of the lattice spacing are governed by the Gaussian fixed point in the UV. On the other hand, when the theory is engineered to be close enough to an IRFP, its scaling behaviour regarding the change of both the lattice spacing and the finite volume must be completely determined by the IR scale invariance. In this situation, the "Symanzik-type" polynomial extrapolation of Eq. (34) is not applicable, and alternative methods have to be adopted. We will discuss this issue in Sec. IV.
It is well known that the continuum extrapolation is the main source of systematic errors in the step-scaling investigation of the coupling constant. To check that this procedure is implemented reliably in our work, we make use of two discretisations, namely the clover and the plaquette, to compute the energy density defined in Eq. (12). These two discretisations contain different lattice artefacts. On the other hand, any result obtained with these methods should extrapolate to the same continuum limit, if the discretisation effects are under control. Figure 3 shows representative plots of the continuum extrapolation at c τ = 0.375 and 0.5. We first notice that all the extrapolations are mild. Even at strong coupling, the change ofḡ 2 latt (g 2 0 , 2L) from our coarsest lattice (L = 8) to the continuum limit is at the level of a few percent. This is partly because we extract the renormalised coupling using a result from lattice perturbation theory for the factorN in Eq. (14). For the case of c τ = 0.375 at strong renormalised coupling, the two discretisations do not lead to compatible results in the continuum limit. This renders the analysis unreliable in the IR regime. We stress that all fits for the continuum extrapolation produce good or acceptable χ 2 /d.o.f. (typically between 0 and 2) at this value of c τ , andḡ 2 latt (g 2 0 , 2L) only weakly depends on (a/L) 2 . Nevertheless, this does not mean that the procedure is under control. We further comment that this observation is made possible because our data are obtained at small enough statistical errors. In the same figure, it is demonstrated that effects of the lattice artefacts are reduced when c τ is increased. This is expected, since the smearing radius of the gauge field grows with c τ . For the case of c τ = 0.5, the clover and the plaquette discretisations produce consistent results in the continuum limit at all values of the coupling investigated in this project. In the current work, this extrapolation is under control, i.e., the continuum-limit results obtained from the two discretisation methods are compatible, up to g 2 GF ∼ 6 when c τ ≥ 0.45. This can also be seen clearly in Fig. 3.

C. Results and discussion
We present results from our main analysis in this section. As discussed in Sec. III A, we perform the bare-coupling interpolation using a non-decreasing polynomial function, Eq. (29), with the perturbation-theory constraint in Eq. (31). As already pointed out in previous similar studies, this interpolation procedure is not introducing significant systematic effects. In this work we observe that this is also true in our analysis, by varying the numbers of parameters reported in Table II. On the other hand, we find that the systematic errors associate with the continuum extrapolation can be significant. Therefore we concentrate on the discussion of this aspect of the error estimation in this section. Figure 4 shows results of r σ = σ(u)/u, as defined in Eq. (23), with the renormalised coupling computed using the clover and the plaquette discretisations. We first observe that the running of the coupling constant is very slow in SU(3) gauge theory with 12 flavours. Doubling the length scale leads to at most 4 to 5% change in the coupling constant in the range that our investigation is performed. Compared to a "fast running" theory, such as QCD in  which the coupling can increase by a few dozen percent when the length scale is doubled, the running is very slow in this theory. In order to provide evidence for the existence of an IRFP, one has to demonstrate that r σ is unity in both UV and IR regimes, while deviates from one in between. Therefore it is desirable to have high-precision data for such study of this theory. In this work we achieve good enough accuracy, and it is clearly discernible that the theory flows out of the vicinity of the UV Gaussian fixed point when the length scale is increased. At low energy, where g 2 GF ∼ 6, the results of r σ indicate that the coupling runs significantly slower than the two-loop perturbative prediction. For the case of c τ = 0.5 presented in Fig. 4, r σ is almost consistent with unity in this regime. This provides evidence that the scaling of the theory may be governed by IR conformality in this region. However, as already pointed out in Sec. III B, the continuum extrapolation for the results in Fig. 4   which may not be applicable near an IRFP. It requires further scaling test to clarify this issue. We will discuss this point in detail in Sec. IV.
Our analysis relies on the interpolation method reported in Sec. III A, in order to efficiently perform the timeconsuming tuning procedure of Eq. (16). Although we have many data points for this interpolation (Table I), it will still introduce correlation amongst r σ computed at different input g 2 GF . That is, values of g 2 GF (2L)/g 2 GF (L) at different g 2 GF (L) presented in Fig. 4 may be correlated. Therefore it is necessary to study the statistical significance of results in these plots. Regarding this purpose, we investigate the likelihood function, where r σ,k is the ratio r σ at the input coupling g 2 GF = u k , andr σ,k is the central value of r σ,k in our numerical computation. The elements of 2 × 2 covariant matrix, Cov ij , can be determined using the bootstrap samples of r σ,i and r σ,j in the numerical calculation. Figure 5 displays the likelihood functions for the study of the correlation between r σ at input g 2 GF = 6.0 and several choices of r σ,k , at c τ = 0.45. It is clear that this ratio computed at input g 2 GF = 6.0 is at least mildly correlated with that extracted at input g 2 GF ≥ 5.3. Notice that we have at least two data points for every lattice volume between these two values of the renormalised coupling. This investigation shows the necessity of having simulations at many choices of the bare coupling for eachL, in order to reduce the correlation amongst r σ computed at different input renormalised couplings.

IV. STRATEGY FOR THE CONTINUUM EXTRAPOLATION AND FINITE-SIZE SCALING
As discussed in Sec. III B, implementation of the continuum extrapolation employing the fit formula of Eq. (34) is inspired by the "Symanzik-type" argument. This approach is applicable when the bare parameters are tuned such that the effects of the lattice spacing are only related to the UV Gaussian fixed point. In the present study, this is reached when g 2 0 is close enough to zero. Under this circumstance, a major origin of scaling violation are the irrelevant operators that can be included in the theory. The classical dimensional analysis is a good approximation in this region, and it leads to simple power-law dependence on the cut-off. For a generic observable, M latt , computed on the lattice, the approach to the continuum limit is governed by the behaviour where Λ i (i = 1, 2, . . . , N IR ) are all the possible IR energy scales that are well below 1/a. Clearly, M 0 is the continuum limit of M latt . Quantum fluctuations in the above equation can be accounted for by using perturbation theory. Because of the Gaussian nature of the fixed point, they introduce logarithmic dependence on the lattice spacing in the coefficients, M n,i . These logarithms are often discernible in numerical analysis only when very high-precision data are available, therefore they are normally not included in the fitting procedure.
To employ Eq. (37) for conducting the continuum extrapolation in search of an IRFP using the step-scaling method, it is essential to make certain that the dimensionfull lattice size, L, is in the long-distance region that is governed by possible IR conformality, while the scaling property of the theory at the lattice spacing is still dominated by the UV Gaussian fixed point. Therefore, to adopt this Symanzik-type continuum extrapolation for distinguishing between theories with IR scale invariance and slow-running behaviour, one may have to perform lattice simulations extremely close to the limitL → ∞. This is particularly crucial for the study of a theory that contains a small β−function, such that the UV and the IR scaling regimes can be separated by many orders of magnitude in the difference of scales. Since one normally works with the lattice size, in current step-scaling investigation of the running coupling, it is challenging to achieve this separation. Therefore, one has to be cautious when utilising the Symanzik-type strategy, Eq. (37), for confirming the existence of an IRFP.
Given the usual choices of the lattice size in Eq. (38), it is plausible that if an IRFP exists in the theory, the bare couplings can be tuned such that the scaling with respect to the change in both a and L is controlled by IR conformality. In fact, in all the contemporary lattice calculations employing the step-scaling for probing IR scale invariance in gauge theories [37,38,54,[74][75][76][77][78][79][80][81][82][83][84][85][86][87][88][89], the values of g 2 0 are often larger than unity. Therefore the continuum extrapolation in these computations (including our present work) may not be guided by the simple polynomial formula as in Eq. (37). Below we examine the alternative scenario in which the continuum limit is reached according to approximate IR conformality.
Near an IRFP at strong coupling, the classical dimensional analysis receives significant corrections from quantum fluctuations, and the cut-off dependence may no longer be as simple as Eq. (37). The anomalous dimensions of the operators in the theory can lead to dependence on fractional powers of a/L. Investigation for details of the scaling laws and the continuum limit near possible strong-coupling fixed points is not new in lattice field theory computations. Recent examples are the studies of the Higgs-Yukawa model in Ref. [90], and the three-dimensional scalar theory in Ref. [91]. Here we will first illustrate this point in the context of this work by examining a generic coupling, g R , renormalised at the length-scale ρ. In the vicinity of a strongly-coupled IRFP, the β−function can be well approximated by the linearised form, where g * is the location of the IRFP, and γ * is the slope of the β−function at this zero. Notice that the value of g * depends on the choice of the renormalisation scheme, while γ * is a universal quantity and takes real positive value. Integrating Eq. (39) between two length scales, l 1 and l 2 , we obtain which clearly indicates the possibility of having dependence on non-integer powers of l 1 and l 2 . For the purpose of our discussion, we introduce another scale, L ref , such that and work with fixed lattice spacing. To proceed, in the following discussion we will present our argument using the GF-scheme renormalised coupling,ḡ 2 latt (g 2 0 ,L), as defined in Eq. (14). Expressing all the length scales in lattice units, and identifying l 1 and l 2 in Eq. (40) with L ref and L, one obtains in the vicinity of the IRFP. This equation can be regarded as a finite-size scaling formula. Confronting it with lattice data enables us to confirm/exclude IR conformality, and it leads to the determination of g * and γ * . In addition to fixing the bare coupling, we can further choose to work at a particular value ofL ref in the analysis. It has to be stressed again that the renormalised couplings,ḡ 2 latt (g 2 0 ,L) andḡ therefore one has to work in a regime where lattice artefacts are small compared with the statistical uncertainties. This can be checked in practice by using different discretisations and/or different lattice sizes to extractḡ 2 latt .
We further notice, from Fig. 2, that in this work the change ofḡ 2 latt (g 2 0 ,L) is small when varyingL between 8 and 24 at fixed lattice spacing. When the coupling is very small, this is due to the effect of the Gaussian UVFP. At intermediate and strong couplings, such behaviour arises from the smallness of the β−function. Therefore, away from the asymptotic-freedom regime, we can fit our data, at a particular choice ofL ref and a, with the formula, where g l and γ are the free parameters, with the definition, Equation (43) can be regarded as the consequence of a "locally linearised" β−function, which is a good approximation only when one works with small variations of the coupling around g ref . This is the reason why g l and γ depend on g ref . Nevertheless, when the theory is tuned to be close to an IRFP, this equation must converge to Eq. (42), and g l and γ will approach constant values, g * and γ * .
For each fit, we specify a value for g 2 0 (hence g ref ), and extract g l and γ. When conducting this procedure in a region without IR conformality, g l and γ will show dependence on the input g ref . On the other hand, when the theory is engineered to be in the neighbourhood of an IRFP by tuning the bare coupling, these two quantities should show a clear trend to converge to g * and γ * . In summary, we can utilise our data and perform the fit to Eq. (43) at fixedL ref = 8, and scan though many values of g 0 (hence g ref ) in the strong-coupling regime. At each Ref. [38] Ref. [  choice of g 0 , we carry out a fit. If our data indicate the existence of an IRFP, both g l and γ should show plateau behaviour when plotted against g ref , and different discretisations forḡ latt will lead to consistent results. Figure 6 is the outcome for γ determined using the above analysis procedure at c τ = 0.5. In this plot, we notice that the largest g ref for the clover discretisation is slightly smaller than the largest input g GF (L) value which is 5.8 for the same c τ in the step-scaling analysis presented in Sec. III B. This is becauseḡ latt (g 2 0 ,L) grows withL at fixed g 2 0 as a general trend in our data, and we follow the principle that no extrapolation in g 2 0 is implemented in this work. For the finite-size scaling test discussed in this section, the largest g ref must be chosen to be the value ofḡ latt (g 2 0 ,L = 8) at the largest g 2 0 where we have data for theL = 24 lattice. Therefore, according to Table I, the maximal bare coupling in the analysis leading to the result in Fig. 6 is g 2 0 = 1.442. On the other hand, the largest bare coupling for computing the input g GF (L) in the step-scaling analysis is the maximal value of g 2 0 for theL = 16 lattice. It can be seen in Table I that this is at g 2 0 = 1.449. This small difference in g 2 0 can produce minor but visible difference in the renormalised coupling, since in this regimeḡ latt increases rapidly with g 2 0 , as demonstrated by the plots in Fig. 2.
From the plot in Fig. 6, it is obvious that results from the clover and the plaquette discretisations are not compatible, and there is no plateau in the strong-coupling region. This is consistent with our previous analysis, namely that the theory is not governed by IR conformality for the values of coupling probed in our simulations.

V. COMPARISON WITH PREVIOUS WORKS
In this section, we compare our result to previous lattice step-scaling investigations of SU(3) gauge theory with twelve flavours. In Refs. [38,54] this was carried out using the SF and TPL schemes, respectively. These two calculations made use of the same lattice actions, namely, the Wilson plaquette gauge action and unimproved staggered fermions. As summarised in Table III , while it was claimed in Refs. [38,54] that the theory can be IR scale-invariant, here we do not see compelling evidence for this conclusion. In comparison with these two previous lattice studies, several aspects of the computation have been improved in our current project. First, the maximal lattice size and the bare coupling in this work are larger than those in Refs. [38,54]. In principle, this allows us to probe the theory at greater length scale. Secondly, in the present calculation, the statistical error for the renormalised coupling is at the subpercentage level, while it is around or bigger than 2% in the two earlier works. Finally, the use of the GF scheme enables us to obtain our result with two different lattice discretisations, making it feasible to estimate the systematic effects arising from the continuum extrapolation. Such procedure is not possible in Refs. [38,54], since there are no alternative discretisations for the observables employed for determining the renormalised couplings. In fact, it was demonstrated that the TPL-scheme coupling computed in the continuum limit in Ref. [54] is unreliable, upon adding lattice data atL = 24 with similar values of the bare coupling [55].
It has to be stressed that in Table III, the results from Refs. [38,54] are both obtained using the Symanzik-type ansätz for the continumm extrapolation. As was already pointed out in Sec. III B, and discussed in detail in Sec. IV, one has to be cautious when using this approach to confirm the existence of an IRFP with lattice simulations.
We notice that the authors of Ref. [35] have performed a lattice computation for the GF-scheme coupling for twelve-flavour SU(3) gauge theory, employing a procedure that is similar to the step-scaling method. In Ref. [35], the fundamental-adjoint plaquette gauge action and the nHYP-smeared [92,93] staggered fermions are used. For the analysis procedure, the clover discretisation has been adopted to extractḡ 2 latt (g 2 0 ,L), as defined in Eq. (14) in the current paper, with the normalisation factor,N , calculated using the continuum perturbation theory. In the strong-coupling regime, these authors search for the intersections of pairs of curves representingḡ 2 latt (g 2 0 ,L) as a function of g 2 0 , atL and sL with s being the step size. Such intersections are interpreted as the consequence of IR scale invariance. The values ofḡ 2 latt at these intersections are then extrapolated to the limit of vanishing a/L with the Symanzik-type ansätz, and the result is regarded as the location of the IRFP in the continuum limit. In Ref. [35], it is claimed that this IRFP is reached at g 2 GF ∼ 7. This procedure cannot be implemented in the present work, because the above pairwise intersections are not observed in our data for c τ ≥ 0.45, where we have the lattice artefacts under control.
The conclusion in Ref. [35] does not contradict the result of this work. Our main observation is that the scaling behaviour of the GF-scheme coupling in SU(3) gauge theory with twelve flavours is not governed by IR conformality at g 2 GF ∼ 6. Of course the effects of the possible IRFP can appear at g 2 GF > 6. This is beyond the scope of this project. On the other hand, it will be interesting to examine the reliability of the procedure in Ref. [35] by carrying out the same analysis using the plaquette discretisation.

VI. CONCLUSION
In this article, we present our step-scaling analysis of the coupling constant in SU(3) gauge theory with 12 massless flavours, using the Gradient-Flow scheme. In this theory the β−function is very small, such that doubling the length scale induces at most 6% variation in the renormalised coupling according to two-loop perturbation theory. Therefore, to make any statistically-meaningful statement regarding possible IR conformality in this theory, it is desirable to have lattice data with error at the subpercentage level for the extracted renormalised coupling. To our knowledge, our work is the first computation that achieves such precision.
It is well known that the continuum extrapolation is the main source of the systematic error in the step-scaling approach. The implementation of the Yang-Mills gradient flow reduces the cut-off effects in our calculation. In this project, we obtain the renormalised coupling via the computation of the gauge field energy density using two different lattice discretisations, namely the plaquette and the clover operators. Such strategy enables the estimate of the systematic error arising from lattice artefacts. We find that at large enough flow time, such that this extrapolation is under control.
Being able to have good control of both statistical and systematic errors, we manage to demonstrate, in a statistically-meaningful manner, that the theory flows out of the vicinity of the UV Gaussian fixed point at g 2 GF ∼ 2, and the running of the coupling begins to be significantly slower than the two-loop perturbative prediction around g 2 GF ∼ 5. At c τ = 0.5, our result indicates that the ratio r σ = g 2 GF (2L)/g 2 GF (L) is almost consistent with unity at g 2 GF (L) ∼ 5.8. This conclusion is reached using the Symanzik-type formula in performing the continuum extrapolation.
In this paper, we discuss the application of the continuum extrapolation ansätz a'la Symanzik in the search for possible IRFP through the step-scaling approach. This Symanzik-type method is based on the scenario that the property of the theory at the cut-off scale is governed by the UV Gaussian fixed point, while its scaling behaviour at the lattice size can be dominated by IR conformality. This is obviously very challenging to achieve in practice. To confirm the existence of the IRFP, we argue that it is essential to examine the theory with the assumption that the scaling of the theory at the lattice spacing is also determined by IR scale invariance. Following this argument, we perform a finite-size scaling test of SU(3) gauge theory with twelve flavours. The result of this test indicates that the behaviour of the theory, as probed using our lattice data, is not governed by possible IR conformality. That is, our result does not support the existence of an IRFP in this theory in the region g 2 GF (L) ≤ 6. In summary, lattice computations for the determination of the conformal windows for various gauge theories have matured significantly in recent years. The importance of controlling errors in these calculations is receiving growing attention. The work presented in this paper is our first attempt in this research avenue with high-accuracy lattice data. To make further progress in this direction, it would be desirable to implement the procedure with improved actions in future lattice simulations.  (10) 0.890(9) 11.15 0.752 (4) 0.770(4) 12.00 0.675 (4) 0.692(4) 13.85 0.557 (2) 0.570 (2) (2) 0.730(2) 14.00 0.569 (2) 0.590(2) 16.00 0.473 (2) 0.491(2) 18.00 0.412 (1) 0.427(2) 20.00 0.366 (2) 0.380 (2) TABLE VI: Raw data for the renormalised couplings extracted using the clover and the plaquette discretisations atL = 12 and 24 with cτ = 0.5.