Unveiling confinement in pure gauge SU(3): flux tubes, fields, and magnetic currents

A characteristic signature of quark confinement is the concentration of the chromoelectric field between a static quark–antiquark pair in a flux tube. However, the structure of this flux tube, and hence of the confining force, has not been completely understood. Here we perform new lattice measurements of field distributions on smeared Monte Carlo ensembles in SU(3) gauge theory. On the basis of these simulations we demonstrate that the confining force can be understood using the analogy with the basic principles of electromagnetism as elucidated by Maxwell. We derive a chromomagnetic Lorentz force density coupling the chromoelectric field to chromomagnetic currents and integrate this force density over the flux tube interior to obtain a Maxwell-like force that squeezes the flux tube in the transverse direction. We show that the strength of this transverse confining force is equal to the value of the string tension calculated numerically from the chromoelectric field on the midplane between the quarks, verifying the consistency of these two complementary pictures of confinement.


Introduction
The confinement of quarks and gluons inside hadrons remains an open problem of Quantum Chromodynamics (QCD).A theoretical explanation of this phenomenon is still missing and our current understanding is based on models of the QCD vacuum (for a review, see Refs.[1,2]) and Monte Carlo numerical simulations of QCD on a space-time lattice.
A great deal of numerical evidence shows that a static quark and antiquark interact via a confining linear potential a e-mail: mbaker4@uw.edub e-mail: chelnokov@itp.uni-frankfurt.dec e-mail: leonardo.cosmai@ba.infn.itd e-mail: cuteri@itp.uni-frankfurt.dee e-mail: alessandro.papa@fis.unical.itfor distances equal to or larger than about 0.5 fm.This linear potential is almost entirely due to the electric 1 field, which is mostly longitudinal, i.e., oriented along the line connecting the static quark and antiquark [3][4][5][6].
Many numerical studies in SU (2) and SU(3) Yang-Mills theories  have characterized the shape of the electric field in the transverse plane at the midpoint of the line connecting the static quark and antiquark.More recently, numerical investigations have extended their reach, achieving a detailed description of all components of the color fields around static sources [34,35], as well as the spatial distribution of the stress energy-momentum tensor [35][36][37] and the flux densities for hybrid static potentials [38,39].This growing numerical phenomenology about color fields near static sources could give new hints in quest for the mechanism of confinement.
The present paper extends the results of Refs.[34,35].The main message in those papers was that the electric field can be viewed as the superposition of a 'perturbative' part, which is the sole contributor to components transverse the quark-antiquark axis, and a 'non-perturbative', longitudinal part, all components of the magnetic field being negligibly small everywhere.Here we focus instead on the distribution of electric and magnetic sources and currents, as they can be inferred from the shape of color fields.In particular, we present for the first time evidence of the solenoidal magnetic current responsible for the longitudinal electric field and study its behavior toward the continuum limit, leading to new understanding of confinement.
The organization of the paper is as follows: Sections 2 to 5 are devoted to the theoretical background, Section 6 describes our numerical setup, Section 7 contains our new results, and conclusions are drawn in Section 8.

Connected correlator and the field strength tensor
In previous papers [34,35], using lattice simulations of SU (3) pure gauge theory we have measured the spatial distributions of the color fields induced by a quark-antiquark pair separated by a range of distances d ranging from 0.37 fm to 1.25 fm.
These distributions were obtained from lattice measurements of the connected correlation function ρ conn W,µν [11] of a plaquette U P = U µν (x) in the µν plane, and a Wilson loop W in the 41 plane (see Figure 1), N = 3 being the number of QCD colors.
The correlator ρ conn W,µν provides a lattice definition of a gauge-invariant field strength tensor F µν q q ≡ F µν carrying a unit of octet charge, while possessing the space-time symmetry properties of the Maxwell field tensor of electrodynamics.
W UP L (Schwinger line) Figure 1 The connected correlator between the plaquette U P and the Wilson loop (subtraction in ρ conn W, µν not explicitly drawn).The longitudinal electric field E x (x t ) relative to the position of the static sources (represented by the white and black circles), for a given value of the transverse distance x t .(x l = x, x t = y.) When the plaquette U P lies in the 41 plane, the measured 41 component of the field tensor determines E x , the component of the electric field along the q q axis E x = F 41 ; i.e., the longitudinal component of the electric field at the position corresponding to the center of the plaquette.
When U P is in the 42 plane, F 42 = E y , a component of the electric field transverse to the q q axis, and when U P is in the 23 plane, F 23 = B x , the longitudinal component of the magnetic field.In the numerical evaluations of Refs.[34,35] all components B i of the magnetic field were equal to zero within statistical errors.
3 Separation of the electric field into 'perturbative' and 'non-perturbative' components The transverse components E y of the simulated electric field E were fit (Ref.[35]) to the transverse components of an effective 'perturbative' Coulomb-like field E C satisfying the condition: Subtracting E C from the simulated field E gives the 'nonperturbative' longitudinal field E NP pointing along the q q axis êx : Throughout this paper all the field derivatives (e.g., ∇ × E and ∇ • E) are calculated numerically on the lattice from the measured electric field E(x l , x t ).See Fig ( 1).

