The confining color field in SU(3) gauge theory

We extend a previous numerical study of SU(3) Yang–Mills theory in which we measured the spatial distribution of all components of the color fields surrounding a static quark–antiquark pair and provided evidence that the simulated gauge invariant chromoelectric field can be separated into a Coulomb-like ‘perturbative’ field and a ‘non-perturbative’ confining field. In this paper we hypothesize that the fluctuating color fields not measured in our simulations do not contribute to the string tension. Under this assumption the string tension is determined by the color fields we measure, which form a field strength tensor pointing in a single direction in color space. We call this the ‘Maxwell picture of confinement’. We provide an additional procedure to isolate the confining field. We then extract the string tension from a stress energy-momentum tensor having the Maxwell form, constructed from the simulated non-perturbative part of the field strength tensor. To test our hypothesis we calculate the string tension for values of the quark–antiquark separation ranging from 0.37 fm to 1.2 fm. We also calculate the spatial distributions of the energy-momentum tensor surrounding static quarks for this range of separations, and we compare with the distributions obtained from direct simulations of the energy-momentum tensor.


Introduction
Quantum chromodynamics (QCD) is universally accepted as the theory of strong interactions. Nobody doubts that the well established phenomenon of confinement of quarks and gluons inside hadrons is encoded into the QCD Lagrangian. Yet, our current understanding does not go beyond that provided by a number of models of the QCD vacuum (for a review, see Refs. [1,2]). In particular, a theoretical a priori explanation of the so called area law in large size Wilson loops, which is closely related to a linear confining potential between a static quark and antiquark at large mutual distances, is still missing.
In such a challenging situation, first-principle Monte Carlo simulations of QCD on a space-time lattice represent an indispensable tool not only for checking (or ruling out) models of confinement, but also for providing new numerical "phenomenology" and possibly stimulating original insights into the mechanism of confinement.
Numerical simulations have established that there is a linear confining potential between a static quark and antiquark for distances equal to or larger than about 0.5 fm. This linear regime extends to infinite distances in SU(3) pure gauge theory, and, in the presence of dynamical quarks to distances of about 1.4 fm, where string breaking should take place [3][4][5][6]. The long-distance linear quark-antiquark potential is naturally associated with a tube-like structure ("flux tube") of the chromoelectric field in the longitudinal direction, i.e. along the line connecting the static quark and antiquark [7][8][9][10].
A wealth of numerical evidence of flux tubes has accumulated in SU(2) and SU(3) Yang-Mills theories . Most of these studies concentrated on the shape of the chromoelectric field in the transverse plane at the midpoint of the line connecting the static quark and antiquark, given that the other two components of the chromoelectric field and all the three components of the chromomagnetic field are suppressed in that plane.
Recent times have witnessed an increasing numerical effort toward a more comprehensive numerical description of the color field around static sources, via the measurement of all components of both chromoelectric and chromomagnetic fields on all transverse planes passing through the line between the quarks [38]; of the spatial distribution of the stress energy-momentum tensor [39,40]; and the flux densities for hybrid static potentials [41,42]. A more complete numerical description of the color field around the sources brings improved visualization, enabling us to grasp features otherwise less visible.
In the numerical study [38] we simulated the spatial distribution in three dimensions of all components of the chromoelectric and chromomagnetic fields generated by a static quark-antiquark pair in pure SU(3) lattice gauge theory. We found that, although the components of the simulated chromoelectric field transverse to the line connecting the pair are smaller than the simulated longitudinal chromoelectric field, these transverse components are large enough to be fit to a Coulomb-like 'perturbative' field produced by two static sources parameterized by effective charges ±Q of the sources (see Eq. (5) below).
The longitudinal component of this Coulomb-like 'perturbative' field accounts for a fraction of the simulated longitudinal chromoelectric field. We then identified the remaining longitudinal chromoelectric field as the confining 'nonperturbative' part of the simulated SU(3) flux tube field.
It is this non-perturbative part of the simulated field which contributes to the coefficient of the linear term in the heavy quark potential, the string tension.
In this paper we extend our simulations to a wider range of quark-antiquark separations. We extract the string tension from these simulations and compare our analysis with the results of recent simulations [39] of the energy-momentum tensor in SU(3) Yang-Mills theory.
We present a new procedure (the curl procedure) to extract a perturbative Coulomb field E C from the transverse components of the numerically simulated chromoelectric field. We avoid the use of a fitting function, directly imposing the condition that E C is irrotational (see Eq. (6) below). This provides a second method for implementing the underlying idea of our previous paper; that is, the chromoelectric field generated by a quark-antiquark pair can be separated into perturbative and non-perturbative components by a direct analysis of lattice data on the color field distributions generated by the pair.
As noted in [38], we can extract the value of the string tension from the non-perturbative field by utilizing the fact that the value of the chromoelectric field at the position of a quark is the force on the quark [43]. However, the Coulomb-like field (Eq. (5)) does not give a good description of the trans-verse components of the chromoelectric field at distances closer than approximately two lattice steps from the sources [38], so that we must use the curl procedure to isolate the confining field in order to extract the string tension directly as the force.
The color fields F μν we measure, defined by the gauge invariant correlation function ρ conn W,μν (see Eqs. (1) and (2), below), point in a single direction in color space, parallel to the direction of the 'source' Wilson loop. In this paper we construct a stress energy-momentum tensor T μν having the Maxwell form from the 'measured' flux tube field tensor F μν , and extract the string tension (see Appendix A). This leads to a picture of a confining flux tube permeated with lines of force of a gauge invariant field tensor F μν carrying color charge along a single direction.
The 'Maxwell' energy-momentum tensor T μν does not account for the contribution to the quark-antiquark force from the fluctuating color fields not measured in our simulations. On the other hand, the complete Yang-Mills energy momentum tensor T YM μν simulated in Ref. [39] includes these fluctuating contributions, so that comparison of T YM μν with the 'Maxwell' energy-momentum tensor T μν constructed from the chromoelectric and chromomagnetic fields measured in our simulations provides a measure of the fluctuating contributions to the stress tensor.
We noted in our previous paper that the Coulomb-like 'perturbative' field (Eq. (5)) generated a stronger long distance Coulomb force between the heavy quarks than the Coulomb force measured in lattice simulations of the heavy quark potential [44][45][46] indicating the importance of fluctuations for the Coulomb contribution. In this paper we reexamine this issue.
The paper is organized as follows: in Sect. 2 we present the theoretical background and the lattice setup. In Sect. 3 we show some results on the spatial distribution of the color field around the two static sources and review the procedure to extract its non-perturbative part by subtraction of the Coulomb-like perturbative part identified by a fit of the transverse components of the chromoelectric field; in Sect. 4 we describe the new curl procedure to isolate the non-perturbative part; in Sect. 5 we show how to determine the string tension and other fundamental parameters describing the (non-perturbative) flux tube; finally, in Sect. 6 we discuss our results and give some ideas for future work and in Sect. 7 we describe the 'Maxwell picture of confinement' and how it emerges from our work.

