Three-quark potentials in an $SU(3)$ effective Polyakov loop model

Three-quark potentials are studied in great details in the three-dimensional $SU(3)$ pure gauge theory at finite temperature, for the cases of static sources in the fundamental and adjoint representations. For this purpose, the corresponding Polyakov loop model in its simplest version is adopted. The potentials in question, as well as the conventional quark--anti-quark potentials, are calculated numerically both in the confinement and deconfinement phases. Results are compared to available analytical predictions at strong coupling and in the limit of large number of colors $N$. The three-quark potential is tested against the expected $\Delta$ and $Y$ laws and the $3q$ string tension entering these laws is compared to the conventional $q\bar{q}$ string tension. As a byproduct of this investigation, essential features of the critical behaviour across the deconfinement transition are elucidated.


Introduction
The interest in studying the interquark potential for a three-quark system is not a recent issue at all. It has instead a long history due to its importance in the spectroscopy of baryons. The first studies date back to the mid '80s [1,2] and after more than a decade a new turn of research has started around the year 2000 which continues till now [3][4][5][6][7][8][9][10][11][12][13][14][15][16][17]. New results are somewhat contradictory, which could be reasonably explained by the difficulty of accurate measurements of the three-quark potential. But from these discussions spanning many years, two main Ansatze emerged to describe the three-quark potential: the ∆ law and the Y law. Denoting by r 1 , r 2 , r 3 the sides of the triangle having the quarks at its vertices, the ∆ law is defined by which describes a potential linearly rising with half the perimeter of the triangle. The Y law describes the three-quark potential as linearly rising with the minimal total length of the flux lines connecting the three quarks, where Y is the sum of the distances of the three quarks from the Fermat-Torricelli point F, which is the point such that this sum is the least possible. When all inner angles of the triangle are less than 2π/3, one has where A is the area of the triangle; if one of the angles is larger than 2π/3, we have instead Y = min (r 1 + r 2 , r 1 + r 3 , r 2 + r 3 ) ≡ Λ , which gives rise to the Λ law, Some earlier [4,8,9] and the most recent studies [13,16] in the SU (3) pure gauge theory seem to support the Y -ansatz, while other simulations [3,[5][6][7] prefer the ∆-ansatz, at least for not too large triangles. An even more complicated picture emerged after simulations of the simpler, Z(3) Potts model in two-dimensions, which is believed to capture the most essential features of the gauge model [7,10]. Namely, it was conjectured that there might be a smooth crossover between the ∆ law and the Y law when the size of triangles grows (see, however, [11] where this scenario has been criticized). Also, the paper [10] proposes a new ansatz in which both the Y law and the Λ law are present.
In this paper we are going to study an SU (3) spin model which is an effective model for Polyakov loops and can be derived from the original gauge theory in the strong coupling region. For simplicity, we consider, following [10], only its twodimensional version. Our primary goal is to get some analytical predictions for the three-point correlation function of the Polyakov loops and compare them with numerical simulations. For that we use the SU (3) spins both in the fundamental and adjoint representations. The main tool of our analytical investigation is the large-N expansion. Within this expansion we demonstrate that the fundamental three-point correlator is described by a sum of the Y and Λ laws. The ∆ contribution is not present. In turn, the connected part of the adjoint three-point correlation follows the ∆ law in the confinement phase. In addition, we study the critical region of the model and confirm that it belongs to the universality class of the two-dimensional Z(3) spin model. This paper is organized as follows. In the next section we introduce our notations, define SU (N ) Polyakov loop model and its dual. Certain analytical predictions for two-and three-point correlation functions are obtained in the strong coupling expansion and in the large-N limit. Moreover, we check the restoration of the rotational symmetry for the 3-quark system. Section 3 outlines some details of our numerical simulations. Here we compare numerical data with the strong coupling expansion and study the critical behaviour of the model using the finitesize scaling analysis. Section 4 presents the results of Monte Carlo simulations for the fundamental and adjoint two-and three-point correlations in the confinement region. Results for the same quantities above critical temperature are described in Section 5. In Section 6 we summarize our results.