The 'Maxwell' picture of the Yang-Mills flux tube
In our previous paper [35] we showed that the string tension σ can be obtained from lattice simulations of the connected correlator ρ conn W,µν , given in Eq. ( 2), defining the Maxwell-like field F µν , as the integral over the midplane x l = d/2 between the quarks of the longitudinal component T xx of the Maxwell stress tensor, constructed from the non-perturbative electric field E NP (x l , x t ): where Independently, the string tension can be determined from the value of E NP at the position of the quark.
In this paper we use the lattice simulations of F µν to evaluate the divergence of the Maxwell stress tensor T αβ and obtain the force density f β in the Yang-Mills flux tube; Identifying the magnetic currents, we can express the force density (8) as the sum of an electric Lorentz force density and a magnetic Lorentz force density corresponding to a flux tube comprised of gauge invariant electric and magnetic currents measured in our simulations.This description emerges from using the divergence of the Maxwell stress tensor to calculate forces, without requiring that the field tensor F µν satisfy Maxwell's equations.
To test this picture we perform new simulations of the electric field E, and numerically evaluate the magnetic current density J mag and the electric charge density ρ el in the flux tube.
If F αβ had been the field tensor of electrodynamics, the second term on the RHS of Eq. ( 8) would vanish because of the homogeneous Maxwell equations, and the first term would be fixed by the inhomogeneous Maxwell equations in terms of the electric current density of the charged matter.The right hand side of Eq. ( 8) then would reduce to the expression for the Lorentz force density.
For the Yang-Mills flux tube, F αβ (see Eq. ( 2)) is the Maxwell-like field tensor measured in our simulations and the derivatives on the right hand side of Eq. ( 8) define the electric current density J el β and the magnetic current density J mag β generated by the field F αβ : • J el 4 and J mag 4 are the electric charge density ρ el and magnetic charge density ρ mag , respectively; • J el i and J mag i , i = 1, 2, 3 are the components of the vector electric current density J el and the vector magnetic current density J mag , respectively.
The current densities expressed in terms of the electric components E k = F 4k and the magnetic components B i = 1 2 ε i jk F jk of the field tensor F µν have the form while Eq.(10) for the Lorentz force density f assumes the form Since the magnetic field B simulated in the static flux tube vanishes within statistical error, Eq. ( 13) simplifies to Furthermore, since ∇ × E C = 0, the perturbative Coulomblike field E C does not contribute to the magnetic current density.Therefore, We will see from our simulations that the electric charge density ρ el = ∇ • E is significantly different from zero only at the positions of the quark sources.
Consequently, the force density f interior to the flux tube directed towards the flux tube axis.We can imagine cutting the flux tube along any plane containing its axis.We then calculate the force F on one half of the flux tube cut by the plane.
Making use of rotational symmetry, we introduce polar coordinates (x t , θ ) in the transverse (y, z) plane.θ is the angle between the y-axis and êx t , the unit vector in the direction of f ; i. e. , y = x t cosθ , z = x t sinθ .
The force F on this half flux tube is then given by: where f (x l , x t ) is the magnitude of the Lorentz force density, Eq. ( 16).The integral cos θ dθ = 2, and the integral sin θ dθ = 0, so that the expression for the force F yields the result directed toward the flux tube axis.
Integrating the force density f over one half of the cut flux tube yields the force F on that half; integrating f over the other half of the tube yields an equal and opposite force on that half, pushing the two halves together.The flux tube is thus confined in the transverse direction by this 'squeezing' force F.
We use our numerical simulations of the connected correlator to evaluate the integrand in Eq. ( 18) determining the magnitude of the magnetic Lorentz force density in the half flux tube.Figure 9 shows the magnitude of this force density generated by magnetic currents for three values of β .
In Section (7.3) we show the confining force F obtained from the integration Eq. ( 18) over the half flux tube is compatible within systematic errors to the value of the string tension σ calculated from the integral Eq. ( 5) of T xx over the midplane between the quarks, checking the consistency of these two representations of confinement.

