Pure gauge QCD flux tubes and their widths at finite temperature

We study the flux tubes produced by static quark-antiquark, quark-quark and quark-gluon charges at finite temperature in pure gauge SU(3) lattice QCD. Our sources are static and our lattice correlators are composed of fundamental and adjoint Polyakov loops. To signal the flux tubes, we compute the square densities of the chromomagnetic and chromoelectric fields with plaquettes, in a gauge invariant framework. We study the existence and non-existence of flux tubes both above and below the deconfinement phase transition temperature Tc. Using the Lagrangian density as a probability distribution, we also compute the widths of the flux tubes and study their widening as a function of the intercharge distance. We determine our results with both statistical and systematic errors. Our computations are performed in NVIDIA GPUs using the CUDA language.


I. INTRODUCTION
The understanding of confinement and deconfinement in QCD remains a central problem of particle physics. A major evidence of QCD confinement is the flux tube arising between quark-antiquark static charges, both from gauge invariant pure gauge lattice QCD simulations [1][2][3][4] and from experimental observations like Regge trajectories [5][6][7][8][9] consistent with linear confining potentials. Even in dynamical QCD where the flux tube breaks due to the creation of another quark and antiquark, a flux tube develops up to moderate quark-antiquark distances. Different flux tubes have also been shown to occur in lattice QCD simulations of different exotic hadrons, typical of SU (3) [10][11][12][13]. It is important for the fundamental understanding of the pure gauge QCD flux tubes to measure the flux tube profile [14][15][16][17][18][19][20] with more quantitative results.
Moreover, in the topic of High Energy Heavy Ion Physics, the interactions between heavy quarks are relevant for the hard probes of the QCD phase diagram [21], and the development of flux tubes between different charges may help to microscopically understand the phenomenological vortex line model for the flow and fragmentation [22]. Thus it is also important to study flux tubes, not only at zero temperature, but also at finite temperature.
In this work, we study whether flux tubes exist or not, between different static charges, at different temperatures above and below pure gauge QCD critical deconfinement temperature T c . In particular, we study the flux tubes created between static quark-antiquark, quarkquark and quark-adjoint charges.
Notice we do not know exactly, neither the theoretical origin of the QCD flux tubes not their effective behaviour, and it is important to explore in more detail their properties in order to, hopefully, one day solve this important * bicudo@tecnico.ulisboa.pt † nuno.cardoso@tecnico.ulisboa.pt ‡ mjdcc@cftp.ist.utl.pt The lattice is represented in black, the Polyakov loops and plaquette are painted in blue and the axis and relevant vectors are painted in red. R is the vector separating the two Polyakov loops,and r is the position of the plaquette. Because the Euclidian space-time is four-dimensional, we represent only one temporal dimension and two spatial dimensions. The plane in green is the projection of the three-dimensional space into the x, z dimensions. The charge axis is the z axis. The mediator plane is the x, y plane, here represented by a green line only since we do not show the y dimension.
problem. For instance, two qualitatively different effective models for the QCD flux tube exist. Already in the 1970's, Nambu [23], 't Hooft [24] and Mandelstam [25] proposed that quark confinement would be physically interpreted using the dual version of the superconductivity [26,27]. The QCD vacuum state would behave like an ordinary superconductor, where Cooperpair condensation leads to the Meissner effect, and the magnetic flux is excluded from the vacuum and squeezed arXiv:1702.03454v3 [hep-lat] 10 Oct 2017 in a quasi-one-dimensional tube, the Abrikosov-Nielsen-Olesen vortex, where the magnetic flux is quantized topologically [28][29][30]. In a superconductor, the fields are approximately classical and the flux tube main parameter is the penetration length λ in the London equation has a direct relation with an effective mass of the interaction particle fields, i. e., the photon. In QCD, the dual gluon mass has been studied by several authors, [31][32][33][34][35][36][37][38], as well as the gluon effective mass [39]. Interestingly, there is also an evidence for a gluon mass in the Landau Gauge [40] and in the multiplicity of particles produced in heavy ion collisions [41]. Recently, the penetration length started to be computed with gauge invariant lattice QCD techniques [42][43][44].
On the other extreme limit, at quark-antiquark distances larger than the penetration length, the flux tube is similar to a quantum string, contrary to the picture of the superconductor with classical electromagnetic fields. Due to its vibration, the quantum string has a Gaussian profile and a finite width, different from the exponential profile of the classical superconductor flux tube [45,46]. Thus a second model of the QCD flux tube is given by the string model, based on the Nambu-Goto Action [47,48].
At zero temperature, the energy of the quantum string with length R and fixed ends, with quantum transverse fluctuations quantum number n, is expressed in the Lüscher term and in the Arvis Potential [49,50], indeed observed in lattice QCD for 4 space-time dimensions [45]. In Eq. (1), D is the dimension of the space time. Note that the Arvis potential is tachyonic at small distances since the argument of the square root becomes negative in this limit, moreover rotational invariance is only achieved for D = 26. Nevertheless, the first two terms in the 1/R expansion, σR + π(n − D−2 24 ) 1 R , are more general than the Arvis potential, since they fit the D = 3 and D = 4 lattice data quite well down to distances much smaller than the Arvis tachyonic distance. The 1/R term is independent of the string tension σ and for the physical D = 3 + 1 has the value − π 12 . This is the Lüscher term [49]. The energy spectrum of a static quark-antiquark and of its flux tube is certainly well defined (not tachyonic) and this was the first evidence of flux tube vibrations found in lattice field theory. Moreover it was shown [45] that the width of the ground state flux tube diverges when R → ∞ with a logarithmic law, where the squared width w 2 is the mean squared radius of the flux tube, computed in its centre. This enhancement of the flux tube transverse radius as R → ∞ is frequently called widening. The widening as been recently extended Volume β T /Tc a √ σ [78] # config.  with two-loop calculations [51]. The flux tube widening has been verified numerically for compact U(1) QED D = 3 + 1 lattices [52], for non-abelian SU(2) D = 2 + 1 lattices [53][54][55][56][57][58][59][60][61][62][63][64][65][66][67][68][69] and, more recently, for the more physical four-dimensional SU (3) case [70,71]. Recently, it has been shown the flux tubes exhibit characteristics of both superconductor and string models, with both a penetration length λ and the quantum widening of w [44]. Moreover, at finite T , close to T c but still in the confining regime of T < T c , it has been predicted in Ref. [72] that widening becomes linear with the inter-charge distance R. This occurs because the flux tube width w also depends on the extent τ of the compactified time distance of the lattice, Since τ is small at finite T , we expect the dominant term to be the linear term in R. This result has been verified for the Ising model [72] and for compact U (1) [52]. Recently, widening for SU (3) has also been studied for Baryons [70,73,74] We intend to test the linear broadening of Eq. (3) as well for SU (3). However, at T > T c , it has been recently claimed by Refs. [75,76] that a flux tube continues to exist. This apparently contradicts previous results on the static quarkantiquark QQ potentials which indicate that linear confinement disappears for T > T c [77]. Thus, we also intend to clarify how the colour fields and possible flux tubes behave above the phase transition.
Here we extend our previous study [44] to finite temperature T . We aim to measure in detail with lattice QCD techniques the profile of SU(3) pure gauge flux tubes in dimensions D = 3+1. We study the colour fields distributions inside the flux tubes formed by Polyakov loops in the static QQ, QQ and QA [79] systems at finite T , both below and above the phase transition temperature T c . We address how the flux tube evolves with the distance between quarks and when the temperature increases beyond the phase transition. Moreover, we compare our results with the T = 0 flux tubes with the static QQ system computed with the Wilson loop [44].
In Section II, we describe the lattice formulation at finite T . We briefly review the Polyakov loops of the different colour charges systems, detail how to compute the colour field and Lagrangian distributions, and discuss the techniques we utilize to increase the signal over noise ratio. In Section III, we show our results for the different squared field densities, both in the charge axis and in the mediator plane, and qualitatively discuss them. In Section IV, we compute and analyse the widening of the QQ flux tube profile in the inter-charge mediator plane when the separation of the charges increase, in particular we analyse the systematic errors of the width and combine them with the statistical errors. Finally, we present our conclusion in Section V.