The model and theoretical expectations 2.1 Partition and correlation functions
We work on a 2d Euclidean lattice Λ = L 2 , with sites x = (x 1 , x 2 ), x n ∈ [0, L − 1], and denote by e n the unit vector in the n-th direction. Periodic boundary conditions (BC) are imposed in all directions. Let W (x) ∈ SU (N ), and TrW be the character of SU (N ) in the fundamental representation. Consider the following partition function on Λ, which describes the interaction of non-Abelian spins: Re TrW (x)TrW * (x + e n ) .
The trace of an SU (N ) matrix can be parameterized with the help of N angles, e.g. by taking W = diag(e iω 1 , · · · , e iω N ), subject to the constraint k e iω k = 1. In this parameterization the action has the form Re TrW (x)TrW * (x + e n ) = N i,j=1 The invariant measure for SU (N ) is given by where and δ(x) is the periodic delta-function. Due to this constraint, the SU (N ) model is invariant only under the global discrete shift ω k (x) → ω k (x) + 2πn N for all k and x. This is just the global Z(N ) symmetry.
The main subjects of this work are the two-and three-point correlation functions for the SU (3) model. In the fundamental representation these correlations are given by while in the adjoint representation the correlations are written as Γ ad 2 (β, R) = χ ad (W (0))χ ad (W (R)) , where we use the relation χ ad (W ) = TrWTrW * − 1.
The partition function (6) can be regarded as the simplest effective model for the Polyakov loops which can be derived in the strong coupling region of 3d lattice gauge theory (LGT) at finite temperature (see, e.g., [18] and references therein). Namely, the integration over the spatial gauge links on the anisotropic (d + 1)dimensional lattice with two couplings β s and β t ≡ β in the limit β s = 0 and for β sufficiently small leads to the d-dimensional spin model (6). It describes the deconfinement phase transition of the pure gauge theory, which is of second order for SU (3) if d = 2. It is widely assumed that the phase transition is in the universality class of the two-dimensional Z(3) (Potts) model. The inverse correlation length (mass gap) is the string tension of the gauge theory. The correlation length diverges when approaching the critical point with the critical index ν = 5/6. Another important critical index η, which is a characteristic of the massless phase, equals 4/15 exactly at the critical point. Thus, on the basis of the universality arguments [19] we expect the same values for these indices also in the effective SU (3) Polyakov loop model. More on the critical behaviour of three-dimensional SU (N ) LGTs can be found in Refs. [20,21].
The model (6) cannot be solved exactly at any finite N and D > 1. Therefore, to get some analytical predictions for the behaviour of the three-point correlation functions we consider the large-N limit of the model. This limit can in turn be solved exactly by using the dual representation which we are going to describe shortly.

Dual representation
In some applications the dual formulation of the Polyakov loop model (6) can be useful. Such formulation for the SU (3) model has been derived in [22]. Here we use the dual representation obtained by some of us in [23,24]. This form of dual theory is valid for all N and can be used both for numerical simulations and for the study of the large-N limit of the theory. For the 2d theory the partition function (6) on the dual lattice takes the form where Q N (s,s) results from the invariant integration over the SU (N ) measure, Here d(λ) is the dimension of the irreducible representation λ of the permutation group S s , s −s = kN and λ is enumerated by the partition λ = (λ 1 , λ 2 , · · · , λ l(λ) ) of s, i.e.
is the length of the partition and λ 1 ≥ λ 2 · · · λ l(λ) > 0. The sum in (15) is taken over all λ's such that l(λ) ≤ N and the convention λ + q N ≡ (λ 1 + q, · · · , λ N + q) has been adopted. For the exact expressions of the different correlation functions we refer the reader to the paper [24].

