The Coulomb flux tube revisited

We perform $SU(2)$ Yang-Mills lattice simulation of the electric field distribution in the Coulomb gauge for different values of $\beta$ to further investigate the nature of the Coulomb flux tube.


I. INTRODUCTION
In the SU (N ) Yang-Mills theory in the Coulomb gauge, there is an instantaneous interaction between color charges, which is analogous to the Coulomb force between electric charges. The key difference, however, is that unlike QED where the Coulomb energy is a function of the distance between sources, in the non-abelian case it is a functional of the gauge field. Therefore in order to obtain the Coulomb potential, it is necessary to specify the state of the gluons. The Coulomb potential is conventionally referred to the state in which external, static sources are suddenly added to a vacuum. In the following, we refer to this as the bare state. This is different from the adiabatic situation, when gluons have time to respond to the presence of the external charges resulting in the true QCD eigenstate. We refer to the later as the minimal energy or Wilson state, since its energy is related to the expectation value of the large Wilson loop.
The Coulomb potential, V C (r), evaluated on a state containing a static quark-antiquark pair is of special interest. In the past few years, it has been extensively studied [1][2][3][4][5][6][7][8][9][10][11][12][13][14] both within continuum and lattice approaches, which contributed to a better understanding of the quark confinement [15][16][17]. In particular, it has been shown that the Coulomb confinement is necessary for the Wilson confinement [4]. This has been confirmed on the lattice [6,7] where it was found that the Coulomb potential rises linearly for large quark-antiquark separation r, with the associated string tension σ C larger by approximately a factor of three compared to the minimal one, σ, obtained from the expectation value of the large Wilson loop.
Since both potentials are confining, it is reasonable to ask how other gauge-field-related observables compare in the two states. Numerical studies with SU (2) and SU (3) lattice Yang-Mills theories have established a picture of a flux tube formation between a quark and an antiquark in the minimal energy state [18][19][20][21][22][23][24][25][26][27][28][29]. Phenomenologically, it was established that both the action and energy densities vanish exponentially in a direction perpendicular to the line joining the quark sources. In principle, these observables can be obtained by measuring a (normalized) * email: sdawid@iu.edu correlation function between a large Wilson loop W (r, t) and appropriately placed plaquette U P , which serves as a chromo-electric (or chromo-magnetic) field probe [27].
The analogous question concerning the bare QQ state was recently addressed by Chung and Greensite in Ref. [30]. There it was found that also the bare state has the flux-tube-like characteristics with an exponentially decaying transverse profile. This is an interesting and unexpected result since analysis of the Coulomb energy density distribution and related observables, e.g. the ghost propagator, in the infinite volume [3,8,31,32] typically predicts a power-law fall-off. Also one could argue that the Gribov-Zwanzinger confinement proposal [15,16] implies long-range Van der Waals forces and thus the absence of flux tubes [8].
The aim of this paper is to shed more light on this result. Specifically, we extend the calculations of [30] that were performed at β = 2.5 to a considerably better precision for data at larger transverse distances, y away from the QQ axis, and we also performed the calculation for β = 2.3 and 2.7. Additionally, to understand the flux tube development, we investigated the Euclidean time evolution of the energy density profile. Finally, we performed an analysis of the data using both power-law and exponential profiles.
The paper is organized in the following way. In Sec. II we give a summary of the lattice setup and describe the measured observables. In Sec. III A we present the key results of our simulations for different β's, and discuss the analytic models. Results of the Euclidean time evolution are presented in Sec. III B followed by summary and conclusions in Sec. IV. The complete data set is given in App. A.