Theoretical background and lattice setup
The lattice operator whose vacuum expectation value gives us access to the components of the color field generated by . 1 a The connected correlator given in Eq. (1) between the plaquette U P and the Wilson loop (subtraction in ρ conn W, μν not explicitly drawn). b The longitudinal chromoelectric 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 a static qq pair is the following connected correlator [15,16,47,48]: Here U P = U μν (x) is the plaquette in the (μ, ν) plane, connected to the Wilson loop W by a Schwinger line L, and N is the number of colors (see Fig. 1).
The correlation function defined in Eq. (1) measures the field strength F μν , since in the naive continuum limit [16] ρ conn where qq denotes the average in the presence of a static qq pair, and 0 is the vacuum average. This relation is a necessary consequence of the gauge-invariance of the operator defined in Eq. (1) and of its linear dependence on the color field in the continuum limit (see Ref. [33]).
The lattice definition of the quark-antiquark field-strength tensor F μν is then obtained by equating the two sides of Eq. (2) for finite lattice spacing. In the particular case when the Wilson loop W lies in the plane withμ =4 andν =1 (see Fig. 1a) and the plaquette U P is placed in the planeŝ 41,42,43,23,31,12, we get, respectively, the color field components E x , E y , E z , B x , B y , B z , at the spatial point corresponding to the position of the center of the plaquette, up to a sign depending on the orientation of the plaquette. Because of the symmetry (Fig. 1), the color fields take on the same values at spatial points connected by rotations around the axis on which the sources are located (the1-or x-axis in the given example).
As far as the color structure of the field F μν is concerned, we note that the source of F μν is the Wilson loop connected to the plaquette in Fig. 1. The role of the Schwinger lines entering Eq. (1) is to realize the color parallel transport between the source loop and the "probe" plaquette. The Wilson loop defines a direction in color space. The color field E that we measure in Eq. (2) points in that direction in the color space, i.e. in the color direction of the source.
There are fluctuations of the color fields in the other color directions. We assume that these fluctuating color fields do not contribute to the string tension, so that the flux tube can be described as lines of force of the simulated field E.
The simulated flux tube field E carries color electric charge and color magnetic current along a single direction in color space. The divergence of E is equal to the color electric charge density ρ el (x) and the curl of E is equal to the color magnetic current density J mag (x). Furthermore, the divergence of E satisfies Gauss' law [23], so that ∇ · E(x) = 0 provided that x does not coincide with a source point. Finally, the confining force is calculated from the divergence of a stress tensor T μν having the Maxwell form Eq. (A.1).
The operator in Eq. (1) undergoes a non-trivial renormalization, which depends on x t , as discussed in a recent work [49]. The procedure outlined in that paper to properly take into account these renormalization effects is prohibitively demanding from the computational point of view for Wilson loops and Schwinger lines with linear dimension of the order of 1 fm, where the interesting physics is expected to take place. For this reason, we adopt here the traditional approach to perform smearing on the Monte Carlo ensemble configurations before taking measurements (see below for details). As shown in the Appendix A of our previous paper [38], smearing behaves as an effective renormalization, effectively pushing the system towards the continuum, where renormalization effects become negligible. The a posteriori validation of the smearing procedure is provided by the observation of continuum scaling: as carefully checked in Ref. [34], fields obtained in the same physical setup, but at different values of the coupling, are in perfect agreement in the range of parameters used in the present work.
We performed all simulations in SU(3) pure gauge theory, with the standard Wilson action as the lattice discretization. A summary of the runs performed is given in Table 1. The error analysis was performed by the jackknife method over bins at different blocking levels.
We set the physical scale for the lattice spacing according to Ref. [44]: for all β values in the range 5.7 ≤ β ≤ 6.92. In this scheme the value of the square root of the string tension √ σ ≈ 0.465 GeV (see Eq. (3.5) in Ref. [44]).
The correspondence between β and the distance d shown in Table 1 was obtained from this parameterization. Note that the distance in lattice units between quark and antiquark, corresponding to the spatial size of the Wilson loop in the connected correlator of Eq. (1), was varied in the range d = 8 a to d = 16 a.
The connected correlator defined in Eq. (1) exhibits large fluctuations at the scale of the lattice spacing, which are responsible for a bad 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. Our setup consisted of (just) one step of HYP smearing [50] on the temporal links, with smearing parameters (α 1 , α 2 , α 3 ) = (1.0, 0.5, 0.5), and N APE steps of APE smearing [51] on the spatial links, with smearing parameter α APE = 0.25. The number of smearing steps quoted in the last column of Table 1 was chosen using the same criterion as in Ref. [38], namely that it was taken large enough that all field components, on all transverse planes and at all values of x t (except possibly the few largest ones) have reached their plateaux. Further details on this criterion and on the comparison between smearing and renormalization are given in the Appendix A of Ref. [38].