Large-N solution
Using the dual representation (14), one can construct an exact solution of the model in the large-N limit [25] and even estimate the first non-trivial corrections specific for the SU (N ) group. As an example, we give here the expression for the most general correlation function and for the partition function in the confinement region in the presence of sources The Gaussian part describes the solution in the large-N limit, while the product over x in the second line presents the first correction due to SU (N ). G x,x is the massive two-dimensional Green function for the scalar field.
This solution, together with a similar one in the deconfinement phase, enables one to calculate both fundamental and adjoint two-and three-point correlations in that limit. Different results are obtained in the small and large β regions separated by the deconfinement phase transition. If we take N = 3, then for the confinement phase we get and in the deconfinement phase we get In the equations above the Green function in the thermodynamical limit is given by where R 2 = x 2 1 + x 2 2 and the functional dependence of the mass m on β is different in the confined and deconfined phases. In the confinement phase the mass m(β) coincides with the qq string tension, while in the deconfinement phase this quantity has the meaning of screening mass. M f (β) and M ad (β) define the fundamental and adjoint magnetizations at a given β, correspondingly. They also depend on the considered phase. For example, M f (β) = 0 in the confined phase. α(β) is another R-independent quantity which appears due to Gaussian integration around the large-N solution. All four quantities -m(β), α(β), M f (β), M ad (β) -are known exactly in the large-N expansion.

3q potential
In what follows our strategy relies on the assumption that the large-N formulae (20)- (27) remain qualitatively valid (up to one correction explained below) at finite N , in particular for N = 3. We expect that the most essential difference between the large-N limit and the N = 3 case exhibits itself in the vicinity of the critical point. Indeed, both our solution [25] and the mean-field solution of Ref. [26] reveal the existence of a third order phase transition at large N . Meanwhile, as described above, the SU (3) Polyakov loop model belongs to the universality class of the two-dimensional Z(3) model. It means, in particular that the critical behaviour of two-and three-point correlation functions is described by a different set of the critical indices ν and η. Therefore, we shall use (20)- (27) as fitting functions, where the quantities m(β), α(β), M f (β), M ad (β) are unknown parameters to be found from fits of numerical data. In most cases, we use the asymptotic expansion for the Green function G(β, x) given on the right-hand side of Eq. (28). As we explained above, only the critical indices appearing in these quantities can vary with N . Also, we introduce here another quantity, namely the index η, in order to describe the power dependence of the correlation function, R −η , on the distance. This could again be important in the vicinity of the critical point. In general, this introduces a correction to the potential of the form and is interpreted as the Coulomb part of the full potential in the two-dimensional theory.
Since the asymptotic behaviour of G(β, x) is known, it follows that we actually know the large-distance behaviour of all two-and three-point functions listed above, but Γ f 3 (β, {x i }). The behaviour of the latter can be analyzed by the saddlepoint method when at least one side of the triangle is large enough. We find two types of the behaviour: 1. All inner angles of the triangle are less than 2π/3. The three-point fundamental correlation function is given by the sum of two terms corresponding to Y and Λ laws where A, B are constants and σ qqq = σ qq . This behaviour resembles the behaviour of the three-point correlation function in the Z(3) spin model [10].
2. One of the angles is larger than 2π/3. In this case the asymptotics is described by the above formula with A = 0. Thus, only the Λ law is present. This again agrees with the Z(3) spin model.
Let us also emphasize that we could not find the ∆ law contribution in our large-N approach. Nevertheless, we attempt to fit numerical data both to Y and ∆ laws in the following. Finally, let us stress that the connected part of the three-point adjoint correlation follows the ∆ law in the confinement phase, as is seen from Eq. (23).

Strong coupling expansion
When β is sufficiently small, one can use the conventional strong coupling expansion to demonstrate the exponential decay of the fundamental two-and threepoint correlation functions. Instead, adjoint correlations stay constant over large distance. To check our codes we have calculated the leading orders of the strongcoupling expansion for the two-point correlator at distance R = 1 and for the three-point correlator in the isosceles-triangle geometry T with base b = 2 and height h = 1. The results read Γ ad 2 (β, 1) = For arbitrary isosceles triangle T with base b and height h one obtains On a cubic lattice Y l = h + b is the minimal sum of the lattice distances from the triangle vertices to an arbitrary lattice point. Then, according to Eq. (2), in the strong coupling region on the finite lattice. Thus, strictly speaking the strong coupling expansion predicts not an exact Y law, as it is often stated in the literature, but rather a Y l law. In general, Y l > Y and we expect that the rotational symmetry will be restored quickly with β and the triangle sides increasing. This should result in the restoration of the genuine Y law. To demonstrate that such a restoration really takes place, we have studied the three-point correlation function for triangles with 2 ≤ b ≤ 10 and 6 ≤ h ≤ 14 at β = 0.41. The fact that the rotational symmetry is already restored at this value of β is shown on Fig. 1, where we compare numerical data with the fits of the form e −σqqqD /D η for D = Y l and D = Y . Clearly, D = Y describes data better than D = Y l .

