Lattice study of infrared behaviour in SU(3) gauge theory with twelve massless flavours

We present details of a lattice study of infrared behaviour in SU(3) gauge theory with twelve massless fermions in the fundamental representation. Using the step-scaling method, we compute the coupling constant in this theory over a large range of scale. The renormalisation scheme in this work is defined by the ratio of Polyakov loops in the directions with different boundary conditions. We closely examine systematic effects, and find that they are dominated by errors arising from the continuum extrapolation. Our investigation suggests that SU(3) gauge theory with twelve flavours contains an infrared fixed point.


I. INTRODUCTION
The origin of electroweak (EW) symmetry breaking is one of the most important research topics in physics. With the progress of experiments at the Large Hadron Collider (LHC), it is urgent for a theoretical understanding for the mechanism of the mass generation and its relation to EW symmetry breaking. One appealing scenario for this mechanism is the technicolour models [1,2]. These models involve new asymptotically-free gauge theories in which the coupling constants become strong at the TeV scale. The strong coupling can induce condensates to generate mass gaps, and asymptotic freedom leads to the absence of the naturalness problem. In order to dynamically suppress the flavour-changing neutral currents (FCNC), and to evade the constraints from precision EW data, it is important that the candidate theories exhibit the "walking" (quasi-conformal) behaviour and contain large anomalous dimension for the technifermion mass term [3][4][5].
In recent years, there has been a significant amount of work in search of gauge theories viable for walkingtechnicolour model building. The most important task in this endeavour is the determination of the critical number of massless fermions, given the gauge group and the fermion representation, above which a theory is conformal in the infrared (IR). For theories involving fermions in the fundamental representation, this is denoted as the critical number of flavours, N cr is the number of flavours above which asymptotic freedom is lost), the theory contains an infrared fixed point (IRFP). A candidate walking-technicolour theory with fundamental fermions is believed to have the number of flavours just below N cr f . This makes the determination of N cr f a task with phenomenological significance, in addition to its importance in field-theoretic studies. Since the couplings must be strong at low energies in these theories, nonperturbative methods, such as the Schwinger-Dyson equation and gauge-gravity duality, have to be employed. Amongst these, lattice gauge theory is the only first-principle tool, and has been applied by many groups in this research avenue   1 .
Of all the theories which have been investigated using the lattice technique, the value of N cr f for SU(3) gauge theories with fundamental-representation fermions remains a controversy. Although several groups [7-11, 26, 33, 46] found evidence that SU(3) gauge theory with N f = 12 is conformal in the IR, authors of Refs. [34,37] argued that chiral symmetry is broken in this theory. In this paper, we report our study of this theory, using the step-scaling method to compute the running coupling constant. We adopt the Twisted Polyakov Loop (TPL) scheme [49][50][51]. This article complements the letter [6] which was released in 2011 with other colleagues on this collaboration, and contains more details of our simulations and improved analysis using more data. In Ref. [6], we concentrated on the analysis with the step size, s, set to 1.5, while here we emphasise the case in which s equals two. Furthermore, in the improved analysis with new data, as presented in this paper, we significantly reduce the correlation between data for the step-scaling functions on different lattice volumes. This makes the continuum extrapolation simpler and better controlled, compared to the analysis published in Ref. [6]. We will discuss this in detail in Sec. V C. Related to this work and Ref. [6], we have also published conference proceedings [52][53][54], as well as for a similar project on SU(2) gauge theory with eight flavours [55].
In addition to computing the running coupling constant, we also obtain the ratio between the step-scaling function and the coupling constant, which becomes one when the β−function is zero. To claim the discovery of the IRFP in an asymptotically-free gauge theory, we have to demonstrate that this ratio is indeed one in the ultraviolet (UV) and the IR, while being obviously different from this value between these two regimes. Our study suggests that SU(3) gauge theory with N f = 12 contains an IRFP around the TPL-scheme coupling constant, Amongst systematic effects that we estimate, errors arising from the continuum extrapolation dominate. We also notice that some of our procedures in performing this extrapolation lead to weaker evidence for the existence of the IRFP. Details of the estimation of systematic errors will be presented in Secs. V and VI.
Our finding for the evidence of the existence of the IRFP agrees with the result of Refs. [7,8], where the Schrödingerfunctional (SF) scheme [56,57] was used in defining the coupling constant, and the calculation was performed using the same gauge and fermion actions. The values of g * are different because of scheme dependence. Here we also stress that the lattice artefacts can be very different in these two schemes. In particular, the SF scheme contains O(a) (a is the lattice spacing) lattice artefacts through the introduction of the boundary terms 2 , while in the TPL scheme the lattice artefacts remain of O(a 2 ), making the continuum extrapolation more reliable.
This paper is organised in the following way. In Sec. II, we review twisted boundary conditions and the Twisted Polyakov Loop scheme. Section III contains the details of our simulation strategy and parameters. We describe our analysis procedure in Sec. V, give our results and discussion in Sec. VI, and conclude in Sec. VII. Appendix A contains the study of the eigenvalue spectrum of the Dirac operator used in this work. Values of plaquette and the raw data for the TPL-scheme coupling constants are presented in App. B.

