Color flux tubes in $SU(3)$ Yang-Mills theory: an investigation with the connected correlator

In this work we perform an investigation of the flux tube between two static color sources in four dimensional $SU(3)$ Yang-Mills theory, using the so called connected correlator. Contrary to most previous studies we do not use any smoothing algorithm to facilitate the evaluation of the correlator, that is performed using only stochastically exact techniques. We first examine the renormalization properties of the connected operator, then we present our numerical data for the longitudinal chromoelectric component of the flux tube, that are used to extract the dual superconductivity parameters. These parameters turn out to be quite different from the ones obtained in similar studies carried out using smoothing techniques.

In this work we perform an investigation of the flux tube between two static color sources in four dimensional SU (3) Yang-Mills theory, using the so called connected correlator. Contrary to most previous studies we do not use any smoothing algorithm to facilitate the evaluation of the correlator, that is performed using only stochastically exact techniques. We first examine the renormalization properties of the connected operator, then we present our numerical data for the longitudinal chromoelectric component of the flux tube, that are used to extract the dual superconductivity parameters. These parameters turn out to be quite different from the ones obtained in similar studies carried out using smoothing techniques.

I. INTRODUCTION
The investigation of the color flux tubes connecting static sources in non-abelian gauge theories has become a standard tool to study color confinement [1][2][3][4][5][6][7][8][9][10][11][12][13]. Indeed, in lattice simulations, color sources are seen to be connected by tube-like structures for all the values of the coupling constant (at zero temperature). This is a strong indication that the mechanism responsible for color confinement is the same at weak and at strong coupling, in which limit color flux tubes naturally emerge [14], and the area-law of the Wilson loops can be analytically proven [15].
To study flux tubes on the lattice we need an observable whose average value will provide us informations about the flux tube details: more specifically the value of this observable has to be related to that of the field strength in the background of a couple of static color sources. A quantity that satisfies this requirement can be built by using the correlator of a Polyakov loops pair with a plaquette: the pair of Polyakov loops represents a couple of static color charges (a Wilson loop is also often used for this purpose) while the plaquette probes the field strength in the background of the static charges.
This general idea is common to all numerical implementations, however in the literature two different ways of defining the basic correlator are present: the first possibility is to use the expression [1] ρ disc = Tr(P r )Tr(P † r ′ )Tr(U p ) Tr(P r )Tr(P † r ′ )) − Tr(U p ) , where P r stands for the Polyakov loop at spatial position r and U p for the plaquette operator; this is known as the "disconnected" correlator. Another possibility is to use the definition [5,6] ρ conn = Tr(P r LU p L † )Tr(P † r ′ ) Tr(P r )Tr(P † where N c is the number of colors and L is the parallel transporter associated to the path shown in Fig. 1; ρ conn is known as the "connected" correlator and L is often called the "Schwinger line" in the literature. Both ρ disc and ρ conn are related to color flux tubes but they are not equivalent; in fact they are associated to different physical observables. This can be readily understood by looking at their naive continuum limit, in which the plaquette operator U p is expanded in powers of the lattice spacing a: ρ disc scales to the continuum as a 4 and gives access to 1 Tr(F 2 µν ) (we are assuming U p oriented in the µν plane), ρ conn scales to the continuum as a 2 and it is linear in F µν (see e.g. [5,6] for more details).
Since ρ disc and ρ conn do not provide equivalent physical informations, the choice of the operator to be used requires some discussion. Two arguments that have been adopted in the past to advocate the use of one operator or the other are the following: on one hand the operator ρ disc is theoretically better understood, since it can be easily shown to be multiplicatively renormalizable (see the discussion in Sec. II), and the renormalization constant needed to cancel logarithmic divergences can be fixed by using a lattice sum rule (see e.g. [7]). On the other hand ρ disc is noisier than ρ conn , and noise reduction was the original motivation for the introduction of the connected correlator in [5,6]: since ρ disc probes the square of the field-strength it is more sensitive to ultraviolet (UV) fluctuations.
Graphical representation of the numerator of the first term of Eq. (2). Ellipses denote the Polyakov loops (at distance d from each other) and the string L connects one Polyakov loop with the plaquette, first reaching the midpoint between the Polyakov loops at d/2, then moving in the transverse direction for a distance dt. The plaquette is drawned parallel to the plane identified by the two Polyakov loops since this is the case that will be studied in this paper, but its orientation can a priori be generic.
The choice of the operator to be used was thus largely based on the importance attributed to fluctuations. Sometimes UV fluctuations have a prominent role in the physical phenomenon to be studied, a prototypical example being the fluctuation-induced broadening of flux tubes [16]. In these cases the choice of the disconnected operator is mandatory, and specific stochastically exact noise-reduction techniques have been typically adopted to measure it [17][18][19][20][21][22].
When fluctuations were not expected to be important for the physical problem studied, the operator ρ conn has been the most common choice [23][24][25][26][27], supplemented by the use of smoothing algorithms to reduce UV noise. In studies performed with dynamical fermions only ρ conn has been used so far [28,29], since the accessible statistics are much lower than in the pure glue case, and stochastically exact error reduction techniques are still not available.
An important point to be noted is that smoothing has always been used to reduce the effect of UV fluctuations in ρ conn , however this standard procedure can a priori also induce systematical errors in the flux tube measure. Indeed in [29] a worrisome dependence of the physical results on the amount of smoothing adopted was noted (see also [24] for the case of Yang-Mills theory).
The aim of this work is to study ρ conn without using any smoothing algorithm, in order to understand if the connected operator can be used in a coherent fieldtheoretical setup to extract physical quantities related to flux tubes. For this purpose we evaluate the connected correlator ρ conn in four dimensional SU (3) Yang-Mills theory using only stochastically exact techniques (i.e. multihit [30] and multilevel [31] algorithms). To phys-ically interpret these data we need to study the renormalization of the connected correlator ρ conn : in our data the singularities related to the continuum limit (that are usually hidden by the use of smoothing) are clearly visible and we need to take care of them. The renormalization of ρ conn is far less trivial than that of ρ disc , however we will show that ρ conn renormalizes multiplicatively, and it can be used to extract physically relevant informations.
The paper is organized as follows: in Sec. II we discuss the issues related to the renormalization of ρ conn , using arguments largely based on [32,33], where the renormalization of cyclic Wilson loops was addressed. In Sec. III we introduce the numerical setup adopted, we present the results obtained for the longitudinal chromoelectric field, and we discuss the physical implications of these results for the dual superconductor model of the vacuum. Finally, in Sec. IV we draw our conclusions.

II. RENORMALIZATION OF ρconn
In order to discuss the renormalization of ρ conn it seems appropriate to start by briefly recalling some general facts about the renomalization of loop operators.
A loop operator is a generalized Wilson loop, which in the continuum can be written as where C is a closed curve and P stands for path-ordered. The systematic study of the divergences associated to these operators in four dimensional gauge theory was initiated in [34], where it was suggested that W C is multiplicatively renormalizable if C is piecewise smooth and not self-intersecting. From a one loop computation with cut-off regularization two sources of divergences were identified in [34]: logarithmic divergences originate from the points at which C is not differentiable, while linear divergences are always present, they exponentiate and globally contribute with a term of the form exp(cL(C)/a), where L(C) is the length of the curve, a is the UV cut-off and c is a constant. If C is smooth it was shown in [35] that W C is finite at all orders of perturbation theory, after the usual charge renormalization is performed. The case of a nonsmooth curve was studied in [36], where it was proven that to each cusp with angle γ a multiplicative renormalization constant has to be associated, whose value depends just on γ. The case of self-intersecting curves was also studied in [36]: in this case operator mixing between operators corresponding to different color contractions at the crossing points has also to be taken into account. The final result is the following: if r intersection points (corresponding to the sets of intersection angles {θ 1 }, . . . , {θ r }) and s cusps (corresponding to the angles γ 1 , . . . , γ s ) are present, then renormalization matrices and renormalization constants exist such that every color contraction W i1,...,ir C can be renormalized by using where the exponentials associated to linear divergences are implied. In [36] the possibility for different color contractions to have different lengths was however not considered, and this could seem to be a source of problems for the renormalization of ρ conn . Let us start by studying the renormalization of ρ conn in a scheme in which no power-law divergences are present (like e.g. minimal subtraction), so that we can use Eq. (4) without worrying of the complications related to linear divergences. In this case Polyakov loops do not need any renormalization, since they are associated to smooth contours, and the denominator Tr(P r )Tr(P † r ′ ) of Eq. (2) is finite once charge renormalization is performed. For the same reason the term Tr(P † r ′ ) in the numerator is also harmless and we can just concentrate on the term The path associated to O 1 is not smooth at three point: the point at which L and L † connects to P r , the point at which they connects to U p and the corner point of L and L † . All the rest of the contour contributes only to linear divergences, that we are neglecting for the moment.
In studying the renormalization of O 1 we have to take into account the mixing with all the operators that can be build from O 1 using different color contractions at the crossing points. Eight color contractions can be built (2 different contractions for each crossing point), however it is easy to realize, using LL † = 1 and analogous relations, that all the contractions that are not equal to O 1 are equal to where the 1/N c factor is needed to keep the same normalization of Eq. (5). From the previous general discussion it follows immediately that O 2 is multiplicatively renormalizable, a fact that will be used soon. Eq. (4) implies that we have to consider in general an 8 × 8 mixing matrix, which is the tensor product of the three basic 2 × 2 mixing matrices, however very stringent constraints are imposed on this 8 × 8 matrix by the fact that only two color contractions are different from each other, and by the fact that O 2 is multiplicatively renormalizable.
Let us discuss explicitly the case d t = 0 (see Fig. 1), in which we have just the tensor product of two 2 × 2 mixing matrices and 4 color contractions. Denoting the two 2 × 2 mixing matrices by Z A and Z B we thus have 2 where the apex (R) stands for "renormalized" and Since O 2 renormalizes multiplicatively, the coefficients Z A 10 and Z B 10 have to vanish, otherwise O (R) 2 would depend also on O 1 . If we use Z A 10 = Z B 10 = 0 we see that the last three lines of Eq. (7) are consistent with each other only if It is then simple to verify that Eq. (7) collapses to where As a consequence we finally have which means that ρ conn renormalizes multiplicatively. The same argument (which is an adaptation of the one used in [32]) can be repeated without changes also when d t > 0, in which case we have to start from the tensor product of three 2 × 2 mixing matrices. Since the renormalization constants only depend on the set of intersection angles, the value of Z 1 is the same for all the positive d t values, but it differs from the one at d t = 0, due to the presence of new logarithmic divergences for d t > 0.
Let us now consider a renormalization scheme in which linear divergences are present. The argument to show that ρ conn is multiplicatively renormalizable also in this case is exactly the same that was used for cyclic Wilson loops in [33], to which we refer for further details and some enlightening examples. Here we will just briefly sketch some basic steps of the proof for the benefit of the reader and to fix the notation. The main ingredient that is needed is the relation between an operator acting in the fundamental representation of SU (N c ) (that we will denote by U ) and the corresponding operator acting in the adjoint representations (that we will denote by U adj ), which is where T a are the SU (N c ) generators, with the normalization Tr(T a T b ) = 1 2 δ ab . By writing P r and U p in the base {1, T a } and using Eq. (13) it is simple to show that We thus see that O 1 − O 2 can be written as a single generalized loop function in which traces in different representations are present: two traces in the fundamental representation (explicitly denoted by Tr) and a trace in the adjoint representation (the summation on a, b). This expression is particularly convenient since it was shown in [33] (using the exponentiation theorems of [37,38]) that linear divergences factorize also in untraced loop operators, and they can be cancelled by multiplicative factors of the form exp(−c r L(C r )/a), where c r is a representation dependent coefficient, a is the UV cut-off, and L(C r ) is the length in physical units of the curve C r associated to the representation r.
Applying this result to the operator where c f and c a are the constants associated to the fundamental and adjoint representations, d 2 +d t is the length of the contour in the adjoint representation (see Fig. 1) and L t is the temporal extent of the lattice in physical units. Z 1 is the renormalization constant associated to logarithmic divergences that was previously introduced in Eq. (11), while the term 4c f is associated to the plaquette, that has length 4a, does not generate linear divergences in the continuum and can be safely neglected. Since to each Polyakov loop a factor e −c f Lt/a is associated, we finally have where Z 1 is independent of d and assumes two different values for d t = 0 or d t = 0.

A. Setup
As anticipated in the introduction, we study ρ conn in four dimensional SU (3) Yang-Mills theory, discretized on the lattice by using the standard Wilson action [39]: The connected correlator defined in Eq. (2) can be estimated by using a combination of multihit [30] and multilevel [31] techniques. The multilevel method can be easily adapted to measure ρ conn , since LU p L † involves variables living in a single time slice of the lattice; multihit can be applied to the links of the Polyakov loops and to a single link of the plaquette. For all the cases studied in this work a single level of the multilevel algorithm was used in simulations, with slices of dimensionless temporal extent ∆ = 4; the number of internal updates of the multilevel algorithm was optimized by using a technique analogous to the one discussed in [31].
Simulations were performed on symmetric lattices (N s = N t ), and in Tab. I we report the simulation details. The plaquette U p which appears in Eq. (2) was always chosen to be parallel to the plane identified by the two Polyakov loops, since we are interested in studying the longitudinal chromoelectric field, which in all previous studies was shown to be the dominant component of the flux tube.
The coupling values were chosen in such a way that d = 4a, 6a and 8a correspond to the same distance in physical units; for this purpose the parametrization of a(β) obtained in [40] was used. Parameters in Tab. I fix the physical distance between the Polyakov loops to about 0.37 fm (using r 0 ≃ 0.5 fm for the Sommer scale [41]), and the dimensionless lattice size was rescaled in order to have (almost) constant physical volume. Test simulations were also performed on the lattice 28 4 with coupling β = 6.0, which excluded the presence of sizable finite size effects in our data. Data points corresponding to different values of d t have been extracted from independent simulations, hence they are statistically independent from each other.

B. Results for the chromoelectric field
From the naive continuum limit of ρ conn in Eq. (2) we can define the longitudinal (due to our choice of the plaquette orientation, see Sec. III A) chromoelectric field by using the expression where in our simulations d is fixed to about 0.37 fm, and d t denotes the transverse distance from the center of the flux tube (see Fig. 1). From the discussion in Sec. II it follows that we can not expect E L defined in this way to have a nontrivial continuum limit. Indeed data shown in Fig. 2 indicate that E L (d, d t ) converges to zero as the continuum limit is approached. To properly define the continuum limit of E L we have to use Eq. (16) and define where Z 1 (a, d t ) is the renormalization constant associated to logarithmic divergences (which is different for d t = 0 and d t = 0), and we introduced the shorthand to denote the multiplicative factor needed to remove linear divergences. To completely define E R L we thus have to fix the three constants c a and Z 1 (a, d t ) (for d t = 0 and d t = 0).
The numerical value of Z 1 (a, d t ) could be computed in perturbation theory, however some interesting physical observables can be studied also without a precise knowledge of this renormalization constant. This is due to the fact that, for d t > 0, Z 1 is a multiplicative factor independent of d t , hence the functional form of E (R) L (d, d t ) for d t > 0 is completely fixed also without any knowledge of Z 1 . However, also to study just the functional form of E (R) L (d, d t ), we need to fix c a . Since c a is a fundamental property of the discretization adopted, independent of the specific adjoint loop function used and of the infrared properties of the theory, we have the freedom of choosing the simplest numerical setup available to fix its value. We decided to extract if from the continuum scaling of the Polyakov loop in the adjoint representation at finite temperature, which is an observable that is easily computed to high precision. For this purpose we performed simulations starting from a 4×16 3 lattice at β = 5.8 (the deconfinement transition on N t = 4 lattices takes place at β c = 5.6925(2), see [42]), then increasing the value of N t keeping the physical temperature constant and the aspect ratio fixed to 4, see Tab. II. The average value of the Polyakov loop in the adjoint representation can be computed by using the relation which is an easy consequence of Eq. (13), and from the discussion in Sec. II if follows that TrP adj scales to the continuum as exp(c a N t ). Numerical results obtained for TrP adj are shown in Fig. 3, in which some deviations from the asymptotic exp(c a N t ) behaviour are also visible. To extract the value of c a fits of the form have been performed, and the stability of the fit under changes of the fit range and of the functional form adopted (i.e. by setting k 2 = 0) has been investigated. As our final estimate we report the value c a = −0.45 (2) .
Using this value for c a we can remove linear divergences from E L (d, d t ), and in Fig. 4 the values of Z L (d, d t , a)E L (d, d t ) are shown as a function of d t (in physical units). The lattice spacing dependence of data in Fig. 4 is much milder than that observed in Fig. 2, however we have to remember that the renormalization factor Z 1 (a, d t ), needed to take care of logarithmic divergences, is still missing. A consequence of this fact is that the scaling to the continuum of points associated to distances d t = 0 and d t > 0 is different, as can be seen in Fig. 4 and will be most clearly evident from the considerations of the next section.

C. Dual superconductivity parameters
According to the dual superconductor model of color confinement [43][44][45] the vacuum of nonabelian gauge theories behaves as a "dual" superconductor, in which condensation of chromomagnetic degrees of freedom produces a "dual" Meissner effect, that squeezes the chromoelectric field lines into flux tubes producing confinement. The characteristic feature of this model is to provide a conceptually simple and physically appealing framework to interpret some nonperturbative aspects of gauge theories.
A generic feature of any superconductor (dual or not) is the presence of two typical lengths in the infrared effective theory: the coherence length ξ and the penetration length λ. The values of these lengths characterize the functional form of the flux tube profile, with the penetration length being associated to the exponential decrease of the field far from the center of the flux tube. In order to determine both ξ and λ starting from the data presented in the previous section, we will follow the approach first adopted in [24] (and then used in [25,26,28,29] Fig. 4, together with 1σ bands obtained from the combined fit described in the text.
In this approach the following parametrization of the longitudinal component of chromoelectric field inside the flux tube is used: where K 0 , K 1 are modified Bessel functions of the second kind, and φ, α and µ are fit parameters. This is the "dual" version of the parametrization introduced in [46] for the longitudinal magnetic field inside a vortex line in type II superconductors. The parameter µ is just the inverse of the penetration length, µ = 1/λ, while the relation between the fit parameters in Eq. (24) and the coherence length is less direct: it can be shown (see [46]) that the Ginzburg-Landau parameter κ = λ/ξ is related to α by the relation Values of κ smaller than 1/ √ 2 correspond to superconductors of type I, while κ > 1/ √ 2 for type II superconductors (see e.g. [47]).
If we try to fit each of the fixed β data sets shown in Fig. 4 by using the parametrization in Eq. (24), we immediately realize that the quality of the fits degrades as the coupling is increased. This is a consequence of the previously noted fact that data points at d t = 0 and at d t > 0 scale to the continuum in different ways, due to the different logarithmic divergences in the two cases. If on the other hand we simply discard the point at d t = 0 from each of our data sets, the precision of our data is not enough for the fit to provide significant information on α, and consequently on the Ginzburg-Landau parameter κ.
We thus decided to perform a combined fit of all our data at d t > 0 keeping three different φ parameters, corresponding to the three lattice spacings used, since φ is sensitive to the multiplicative renormalization Z 1 (a, d). Using c a = −0.45 we obtain the best fit shown in Fig. 5 and the final results for the superconductivity parameters are while the values obtained for the φs parameters shown in Fig. 6. The χ 2 test for this fit gives χ 2 /dof = 16/10 that is somehow large but still acceptable, and the final results are almost unchanged also for c a = −0.47 and c a = −0.43, which means that the uncertainty in c a is not the main source of error in our final results.

IV. CONCLUSIONS
In this paper we presented the results of our study of the longitudinal chromoelectric component of the color flux tube, performed by using the connected correlator. Measures were carried out by means of stochastically exact techniques, without any smoothing, in order to investigate the possibility of using the connected correlator ρ conn in a coherent field-theoretical setup.
We first investigated the renormalization properties of ρ conn , showing that it is multiplicatively renormalizable, and reducing the problem of its renormalization to the determination of three renormalization constants. One of these constants (denoted by c a ) is related to linear divergences, while the other two take care of the logarithmic divergences in the cases d t = 0 and d t = 0 respectively.
We then fixed the value of c a , by studying the β dependence (at fixed temperature) of the Polyakov loop in the adjoint representation. Using the value of c a obtained in this way (c a = −0.45(2)), we removed linear divergences from ρ conn , obtaining the results shown in Fig. 4. While it is important to stress that these data are still not renormalized (since logarithmic divergences have not been removed), it is possible to extract from them quantities of direct physical interest.
In particular, starting from the functional form of the longitudinal chromoelectric field, we evaluated the coherence length ξ and the penetration length λ of the dual superconductor model. The numerical values of these quantities, extracted using a flux tube of length d ≃ 0.37 fm, are reported in Eq. (26). Our estimate of the penetration length λ is in good agreement with previous determinations. Our Ginzburg-Landau parameter (and consequently our coherence length) is instead quite different from the one obtained in similar studies carried out by using smoothing [24,25] (where κ ≈ 0.2 was found), being closer to older results suggesting the SU (2) vacuum to be at the boundary between type I and type II superconductivity [48][49][50][51][52][53].
This result suggests that smoothing could introduce some systematics in the determination of the flux tube profile, and we get the following intuitive picture: ξ is the typical scale of the bulk of the flux tube, which broadens under smoothing, while λ is related to the large distance behaviour of the tails, which is almost unaffected by smoothing. As a consequence we expect smoothing to leave almost unaltered λ and to decrease κ = λ/ξ.
While this picture seems appealing we also have to keep in mind the limitations of our computation: first of all our determination of κ has a 30% relative error, so this effect could just be a statistical fluctuation. Moreover the length of the Polyakov loop used to extract ξ and λ was about 0.37 fm, which is surely not asymptotically large; as a consequence a contamination from the Coulomb component of the flux tube is possible (see [27] for a discussion on this point). To improve on these limitations is surely matter for further studies, just as the determination of the renormalization constants associated to the logarithmic diverges of ρ conn .

ACKNOWLEDGMENTS
It is a pleasure to thank Massimo D'Elia and Alessandro Papa for useful comments. Numerical simulations have been performed on the CSN4 cluster of the Scientific Computing Center at INFN-PISA, and on the MAR-CONI machine at CINECA, based on the agreement between INFN and CINECA (under project INF19 npqcd).