II. ELECTRIC FIELD DISTRIBUTION IN THE PRESENCE OF STATIC QUARKS
In Ref. [6] it was shown that on the lattice in the Coulomb gauge, both the Coulomb and Wilson energies can be calculated from the expectation value of two Wilson lines: where a is the lattice spacing, and L t (x) is a time-like Wilson line of length t starting at position (0, x). In the (Euclidean time) limit t → ∞, potential V (r, t) becomes the Wilson eigenenergy V min (r). In the limit t → 0 this quantity, up to an additive, r-independent constant, approaches the lattice version of V C (r), defined as a correlation of short time-like links where U µ (x) is a link variable at position x = (x 0 , x) in the direction of µ ∈ [0, 3]. Four-vectors (0, 0) and (0, r) represent positions of the quark and the antiquark, respectively. One can understand both limits as a starting and ending point of an equilibration process, which takes a set of gauge fields unperturbed by the presence of the QQ pair and thermalizes it to the true ground state of the theory. As mentioned earlier, in the SU (2) Yang-Mills theory, it was found that at fixed value of β, the string tension σ C computed from V (r, t) as a function of r decreases with increasing t and approaches that of the Wilson energy at large times. In the Coulomb limit t → 0, it is found [6] to be larger by approximately a factor of three, see Fig. 6.
In the Coulomb gauge, the longitudinal component of the chromo-electric field E L is determined by charge distribution (as in the classic theory) via the Gauss' law. Thus one can determine distribution of this field in a state at any time, as the expectation value of TrE 2 L . Since it is expected that the main contribution to the energy density comes from the field component parallel to the QQ axis, which we assume to lie in the x direction, here we calculate this one component contribution as in Ref. [30] using where we have switched to dimensionless distances R = r/a and T = t/a. Plaquette U P (p, T ) is defined analogously to U P (p, 0) in Eqs. (25) and (26) of Ref. [30], i.e. it is oriented in xt-plane and placed at position but at a different time slice: T /2 for even T , or (T − 1)/2 for odd T . Here y = y ⊥ /a is a distance from the QQ axis. The plaquette acts as a probe of the x-component of the longitudinal chromo-electric field at point p, therefore allows to see how Q T changes with the transverse distance y. The key quantity, which is the numerator of the first term in Eq. (3), is depicted schematically in Fig. 1.
The observable Q T (R, y) is a generalized version of Q(R, y) introduced in Eq. (24) of Ref. [30], and reduces to the latter for T = 1, i.e. when the Wilson line is equal to one temporal link L T =1 (x) = U 0 (0, x). From Q(R, y), in Ref. [30], the energy distribution in the Coulomb state was obtained. The generalization given by Eq. (3) allows us to make a connection with the minimal energy flux tube measured as a normalized correlation function between a large Wilson loop and a plaquette [27]. Since in the Coulomb gauge spatial links become very close to the identity matrix, the Wilson loop can be approximated as a product of two temporal Wilson lines, and thus one should be able to observe the convergence of the bare state field distribution to the Wilson state field distribution as increasing T is considered.

A. Lattice framework
We performed Monte Carlo simulation of the pure SU (2) Yang-Mills theory using Wilson's action [33]. We used lattice of size V = 32 4 with periodic boundary conditions, at couplings β = 2.3, 2.5, 2.7, which correspond to lattice spacings a = 0.165, 0.085, 0.045 fm, respectively [34]. Field configurations were generated using the heat-bath algorithm [35] and we considered the lattice equilibrated after initial 10 000 sweeps for each β. Lattice configurations used for data extraction were separated by 300 sweeps to minimize the impact of autocorrelations. The Coulomb gauge was assumed to be fixed is a value of the functional to be minimized, after the i-th gauge fixing iteration. We found that the more strict condition ∆F < 10 −8 did not affect the results noticeably.
To speed up the gauge fixing procedure we implemented an overrelaxation method [36,37] with ω = 1.75. As a check, we repeated the simulation of Ref. [30] with increased statistics by generating 30 000 lattice configurations at β = 2.5, to obtain Q T =1 (R, y). In addition we used approximately 11 000 configurations to compute Q for β = 2.3 and β = 2.7. The Euclidean time dependence was obtained also from the same number of configurations at β = 2.5 and for T = 1, 2, 3, 4, 5. For each configuration, we averaged the value of the observables over four possible translations and three 90 • spatial rotations. Expectation values of observables and statistical errors were obtained from the jackknife method. We do not investigate here the systematic errors due to e.g. finite lattice size or gauge fixing quality and assume they do not affect the main conclusion of our work.