II. TWISTED POLYAKOV LOOP SCHEME
In this section, we give the details of our definition of the renormalised coupling constant in the twisted-Polyakovloop (TPL) scheme [50,51]. This scheme makes use of twisted boundary condition (TBC) [58], which is implemented on the link variables, U µ (n) (µ = x, y, z, t is the Lorentz index andn is the position of a lattice site), through where L ν is the (dimensionful) box size in the ν direction (withν denoting the unit vector), and a is the lattice spacing. The "twisting matrices", Ω ν , act in the colour space. In this work, we apply TBC for ν = x, y, while maintaining periodic boundary condition (PBC) for the other two directions. This means For SU(3), the twisting matrices, Ω x,y , satisfy In this work, we explicitly implement [59] The inclusion of fermions is not straightforward when TBC, Eq. (2), is imposed on gauge fields. In order to maintain gauge invariance and single-valuedness of the fermion field, ψ(n +xL x /a +ŷL y /a), under the application of two boundary twistings (inx andŷ directions) with different orderings, it is necessary to introduce the "smell" degrees of freedom [60]. This quantum number is carried by fermions. The number of smells, N s , is equal to the number of colours, N c . Twisted boundary condition on fermion fields is given by, where a and b are colour indices, and the twisting matrices, Ω µ , have been generalised to act on the smell degrees of freedom (indices α and β). The factor e iπ/3 is introduced only for ν = x, y, to remove the zero-momentum modes in these directions. For ν = z, t directions, we implement ordinary PBC, ψ(n +νL ν /a) = ψ(n). Since the smell quantum number is not carried by the gauge fields, it can be considered as additional flavours. Therefore the number of flavours in simulations involving dynamical fermions with TBC has to be a multiple of N s (=N c ).
The Polyakov loops in the twisted directions, ν = x, y, are P x (n x ,n y ,n t ) = Tr  Ω x e i2πny a/(3Ly)   , P y (n x ,n z ,n t ) = Tr The extra factors outside the square brackets are introduced to maintain gauge and translation invariance. The renormalised coupling constant can be defined via the ratio between correlators of Polyakov loops in the twisted and periodic directions, where P z,t are ordinary Polyakov loops in the directions with PBC. In this study, we always use hypercubic lattice L x /a = L y /a = L z /a = L t /a = L/a. The proportionality factor k can be extracted by computing the above ratio in perturbation theory to O(g 2 0 ), where g 0 is the bare coupling constant. Using lattice perturbation theory, one obtains the lattice version of this factor [6], The coupling,ḡ latt , defined in Eq. (8) contains lattice artefacts, therefore depends on the lattice spacing as well as the volume. Its continuum-limit counterpart at fixed physical volume is defined as, The TPL scheme, as defined in Eq. (8), contains the feature that the renormalised coupling constant has the fixed value 1/k ∼ 5.6 in the IR limit (L → ∞). Therefore, in order to firmly establish the existence of the IR fixed point, we have to show thatḡ c is significantly different from this value at the fixed point.
Contrary to the SF scheme, the O(a) lattice artefacts are absent in the TPL scheme. As explained in the following sections, it is important to control the continuum extrapolation in the step-scaling study of the running coupling constant. This makes the use of the TPL scheme very desirable. In Sec. V, we will show the lattice-spacing dependence of the TPL-scheme coupling constant.

III. SIMULATION SETTING
We give the details of our lattice simulation in this section. As discussed in Sec. II, the number of flavours in our calculation must be a multiple of N s = N c = 3. Since we are using staggered fermions, it also has to be proportional to the number of tastes, N t = 4. In this work, we investigate the SU(3) gauge theory coupled to twelve flavours, which is allowed by these constraints.