Spatial distribution of the color fields
Using Monte Carlo evaluations of the expectation value of the operator ρ W, μν over smeared ensembles, we have determined the six components of the color fields on all twodimensional planes transverse to the line joining the color sources allowed by the lattice discretization. These measurements were carried out for several values of the distance d between the static sources, in the range 0.37 fm to 1.19 fm, at values of β lying inside the continuum scaling region, as determined in Ref. [34].
We found that the chromomagnetic field is everywhere much smaller than the longitudinal chromoelectric field and is compatible with zero within statistical errors (see, e.g., Fig. 3 of Ref. [38]). As expected, the dominant component of the chromoelectric field is longitudinal, as is seen in Fig. 2, where we plot the components of the simulated chromoelectric field E at β = 6.370 as functions of their longitudinal displacement from one of the quarks, x l , and their transverse distance from the axis, x t .
While the transverse components of the chromoelectric field are also smaller than the longitudinal component, they are larger than the statistical errors in a region wide enough that we can match them to the transverse components of an effective Coulomb-like field E C produced by two static sources. For points which are not very close to the quarks, this matching can be carried out with a single fitting parameter Q, the effective charge of static quark and antiquark sources determining E C .
To the extent that we can fit the transverse components of the simulated field E to those of E C with an appropriate choice of Q, the non-perturbative difference E NP between the simulated chromoelectric field E and the effective Coulomb field E C , will be purely longitudinal. We then identify E NP as the confining field of the QCD flux tube.
To illustrate this idea, let us fix, for the sake of definiteness, β = 6.240 and put the two sources at a distance d = 16a = 1.02 fm. We then consider the plane, transverse to the longitudinal x-axis connecting the two sources, at a distance x l = 4a from one of them, and evaluate the components E x , E y and E z of the chromoelectric field in 0 0.5