Details of numerical simulations
To calculate the correlation functions from numerical simulations we used two different approaches. The first is the simulation of the model in terms of the eigenvalues ω i (x) of the SU (3) spins, described in more detail in [21]. In this approach (denoted as standard in the following), an updating sweep consisted in the combination of a local Metropolis update of each lattice site, followed by two updates by the Wolff algorithm, consisting in Z(3) reflections of the clusters. An alternative approach is the simulation of the dual model (14), using the heatbath update for the link variables k, q and the dual site variables r. In this case, we can measure only observables invariant under the global Z(3) symmetry. In both approaches we measured two-and three-point correlation functions in the fundamental and adjoint representations, taking for the two point correlation function pairs of points separated by R in one of the two lattice directions, with R = 2, 4, . . . , L/2. For the three-point correlation functions two geometries were studied: isosceles triangles with base b and height h, and right-angled triangles with the catheti (of lengths a 1 and a 2 ) along the two lattice directions. In both cases, b and h, and a 1 and a 2 , took independent values in the set {2, 4, . . . , L/2}.
In addition to the two-and three-point correlations (10)-(13), the magnetizations and their susceptibilities were measured: The x-dependent values are averaged over all sites of the lattice. For each simulation we performed 10 4 thermalization updates, and then made measurements every ten whole lattice updates (sweeps), collecting a statistics of 10 5 -10 6 . To estimate statistical errors a jackknife analysis was performed at different blocking over bins with size varying from 500 to 10000.
A comparison of the two simulation methods showed that the dual code performs better at small values of β, while giving much larger fluctuations than the standard one when β is close to its critical value. What is more important -at larger values of β for the fundamental correlation function the fluctuations rapidly increase with the distance between the points. Due to this, most of the results presented here have been obtained in the standard approach, and the dual code was used only for cross-check purposes.

Comparison with strong coupling
To test our algorithms we performed a set of simulations at small values of β (β < 0.15), and compared the obtained values of Γ f 2 , Γ f 3 , Γ ad 2 and Γ ad 3 with the corresponding determinations in the strong coupling expansion (Eqs. (31)-(34)). The results of the comparison are shown in Fig. 2 for the correlations in the fundamental representation, and in Fig. 3 for the correlations in the adjoint representation. It can be seen that the two-point correlation, both in the fundamental and adjoint representation, is in good agreement with the strong coupling expansion. For the three-point correlation, due to its small absolute value, statistical errors in the standard simulation are too large to make any statement about agreement with the strong coupling prediction. The results for the adjoint correlation from the dual code are compatible with the strong coupling expansion up to β = 0.1. Since the results of the two simulation codes agree in the region around β c , where most of our simulations were carried out, we are confident in the reliability of our measurements.