A. Step scaling
Our goal is to measure the evolution of the running coupling constant over a wide range of scale. Given that the lattice imposes infrared (the volume) and ultraviolet (the lattice spacing) scales, the most convenient way to achieve this goal is the step-scaling technique. In this approach, we first measure the renormalised coupling constant,ḡ latt , on the lattice in the scheme defined in Eq. (8). Since we perform computation at vanishing fermion mass,ḡ latt only depends on the lattice spacing and the lattice volume, L/a. Choosing a few values of L/a, we then simulate at a wide range of β ≡ 6/g 2 0 , where g 0 is the lattice bare coupling constant. This enables us to tune β (lattice spacing) to obtain the renormalised coupling in the continuum limit, g c (L) =ḡ latt (β 1 , L/a 1 ) =ḡ latt (β 2 , L/a 2 ) = . . . =ḡ latt (β n0 , L/a n0 ) , where n 0 is the number of choices of L/a. Sinceḡ c is independent of the lattice spacing, it is renormalised at the length scale L. In this work, we perform lattice simulations at L/a = 6, 8, 10.
Using the combinations of (β, L/a) which lead to the sameḡ c (L) (or u =ḡ 2 c ), we compute the lattice step-scaling function, where i = 1, 2, . . . , n 0 as in Eq. (11), and s is the step size. Since we can obtain n 0 results for Σ at the same physical volume, L, with different lattice spacings, this allows us to determine the continuum-limit step-scaling function, In this work, we choose the step size s = 2, leading to the need for simulations performed on the lattice volumes, sL/a = 12, 16, 20, with s = 2.
For convenience, we define The step-scaling function is a scheme-dependent quantity, since it is simply the renormalised coupling constant computed at a certain scale. To facilitate a better method in demonstrating the existence of the IRFP, we compute the ratio This ratio becomes one at the zeros of the β-function. The existence of such zeros is independent of the renormalisation scheme used in the calculation. In order to show that the gauge theory under investigation does contain an IR fixed point, we have to verify that r σ (u) is one at both UV and IR regimes, while deviating from this value in between.
A major source of systematic errors in the step-scaling method is the continuum extrapolation. It is a challenging task to properly address this issue. In order to have more information regarding this extrapolation and its possible systematic effects, we also perform simulation with L/a = 14, and resort to an interpolation procedure to obtain data for the TPL-scheme renormalised coupling on the lattice size L/a = 7. This enables us to carry out the investigation with (L/a = 6, 7, 8, 10) −→ (2L/a = 12, 14, 16, 20) .
Here we stress that staggered fermions are used in this work, therefore it is not possible to have data directly on the L/a = 7 lattice. The interpolation procedure for obtaining such data is explained in detail in Sec. V C. This interpolation in volume can introduce systematic effects, although it may result in more information regarding the continuum limit. Therefore, we only use the 4-point step-scaling analysis in Eq. (18) as a means to estimate errors in the continuum extrapolation.