Direct simulation of the magnetic current density
The possible presence of magnetic currents in SU(3) lattice gauge theory theory was pointed out in Ref. [19], where it was noted that, in contrast to the magnetic monopoles in U(1) lattice gauge theory, the magnetic currents in non-Abelian lattice gauge theory need not be quantized.
In this paper we evaluate the magnetic current density numerically, using our lattice measurements of the connected correlator.
To evaluate numerically, we replace the field derivatives in Eq. ( 19) by the corresponding differences of expectation values of fields measured in our lattice simulations of the connected correlator; i.e. (see Section 2, Figure 1 and Eqs. ( 1) and ( 2)), Replacing Eq. ( 20) and Eq. ( 21) in Eq. ( 19) yields the magnetic current distributions J mag z corresponding to each of the q q separations for which we carried out simulations of the color field distributions.
Before presenting the numerical results of these measurements, we now show that they can be identified as measurements of expectation values of loops constructed from plaquettes lying on opposite faces of an (x, y, x 4 ) cube, as depicted in Figure 2, and consequently are measurements of the flux of magnetic current out of the cube, as proposed in Ref. [19].
The black dot in Figure 2 stands for the Wilson loop W and connecting Schwinger lines L and L + attached to the plaquette U 41 (x l , x t ) in Figure 1.The loop on the left in Figure 2 is then a condensed representation of Figure 1 measuring the longitudinal electric field E x (x l , x t ) q q.In the loop on the right in Figure 2, the plaquette U 41 (x l , x t + 1), translated by one unit in the y direction and measuring E x (x l , x t + 1) q q, is attached to Schwinger lines extending L and L + by single links U 2 (x l , x t ) and U + 2 (x l , x t ) in the y direction.
Eq. ( 20) for ∂ E x ∂ x is then the expectation value of the difference of two loops constructed from plaquettes U 41 lying on opposite xx 4 faces of the (x, y, x 4 ) cube in Figure 2.
Likewise, Eq. ( 21) expresses in Eq. ( 12) gives the contribution to J mag z from plaquettes U 12 on the xy faces of the 'magnetic current' (x, y, x 4 ) cube.