Coulomb fit to
The three components of the chromoelectric field measured at β = 6.240, d = 16a = 1.02 fm, for x l = 4a and the three components of the perturbative Coulomb field, obtained from fitting the transverse E y field component to the form (5) this transverse plane. The lattice determinations of E y on this plane can be fitted by the y-component of an effective Coulomb field where r Q and r −Q are the positions of the two static color sources and R 0 is the effective radius of the color source, introduced to explain, at least partially, the decrease of the field close to the sources. This fit is shown in Fig. 3 -see black dots and black solid line. Using the values of the fit parameters Q and R 0 obtained by the fit of E y , one can construct E C z and E C x and compare them to lattice data. Furthermore, the Coulomb-like content of E z fully accounts for the z-component of the chromoelectric field (see red dots and red solid line in Fig. 3), but E C x accounts for only a fraction of the longitudinal component of the chromoelectric field (see blue dots and blue solid line in Fig. 3). This strongly suggests that the non-perturbative component of the chromoelectric field is almost completely oriented along the longitudinal direction. It can be isolated once the parameters of the Coulomb-like component are determined by a fit to the y-and/or z-components of the lattice determination of the chromoelectric field.
The procedure we have just illustrated in a specific case, can be carried out in a systematic manner. We observe that in making the fit we must take into account that the color fields are probed by a plaquette, so that the measured field value should be assigned to the center of the plaquette. This Table 2 Values of the fit parameters Q and R 0 extracted from Coulomb fits of the transverse components of the chromoelectric field and values of the longitudinal chromoelectric fields at (d/2, 0), the midpoint between the sources and transverse distance zero, for several values of distance d. E x (d/2, 0) is the unsubtracted simulated field and E NP x (d/2, 0) is the non-perturbative chromoelectric field. For comparison, in the last column the non-perturbative chromoelectric field E NP x curl (d/2, 0) obtained using the irrotational property of the pertur-bative field (see Sect. 4) is given. For the parameters of the Coulomb fit we quote, along with the statistical error, a systematic uncertainty that accounts for the variability in the values of the fit parameters extracted from all acceptable fits to E y and E z at different x l values (for more details, see Appendix B of Ref. [38]). As the distance between the sources is made smaller and smaller the quality of the Coulomb fits deteriorates and Q and R 0 cannot be reliably extracted for d ≤ 0.51 fm  (4) also means that the z-component of the field is probed at a distance of 1/2 lattice spacing from the x y plane, where the z-component of the Coulomb field E C z is non-zero and can be matched with the measured value E z for the same value of Q. For further details about the fitting procedure and the extraction of the fit parameters we refer to Appendix B of Ref. [38].
In Table 2, we list the values of the effective charge Q obtained from the lattice measurements of E z and E y at the values of d, the quark-antiquark separation, considered in this work.
The statistical uncertainties in the quoted Q values result from the comparisons among Coulomb fits of E y and E z at the values of x l , for which we were able to get meaningful results for the fit. The values of R 0 in physical units grow with the lattice step a, while in lattice units they show more stability. This suggests that the effective size of a color charge in our case is mainly explained by lattice discretization artifacts and the smearing procedure, and is not a physical quantity (see Appendix B of Ref. [38]).
Evaluating the contribution of the field of the quark to E C in Eq. (5) at the position r -Q of the antiquark and multiplying by the charge −4π Q of the antiquark yields a Coulomb force between the quark and antiquark with coefficient −4π Q 2 . By comparison, the standard string picture of the color flux tube gives a Coulomb correction of strength −π/12 to the long distance linear potential (the universal Lüscher term arising from the long wave length transverse fluctuations of the flux tube [52]). In addition, the strength π 12 of the Luscher term is approximately equal to the strength of the Coulomb force determined from the analysis of lattice simulations of the heavy quark potential at distances down to about 0.4 fm [44,46].
By contrast, the strength −4π Q 2 of the Coulomb force generated by E C is roughly 4 times larger than π 12 for the values of the effective charge Q listed in Table 2 and determined from our simulations of ρ conn W,μ,ν . Therefore the fluctuating color fields not measured in our simulations must be taken into account in calculating the Coulomb correction to the long distance heavy quark potential.
In Fig. 4 we plot the longitudinal component E NP x of the non-perturbative field in Eq. (4) as a function of the longitudinal and transverse displacements x l , x t at β = 6.370. As expected, E NP x is almost uniform along the flux tube at distances not too close to the static color sources. This feature is better seen in Fig. 5, where transverse cross sections of the field E NP x (x l , x t ), plotted in Fig. 4, are shown for the values of x l specified in Fig. 5. For these values of x l the shape of the non-perturbative longitudinal field is basically constant all along the axis. A similar scenario holds in the other lattice setups listed in Table 1.
In Table 2 we also compare the values of the measured longitudinal chromoelectric field E x with those of the nonperturbative field E NP x on the axis at the midpoint between the quark and antiquark, for all ten values of their separation d. Given that E NP x is almost uniform along the axis,    While the Coulomb field (5) gives a good description of the transverse components of the chromoelectric field when the distance from the sources is not too small, it does not give a good description at smaller distances, approximately two lattice spacings from the sources. This can be either the result of the non-spherical form of the effective charges, or an effect introduced by the discrete lattice.
To extract the confining part of the chromoelectric field in the data it is then preferable to have a procedure which avoids the use of an explicit fitting function, and which can work close to the quark sources. With this aim in mind we use the following two steps to separate the field into 'perturbative' and 'non-perturbative' components.
1. We identify the transverse component E y of the simulated field with the transverse component E C y of the perturbative field, E C y ≡ E y . 2. We impose the condition that the perturbative field is irrotational, curl E C = 0.
Condition (1) implies that the nonperturbative field is purely longitudinal, E NP y = 0. Condition (2) will then fix the longitudinal component E C x of the perturbative field as well as the longitudinal component E NP x of the non-perturbative field.
One possible way for testing how sound the curl procedure is would consist in considering the spatial distribution of the y and z field components in a plane orthogonal to the axis connecting the static sources. This kind of test, on the other hand, would be automatically satisfied for a field that is symmetric under rotations around the x axis and under reflections with respect to the x y plane. The rotational invariance tests have already been performed [27,29] for the longitudinal component of the field in the same setup as in the present paper. Also, the invariance of all components of the field under reflections and 90-degree rotations was checked when obtaining data for the current study.
To implement the irrotational condition (2), taking into account that the fields are measured at discrete lattice points, the sum of the measured fields along any closed lattice path is zero. For example, on a plaquette this amounts to One can easily solve this equation for E C x obtaining This of course, leaves one unknown on each transverse slice of the field -the value of E C x (x, y max + 1), but if the value of y max is large enough, the perturbative field at that distance should already be small, so in our analysis we just put E C x (x, y max + 1) = 0. To check that this indeed makes a little change to our results, we have used a separate procedure in which we fixed E C x (x, y max + 1) = E x , in practice making E NP x = 0 at the largest transverse distance. This procedure gave similar results.
After the estimation of the perturbative longitudinal field E C x one can subtract it from the total field, obtaining the nonperturbative component (see Fig. 6). One can see that the nonperturbative part of the flux tube exhibits very little change along the line connecting the quark-antiquark pair; even at the smallest distances from the sources the non-perturbative field remains smooth (This is seen more clearly in Fig. 7).