A. Results for different lattice couplings
In the case of the Wilson state, it was shown in [24] that the profile of the Wilson flux tube, for small separations R, can be calculated from perturbation theory 1 and to the leading order in α s , falls off as 1/y 6 . For large quark separations it is observed that the Wilson energy density   profile changes from power-law to exponential [27]. In this subsection we discuss the energy distribution in the Coulomb state, Q T =1 (R, y) calculated at different values of the coupling β. To understand dependence on the transverse distance y, we first try two simple models, one is the power law, (PL) Q PL , motivated by [24] and the other exponential (CG) Q CG used in Ref. [30] While b = 3 is predicted for the Wilson state, analytical calculations in the Coulomb gauge predict b ≈ 2, independent of R [8]. Thus in the Coulomb gauge one might expect power-law behavior in y for small values of R and as a consequence, it is the y-behavior of the energy profile at large R that should be examined to discriminate between an exponential and power-law decay. Our most precise measurement was performed at β = 2.5. Sample plots of the energy transverse profile for ("small") R = 2 and ("large") R = 7 are shown in Fig.  2, and for all other quark separations are summarized in Fig. 9 in App. A. The fit parameters, fit intervals and corresponding values of χ 2 /d.o.f. for CG and PL models are presented in Tabs. I and II. Details of the fitting procedure are given in the App. A. It appears that the data favor the PL model over the simple exponential (CG), however, it does so with varying quality depending on R. Specifically, for large R, we could not reach conclusive results and both models appear to be deficient. This may indicate that at β = 2.5, Q T =1 (R, y) does not properly reproduce the Coulomb energy density distribution. Since the Coulomb energy is obtained in the continuum limit, or at least for large values of β, it would appear that β = 2.5 is not large enough.
We thus hypothesize that Q T =1 (R, y) evolves from a power-law behavior for large β to an exponential (small β) and consequently performed a calculation at β = 2.3 and β = 2.7. These values correspond to approximately 3.89 (7) 1.89 (7)  We followed the fitting procedure from Ref. [30] and deviated from it only for R = 6, 7, 8, for which fit intervals are [2,7] compared to [1,7] in Ref. [30]. For R 5 our values agree with Tab. I in Ref. [30]. constant ratios of the lattice spacing a β=2.3 /a β=2.5 ≈ a β=2.5 /a β=2.7 ≈ 2. For example, the energy distribution profile at β = 2.5 for R = 2 can be compared with that at β = 2.3 for R = 1, and at β = 2.7 for R = 4, and all of them correspond to r = aR ≈ 0.17 fm. In Fig.  3 we show the energy profiles for the three values of β.
One can see that indeed for larger β's the Coulomb flux tube profile appears to follow a power-law, while at the lowest value, β = 2.3, already starting at low values of R, R = 2, it appears much closer to an exponential -as discussed above for small R the y-profile is expected to follow a power-law. Taking a closer look at the β = 2.3 data, shown in Fig. 8 in App. A, we see a clear exponential decay for large R. To see if the energy distribution for this value of the coupling β can be described by the Wilson flux tube profile, we employed an "improved" exponential model (CCB): proposed 2 in Ref. [27]. Here ν and λ are free parameters and their values for our data set can be found in Tab. III. The CG and CCB models are similar, as they both describe the exponential decay for large y. However CCB model also includes flattening of the energy profile close to the QQ axis. From the plots in Fig. 8 it can be seen that the fit accuracy is improved as R increases. For R = 1, were a power law is expected we fitted the β = 2.3 data with PL model in an interval y ∈ [3,8], obtaining b = 2.83 (7), which is close to the perturbative prediction of b = 3 for the minimal state energy density shape. Even though it might look that the plots for R = 2 and even R = 3 might follow a power law, we were not able to obtain a good fit with such models. For β = 2.7 the energy density transverse profile appears to follow a power law for all values of R and neither CG nor CCB model gives a comparable description, see Fig. 4 and Fig. 10 in App. A. The best PL fit parameters can be found in Tab. IV. The fit indicates that the profile falls off approximately as 1/y 4 with a transverse distance y, in agreement with Ref. [8]. The exact power depends on an interval used for fitting, converging to 2 when more of the low-y points are excluded from the fit. The profile of the bare state energy density close to the QQ axis could not be satisfactorily reconstructed by the straightforward PL model, which should be improved e.g. by inclusion of a factor correcting the small-R and small-y dependence. Such an factor should also mimic the change of the power b ≈ 2.7 for small R to b ≈ 2 for large R. We have tried to employ a perturbative prediction from Eq. (28) of Ref. [8] obtained from the Dyson-Schwinger equations in the Coulomb gauge. It did not lead to a con-  clusive results, predicting even more sever flattening of the small-y shape at any value of R, than the PL model. In particular, it predicts a fall-off of the energy density at the middle point between quarks, y = 0, as Q ∝ R −2 , while our data set favors a higher power, see Fig. 5.