○ ○ x y x4
Figure 2 A pair of paths of closed loops defining a contribution to the magnetic current in the z direction, according to [19].
We can generalize Eq. ( 19), Eq. ( 20) and Eq. ( 21), using the 4-dimensional form Eq. ( 9) for the magnetic current density J mag α .Making the replacements in Eq. ( 9) yields the magnetic current density J mag α q q measuring the flux of magnetic current out of a ( β , μ, λ ) magnetic current cube (the (y, x 4 , x) cube in Figure 2.
The translated plaquette U λ µ (x+a β ) is attached to single links U a ( x) ≡ e iaA β and U + β (x) ≡ e −iaA β , accounting for the extensions of the attached Schwinger lines, so that: Inserting Eq. ( 22) in Eq. ( 9) yields The sum over β , μ, λ is a sum over the six faces of the β , μ, λ magnetic current cube contributing to the α component of the magnetic current density.
If we were to first take the continuum limit a → 0 for fixed q q separation in lattice units (naïve continuum limit), Eqs. ( 23) and ( 24) would reduce to the Bianchi identity: where Instead, we simulated J mag α (x) q q, reducing the lattice spacing, keeping the quark-antiquark separation d fixed in physical units.We have carried out numerical simulations for pure SU(3) gauge theory on a 48 4 lattice for three values of the gauge coupling β , all of which correspond to the quarkantiquark separation d = 0.511 fm (See Table 1).We measured the color fields, as defined in Eq. ( 1), generated by a quark-antiquark pair separated by distance d.For each value of the gauge coupling β in Table 1 measurements were performed every 25 upgrades of the gauge configuration.The value of β was determined for each value of d[lattice units] such that d ≈ 0.511 fm in physical units, where the physical scale for the lattice spacing a(β ) was set according to Ref. [40]: for all β values in the range 5.7 ≤ β ≤ 6.92.In this scheme (see Eq. (3.5) in Ref. [40]) we have, for the value of the square root of the string tension, √ σ ≈ 0.464 GeV.The connected correlator defined in Eq. ( 1) exhibits large fluctuations at the scale of the lattice spacing, which are responsible for small signal-to-noise ratio.To extract the physical information carried by fluctuations at the physical scale (and, therefore, at large distances in lattice units) we smoothed out configurations by a smearing procedure.
Figure 3(a) shows that, at x t = x l = d/2, both the full and the non-perturbative longitudinal electric fields display plateaus as a function of the amount of smearing -quantified as the number of smearing steps times the squared lattice spacing -at all considered values of β (and the same holds true at all values of x t ).The value of the longitudinal component E x for the full field shows no sign of degradation of the signal within the explored range of smearing steps, while for the longitudinal component E NP x there is a very wide maximum, followed by a decrease of the smeared value when the number of smearing steps is further increased.We extract x land x t -dependent optimal amounts of smearing for the field based on the position of the above mentioned maximum.
The smearing procedure can also be validated a posteriori by the observation of continuum scaling, that is by checking that fields obtained in the same physical setup, but at different values of the coupling and of the optimal number of smearing steps, are in good agreement in the range of gauge couplings used.This is seen, for the non-perturbative field, in Figure 4, where E NP x (x t ) at a transverse plane x l = d/4 is plotted for three values of β .When it comes to the behavior versus the number of spatial smearing steps of the magnetic current density (see Figure 3(c)) the situation is somewhat different: We observe a significant degradation of the signal after some small amount of smearing.Smearing still plays a role in improving the signal-to-noise ratio and eliminating short-distance UV fluctuations, but a small amount of smearing seems sufficient.Also in this case, our choice of identifying the optimal amount of smearing by the maximum in ( ∇ × E) z as a function of the amount of smearing is validated by how nicely results obtained at different values of β scale at the maximum, as opposed to at larger amounts of smearing (see Figure 5).
Finally, the behavior versus smearing of the full and non-perturbative electric charge density (see Figure 3(b)) is indicative of the fact that these quantities will vanish when going to the continuum limit (see also Figure 6).In this case we fix the amount of smearing by the plateau displayed by ∇ • E at our smallest β and then we scale that amount of smearing on our finer lattices by taking into account the diffusive nature of the smearing process.

Numerical Evaluation of ρ el
As shown in section ( 6), the behavior under smearing of the divergence of both full and non-perturbative fields does not show a maximum and a further decay, but rather approaches a plateau, which in most cases is close to zero.
Additionally, the scaling of the smeared divergences is worse than for other observables, suggesting that the values we obtain have large contributions from discretization errors.In this case we have to use another approach to the choice of the optimal smearing step, namely fixing it at a value that lies in the plateau region for all lattice points.
Figure 7 shows ∇ • E, the divergence of the full field E for d = 0.51 fm at β = 6.769, after 96 smearing steps.As one can clearly see, the values of the divergence are significantly different from zero only at distances of about two lattice steps away from the sources.
The visible peaks of the divergence of the full field are removed by the subtraction of its perturbative part, so the non-perturbative charge density becomes negligible.

Numerical Evaluation of J mag
As seen in Figure 8 the magnetic current density J mag (x), in contrast to the electric charge density ρ el (x), is manifest throughout the whole length of the flux tube.J mag (x) points in the direction êθ , circulating about the flux tube axis.We use the procedure defined in the previous section, to extract the amount of smearing resulting in the local maximum of J mag (x) at each point (x l , x t ). Figure 8 shows the dependence of the magnetic current density on x l and x t for d = 0.51 fm and two values of β .
The continuum scaling behavior of the magnetic current density is shown in Figure 5.One can see that all three scales give compatible values of the current density except for the smallest β at the points close to x t = 0.
Note that, since the perturbative field is defined explicitly as that for which the curl is equal to zero, the non-perturbative current density, as opposed to the charge density, is exactly equal to the full current density, and thus is not dependent on the conjecture that the non-perturbative field is purely longitudinal.

Calculation of the confining force F and the string tension σ
Here we evaluate the magnitude of the force F and of the string tension σ by numerically calculating the integrals in Eq. ( 5) and Eq. ( 18).First, we determine the integrand of Eq. ( 18) over a grid of (x l , x t ) points, using the averages over ensembles of J mag z and E NP x .The behavior under smearing J mag z shows a fast growth up to a broad maximum, followed by a slow decay; under smearing E (NP) x behaves similarly, except that the maximum is much less pronounced and is typically followed by an extended plateau.In both cases we chose as optimal smearing step as the one at the local maximum.Then, we spline-interpolate the discrete integrand by a smooth function using the tool Interpolation of Wolfram Mathematica.
In Figure 9 we show the behavior of the integrand of Eq. ( 18) on the (x l , x t )-plane for three values of β .The numerical integration was finally performed using NIntegrate in Mathematica.
Our results for √ F are summarized in Table 2.The statistical errors on the results were calculated by replicating the numerical integration 100 times after a flat reshuffling of the values of J mag z

and E NP
x entering the integrand of Eq. ( 18) within their respective 3σ uncertainties, neglecting correlations between observables.
Table 2 Summary of the numerical determination of the (square root of the) force, Eq. ( 18) and the (square root of the) stress tensor Eq. ( 28); the last four columns give the edges of the integration domains.See the text for how the assessment of systematics on √ F is performed.[fm] [fm] Figure 9 The integrand of Eq. ( 18), 2x t f (x l , x t ) = 2x t times the magnetic force density, for three values of β , all corresponding to the same quark-antiquark separation in physical units, d = 0.511 f m.The force density falls off for large x t for all values of x l , manifesting confinement.
We also estimated systematic errors on √ F. Estimates of the positive systematic errors were obtained by integrating over the variable x l up to the midpoint and then doubling the result assuming symmetry.This choice is motivated by our observation of an asymmetry, in the x l direction w.r.t. the midpoint, of the integrand of Eq. ( 18) that is, in fact, systematically larger correspondingly to shorter Schwinger lines (x l < d/2), as a result of an analogous asymmetry both in the magnetic current and in the non-perturbative field, which we ascribe to an interplay of renormalization and noise reduction through smearing.Estimates of the negative systematic error (available only for the two larger values of β ) were, instead, obtained by performing the integration in x l over the smallest domain, i.e. the region of integration in the case of β = 6.240.By considering the smallest available domain, while clearly underestimating the actual result, we eliminate the system-atic discrepancy between results at different lattice spacings coming from the fact that shorter and shorter distances from the sources can be attained on finer and finer lattices.
Table 2 also shows the calculation of √ σ from the nonperturbative field at the midpoint of the flux tube (See Eq. ( 5)): In Eq. ( 28), after a spline interpolation of the numerical data for E NP x (x midpoint l , x t )), the integral has been numerically evaluated.
We observe that the resulting value for √ F remains fairly stable approaching the continuum and is only slightly larger than √ σ , which is estimated to be 0.464 GeV [40].
In the Appendix we compare our results for the string tension with the value of the non-perturbative field E NP at the position of the quark sources.