The string tension and the width of the flux tube
The forces between charged particles in electrodynamics are determined by a stress tensor T μν constructed from fields F μν satisfying Maxwell's equations (see Eq. (12.113) in Ref. [53]). Similarly, the force between quarks and antiquarks in Yang-Mills theory is determined by the stress tensor T μν , Eq. (A.1), constructed from the field tensor F μν obtained from our simulations.
The quark-antiquark force F is then the integral of the longitudinal component T x x = (E x (x l = d/2, x t )) 2 /2 of the stress tensor over the median plane x = d/2 bisecting the line connecting the quarks, Eq. (A.8). The non-perturbative quark-antiquark force F NP = −ê x σ determining the string tension σ has the corresponding expression in terms of the non-perturbative longitudinal component of the stress tensor  The square root of the string tension is then equal to We have evaluated the integral (Eq. (8)) in two ways:

by direct numerical integration, using the values of E NP
x determined by our simulations, and 2. analytically, by fitting the numerical data for the transverse distribution of E NP x (x t ) as in [26][27][28][29][30] to the Clem parameterization of the field surrounding a magnetic vortex in a superconductor [54].
where φ, μ and α are fitting parameters. In the dual superconducting model [55][56][57][58] λ = 1 μ is the penetration depth and is the Landau-Ginzburg parameter characterizing the type of superconductor. Figure 8 shows an example of the fit of the data to the Clem functional form Eq. (9) for the transverse distribution of E NP x (x t ), obtained using the curl procedure. We can obtain a second expression for the string tension by utilizing the result [43,59] that the force on a quark is equal to the value of the chromoelectric field at the position of the quark. The string tension is then equal to the corresponding value of the confining part of the chromoelectric field Equations (8) and (11) provide two independent ways to extract the string tension from simulations. As mentioned earlier, we must use the curl procedure to isolate the confining field to extract the string tension by Eq. (11).
To obtain additional information about the structure of the chromoelectric flux tube we have calculated the mean square root width: Just as we have evaluated the integral (8) for the string tension, we have evaluated the integral (12) for the mean square root width both numerically, using the data for E NP x (x t ), and analytically, in terms of Clem parameters, fitting the longitudinal component of E NP x (x t ) in the median plane to the Clem parametrization (9). The results of that fit are given in Table 3. In most cases the parametrization in Eq. (9) gives a good description of the field, shown by the values of χ 2 r , though the parameters themselves are somewhat unstable, which reflects the strong correlation between the parameter estimates.
We compared two different methods for calculating the integrals in Eqs. (8) and (12). First, we carried out the numeric integration in Eqs. (8) and (12), respectively, postulating the rotational symmetry of the field. This approach was repeated for both the non-perturbative field obtained using the "curl procedure" (resulting in √ σ int and w 2 int ) and the field obtained using Coulomb subtraction (resulting in √ σ Coulomb and w 2 Coulomb ). Next we calculated the values Table 3 The Clem parameters describing the non-perturbative field transverse section going through the midpoint between the quark and antiquark positions. The data for the fit is obtained using the curl subtraction procedure, taking the perturbative field at y max+1 equal to zero  Table 4 The string tension estimated using the non-perturbative field from the curl procedure by employing different methods (from left to right: numerical integration of the field, analytical integration of the Clem function with parameters given in Table 3, estimation of fields at sources). In the last column we report also the value of the string tension obtained by numerically integrating Eq. (8) and using the nonperturbative field from the Coulomb subtraction E Coulomb  (5) 0.593 (28) of the string tension and the mean square root width using the Clem parameters given in Table 3 to get the values denoted as √ σ Clem and w 2 Clem . (One remark should be made for the width -while we know that the value of E C x (x, y max + 1) in Eq. (7) is small (O(y −2 max )), in the numerator of Eq. 12 this small constant will be multiplied by y 2 (y 3 after the integration over polar angle), which will cause the error introduced to increase with y max . Indeed, the comparison with the analysis done taking E NP x (x, y max ) = 0 shows large discrepancies in this case.) Finally, we evaluated Eq. (11) for the string tension, σ = E NP x (x l = 0, x t = 0), using the curl procedure to determine the magnitude of the non-perturbative field at the sources.
Our results are gathered in Tables 4 and 5 , where we use the notation σ 0 ≡ E N P x (x l = 0, x t = 0) and in Figs. 9, 10. The data shown in Fig. 9 give a consistent value of √ σ for all values of the separation d, with scatter that increases with d as the resolution diminishes. The values of √ σ lie close to 0.465 GeV, the value used in the parameterization Ref. [44] Let us review the basis of our calculations. Our hypothesis is that the string tension is determined by the field E we measured (the 'Maxwell picture of confinement'). We have determined σ from both the transverse structure of the flux tube (Eq. (8)) and its longitudinal structure (Eq. (11)) as shown in Fig. 6c, in which the non-perturbative field has been isolated.
We emphasize that, as discussed in Sect. 3, the 'Maxwell picture of confinement' cannot be used to obtain the Coulomb correction to the string tension. This implies that the fluctuating fields not measured in our simulations must contribute to the Coulomb force.
On the other hand, the Coulomb correction has been obtained by recent direct simulations of the stress energymomentum tensor in Yang-Mills theory [39]. The Yang-Mills stress tensor accounts for the contributions of fluctuating fields but cannot be directly related to measured fields, in contrast to the Maxwell stress tensor T x x = (1/2)E 2 x , determining the string tension (see Eq. (A.9)). Table 5 The flux tube width estimated using the non-perturbative field from the curl procedure by employing different methods (from left to right: numerical integration of the field, analytical integration of the Clem function with parameters given in Table 3). In the last column we report also the value of the width obtained by numerically integrating Eq. (12) and using the non-perturbative field from the Coulomb subtraction E Coulomb