II. SU (3) LATTICE QCD FRAMEWORK
We aim to measure the SU (3) colour flux tube produced by static charges, in a lattice QCD framework. We utilize the quenched QCD configurations detailed in Table I.
Our two charges are separated by a lattice vector (fourdimensional) R with spacial components only. We choose our frame with the charge axis in the z direction and the median point of the charges as the origin. The mediator plane is the x, y plane. Moreover we have the Euclidian time axis in the fourth dimension t. This is illustrated in Fig. 1.
The relevant observables of the flux-tube system can be extracted from the correlation of the plaquette µν and charge operators O. The plaquette measures the fields and is computed with four gauge links U , (4) where r is the four-dimensional position of the plaquette, see Fig. 1.
We aim to compare the field density inside the flux tube to the field density in the vacuum. For the vacuum we utilize a reference point r ref sufficiently far from the centre of the flux tube. We measure the following correlator [80], Our operators O are combinations of fundamental representation Polyakov loops L, where A stands for a static charge in the adjoint representation of SU (3), and is the fundamental Polyakov loop and N t is the number of time slices of the lattice. We also use the periodicity in the time direction for the Polyakov loops to average the plaquette over the time direction, Therefore, using the plaquette orientations (µ, ν) = (2, 3), (3, 1), (1,2), (1,4), (2,4), (3,4), we can relate the six components in Eq. (5) to the components of the chromoelectric and chromomagnetic fields, and also calculate the total action (Lagrangian) density, In order to improve the signal over noise ratio in the QQ and QQ systems, we use the extended multihit technique detailed in Ref. [44], which is an extended version of the multihit technique [81,82]. The technique consists in replacing each temporal link by its thermal average with the first N th order neighbours fixed, whereas the simple multihit would just take the thermal average of a temporal link with the first neighbours. We apply the heat-bath algorithm to all the links inside, averaging the central link,  By using N = 2 we are able to greatly improve the signal, when compared with the error reduction achieved with the simple multihit. Of course, this technique is more computer intensive than simple multihit, while being simpler to implement than multilevel [83]. The only restriction is R > 2N for this technique to be valid. Moreover, just below the phase transition, we need to make sure that we don't have contaminated configurations as already mentioned in [84]. By plotting the histogram of Polyakov loop history for β = 6.055 we are able to identify a second peak, see Fig. 2. Thus we remove all the configurations that lie on the second peak [84]. Therefore, in Table I the value with asterisk corresponds to the configurations after removing these contaminated configurations.