Conclusions
Using the connected correlator function ρ conn W,µν to provide a lattice definition of a gauge invariant field strength tensor F µν in SU(3) gauge theory, we have isolated a non-perturbative gauge invariant longitudinal electric field, E NP , and magnetic currents, J mag , that circulate about the axis of the flux tube.The magnetic current density J mag (x), in contrast to the electric charge density ρ el (x), is manifest throughout the whole length of the flux tube.The behavior of the expectation value of magnetic currents over smeared Monte Carlo ensembles suggests that they can be safely extrapolated to the continuum and have therefore the status of physical observables.
The interaction of the magnetic currents with the nonperturbative longitudinal electric field gives rise to a magnetic Lorentz force density f = J mag × E NP filling the flux tube, directed toward the flux tube axis.Figure 9 shows the calculated magnitude of this force density rises to a maximum at x t ≈ 0.1 fm all along the flux tube length, and then falls off for large x t , manifesting confinement.
Integrating the force density f over one half of a flux tube cut by a plane through its axis produces a force F on that half; integrating f over the other half of the tube produces an equal and opposite force on the other half, pushing the two halves together.Thus the flux tube is confined by these Maxwell-like forces on magnetic currents; the flux inside cannot spread out.
We find that the magnitude of the confining force F is compatible, within systematic errors, to the magnitude of the string tension σ , which is in turn determined by the distribution of the stress tensor on the midplane perpendicular to the flux tube axis Eq. ( 28), checking the consistency of our numerical calculations.
On the basis of our numerical simulations of field distributions in SU(3) gauge theory and using the analogy with the basic principles of electromagnetism, we have thus obtained a consistent physical picture of the interior of the flux tube, thereby unveiling confinement.This investigation was in part based on the MILC collaboration's public lattice gauge theory code (https://github.com/milc-qcd/).Numerical calculations have been made possible through a CINECA-INFN agreement, providing access to HPC resources at CINECA.LC and AP acknowledge support from INFN/NPQCD project.FC and VC acknowledge support by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) through the CRC-TR 211 "Strong-interaction matter under extreme conditions"project number 315477589 -TRR 211.FC acknowledges the support by the State of Hesse within the Research Cluster ELEMENTS (Project ID 500/10.006).
Appendix A: The force and the field at the position of the source We use Figure 1 and Eqs. ( 1) and (2) of our paper to evaluate the field E x = F 41 as follows.Let U P = U 41 (x 4 = 0, x 1 = x l = a/2, x t = 0), i.e., the center of the plaquette is the point x 4 = 0, x 1 = a/2.Then connect this plaquette to the center of the link U * 4 (x 4 = a/2, x 1 = 0) of the Wilson loop W in Figure 1 with the following choice of the Wilson Lines L and L * in Figure 1: • The Schwinger line L in Figure 1  where C is the loop bounding the area of a loop in which the area of the plaquette U P has been removed from the area of the original loop depicted in Figure 1.Eq. (A.1) follows from the fact that in the product W LU P the matrix U * connecting the point x 4 = a/2 with the point x 4 = 0 in W eliminates a corresponding matrix U connecting these two points in U P .(UU * = 1).Similarly the product U P L * W eliminates matrices connecting the points x 4 = −a/2 and x 4 = 0 in U P and W .As a consequence in the tr (W LU P L) the lines connecting x 4 = a/2 and x 4 = a/2 in both W and in U P have been eliminated.There remain direct connections between W and U P at x 4 = a/2 and at x 4 = −a/2, yielding the result tr W (C ), where C is a loop traversing the remaining three lines in U P and then connecting to W at these points.This is result asserted in Eq. (A.1).
Next note that to order a 2 the second term in Eq. (1) of our paper equals one since tr U P = N to that order.Thus if we multiply Eq. ( 1) by tr (W ) and use Eq.These values are larger than our previous evaluations (Table 2) and do not exhibit scaling, reflecting the uncertainty in our calculations of E NP close to the sources.