Conclusions and outlook
In this paper we have determined the spatial distribution in three dimensions of all components of the color fields generated by a static quark-antiquark pair. We have found that the dominant component of the color field is the chromoelectric one in the longitudinal direction, i.e. in the direction along the axis connecting the two quark sources. This fea- ture of the field distribution has been known for a long time. However, the accuracy of our numerical results allowed us to go far beyond this observation. First, we could confirm that, as observed in [38], all the chromomagnetic components of the color field are compatible with zero within the statistical uncertainties. Second, the chromoelectric components of the color fields in the directions transverse to the axis connecting the two sources, though strongly suppressed with respect to the longitudinal component, are sufficiently greater than the statistical uncertainties that they can be nicely reproduced by a Coulomb-like field generated by two sources with opposite charge (everywhere except in a small region around the sources). In Ref. [38] we subtracted this Coulomb-like field from the simulated chromoelectric field to obtain a non-perturbative field E NP according to Eq. (4) and found that the dependence of the resulting longitudinal component of E NP on the distance x t from the axis is independent of the position x l along the axis, except near the sources, thus suggesting that the non-perturbative field found in this way from lattice simulations can be identified as the confining field of the QCD flux tube.
In this work we have improved the approach of Ref. [38] by presenting a new procedure to subtract the Coulomb-like field, which does not rely on any preconception about its analytic form, but is based only on the requirement that its curl is equal to zero.
Moreover, we have carefully analyzed the spatial distribution of the subtracted, non-perturbative part of the longitudinal chromoelectric field to extract from it some relevant parameters of the flux tube, such as the mean width and the string tension, both by means of a fully numerical, modelindependent procedure and by a prior interpolation with the dual version of the Clem function for the magnetic field in a superconductor.
We have also used our determinations of the color field components to construct the Maxwell stress tensor. Details about its determination and a comparison with the recent literature about this topic [39] are presented in Appendix A.
In conclusion, we have shown that the separation of the chromoelectric field into perturbative and non-perturbative components can be obtained by directly analyzing lattice data on color field distributions between static quark sources, with no need of model assumptions. To the best of our knowledge, this separation between perturbative and non-perturbative components has not been carried out previously. It provides new understanding of the chromoelectric field surrounding the quarks. We have used the non-perturbative field to calculate the string tension and the spatial distribution of the energy-momentum tensor surrounding the static quarks, under the assumption that the fluctuating color fields not measured in our simulations do not contribute to the string tension. The extension of our approach to the case of QCD with dynamical fermions with physical masses and at non-zero temperature and baryon density is straightforward [33].