B. The Euclidean time dependence
As already discussed, the time development of the Coulomb potential and its evolution towards the minimal potential in principle can be analyzed by considering the large t limit in Eq. (1). We extracted the behavior of the string tension as the Euclidean time progresses and observe its convergence to the minimal string tension, see Fig. 6, which agrees with Ref. [6].
We have already seen that, as one takes relatively small value of β, it is possible to obtain an exponential fall-off of the bare state energy density transverse profile for large values of R. An interesting question arises if it converges to the minimal flux tube when the limit T → ∞ is taken in Eq. (3) for larger values of β, for which a power-law behavior was found instead. As already discussed, in the  Coulomb gauge the spatial links can be approximated by the identity matrix. If we thus made an assumption that the Wilson loop W (R, T ) ≈ Tr L T (0)L † T (R) , then for large values of T Eq. (3) would agree with Eq. (7) from Ref. [27], i.e. the formula for the chromo-electric contribution to the energy density profile of the minimal flux tube. To answer this question we investigated Q T (R, y) for T from T = 1 up to T = 5, and the results are presented in Fig. 7 for R = 2 and R = 7, and Fig. 11 in App. A for all R's. From the plots it seems that the quantitative behavior of the profile does not change with the growing T ; it increases globally (i.e. for all values of y) but the functional form of the profile appears unaltered. This is quite surprising. It might indicate that the flux tube described by Q T =1 at β = 2.5 is already "equilibrated" in a sense that it represents closely the minimal flux tube rather than the bare state energy density distribution, and one should move closer towards the continuum limit (larger β) to see a noticeable difference. This requires further investigation, with a better statistics and in a larger volume, where larger values of R, y and T can be considered. and a = 0.20 (2). The R = 1 point was excluded from the fit, because it does not satisfy µR 1 condition from Ref. [8] for µ = 0.63. The dotted, orange line represents theoretical prediction βQ = 32σC/π 3 R 2 . The Coulomb string tension σC at β = 2.7 was obtained as a by-product of the energy density profile calculation from Eq. (2). It was extracted from the fit VC(R) = σCR + β/R + γ, with free parameters σC, γ, β, and is equal to σC = 0.0409(3). String tension T asymptotic value tension Figure 6. Dependence of the string tension on T . It was obtained following the procedure of Ref. [6], for β = 2.5 and around 1700 lattice configurations. The horizontal line is the asymptotic value from Ref. [24] IV. CONCLUSIONS In this paper we have extended the calculation of [30] to other β values. We found that our numerical results for Q T =1 at β = 2.5 agree with those presented in Ref. [30] with increased statistics. The concussion from [30] was that the Coulomb flux tube vanishes exponentially with the transverse distance, and has a width larger than the minimal flux tube. This is difficult to reconcile with theoretical prediction. We performed a somewhat different analysis of the same quantity by emphasizing large-y values. Furthermore by considering other values of β we showed that with increasing β the genuine, it is likely that Coulomb flux tube with power-law fall-off develops. As a consequence we were concluded that it is possible that the energy density profile evolves from the Wilsonlike to the Coulomb one in the continuum limit. This was supported by our study of the Euclidean time development of the profile, however, quantitative analysis requires better statistics and in larger volumes.
The power-law model was not able to describe the small-y data properly, predicting too small values of the energy density near the QQ axis. The same happened with the theoretical prediction of Ref. [8]. This might indicate that on the axis the Coulomb energy density contains a significant non-perturbative contribution which should be explained.