Critical behaviour
A clear indication of the two-phase structure of the model is provided by the scatter plots of the complex magnetization at different values of β, shown in Fig. 4.
To precisely locate the β c at which the phase transition occurs, we have studied the magnetization susceptibility for different lattice sizes L, extracting the value of β pc (L) from a fit of the peak of the susceptibility with a Lorentzian function.   The value for the critical index ν is in agreement with the critical index ν = 5/6 of the two-dimensional three-state Potts model, to whose universality class our model is believed to belong. A direct extraction of the critical exponent ν, performed in the subsection 4.3, gives a compatible result, which is sensitive to the choice of the region of β values where critical scaling is supposed to hold. As a second check of the order of the phase transition and of the universality class, we studied the dependence of the peak value of the magnetic susceptibility for different lattice sizes using the scaling law (see the right panel of Fig. 5) We found B = 0.0282 (27), γ/ν = 1.737 (17), χ 2 r = 0.30 . The potential parameter σ qq is extracted from the measurements of the observable Γ f 2 . Following Eq. (20) and the explanation in subsection 2.4, we expect One can extract σ qq from the following ratio: We have compared our Monte Carlo data for σ eff qq (R) with the formula σ qq + η/2 ln [R/(R − 2)]. The interval of β values that we considered for the extraction of σ qq was [0.41, 0.42], since for β < 0.41 the two-point correlation drops too fast to be of any use. The values of σ eff qq obtained in the selected range of β do not show any significant difference when moving from a 64 × 64 to a 128 × 128 lattice, thus making unnecessary to perform simulation on even larger lattices.
As an alternative method for extracting σ qq , we measured the wall-wall correlation function, Tr W (x, y 1 ) Tr W † (x + R, y 2 ) , which is known to obey the exponential decay law, with no power corrections, Introducing, similarly to (44), we found that σ eff ww (R) exhibits a long plateau at each considered β value; we took as plateau value σ ww the value of σ eff ww (R) at the smallest value of R after which all values of σ eff ww (R) agree within statistical uncertainties. Results for the 128 × 128 lattice are summarized in Table 1: we can see that results of σ qq obtained from the fitting of σ eff qq (R) according to (44) are in good agreement with the plateau values of σ ww . We ascribe the difference between the σ qq values obtained for different choices of R min to possible systematic effects arising from the difficulty to extract the parameters of the exponential decay corrected by a power law and to our treating of the correlation function errors as independent values. These systematic errors seem also to result in a value of σ ww being in most cases slightly higher than the estimates of σ qq .

Extraction of
First we studied the dependence of Γ f 3 on the geometry, considering ∆ and Y laws. In Fig. 6 we see that if we consider σ to be proportional to σ ww we get a reasonable collapse for all β values except the largest one (β = 0.424), which might be too close to the critical point for our lattice size, L = 128. Still this does not allow us to discriminate between the two laws.
We turned therefore to fits with the two laws of the three-point correlation function for small (σ ww Y < 2) triangles, excluding those having an angle larger than 2π/3. The results of these fits for β = 0.423 are shown in Fig. 7. In this case the fit with the ∆ law gives A = 2.319 (14), σ = 0.0631(4), η = 0.503(4), χ 2 r =  To extract σ qqq we followed a procedure similar to the one used for σ qq . We supposed, according to (21) and (30), the decay law In this subsection R is a distance parameter depending on the geometry, that could be ∆ or Y , respectively for the ∆ and Y laws, given in Eqs. (1) and (2). In this proposed fitting function we excluded the second term, corresponding to the Λ law, since, for the size of triangles considered in the fitting procedure, we could hardly distinguish it from the first one, corresponding to the Y law 1 We calculated the following ratio for the pairs of triangles (b + 2, h + 2) and (b, h) for the isosceles geometry and (a 1 + 2, a 2 + 2) and (a 1 , a 2 ) for the right-angled geometry, where R 1 and R 2 are the distance parameters of the two triangles. The ratio is equal to where, for sufficiently large distances we can assume R 1 ∼ R 2 ∼ R and get It turned out that for part of the triangle pairs σ eff qqq and its jackknife error estimate could not be extracted reliably, due to one of the correlation functions being too close to zero. We removed from the study all the triangle pairs in which for at least one of the triangles, at least one of the jackknife samples gave negative correlation. The actual number of the triangle pairs for which the extraction of σ eff qqq was possible, strongly depends on the value of β; for example, for the isosceles geometry we had 274 triangle pairs for β = 0.41, 558 for β = 0.42 and 783 for β = 0.423.
After extracting σ eff qqq , we plotted it directly against the half-perimeter ∆ and against the sum Y of the distances of the triangle vertices from the Fermat-Torricelli point. We overlapped these plots with the plots of σ eff qq and σ eff ww versus R (Figs. 8-10). We see that on the plots for the ∆ law the values of σ eff qqq fail to collapse into a single line, while the collapse is much better for the Y law for all β values we studied. The residual spread of the points can be at least partially explained by different triangle pairs having different R 1 − R 2 values, which are not distinguished on these plots. Another observation that supports the Y law is that the collapse line for the σ eff qqq closely matches the line of σ eff qq , which suggests that not only the sigma values entering the two-and the three-point correlation are the same if we consider the Y law, but also that the parameters η are similar in these cases.
It is worth mentioning that in our study the results of the extraction of σ eff qqq are compared for triangles that have strongly different geometries: there are triangles that have similar Y distance, but some of them can have a small base and a large height, while others can have small heights and large bases. In particular, we did not exclude triangles with angles larger than 2π/3 from the study of σ eff qqq , despite the fact that for them the Fermat-Torricelli point coincides with one of the vertices, thus leading to a different dependence of Y on h and b. The fact that even these "extreme" triangles obey the Y law is explicitly demonstrated in Fig. 11 (instead, the most outlying data points turn out to be the ones with smallest triangle base), where the caption σ qqq (Y → Λ) implies that for such triangles the Fermat-Torricelli point coincides with one of the vertices turning the Y law into the Λ law. As can be seen from Figs. 8-10 these differences in geometry give negligible corrections to the values of σ eff qqq up to β = 0.423, providing us with an additional point in support of the Y law.
based on the use of hyperspherical three-body variables [17].