A closing thought
As a closing thought we review the two hypotheses from which we concluded that the string tension can be extracted from the measured gauge invariant field tensor F μν having the space-time symmetry properties of the Maxwell field tensor of electrodynamics.
1. We assume that the fluctuations of the color fields in color directions other than the one we measure do not contribute to the string tension. Then the quark-antiquark force is determined by the Maxwell stress tensor T μν , constructed from the components of the gauge invariant field tensor F μν measured in our simulations. The measured chromoelectric field E can be visualized by lines of force which, in Yang-Mills theory, correspond to the Coulomb lines of force between electrically charged particles in electrodynamics. 2. In order to determine the string tension from the stress tensor T μν , the Coulomb-like perturbative color field E C of the sources must be separated from the measured chromoelectric field E, thereby isolating its confining part E NP . We determine E C by assuming that the components of the perturbative field E C transverse to the flux tube axis are equal to those of the measured field E. We then implement this assumption with a parameterization-free procedure (the curl procedure) to determine the longitudinal components of E C .
Under assumptions (1) and (2) the confining quark-antiquark force is calculated from the Maxwell stress tensor T μν constructed from the components of the non-perturbative field tensor F NP μν , obtained by subtracting the contribution of the quark sources from the measured field tensor F μν . We call this the 'Maxwell picture of confinement'.
The values of the string tension calculated in this manner from lattice data for a range of quark-antiquark separations (shown in Fig. 9) thus provide a test, of both our assumption that the fluctuating color fields not measured in our calculations do not contribute to the string tension, and of our assumption used to isolate the confining part of the simulated field.
To the extent that our assumptions are realized, the 'Maxwell picture of confinement' emerges from our simulations. Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecomm ons.org/licenses/by/4.0/. Funded by SCOAP 3 .