∂
E y ∂ x as the expectation value of the difference of a corresponding pair of loops constructed from plaquettes U 42 lying on opposite yx 4 faces of the cube, so that J mag z can be expressed as the sum of contributions from the yx 4 and xx 4 faces of the cube.The term − ∂ B z ∂ x 4

2 [Figure 3
Figure 3 Full and non-perturbative longitudinal electric field at x l = x t = d/2 (3(a)), full and non-perturbative electric charge density at x l = x t = d/4 (3(b)), and magnetic current density at x l = d/4, x t = d/2 (3(c)) as functions of smearing.

Figure 4 Figure 5
Figure 4 Continuum scaling of the longitudinal component of the nonperturbative field E NP x at a transverse plane x l = d/4 ≈ 0.128 fm.

Figure 6
Figure 6 Scaling of the divergence of full field at a transverse plane x l = d/4 ≈ 0.128 fm

Figure 7
Figure 7 Divergence of the full field produced by a quark-antiquark pair separated by d = 0.51 fm, corresponding to β = 6.769 after 96 smearing steps.
connecting the point x 4 = 0, x 1 = 0 on W to the point x 4 = a/2, x 1 = 0 on U 41 ; • The Schwinger line L * in Figure 1 connecting the point x 4 = −a/2 , x 1 = 0 on U 41 to the point x 4 = 0, on W .With this choice of U P , L and L * , it then follows that tr (W LU P L * ) = tr (W (C )) , (A.1)

a 2 =
(A.1), we obtain: tr (W ) ρ conn W,41 = tr (W (C )) − tr (W ) .(A.2) The left hand side of Eq. (A.2) determines the field at the position of the quark.The right hand hand side of Eq. (A.2) determines the change in the heavy quark potential when the piece of the Wilson loop between x 4 = −a/2 and x 4 = a/2 is displaced by a distance a in the x direction.Then, making use of Eq. (A.2), we can calculate the derivative of the heavy quark potential, and, therefore, the string tension from the field at the position of the quark [43, gE x (0) , (A.3) where ∂ A denotes the derivative with respect to the area A. This provides us with an alternative method to evaluate the string tension, giving the following results: β = 6.240 , √ σ = 0.56353(81) GeV , β = 6.544 , √ σ = 0.5962(38) GeV , β = 6.769 , √ σ = 0.617(16) GeV .

Table 1
Summary of the numerical simulations.