III. RESULTS FOR THE SQUARED FIELD DENSITIES IN THE CHARGE AXIS AND IN THE MEDIATOR PLANE
In this section, we present the results for different β values using a fixed lattice size of 48 3 × 8. Our diferent ensembles are detailed in Table I. We compute the lattice spacing, in units of the string tension at zero temperature, using the parametrization of Ref. [78]. All our computations are fully performed in NVIDIA GPUs using our CUDA language.
Our results consist in the squared field densities in two spacial subspaces, the charge axis (z axis) including the two charges, and the mediator plane (plane x, y) of the two charges. As a first qualitative analysis, since confinement should become weaker for higher temperatures, the flux tubes should become less squeezed. We expect the flux tubes to be less dense at higher temperatures. Indeed in the right panels of Fig. 3, it is clear the intensity of the fields does decrease with the temperature.

B. QQ field densities at T > Tc
In Fig. 4, we show the results for the QQ system at temperatures at T > T c . As in Fig. 3  What is clear now is that the intensity of the flux tube does decrease while the inter-charge distance increases. This is visible in the right pannels of Fig. 4, in a behaviour different from the equivalent sub-figures of Fig. 3. This suggests the flux tubes no longer exist above the deconfinement temperature T c .

C. QQ field densities at T > Tc
In Fig. 5, we show the results for the QQ system at T > T c (below T c the Polyakov loop of non colour-singlet systems vanish). For an easier comparison with the QQ system, the six different sub-figures are ordered exactly as in Fig. 4.
It is remarkable that the field densities of the QQ system in Fig. 5 are apparently identical, modulo statistic errors, to the ones of the QQ system in Fig. 4. This similarity was not anticipated, and it may be relevant for the various QCD models based in Polyakov loops.
To check in more detail the difference, we plot in Fig. 6 the difference between the squared field densities of the QQ and of the QQ system at our three temperatures T > T c . The difference is consistent with zero modulo the statistical error bars, i e the two field densities are essentially identical.
This may possibly be interpreted as an evidence for the uncorrelation of the different Polyakov loops at temperatures T > T c , and thus for the non-existence of a flux tube, which should be intrinsically non-linear. It is well known [85] that, in the confined phase below T c , the Polyakov loop has Z3 symmetry, in the sense it tends to take values close to the three cubic roots of 1, at 1, −1+ ; and its average value vanishes.
When we have a QQ system, the fields correlate in a non-linear flux tube, and the average no longer vanishes for T < T c , whereas the Polyakov loops of a QQ system  vanish. Now, in the deconfined phase above Tc, the Z3 symmetry is broken, the Polyakov loop gets closer to 1, and the Polyakov loop average over all configurations is real. This is for instance utilized in matrix models for the deconfinement phase transition [86][87][88][89][90]. In this sense the Polyakov loop of a Q is on average identical to its complex conjugate (see Eq. (6), the Polyakov loop of aQ. Thus, in case there are no non-linear correlations between the two Polyakov loops present on the lattice, it is plausible the field densities for the QQ and QQ are identical.

D. QA field densities at T > Tc
Moreover, in Fig. 7 we study the effect of including a static adjoint source in the system with an adjoint Polyakov loop, as detailed in Eq. 6.
In Fig. 7 we show the squared field densities E 2 , − B 2 , L for three systems, a single adjoint source A system, the adjoint source-quark QA system and the adjoint source-adjoint source AA system, all for our highest temperature with β = 6.4249, T = 1.690T c . This is just a first study, possibly interesting for the effective models of QCD with Polyakov loops, and we do not perform an analysis as detailed as in the QQ and QQ systems. Nevertheless, as in the case with quark sources QQ and QQ, the plots suggest the total square fields of the QA system are approximately similar to a simple linear sum of the square fields produced by two charges. Again, we find this linear superposition contradicts the existence of flux tubes at T > T c , since flux tubes are clearly nonlinear objects.

IV. ANALYSIS OF TUBE WIDENING, INCLUDING SYSTEMATIC ERRORS
In this section, we analyse the flux tube profiles in the mediator plane, equidistant between the charges. Examples of profiles are shown in Fig. 8, where we compare the QQ and QQ profiles.
At T below T c only the colour singlet QQ system produces finite Polyakov loops. Moreover, at T above T c the profiles seem to be additive, in the sense the QQ profile is almost identical to the QQ profile, as discussed in Section III. Thus, in this Section, we specialize in the profiles of the QQ system only.
Moreover, we combine the squared field densities in the Lagrangian density to get a clearer signal, with smaller statistical errors. We also make use of the axial discrete symmetry to increase the statistics of points with equal distance r = x 2 + y 2 to the axial charge axis z. Our main concern is to compute quantitative results from the profiles for different temperatures T and inter-charge distances R.

A. Ansatz for the field density profile
We first fit the flux tube profile of our squared field density F 2 with the ansatz proposed in Ref. [44], with three physical parameters: the axis central intensity of the flux tube F 0 2 , the penetration length at large  distances from the axis λ and the parameter ν related to the second derivative −2F 0 2 /(λν). We also have the unphysical parameter K which is due to the statistical fluctuations of the fields at the reference point r ref of Eq. (5) and due to its finite distance from the centre of the flux tube. The parameter K accounts for the error from the (arbitrary) choice of the reference point, it is small since it is vanishing for high statistics and the profile decreases at least exponentially with |r − r ref |.
Moreover, with our fit we also compute another quantitative parameter [44], considering the normalized F 2 (r) − K as a profile density, the root mean square width, w = r 2 , of the flux tube profile,

B. Computation of the systematic errors
To compute the width, w, which is the main quantitative result of this paper, we first must choose what components E i 2 and B i 2 we adopt as profile density F 2 . Note all components have a similar behaviour, and thus we can choose their most favourable linear combination. We opt for the Lagrangian density L in Eq. (10), which has a better signal-to-noise ratio.  Moreover there are systematic errors present in our fit of F 2 (r) and in the computation of the width w. We must select the interval [0, r max ] in the distance r where we fit the profile of the flux tube. This leads to a systematic error, as in Ref. [91][92][93]. Because the profile vanishes at least exponentially with distance, the points with larger r correspond to a vanishing profile, with no relevant physical information. Moreover, the statistical error increases with distance. We estimate the most significant intervals are the ones where r max corresponds to a value for the Lagrangian density between L 0 /50 and L 0 /500. For instance, in the case of a Gaussian distribution, this would correspond to a fraction of distribution included between 99.5 % and 99.96 %.
To extract the width from the Lagrangian density, as illustrated in Fig. 9, we proceed as follows. We crudely estimate L 0 from the point in the charge axis, at r = 0. We fit the Lagrangian density with our ansatz in Eq. (12), with r max in the interval between L 0 /50 and L 0 /500 crudely estimated. With our fit we then get a correct estimate of L 0 , L 0 /50 and L 0 /500, and redo our fits with the correct r max . Then, for all possible different intervals with r max ∈ [L 0 /50, L 0 /500], and respective fits of L, we determine the parameters F 0 2 = L 0 , λ, ν and w 2 .
The parameters come with a statistical error from each fit range [0, r max ]. Moreover, combining all fit ranges, we obtain a systematic error [91][92][93]. We obtain a barycentre for the systematic error value, and two different upper and lower error bars considering the maximum differences to the barycentre. Finally, for the total error bars, the systematic upper and lower error are averaged, and combined with the statistical error to provide a total error. In Table II we detail the statistical, systematic and total error bars in the case of β = 5.96.

C. Results for the width and central density
We find the central value parameter L 0 and the width w have small statistical and systematic error bars. However, our data is not precise enough to determine both the parameters λ and ν with small error bars. In some of our cases, these two parameters have large error bars, due to redundancy. To remove this redundancy, we would need more data, in order to reduce the statistical error bars.
We show all our results for the central value parameter     tube density L 0 has a large step downwards when T c is crossed. A gap is clearly visible between the T < T c and the T > T c data in Fig. 10.
Moreover we show all our results, for all distances and temperatures, for the squared width w 2 in Fig. 11. While at temperatures below T c the width is clearly growing with distance, there is apparently no evidence for widening at temperatures above T c since the square widths are apparently constant.
Finally, we show a detailed analysis of the widening of the flux tubes, expected to occur only at temperatures T < T c . We plot the square width w 2 , in three separated plots for our three temperatures T < T c , as a function of the QQ inter-charge distance R, in Figs. 12 and 13.
For the finite 0 < T < T c data computed here, we observe in Fig. 12 for the first time a SU (3) behaviour previously studied for instance in compact U(1) [52]. Indeed our data is consistent with a linear fit, as predicted by Ref. [72].
We also re-analyse in Fig. 13 the data of Ref. [44] for T = 0 with the present technique to compute the systematic errors, and we confirm the logarithmic behaviour of the width.

V. CONCLUSIONS
We compute the square densities of the chromomagnetic and chromoelectric fields produced by different Polyakov loop sources, above and below the phase transi-  tion. We fit the profile of the QQ flux tubes and compute physical parameters, including the flux tube width, with statistical and systematic errors. As the distance increase between the sources, the fields square densities decrease. Below the deconfinement critical temperature, this decrease is moderate and is consistent with the widening of the flux tube as already seen in studies at zero temperature [44], moreover the field intensity clearly decreases when the temperature increases, as expected from the critical curve for the string tension [84].
Above the deconfinement critical temperature, at T > T c , the fields rapidly decrease to zero as the quarks are pulled apart, qualitatively consistent with screened Coulomb-like fields. While the width of the flux tube below the phase transition temperature increases with the separation between the quark-antiquark, above the phase transition we find no evidence for widening. Moreover, the squared field densities are additive, in the sense the fields produced by a quark Q, an antiquarkQ and a colour adjoint source A approximately add up together when these sources coexist. In the same perspective, the QQ and the QQ square fields are essentially similar. This is in contradiction with the squeezing of the fields into a flux tube which should be non-linear. Thus we find evidence for the non-existence of flux tubes at temperature above the deconfinement temperature, T > T c .
As an outlook, it would be interesting to complete the present study with further tests of the additive nature of squared field densities above T c . We also would like to produce results with smaller error bars, in order to be able to measure precisely the penetration length parameter λ at finite T , as we did for vanishing T in Ref. [44]. We also plan to produce the different Polyakov loop -Polyakov loop potentials, relevant for modelling the deconfinement transition [89,90,94]. It would also be interesting to observe the cross-over between a logarithmic widening at small T and a linear widening at larger T < T c as in Eq. (3). It will be necessary to develop a new technique [95] to match the T = 0 Wilson loops [44] with the higher T Polyakov loops computed here.