Extraction of critical index ν from the scaling of the two-point string tension
We used the values of the string tension obtained from the wall-wall correlation function close to the critical point to extract the critical index ν.
The values of σ ww , as well as the result of the fit in the region 0.422 < β < β c with the scaling function A = 12.5(2.8), β c = 0.424255(34), ν = 0.806(37), χ 2 r ≈ 3.69 , are shown in the Fig. 12. The value of the critical index ν obtained in this way is compatible with both the critical index ν = 5/6 of the three-state Potts model and with our previous estimate in (41). We note, however, that the scaling region in this case is extremely narrow, and the value of ν is quite sensitive to the inclusion of the points outside this region.

Adjoint correlations in the confinement phase
We have performed measurements of the two-and three-point correlation functions in the adjoint representation, defined in Eqs. (12) and (13), at some values of β below the critical one.
Following formulae (22) and (23) and replacing in them the massive Green function G(β, R) with its asymptotic behaviour, we got the following models: Note that Γ ad 3 ({x i }) exhibits the ∆-law decay after subtracting terms proportional to (powers of) the magnetization.
The results of the fitting of the adjoint correlations to the models in Eqs. (53) and (54) are given in Table 2 (see also Figs. 13 and 14).
Unusually low χ 2 r values in the fits arise due to treating the measurements of the correlations at different distances as independent despite being obtained from the same set of measurements. This fact makes the error estimates of the fit parameters unreliable. Also, the estimation of the η value is inaccurate, since it describes short-range corrections to the exponential decay, on which only a few points in the fitting range have impact. This fact is especially visible from the covariance between η and σ which is very close to −1.
Despite that, the fact that the parameters σ and M show some degree of stability when going from the description of the two-point correlation to the threepoint one, and also the compatibility of the results for σ with the σ qq and σ ww values given in Table 1, support the validity of the suggested descriptions for the adjoint correlations.
We would like to stress that the values of M ad extracted from the fits do agree with direct measurements of this quantity as defined in Eq. (38).