Appendix A: The 'Maxwell' stress tensor
In this Appendix we consider the "Maxwell" energy-momentum tensor T μν as a function of the field tensor F μν characterizing the SU(3) flux tube, which is in turn defined in terms of the gauge invariant correlation function ρ conn W,μν of Eqs. (1) and (2) and points in a single color direction parallel to the color direction of the source (which is determined dynamically). Its six tensor components (the electric and magnetic fields E and B) correspond to the six orientations of the plaquette U P relative to that of the Wilson loop (see Fig. 1a).
The simulated fields E and B have the space-time symmetries of the Maxwell fields of electrodynamics, while carrying color charge in a single direction in color space. The energymomentum tensor T μν lies in the same direction in color space as the simulated fields E and B and has the (Euclidean) Maxwell form: Its spatial components μ = i, ν = j, with i, j = 1, 2, 3 determine the Maxwell stress tensor: Taking μ = i and ν = j = i in Eq. (A.1) gives while the diagonal time component −T 44 of T μν determines the energy density, We use cylindrical coordinates (x, r, θ), r ≡ y 2 + z 2 , tan θ ≡ z/y, and the corresponding unit vectorsê r ,ê θ : e r =ê y cos θ +ê z sin θ, (A.5) e θ = −ê y sin θ +ê z cos θ (A.6) (x is the longitudinal direction of the flux tube, i.e. the axis along which the static sources are located). The force exerted by the antiquark on the quark can be expressed, by means of the stress tensor, as a surface force F acting on the infinite plane x = d 2 bisecting the line connecting the pair: , and taking into account that the measured magnetic field B is compatible with zero (A.10) The angular average over the radial vectorê r in Eq. (A.10) vanishes. Furthermore by symmetry the transverse field E r on the mid-plane x = d 2 vanishes, so that the quark-antiquark force in Eq. (A.10) becomes Replacing E x (r ) by the non-perturbative field E NP x (r ) in Eq. (A.11) gives the non-perturbative quark-antiquark force F NP , Eq. (A.12) determines the string tension σ in terms of the longitudinal component of the non-perturbative field E NP x (r ), the confining component of the SU(3) flux tube. We have already presented it, in a slightly different notation, in Eq. (8).
Using Eqs. (A.5) and (A.6) in Eq. (A.3) we can obtain the components of the Maxwell stress tensor in cylindrical coordinates: (A.16)  (27) The remaining non-vanishing component of T μν is the energy density T 44 , Further, we note that the trace of the stress tensor T μν evaluated from Eqs. (A.13)-(A.17) vanishes independent of the simulated flux-tube fields E x (x, r ) and E r (x, r ): We have calculated the non-perturbative content of T rr on the symmetry plane (where T rr = T θθ = −T 44 = −T x x ) versus r for three different values of the quark-antiquark distance: d = 0.51 fm (at β = 6.240), d = 0.69 fm (at β = 6.539) and d = 0.95 fm (at β = 6.299). Results are presented in Fig. 11, where also the full (non-perturbative plus Coulomb) content of T rr is shown.
The width of the energy density distribution T NP 44 can be obtained through Eq. (12), with E NP x replaced by T NP 44 as given in Eq. (A.18); results are presented in Table 6. Since T NP 44 is proportional to E NP x 2 the width of the T NP 44 component of the Maxwell stress tensor obtained from the nonperturbative field given in Table 6 is systematically smaller than the width of the nonperturbative part of the longitudinal chromoelectric field component E NP x given in Table 5. (The Fig. 11 The diagonal components of the Maxwell stress tensor recovered from the full field E x (filled circles) and non-perturbative field E square of the field decreases more rapidly with distance than the field itself.) We now compare the above results obtained using our measured flux tube fields to evaluate the 'Maxwell' energymomentum tensor T μν with the corresponding results of recent direct simulations [39] of the expectation value of the energy momentum tensor T YM μν in the presence of a quarkantiquark pair. The latter simulations, which measure the energy and stresses in all color directions directly, were carried out in the plane midway between the quark and the antiquark, for three values of their separation.
The tensor T YM μν has the form [60] T μν = F a μα F a αν − g μν F a αβ F a αβ /4 , (A. 21) where F a μν is the Yang-Mills field tensor in the adjoint representation of SU (3), where f abc are the structure constants of the SU(3) algebra. In our definition the field is squared after color projection, whereas in Eq. (A.21) the sum over color components is taken after squaring. Moreover, the stress tensor in Ref. [39] is renormalized (this motivates the superscript R in the formulas below). In [39] the expectation value of the energy-momentum tensor in the background of a quark-antiquark pair is denoted by T R i j (r ) QQ . We will use this notation in comparing our results (A.18) and (A.19) with that work (Fig. 10).
In [39] the r dependence of the components of T R μν (r ) QQ was plotted for the 3 values of the quark-antiquark separation for which simulations were made, and these 'noticeable' features of the results were pointed out: 1. Approximate degeneracies between temporal and longitudinal components and between radial and angular components are found for a wide range of r ; We emphasize that the two inequalities in Eq. (A.23) are general consequences of (A.18) and (A.19), independent of the values of the simulated field E x (x, r ). In contrast, a recent study [40] of the stress tensor distribution in the Abelian Higgs model found that these relations could only be satisfied within a very narrow range of the model parameters. is partially restored inside the flux tube. 3. Each component of the energy-momentum tensor at r = 0 decreases as the separation becomes larger, while the transverse radius of the flux tube, typically about 0.2 fm, seems to increase for large separations [31,61,62], although the statistics are not sufficient to discuss the radius quantitatively. We see some indication of the increase in the width of the distributions of the diagonal components of the Maxwell stress tensor in Fig. 12 and Table 6 for all ten values of the quark-antiquark separation. However, this width is greater than 0.2 fm, the transverse radius of the flux tube estimated by [39].
Combining Eq. (A.23) with Eq. (A.24), we obtain which can be clearly seen from Fig. 3 of Ref. [39], where the components of T R i j (r ) QQ were plotted. The 'Maxwell' stress tensor does not include the contributions to Eq. (A.21) of the fluctuating color fields not measured in our simulations. Comparison of the spatial distributions of the diagonal components of the Yang-Mills stress tensor with the corresponding distributions of the 'Maxwell' stress tensor then provides a measure of the contributions of the fluctuating color fields.