B. Details of simulation parameters
Our calculation is performed using the Wilson plaquette action for the gauge fields, and unimproved staggered fermions. We implement the standard Hybrid Monte Carlo (HMC) algorithm using the Omelyan integrator with multi-time steps [61,62]. To compute the inversion of the lattice fermion operator, biCGstab solver with convergence condition that the residue is smaller than 10 −16 for molecular dynamics, and the accuracy of 10 −24 for the Metropolis tests, are used. A significant fraction of of our simulations were carried out on Graphics Processing Units (GPU's), where a mixed-precision solver with defect correction was implemented. The GPU codes were developed with CUDA [63].
To thermalise configurations in the Markov chains, we have used two procedures. In the first procedure, we start a simulation from a trivial gauge-field configuration U µ (n x ,n y ,n z ,n t ) = 1, with fermion mass Then, we gradually decrease the mass to zero. In this process, we monitor the Polyakov loops in the untwisted directions, and make certain that the imaginary parts are non-vanishing. This ensures that the Markov chains progress mostly near the true vacua [6]. In the second procedure, we start with a configuration, U µ (n x ,n y ,n z ,n t ) = 1 elsewhere, which always results in non-zero imaginary parts in the Polyakov loops in the untwisted directions. It also produces the largest gap in the vicinity of zero in the fermion matrix. In this case, we can start the simulation directly with zero fermion mass, making this procedure significantly more efficient than the one implemented with the initial conditions of Eqs. (19) and (20). In both cases, we observe that the simulations always stay near the true vacua, and tunnelling amongst these vacua occur occasionally. We will discuss this issue in more detail in Sec. IV B.
In order to implement the step-scaling investigation of the running coupling constant as discussed in Sec. III A, we carry out simulations at the lattice volumes, L/a = 6, 8, 10, 12, 14, 16, 20.
For each volume, we simulate at several β values between 4 and 99, in the gauge action. Since the running is expected to be slow in SU(3) gauge theory with twelve flavours, this large range of β is necessary to trace the coupling constant from the UV to the IR regimes. We aim at determining the Polyakov-loop correlators with statistical errors around 2.5% or smaller. For this purpose, a significant amount of gauge-field ensembles have to be generated. The raw data for the TPL-scheme renormalised coupling, as defined in Eq. (8), are given in App. B.

IV. PLAQUETTE, POLYAKOV LOOP AND THE VACUUM STRUCTURE
In this work, we perform several detailed checks on the simulations, in order to ensure that we are estimating the autocorrelation and performing the continuum extrapolations reliably. These checks include the lowest-lying eigenvalue spectrum of the Dirac operator, the plaquette values, and the phases of the Polyakov loops. The computation of the lowest-lying eigenvalues is presented in App. A, while in this section we address the other two topics.

A. Plaquette
As shown in App. B, some of our simulations are performed at small β values (coarse lattice spacings). It is necessary to check that these simulations are still in the weak-coupling phase, in order to make certain that at these β values, the theory is still in the same universality class as that with high−β (fine lattice spacings). This is essential in order to ensure that the continuum limit can be reliably taken in our calculations. For the above purpose, we examine the expectation values of the plaquette for many of our HMC simulations. The results are summarised in Table III in App. B. These expectation values are plotted in the Upper panel of Fig. 1, where we also show the predictions from the weak coupling expansion for pure Yang-Mills theory, By comparing our data with this function, it is evident that all our simulations are in the weak-coupling phase, and are safe from being in the novel phase observed in Ref. [19] 3 . We have also studied the volume dependence of the plaquette, by first fitting the data obtained on the largest lattice, L/a = 20, to a weak-coupling expansion formula (p i are the fit parameters), then computing the difference between the data points to this curve. The result of this investigation is shown in the lower panel of Fig. 1. This shows that finite-size effects are minor in the computation of the plaquette in this work.

B. Polyakov loops and vacuum tunneling
The study of the plaquettes in the last section confirms that our simulations have been carried out in the weak coupling phase. In this phase, as pointed out in Ref. [6], the true vacua in SU(3) gauge theory with fermions are always those in which the vacuum expectation values of the Polyakov loops in the untwisted directions are non-vanishing and complex. On the other hand, in the vicinity of the false vacua, the untwisted Polyakov loops are real. As for the Polyakov loops in the twisted directions, we expect that they will scatter around zero, configuration by configuration.
Markov chains in our simulations can be trapped in the false vacua. However, by using the above property of the untwisted Polyakov loops, we can monitor the simulations and ensure that they are mostly progressing near the true vacua.
Investigating the Polyakov loops trajectory by trajectory, we first confirm that in the twisted directions, they are fluctuating around zero for all simulations in this work. This is shown for two typical cases in the upper panels of Figs. 2 and 3.
Next, we study the Polyakov loops in the untwisted directions. In all our simulations, their values are non-vanishing and complex in all trajectories. The complex phase fluctuates around ±2π/3, indicating that the Markov chains are progressing near the true vacua. The lower panels of Fig. 2 demonstrate a case (L/a = 16, β = 11.15) in which the simulation stays near the vacuum with the phase of Polyakov loop being −2π/3. For simulations performed at smaller L/a (fewer total degrees of freedom) and larger β (stronger coupling), tunnelling between the two true vacua may occur. One of such cases is shown in the lower panels of Fig. 3. Every time this takes place, we then investigate the Polyakov loop correlators trajectory by trajectory, ensuring that these correlators do not exhibit any "discontinuous" behaviour when the tunnelling happens. In Fig. 4, we show the result of this study for the corresponding simulation presented in Fig. 3. From these plots for the Polyakov loop correlators in the twisted and untwisted directions, we conclude that tunnelling between the true vacua does not result in artefacts which complicate the estimation of autocorrelation time.

V. ANALYSIS DETAILS
In this section, we explain the details of our analysis. The statistical analysis in this work is performed using the bootstrap procedure, in which 1000 bootstrap samples are generated for each (L/a, β).

A. Autocorrelation and data binning
As presented in App. B, we perform our calculations with a large number of HMC trajectories. The first step in our analysis is the binning of the raw data. In order to make certain that the binning procedure is reasonable, we study the autocorrelation of the ratio, appearing in the left-hand side of Eq. (8), between the Polyakov loop correlators.
To describe our investigation, we start from the autocorrelation function of primary quantities, [64][65][66] Here,α andβ label the types of primary quantities. In our case, O 1 (i) and O 2 (i) are Polyakov loop correlators of the i-th sample in the twisted and in the periodic directions, respectively. The quantityŌα is the average of Oα(i),Ōα = (1/N ) N i Oα(i). By using Γαβ, the autocorrelation function of the Polyakov loop ratio, as in the left-hand side of Eq. (8), can be written as, with, We define the normalised autocorrelation function, which is normally assumed to behave as, The quantity τ A is the autocorrelation time of single exponential autocorrelation.
Since the integrated autocorrelation function is less noisy than ρ(τ ), we use it to estimate the autocorrelation time between Polyakov-loop-correlator ratios. Upon integrating over τ , we obtain The single-exponential form in Eq. (29) is often a poor approximation to ρ(τ ), when the system contains degrees of freedom that are characterised by very different autocorrelation times. In general, the autocorrelation function can be multi-exponential, The integrated autocorrelation is This function reaches a plateau k τ A for all k. We use this criteria for the estimation of autocorrelation without explicitly determining τ (k) A and a k . A more detailed study of autocorrelation times for conformal field theories will be reported in a separate paper [67].
In our numerical calculation, the integrated autocorrelation is defined as, To estimate error in Θ(τ ), we apply the Madras-Sokal formula [68], Figure 5 shows Θ(τ ) for the representative cases in this work. The separation between two decorrelated trajectories can be estimated by investigating the plateau of Θ(τ ). As demonstrated in Fig. 5, this separation depends on the physical volume, L. It is around 20 on the smallest volumes, and about a few hundred to 1000 on the largest volumes.
In App. B, we show the details for the numbers of HMC trajectories in our simulations. For each choice of (L/a, β), we divide the trajectories evenly into ∼ 200 bins by averaging over them in each bin. These bins are then used to create 1000 bootstrap samples. From the result presented in this section, it is evident that our bin sizes are large enough compared to the autocorrelation times. This ensures that the data amongst these bins are decorrelated.
We have also confirmed this with the Jackknife analysis using these and larger bin sizes. The statistical errors in this approach are almost the same as those in our bootstrap analysis, and they are stable against the change of the bin sizes. group is the result of using 200 bins, as chosen in our bootstrap procedure. The β values for the two blue points are slightly shifted for the purpose of presentation. They correspond to choosing the numbers of bins to be 100 (left) and 500 (right). From these plots, it is apparent that having 200 bins leads to enough trajectories in each bin, in order to correctly estimate statistical errors.

B. Interpolation in β (bare coupling constant)
In the step-scaling study of the running coupling constant, we first have to perform the tuning of the β values in Eq. (11), for the lattice volumes L/a = 6, 8, 10. In principle, this can be achieved by repeatedly adjusting β and carrying out new simulations, until Eq. (11) is satisfied to high accuracy. As discussed at the end of Sec. III A, we also want to obtain the TPL-scheme renormalised coupling on the L/a = 7 lattice through interpolation in volume, in order to estimate systematic errors in the continuum extrapolation. For this purpose, one has to tune a different set of β values for L/a = 6, 8, 10 and interpolate to L/a = 7 at each step of this tuning.
The above procedure is very time-consuming, and becomes impractical for studies in which one has to trace the coupling constant across a large range of length scale. This is the case in the current work. Therefore we resort to a variation of the above method. That is, we simulate at many β values for each L/a, and perform interpolations in β for the renormalised coupling constant, volume by volume. The choices of these β values are presented in App. B. The use of this interpolation method inevitably introduces systematic effects in our calculation. We will address this issue in this section.
Since we are simulating at a large range of bare coupling constant, it is a challenging task to have a well-inspired interpolation function in β. One reasonable way to proceed is to note that in the large−β (small bare-coupling) regime, one-loop perturbation theory has to be valid, and therefore at fixed L/a, u latt ≡ḡ 2 latt (β, L/a) ≈ where g 0 is the bare gauge coupling. This motivates the use of polynomial functions in 1/β to perform the interpolation. Since we have data for many β values (see App. B) for each L/a, it is in principle possible to have high degrees of polynomials for these fits. Such high-degree polynomials will generally fit all the data points. On the other hand, the Runge phenomenon may occur in this procedure, resulting in artificial oscillatory behaviour of the fit functions. In order to avoid this artefact in the β−interpolation, we note that the renormalised coupling should always be non-decreasing with growing lattice spacing (i.e., decreasing β) at fixed L/a, otherwise the theory will be in the strong-coupling phase and the continuum extrapolation cannot be reliably performed. From our study of the plaquettes in Sec. IV A, it is evident that our simulations are all carried out in the weakcoupling regime. This is also reflected on the data points plotted in Fig. 7, in which we see that all our renormalised couplings, u latt , are non-decreasing when β decreases. This leads to the use of the non-decreasing polynomial, in the β−interpolation procedure at fixed L/a. We implement the constraint from perturbation theory, Eq. (35), which results in, This constraint leads to the number of fit parameters, where N deg and N h are defined in Eq. (36). The use of the non-decreasing polynomial ansatz makes the Runge phenomenon milder compared to the simple polynomial fits. The inverse of the fit function in Eq. (36) is also single-valued. This is essential in the step-scaling method. The results of applying this (uncorrelated) fitting procedure in the β interpolation are shown in Fig. 7. The optimal choices of N param , leading to the best (smallest) χ 2 /d.o.f., are listed in Table I.
In order to estimate systematic error resulting from the interpolation in β, we change the fit function from Eq. (36) to a simple polynomial function, with the constraint,c from the validity of perturbation theory at high−β. This constraint results in the number of fit parameters, The values ofÑ param for the best χ 2 /d.o.f. are presented in Table I.
As indicated at the end of Sec. III A, it is desirable to gain more information regarding systematic errors in the continuum extrapolation. In view of the fact that our reference input renormalised couplings are computed on L/a = 6, 8, 10, a practical way to proceed is to have data for L/a = 7. This enables us to attempt the step-scaling study, (L/a = 6, 7, 8, 10) −→ (sL/a = 12, 14, 16, 20), where s = 2, without having to perform simulations on large lattices, such as L/a = 24.
Since staggered fermions are used in this work, we have to use an interpolating procedure to obtain u latt for L/a = 7. To have a well-motivated method for this interpolation, we resort to the β−function of the theory. It is well-established that the coupling constant in SU(3) gauge theory with twelve flavours runs slowly compared to, e.g., QCD. This is reflected on the fact that a small change in the renormalised coupling has to result from a significant variation of the scale. As shown in Fig. 7, this is indeed the case. Namely, enlarging the box size by a factor of two induces very little changes in u latt , and one can locally approximate the β−function using a linear form where a l and b l are unknown parameters. We stress that this approximated form is not based on perturbation theory, and is only valid within a small range of u latt . That is, in different ranges of u latt , the parameters, a l and b l , have different values.
To determine u latt on the L/a = 7 lattice, we use our data on the L/a = 6, 8, 10, 12 lattices, and interpolate with the function,  at fixed lattice spacing. The unknown coefficients, A L , B L , and C L are related to a l and b l , and the integration constant in solving Eq. (43). Figure 8 shows two representative plots for the interpolation using Eq. (44). It is obvious that the interpolation is smooth, and the values of the coefficients, A L , B L , and C L , can vary significantly in different ranges of u latt . The fits presented in Fig. 8 are performed on the L/a = 6, 8, 10, 12 data without β−interpolation.
Equations (43) and (44) are used to motivate an interpolation function in L/a at fixed a (β value). However, the the effects of the lattice spacing can appear as powers of (a/L) 2 in our simulations. This means the data points used in each of this volume interpolation may have different lattice artefacts, leading to systematic effects introduced in this procedure. In view of this, we do not include the L/a = 7 data in our central analysis procedure, and only use them to perform the step-scaling investigation in Eq. (42) as a means to estimate errors in the continuum extrapolation.
Another issue in this volume-interpolation method for obtaining the L/a = 7 data is statistical correlation. The procedure is carried out using uncorrelated fits in this work. However, it is natural to expect that there will be correlation between the L/a = 7 (interpolated) data and those extracted directly from independent simulations on L/a = 6, 8, 10, 12. This correlation has to be closely examined, since all these data are used in the investigation of the continuum extrapolation, as discussed in Sec. V D. For this purpose, we study the likelihood function, where u i denotes u latt computed on the lattice volume L/a = i. Here u i is kept as a variable, andū i is its central value for this quantity from our simulation. The symbol Cov is the covariance matrix which can be computed from the bootstrap samples of u i and u j obtained from numerical calculations.
Our investigation shows that, although the coupling constant on the L/a = 7 lattice is interpolated using those on the L/a = 6, 8, 10, 12 lattices, it only shows significant correlation with that on the L/a = 8 lattice. In Fig. 9, we display an example of this likelihood-function study performed for β = 5.53. It is obvious from this figure that coupling constants on L/a = 7 and L/a = 6, 10, 12 exhibit very small correlation, while it is the opposite between L/a = 7 and L/a = 8. The corresponding covariance matrices for the example in Fig. 9 are The volume, L/a = 7, is one of the "small" lattices, on which we compute the reference coupling instead of the step-scaling function. The importance of the above study is the demonstration that there is negligible correlation between data on this lattice and that on the "large" lattice, L/a = 12, from which we compute the step-scaling function. Furthermore, the statistical errors of the data obtained on all our small lattices are small. In view of this, it is reasonable to expect that this correlation between the L/a = 7 and L/a = 8 TPL coupling constants does not necessitate correlated fits in the continuum extrapolation.
The above study of the data correlation also leads to the conclusion that one has to be very cautious about interpolating u latt in L/a. In certain analysis procedures, such as the one we adopted in Ref. [6] by setting the step size to 1.5, large correlation amongst data used in the continuum extrapolation can occur.

D. Continuum extrapolation for the step-scaling function
The last step in our analysis is the continuum extrapolation for the step-scaling function, σ(u), defined in Eqs. (14) and (16). Since unimproved staggered fermions and the Wilson plaquette action are used in this work, we will investigate (a/L) 2 dependence in the lattice step-scaling function, Σ(β, L/a, u, s = 2). In Fig. 10, this dependence is displayed at representative values of u in the regimes of weak, intermediate and strong coupling. From this figure, it is obvious that effects of the lattice artefacts grow with increasing u, as expected. In the region u < 0.8, we see that the step-scaling functions show insignificant dependence on the lattice spacing, and are almost consistent with the input reference coupling. On the other hand, in the strong-coupling regime, the a−dependence in Σ becomes noticeable, necessitating good control of the continuum extrapolation in the investigation of the existence of the IRFP. It is worth noting that the lattice artefacts tend to make the step-scaling function larger than its continuumlimit counterpart, especially in the strong-coupling regime. This feature is different from what was discovered in the Schödinger-functional scheme [7,8].
In performing the continuum extrapolation for our central analysis procedure, we use our simulation results for Σ(β, L/a, u, s = 2), obtained at sL/a = 12, 16, 20, and carry out the linear fit (σ l (u) and A l are the fit parameters), with the β−values for various L/a determined by tuning the coupling, u, to be the same on the corresponding small lattices (L/a = 6, 8, 10). This procedure does not include the L/a = 7 data which are extracted with an additional volume-interpolation, as detailed in Sec. V C.
To estimate systematic errors in the continuum extrapolation, we include the volume-interpolated, L/a = 7 data, as well as the step-scaling functions computed on the lattice, sL/a = 14. We first perform the quadratic fit (σ q (u), A q and B q are the fit parameters), to implement the 4-point step-scaling method in Eq. (18).
In order to further account for systematic effects arising from the continuum extrapolation, we perform two additional linear fits: 1. Using the data for the step-scaling functions from sL/a = 14, 16, 20 (L/a = 7, 8, 10).
2. Using the data for the step-scaling functions from sL/a = 12, 14, 16, 20 (L/a = 6, 7, 8, 10). Figure 11 shows representative plots of the continuum extrapolation using the above procedures (quadratic fit and the three linear fits). From these plots, we observe that σ l and σ q are well consistent with each other at intermediate and strong couplings. In the weak-coupling regime (top row of Figure 11), we notice that the quadratic fit, and the 3-point linear fit using the L/a = 7, 8, 10 data are not consistent with the other two procedures. They result in σ(u) smaller than u after the continuum extrapolation. However, we stress that in this regime, the lattice step-scaling function, Σ, demonstrates very mild lattice-spacing dependence, and is almost consistent with the input reference u. This is the consequence of asymptotic freedom. Furthermore, our data do not show significant O(a 4 ) contributions in the continuum extrapolation at strong and intermediate couplings (center and bottom rows of Figure 11), where the lattice artefacts are expected to be larger compared to the small−u region. In view of this, we conclude that the quadratic fit, and the 3-point linear fit using the L/a = 7, 8, 10 data can be artificially amplifying statistical fluctuations and leading to unreliable results in the weak-coupling regime. In order to properly address this issue, one has to generate data with very high statistical accuracy (e.g., < 0.5%) at large β values. This is beyond the scope of this work, since our main focus is on the existence of the IRFP in the strong-coupling regime.

VI. FINAL RESULTS AND DISCUSSION
In this section, we present the final results of our analysis, and discuss the estimation of systematic errors. We begin by showing the result from our central analysis procedure for the ratio r σ = σ(u)/u, defined in Eq. (17). In performing this central-procedure analysis, we first interpolate in the bare coupling, β, for simulation data obtained at L/a = 6, 8, 10, 12, 16, 20, using the non-decreasing polynomial function in Eq. (36) with the constraint from Eq. (37), and the polynomial degrees and the numbers of fit parameters presented in Table I. We then carry out the step-scaling of L/a = (6, 8, 10) −→ 2L/a = (12,16,20), (49) by extrapolating the step-scaling function to the continuum limit with the linear form in (a/L) 2 , Eq. (47). Result of this central analysis is shown in Fig. 12, which demonstrates evidence for the existence of an IRFP.   Next, we discuss the estimation of systematic effects arising from the β−value (bare-coupling) interpolation and the continuum extrapolation. For this purpose, we perform the changes in the central procedure. These changes are carried out independently, i.e., we vary one component in the central procedure, while keeping the other fixed.
We begin by varying the β−interpolation in the central procedure. This is carried out by changing the nondecreasing fit function in Eq. (36), to the simple polynomial form in Eq. (39)   the numbers of parameters reported in Table I. The result of this procedure is shown in the top-left panel of Fig. 13.
In order to estimate systematic errors associate with the continuum extrapolation, we perform various fits with the inclusion of the L/a = 7 and L/a = 14 data. First, we perform the quadratic fit using Eq. (48). This leads to the result for r σ (u) as depicted in the top-right panel of Fig. 13. As expected, this extrapolation strategy results in large statistical errors. In addition to the quadratic fit, we also carry out the two linear continuum extrapolations discussed in Sec. V D, L/a = (7, 8, 10) −→ 2L/a = (14,16,20), (50) L/a = (6, 7, 8, 10) −→ 2L/a = (12,14,16,20).
The result from the first these procedures is presented in the bottom-left panel of Fig. 13, while that from the second one is shown in the bottom-right panel of Fig. 13. As discussed at the end of Sec. V D, the quadratic fit, and the 3-point linear fit using the L/a = 7, 8, 10 data can lead to unreliable continuum extrapolations in the weak-coupling regime. Therefore, in Fig. 13 we only show results at intermediate and large u for these two procedures.
In Fig. 12, and in the top-left and the bottom-right plots in Fig. 13, it is observed that these procedures lead to r σ (u) consistent with one in the UV and the IR, while statistically different from this value between these two regimes. This suggests that there exists an IRFP in SU(3) gauge theory with twelve flavours. However, the continuum extrapolations using 4-point quadratic fit (the top-right plot in Fig. 13) and 3-point linear fit without the L/a = 6 data (the bottom-left plot of Fig. 13) lead to weaker evidence for the IR conformal behaviour. For these two procedures, in addition to the difficulty in the continuum extrapolations in the weak-coupling regime (discussed at the end of Sec. V D), we also observe large errors in the IR regime, leading to no apparent feature that r σ (u) crosses one. This phenomenon is actually the consequence of the "double crossing" behaviour in some bootstrap samples. Namely, in these samples, r σ (u) crosses one from above, and then turns around to cross the same value from below in a slightly larger u. We stress that out of 1000 bootstrap samples we have created in this work, r σ (u) in more than 680 (1σ) of them cross the unity from above in the IR regime, when the continuum extrapolations are performed with the quadratic fit or the 3-point linear fit without the L/a = 6 data. This leads to hints of the existence of an IRFP using these analysis procedures. In order to illustrate this point, in Fig. 14 we plot 100 bootstrap samples in the intermediate− and strong−u regions in these two procedures.
In Table II, we summarise the values of g 2 * obtained from the above procedures. Because it is challenging to precisely estimate systematic effects, as discussed above, we take a conservative approach to conclude that in SU(3) gauge theory with twelve flavours, our data suggest the existence of an IRFP around, This result is similar to what we obtained with other collaborators using a different analysis procedure [6] by setting the step size to be 1.5. Here we stress that it is more challenging to control systematic effects and the correlation amongst data points in the procedure in Ref. [6], because of the need for many interpolations in lattice volumes when computing the lattice step-scaling function, Σ.  The result in Eq. (51) is much smaller than that obtained in the SF scheme [7,8], The significant difference clearly indicates that the two schemes are very different. It should also be noted that in the TPL scheme, there is an upper bound for the renormalised coupling constant, as discussed in Sec. II. This may result in slower running behaviour compared to the SF scheme.

VII. CONCLUSION
In this paper, we present our work on the lattice study of IR behaviour in SU(3) gauge theory with twelve flavours. We use the step-scaling method to investigate the running coupling constant over a large range of scale. Our renormalisation scheme is defined via the ratio of Polyakov loop correlators in the twisted and untwisted directions.
In particular, we compute the ratio, r σ (u) defined in Eq. (17), between the step-scaling function and the input renormalised coupling. In our central analysis procedure, we perform the continuum extrapolation using the 3-point linear fit with the reference coupling computed on the L/a = 6, 8, 10 lattices. Data on these lattices are free of volume interpolation. Using this procedure, we find that this theory contains an IRFP at around g 2 * ∼ 2. In this work, we have investigated systematic errors in the bare coupling interpolation for each lattice volume, and the continuum extrapolation. We have performed reasonable variations on these interpolation and extrapolation, and carefully examined possible correlation amongst data points used in the continuum extrapolation. We find that the dominant systematic effect arises from the continuum extrapolation. To gain information about possible errors in this extrapolation, we compute the step-scaling function on the L/a = 14 lattice, obtain the reference input renormalised coupling for the L/a = 7 lattice using an interpolation procedure, and then study the continuum limit using the 4-point linear and quadratic fits, as well as the 3-point linear fit without the L/a = 6 data. We find that all our analysis procedures result in evidence for the existence of an IRFP, although the latter two continuumextrapolation methods result in significant errors. In view of this, the result of our work suggests that SU(3) gauge theory with twelve fermions in the fundamental representation contains an IRFP.
Our finding shows that the conformal window for SU(3) gauge theories with fundamental fermions may lie below N f = 12. Although this conclusion agrees with most other studies [7-11, 26, 33, 46], the result in Ref. [34,37] leads to the opposite conclusion. Combining this information with the recent result from the N f = 10 calculation [12], this can indicate that the N f = 12 is already very close to the lower bound of the conformal window.  In this appendix, we present details for the plaquette values and the TPL scheme renormalised coupling constants obtained at various lattice volumes, L/a, and bare couplings, β. We also give the numbers of HMC trajectories for the computation of the TPL coupling.