Correlation functions in deconfinement phase
Using the same approach adopted in the subsection 4.4 for the description of adjoint correlations in the confinement phase, we get the following models for the where m can be interpreted as the chromoelectric screening mass. We used these formulae to describe the correlation functions measured above β c . The remarks given at the end of subsection 4.4 apply here as well. The results of the fits are gathered in Table 3 for the fundamental correlations (see also Figs. 15 and 16) and in Table 4 for the adjoint ones (see also Figs. 17 and 18). We see that the values of M and m extracted from the fits for two-point and three-point correlations are compatible. Moreover, the m values extracted from the fits of the fundamental and adjoint correlations are also compatible between themselves (the same should not apply to the M values, since they represent the average magnetization in two different representations). This supports the validity of the formulae (55)-(58) for the description of the correlation functions in the deconfinement phase.
Also in this case, we found that the values of M ad extracted from the fits do agree with direct measurements of this quantity as defined in Eq. (38).
In the deconfinement phase the value of m becomes the inverse correlation length for the connected part of the correlation, as can be seen from the Taylor  Table 2: Parameters extracted from the fits of the Γ ad 2 and Γ ad 3 at some given β < β c . For each value of β the first line contains the result of the fit of Γ ad 2 to (53), and the next lines contain the result of the fit of the values of Γ ad 3 obtained for the isosceles triangles with fixed base b to (54). expansion of the outer exponent in Eq. (55). When we approach the critical point from above the value of m should vanish as The ratio A dec /A, where A is the amplitude for the scaling of σ in the deconfinement phase given by Eq. (52), equals 2.657 for 2d Z(3) universality class [27].

Summary
In this study we performed an extensive analysis of a two-dimensional effective SU (3) Polyakov loop model. Differently from other approaches of the same kind, in our effective model the basic degrees of freedom are traces of SU (3) matrices and not Z(3) spins. The partition function is therefore integrated with a groupinvariant measure. The motivation for this choice is that it can help to catch some important features of the (2 + 1)-dimensional SU (3) lattice gauge theory at finite temperature that escape approaches based on the center degrees of freedom. For example, one cannot define adjoint correlation functions in Z(3) models. Our main goal was to examine the behaviour of the three-quark potential below the critical point and to distinguish between possible scenarios for the three-point correlation function decay: the ∆ law and the Y law. Considering the triangle whose vertices are the positions of the three sources, the ∆ potential depends on the perimeter of this triangle, while the Y potential depends on the sum Y of the distances of its vertices from the Fermat-Torricelli point. Other studies accomplished in the paper include investigation of the two-and three-point correlation functions in the adjoint representation both below and above critical point and calculation of critical indices in the vicinity of the deconfinement phase transition.
Our main findings can be summarized as follows: • Similarly to the pure gauge SU (3) LGT and Z(3) spin model, the leading contribution to the three-point fundamental correlation function in the strong coupling region is described by the Y l law, as explained in section 2.5. Exact Y law is restored as soon as the rotational symmetry is also restored.
• From the study of the large-N limit of the model we also obtained the general form for the two-and three-point correlations in fundamental and adjoint representations both above and below the critical point. The analytical results suggest that the fundamental three-point correlation behaves as in Eq. (30), i.e. it is described by a combination of Y and Λ laws if all angles of the triangle are less than 2π/3 or by only the Λ law if any of the angles is larger than 2π/3. We have not found analytical support in favour of ∆ law (one possibility is that this law is suppressed in the large-N limit). The comparison with the results of numerical simulations shows that these forms can indeed be used to describe the behaviour of the corresponding correlations.
• The critical behaviour across the deconfinement transition supports the universality conjecture that this model is in the universality class of the twodimensional Z(3) spin model. In particular, we have determined the critical indices ν and η from finite size scaling (the index ν has been also evaluated   • The fact that the assumption of Y law gives much better collapse of the effective string tension to a single curve, when different locations of the three sources are taken with the same value of Y , indicates that in this effective model the Y law is preferred. We have found the agreement of the effective string tensions for the two-and three-point correlation functions. Moreover, all string tensions appearing in the two-quark potential and in the threequark potential with Y or ∆ law agree up to uncertainties. This result is also supported by the analytic study of the model in the large-N limit. Since our results are supportive of the Y law even for small triangles, we do not observe a smooth crossover from ∆ to Y law as conjectured in [6,10]. For triangles with one of the angles larger than 2π/3 the three-point fundamental correlation function follows the Λ law.
• In the deconfinement phase the screening masses for fundamental two-and three-point correlations also coincide up to numerical errors.
• Adjoint correlations share a similar pattern in the deconfinement phase, namely the corresponding screening masses are consistent with each other