V. ACKNOWLEDGEMENTS
We would like to thank K. Chung  On the next pages, we present our data set of Q T =1 (R, y) for all values of β, for R, y ∈ [1,8], together with the fits. Also, in Fig. 11 we show the plots of Q T (R, y) at β = 2.5 for different times T .
Apart from the CG model of Eq. (7) at β = 2.5, the fits for each R were performed in an interval y ∈ [y min , 8], with y min = 1, 2 or 3 depending on R. We claim it is important to include large values of y (as long as they are not affected by the finite volume effects) to discriminate between power-law and exponential behavior. In general only for large transverse distances y the difference between various exponential and power-law models might become significant. Moreover, examining Fig. 2(b), one can see that for small y and large R the lines have curvature and are flatter compared to a simple exponential. This indicates that close to the quark-antiquark axis the exponential model (CG) may not be accurate. Thus the small y should be excluded when fitting with this model. In practice, none of the models we used was successfully in describing the data y-dependence close to the QQ axis, and as a consequence, to obtain satisfactory values of χ 2 /d.o.f. We had to remove points y = 0, 1, and sometimes y = 2, from fit intervals. These points usually come with a small relative errors, thus they affect the value of χ 2 significantly. The fit parameters are summarized in Tabs. I, II, III, IV.
As an illustration of how a choice of the fitting interval affects the quality of the fit, for the exponential model Q CG we followed the fit procedure of Ref. [30]. Specifically, we used the same intervals as in [30] for small R, up to R = 5. For larger values of R, one may argue that as function of y, Q(R, y) falls off exponentially for y 2 and we excluded the points y = 0, 1 from the fit. Even with such an "optimized" data set, we find the resulting values of χ 2 /d.o.f to be quite large as shown in Tab. I. The red line is the best CCB fit to the data. It can be seen, how the results start agreeing with this model as quark and anti-quark are separated further away. For R = 1 we used the PL fit in the interval [3,8], obtaining a = 0.23(5), b = 2.83 (7) and χ 2 /d.o.f.= 0.71. The PL model was not able to describe the behavior of the transverse profile for R 2 and thus is not included on the graphs. The fit parameters are given in Tab Figure 9. Results for Q1(R, y) for different quark sperations R, obtained for β = 2.5 and 30 000 lattice configurations. The red, solid line is the best exponential (CG) fit to the data, performed following the procedure of Ref. [30] for R 5. The green, dashed line is the best PL fit. Fit parameters for both fits are given in Tabs. I and II.  (h) R = 8 Figure 11. Results for QT (R, y) for different quark sperations R and different T 's, obtained for β = 2.5 and around 11 000 lattice configurations. Lines joining the data points are included to guide the eye. The data point for R = 5, y = 8 and T = 4 was negative, so is not visible on the semi-log scale.