Double-winding Wilson loops in SU(N) lattice Yang-Mills gauge theory

We study double-winding Wilson loops in $SU(N)$ lattice Yang-Mills gauge theory by using both strong coupling expansions and numerical simulations. First, we examine how the area law falloff of a ``coplanar'' double-winding Wilson loop average depends on the number of color $N$. Indeed, we find that a coplanar double-winding Wilson loop average obeys a novel ``max-of-areas law'' for $N=3$ and the sum-of-areas law for $N\geq 4$, although we reconfirm the difference-of-areas law for $N=2$. Second, we examine a ``shifted'' double-winding Wilson loop, where the two constituent loops are displaced from one another in a transverse direction. We evaluate its average by changing the distance of a transverse direction and we find that the long distance behavior does not depend on the number of color $N$, while the short distance behavior depends strongly on $N$.


I. INTRODUCTION
What is the true mechanism for quark confinement is not yet confirmed and still under the debate, although more than 50 years have passed since quark model was proposed by Gell-Mann [1] in the beginning of 1960s. In the 1970s, however, the dual superconductor picture was already proposed by Nambu, 't Hooft and Mandelstam [2] as a mechanism for quark confinement. In fact, validity of the dual superconductor picture was confirmed for U (1) pure gauge theory [3], Georgi-Glashow model [4] and N = 2 supersymmetric Yang-Mills theory [5], although it is not yet confirmed for the ordinary nonsupersymmetric Yang-Mills theory [6] and quantum chromodynamics (QCD). Therefore, the dual superconductor picture is now regarded as one of the most promising scenarios for quark confinement, although this does not deny the existence of the other mechanics for quark confinement. See e.g., [7][8][9] for reviews.
In order to establish the dual superconductor scenario, the most difficult issue to be resolved first of all is to guarantee the existence of magnetic monopoles in the pure non-Abelian Yang-Mills gauge theory, which is different from the 't Hooft-Polyakov magnetic monopole [10] in the gauge-scalar model. This issue was circumvented by using the method called the Abelian projection proposed by 't Hooft [11]. The Abelian projection is a gauge fixing which explicitly breaks the original gauge group into its maximal torus subgroup where color symmetry is also broken. By the Abelian projection, magnetic monopoles of the Abelian type [12,13] are indeed realized, but the resulting theory is distinct from the original gauge theory with the non-Abelian gauge group. To avoid the gauge artifact, we must find a procedure which enables one * Electronic address: skato@oyama-ct.ac.jp † Electronic address: akihiro.shibata@kek.jp ‡ Electronic address: kondok@faculty.chiba-u.jp to define magnetic monopoles in a gauge-invariant way. This issue was solved recently for the Yang-Mills theory with the gauge group SU (N ) and any semi-simple compact gauge group [14], by using the non-Abelian Stokes theorem for the Wilson loop operator and the new reformulation of the Yang-Mills theory based on the new field variables obtained by change of variables through the gauge covariant field decomposition of the Cho-Duan-Ge-Faddeev-Niemi-Shabanov [15][16][17][18][19][20][21][22]. See [9] for a recent review. However, these achievements do not necessarily means that the dual superconductivity is the unique scenario for understanding quark confinement. Recently, Greensite and Höllwieser [23] introduced a "double-winding" Wilson loop operator in lattice gauge theory [24] to examine possible mechanisms for quark confinement. The doublewinding Wilson loop operator W (C = C 1 ×C 2 ) is a pathordered product of (gauge) link variables U ∈ SU (N ) along a closed contour C which is composed of two loops C 1 and C 2 , See Fig.1. A more general "shifted" double-winding loop is introduced in such a way that the two loops C 1 and C 2 lie in planes parallel to the x − t plane, but are displaced from one another in the transverse direction, e.g., z by distance R, and are connected by lines running parallel to the z-axis. In the non-shifted case R = 0, the two loops C 1 and C 2 lie in the same plane, which we call coplanar. We denote by S 1 and S 2 the minimal areas bounded by loops C 1 and C 2 , respectively. Note that the double-winding Wilson loop operator is defined in a gauge invariant manner, irrespective of shifted R = 0 or coplanar R = 0.
In [23], they investigated the area (S 1 and S 2 ) dependence of the expectation value W (C = C 1 × C 2 ) of a double-winding Wilson loop operator W (C = C 1 × C 2 ) for the SU (2) gauge group. Consequently, it has been shown in a numerical way that both the original SU (2) arXiv:2008.03684v1 [hep-lat] 9 Aug 2020 lattice gauge theory and center vortex model obey the difference-of-areas (S 1 − S 2 ) law, while the Abelianprojected model obeys the sum-of-areas (S 1 + S 2 ) law. In the coplanar case R = 0, a double-winding loop has been set up as given in Fig.2. In order to discriminate difference-of-areas and sum-of-areas laws, it is efficient to measure the L 1 -dependence of a coplanar double-winding Wilson loop average W (C = C 1 × C 2 ) , with the other lengths L, L 2 , and δL being fixed. For simplicity, we set δL = 0. Then S 1 (= L × L 2 ) and S 2 (= L 1 × L 2 ) are the minimal areas of rectangular loops C 1 and C 2 , respectively. We assume S 1 ≥ S 2 for definiteness hereafter. If W (C 1 × C 2 ) obeys the difference-of-areas law: then ln W (C 1 × C 2 ) must linearly increase in L 1 as L 1 increases. On the other hand, if W (C 1 × C 2 ) obeys the sum-of-areas law: then ln W (C 1 × C 2 ) must linearly decrease in L 1 as L 1 increases. The numerical evidences were obtained as given in Fig.3 which summarizes their results for L 1 dependence of ln W (C 1 × C 2 ) with the other lengths being fixed, e.g., L = 10, L 2 = 1, δL = 0, based on numerical simulations performed on a lattice of size 20 4 at β = 2.4. These results certainly show both the original SU (2) gauge field and center vortex lead to the difference-ofareas law, while Abelian-projected configurations lead to the sum-of-areas law.
From a physical point of view, a double-winding Wil- son loop can be interpreted as a probe for studying interactions between two pairs of a particle and an antiparticle. Then differences among three cases are understood as follows. In the Abelian model, a particle and an antiparticle in a pair are respectively connected by the electric flux with the length of L and L 1 , as indicated in the top panel of Fig.4. The total energy of flux tubes shifted by R > 0 becomes σ (L+L 1 ), where σ is a string tension, if the flux-flux interactions are neglected. This argument will give a reason why the Abelian model gives the sumof-areas law. Moreover, they argue that even in the limit R → 0 the sum-of-areas law remains unchanged in the Abelian model, because electric flux tubes tend to repel each other and they can not coincide in the type II dual superconductor.
For the SU (2) gauge theory, they argue that the "W bosons" play the crucial role, since they are off-diagonal components of the SU (2) gauge field which are not included in the Abelian model. W bosons have charged components W −− and W ++ with respect to the Abelian U (1) group. They explain that charged off-diagonal components W −− and W ++ of the SU (2) gauge field neutralize respectively positive and negative static charges. Consequently, flux tubes exist only for connecting two positive charges and two negative static charges, which leads to difference-of-areas law. See the bottom panel of Fig.4.
In the vortex picture, if a vortex pierces the minimal area of a loop, it will multiply the holonomy around the loop by a factor −1. Therefore, if a vortex pierces two loops C 1 and C 2 simultaneously, it gives a trivial effect. The non-trivial result is obtained only if a vortex pieces the non-overlapping region S 1 − S 2 . This leads to difference-of-areas law.
Quite recently, Matsudo and Kondo [25] have investigated a double-winding, a triple-winding, and general multiple-winding Wilson loops in the continuum SU (N ) Yang-Mills theory. They have found that a coplanar double-winding SU (3) Wilson loop average follows a novel area law which is neither difference-of-areas law nor sum-of-areas law, and that sum-of-areas law is allowed for SU (N ) (N ≥ 4), if the string tension is assumed to obey the Casimir scaling for quarks in the higher representations.
In this way, the study of double-winding Wilson loops itself is interesting because it can be used to test the confinement mechanism in QCD. Moreover, it is worth considering the interactions between two color flux tubes. In this paper, we investigate both "coplanar" and "shifted" double-winding Wilson loops in SU (N ) lattice Yang-Mills gauge theory by using both strong coupling expansion and numerical simulations.
In this paper, we show that the "coplanar" doublewinding Wilson loop average has the N dependent area law falloff: "max-of-areas law" for N = 3 and sum-ofareas law for N ≥ 4, which add a new result to the known difference-of-areas law for an N = 2 "coplanar" double-winding Wilson loop average. Moreover, we investigate the behavior of a "shifted" double-winding Wilson loop average as a function of the distance in a transverse direction and find that the long distance behavior does not depend on the number of color N , while the short distance behavior depends on N .
This article is organized as follows. In section II, we examine how the area law falloff of a "coplanar" doublewinding Wilson loop average depends on the number of color N . In section III, we examine a "shifted" doublewinding Wilson loop, where the two constituent loops are displaced from one another in a transverse direction, especially evaluate its average by changing the distance of a transverse direction. The final section IV is devoted to conclusion and discussion. We also discuss the validity of the Abelian operator studied in [23]. Recently, there are numerical evidences that the dual superconductor for SU (2) and SU (3) lattice Yang-Mills theory is type I [26], although they explain sum-of-areas law on the basis of type II superconductor. We should study the interaction between two flux tubes in the limit R → 0, in case of type I superconductor.

II. A "COPLANAR" DOUBLE-WINDING WILSON LOOP
First of all, we consider the coplanar case R = 0 of a double-winding Wilson loop in the SU (N ) lattice Yang-Mills gauge theory, as indicated in Fig.2. For simplicity, we set δL = 0. Let S 1 (= L × L 2 ) and S 2 (= L 1 × L 2 ) be the minimal areas of rectangular loops C 1 and C 2 , respectively. We assume S 1 ≥ S 2 for definiteness hereafter.

A. strong coupling expansion
Let S g be a plaquette action for the SU (N ) lattice Yang-Mills theory: where the link field U n,µ satisfies U n+μ,−µ = U † n,µ . This action reproduces the ordinary Yang-Mills action − d D x µ<ν tr(F 2 µν ) up to constant in the naive continuum limit (lattice spacing → 0). The diagrammatic expressions of a plaquette variable U n,µν and the plaquette action are given in Fig.5.
Note that the standard Wilson action S W is defined by see e.g., [29]. The difference of the constant term in the action is physically insignificant and we drop it in the strong coupling analysis. By comparing S g and S W , we can find We define a partition function Z by where dU n,µ is the invariant integration measure of SU (N ). Then the expectation value W (C) of an operator W (C) is defined by In order to evaluate the expectation value in eq.(8), we perform the strong coupling expansion: For the large bare coupling constant g, we can expand the weight e Sg into the power-series of 1/g 2 , and perform the group integration over each link variable U n,µ according to the measure dU n,µ . In Appendix A, we summarize the formulas needed for the strong coupling expansion and for the SU (N ) group integration.

SU (2)
First, we study the case of SU (2) gauge group. For a coplanar double-winding Wilson loop, there is a single link variable U for a link ∈ C 1 − C 2 and there is a double link variable U U for a link ∈ C 2 , as shown in the top diagram of Fig.6.
We list some of explicit SU (2) group integration formula as For a single link variable U (resp. U † ) for ∈ C 1 −C 2 , we need at least one additional link variable with an opposite direction U † (resp. U ) to obtain non-vanishing result after integration in eq.(8) according to the integration formulas (10c) for the SU (2) group integrations. Such link variables are supplied from the expansion eq.(9) of e Sg . Since the number of plaquettes which are brought down from e Sg must be equal to the power of 1/g 2 in the expansion eq.(9), the leading contribution to W (C 1 × C 2 ) comes from a set of plaquettes tiling the minimal area S 1 −S 2 with the least number of plaquettes. See the top diagram of Fig.6. For double link variables U U for ∈ C 2 , on the other hand, we do not need additional link variables coming from the expansion of e Sg to obtain the non-vanishing result due to the integration (10d), giving the g-independent contribution.
For the SU (2) gauge group, therefore, the leading contribution to W (C 1 × C 2 ) in the strong coupling expansion comes from the term in which a set of plaquettes tiles the surface with the area S 1 − S 2 , as shown in the top diagram of Fig.6. Therefore, group integrations give the result where σ = log(2g 2 ). This result was first obtained by Greensite and Höllwieser in [23]. We reconfirmed the difference-of-areas law of coplanar double-winding Wilson loops for SU (2). The bottom diagram of Fig.6 shows one of higher-order contributions in the strong coupling expansion for SU (2). This diagram gives non-vanishing contribution due to the integration formula (10f).

SU (N ), (N ≥ 3)
Next, we study the case of SU (N ) (N ≥ 3) gauge groups. We list some of explicit SU (N ) (N ≥ 3) group integration formula as Notice that the SU (N ) case is different from the SU (2) case. For a double link variable U U for a link ∈ C 2 , we need additional N − 2 link variables (U ) N −2 with the same direction to be brought down from the expansion of e Sg in eq.(8) to obtain the non-vanishing result after the integration according to the integration formulas (12e) for the SU (N ) group integrations. See the top diagram of Fig.7. For a single link variable U (resp. U † ) for a link ∈ C 1 − C 2 , on the other hand, we need at least one additional link variable with the opposite direction U † (resp. U ) to obtain non-vanishing result after integration in eq.(8) according to the integration formulas (12c) for the SU (2) group integrations. Therefore, the contribution from the top diagram of Fig.7 is given by where the coefficient p N is calculated by collecting the numerical factors coming from link integrations and the power-series expansions of e Sg . We have another contribution from the bottom diagram of Fig.7. For a double link variable U U with the same direction for a link ∈ C 2 , we have additional 2 link variables (U † )(U † ) with the opposite directions to be brought down from the expansion of e Sg in eq.(8) to obtain the non-vanishing result after the integration according to the integration formulas (12f) for the SU (N ) group integrations. For a single link variable U (resp. U † ) for a link ∈ C 1 − C 2 , on the other hand, we need at least one additional link variable with an opposite direction U † (resp. U ) to obtain non-vanishing result after integration in eq.(8) according to the integration formulas (12c) for the SU (N ) group integrations. Therefore, the contribution from the bottom diagram of Fig.7 is FIG. 7: A set of plaquettes tiling the areas S1 and S2 which gives the leading contribution to a coplanar double-winding Here S1(= L × L2) and S2(= L1 × L2) are respectively the minimal areas bounded by rectangular loops C1 and C2 with S1 ≥ S2. given by where the coefficient q N is calculated in the similar way to p N . For the SU (N ) (N ≥ 3), the leading contribution in the strong coupling expansion may come from one of the two diagrams shown in Fig.7. Since the number of plaquettes brought down from e Sg is equal to the power of 1/g 2 , these two contributions can be written as where coefficients p N , q N are determined by expansion coefficients of the power series expansion of e Sg and SU (N ) group integrations for link variables. Which contribution becomes dominant is naively determined by comparing the power index of 1 g 2 N , which depends on the number of color N .
For N ≥ 4, we find that the second term in eq.(15) gives the dominant contribution in the strong coupling expansion for W (C 1 × C 2 ) , since the inequality holds, Thus we conclude that the sum-of-areas law of a coplanar doublewinding Wilson loop is allowed for N ≥ 4. This result is consistent with the result obtained by Matsudo and Kondo in [25].
From the top panel of Fig.7, we can easily find that the coefficient p N should be calculated for each number of color N , because type of diagrams are different with the number of color N . On the other hands, we can obtain general formula for the coefficient q N , since the diagram of the bottom panel of Fig.7 is common to all numbers of color N . The result is for S 2 ≥ 1 in lattice units. See Appendix B for the detail.
In the following, we show the results for SU (2), SU (3) and SU (4) in more detail.
SU (2) For the number of color N = 2, eq. (15) reduces to where The factor 2 in front of p 2 and q 2 arises from the nonoriented nature of the plaquettes for SU (2), which is to be compared with (11).
where From this result, we find that the first term in eq.(20) gives the dominant contribution to W (C 1 × C 2 ) for sufficiently large areas S 1 and S 2 , which is neither difference-of-areas law nor sum-of-areas law for the arealaw falloff of the coplanar double-winding Wilson loop average. We call this area-law falloff "max-of-areas law" (or max(S 1 , S 2 ) law). This result is also consistent with the result obtained by Matsudo and Kondo in [25].
SU (4) For the number of color N = 4, eq. (15) reduces to where In this case, both terms in eq.(23) behave as sum-of-areas law.

3.
L1 dependence of the W (C1 × C2) From the above discussions, we can understand the L 1 dependence of the coplanar double-winding Wilson loop average W (C 1 ×C 2 ) in SU (N ) lattice Yang-Mills gauge theory for fixed L, L 2 , and gauge coupling g.
For SU (2) gauge group, we plot eq.(17) in Fig.8, which shows the difference-of-areas law behavior of a coplanar double-winding Wilson loop for N = 2. On the other hand, we plot eq.(20) in Fig.9. For SU (3) gauge group, as the coplanar double-winding Wilson loop average follows the max-of-areas law, it is expected that there are no L 1 -dependence of W (C 1 ×C 2 ) for efficiently large areas S 1 and S 2 . In fact, we can see that the plots flatten at L 1 ∼ 4 (resp. L 1 ∼ 1) in top (resp. bottom) panel in Fig.9.

B. Numerical simulation
We examine the L 1 -dependence of W (C 1 × C 2 ) that we discussed above.
SU (2): We generate the configurations of SU (2) link variables {U n,µ }, using the (pseudo-)heat-bath method for the standard Wilson action. The numerical simulations are performed on the 24 4 lattice at β(= 2N/g 2 ) = 2.5. We thermalize 3000 sweeps, and in particular, we have used 100 configurations for calculating the expectation value of coplanar double-winding Wilson loops W (C 1 × C 2 ) . confirm the difference-of-areas law for SU (2). Note that we can also confirm W (C 1 × C 2 ) −1/2 for S 1 = S 2 from Fig.8.   Fig.9. For example, we can see that the plots flatten at L 1 ∼ 4 for L 2 = 8, which means that there are no L 1 -dependence of W (C 1 × C 2 ) . Thus, we numerically confirm the max-of-areas law for SU (3).

III. A "SHIFTED" DOUBLE-WINDING WILSON LOOPS
Finally, we consider the shifted case R = 0 of a doublewinding Wilson loop in the SU (N ) lattice Yang-Mills gauge theory, as indicated in Fig.12. Contours C 1 and C 2 lie in planes parallel to the x-t plane, but are displaced from one another in the z direction by distance R. Just like the previous section, for simplicity, let C 1 (C 2 ) be a rectangular loop of length L, L 2 (L 1 , L 2 ), and S 1 (≡ L × L 2 ), S 2 (≡ L 1 × L 2 ) be the minimal areas of contour C 1 , C 2 respectively.

A. strong coupling expansion
First, we study the shifted double-winding Wilson loop based on the strong coupling expansion.
One of the diagrams which gives a leading contribution in the strong coupling expansion is given by a set of plaquettes tiling the two minimal surfaces S 1 and S 2 , as shown in Fig.13. The results of a group integration for the links U 's on both surfaces become N (1/g 2 N ) S1+S2  for N ≥ 3, and 2N (1/g 2 N ) S1+S2 for N = 2, respectively. The difference of factor 2 in front of N for N = 2 arises from the non-oriented nature of the plaquettes to conclude the N = 2 result: Another type of diagram which also gives a leading contribution in the strong coupling expansion is given by a set of plaquettes tiling the minimal surface S 1 − S 2 and the four sides with the area 2R(L 1 + L 2 ) of a cuboid with a height R, whose bottom is a rectangular of size To summarize the above discussion, the expectation value of the shifted double-winding loop W (C 1 × C 2 ) R =0 from diagrams as shown in Fig.13 and Fig.14 becomes for N = 2, Note that the R → 0 limit of eq.(28) does not agree with the coplanar result eq.(17), although the sum of the second and third terms in eq.(28) from the diagram of Fig.14 reproduce the coplanar result eq.(17) in the limit R → 0. This is because the first term in eq.(28) coming from the diagram of Fig.13 does not have in the limit R → 0 the counterpart of the strong coupling expansion in the coplanar case and hence contributes only to the shifted case with R = 0. For SU (2) gauge group, especially, we perform the detailed study on the R-dependence of a shifted double- In what follows, we rewrite L 2 into T , Let us imagine T direction be time t-axis, L and L 1 direction be spatial x-axis, and R direction be also space z-axis as seen in top side in Fig.15. As is explained in [23], the shifted double-winding Wilson loop at a fixed time can be interpreted as a tetra-quark system consisting of two static quarks and two static antiquarks. The pairs of quark-antiquarks are connected by a pair of color flux tubes, as seen in the bottom side in Fig.15. We study how interactions between the two color flux tubes change, when the distance R is varied. We find that the second term in eq.(28) dominates for R < R C := L1 1+L1/T , and the first term in eq.(28) dominates for R > R C , because the comparison of the two exponents of these terms for S 1 = LT and S 2 = L 1 T reads where we have neglected the third (higher order) term in eq.(28) for the naive estimate of R C . This means that the left diagram of Fig.16 dominates for R < R C , and the right diagram of Fig.16 dominates for R > R C . Therefore, the dominant diagram switches from left to right at a certain value R C of R as R increases, just like the minimal surface spanned by a soap film. In Fig.17, we plot the R-dependence eq.(28) of a shifted double-winding Wilson loop average W (C 1 × C 2 ) for fixed L, L 1 , and L 2 in the SU (2) lattice gauge theory. The second and third terms in eq.(28) have Rdependence, but the first term in eq.(28) does not depend on R. Therefore, the plot gets flattened for R ≥ R C ∼ 1, which is consistent with Fig.16. This behavior does not depend on the number of color N . In fact, SU (3) and SU (4) cases are given as follows.
SU (4) : In Fig.18, we also plot the R-dependence eq.(31) of a shifted double-winding Wilson loop average W (C 1 ×C 2 ) for fixed L, L 1 , and L 2 in the SU (3) lattice gauge theory.

B. Numerical simulation
Next, we examine the R-dependence of W (C 1 × C 2 ) based on numerical simulations on a lattice.
SU (2): In order to calculate the shifted doublewinding Wilson loop average, we use the same gauge field configurations as those used in calculating the coplanar double-winding Wilson loop. However, we have used APE smearing method (N = 5, α = 0.1) as a noise reduction technique. Fig.19 gives the plots obtained for the W (C 1 ×C 2 ) for various values of R where we have fixed L = 5, T (= L 2 ) = 2, L 1 = 3. We see that the behavior of data in Fig.19 is consistent with the analytical result given in Fig.17. 1 ∼ 6. We see that the data in Fig.20 also consistent with the analytical result given in Fig.18 for sufficiently large areas S 1 and S 2 .

IV. CONCLUSION AND DISCUSSION
In this paper, we have studied the double-winding Wilson loops in SU (N ) lattice Yang-Mills gauge theory by using both strong coupling expansion and numerical simulation.
First of all, we have examined how the area law falloff of a "coplanar" double-winding Wilson loop average depends on the number of color N , by changing the size of minimal area S 2 of loop C 2 . We have reconfirmed the difference-of-areas law for N = 2, and have found new results that "max-of-areas law" for N = 3 and sum-ofareas law for N ≥ 4.
Moreover, we have considered a "shifted" doublewinding Wilson loop, where two contours are displaced from one another in a transverse direction. We have evaluated its average by changing the distance of a transverse direction, and have found that their long distance behavior doesn't depend on the number of color N , but the short distance behavior depends on N .
It should be remarked that this "shifted" doublewinding Wilson loop may contain an information about interactions between two color flux tubes. For this purpose, we need to accumulate more data on the fine lattices with more larger size.
Originally, one of reasons why Greensite and Höllwieser considered the double-winding Wilson loops seems to be that they want to evaluate monopole confinement mechanism in lattice SU (2) gauge theory. They have considered an operator which simply replaces SU (2) link variable U n,µ with the Abelian variable u n,µ as an "Abelian" double-winding Wilson loop, and have shown that the expectation value of such a naive operator obeys the sum-of-areas law. But, it is known that such naive operator should work only for a single-winding Wilson loop in the fundamental representation. Recently, Matsudo and his collaborators [28] have given the explicit expression for the Abelian operator which reproduces the full Wilson loop average in higher representations, which is suggested by the gauge-covariant field decomposition and the non-Abelian Stokes theorem (NAST) for the Wilson loop operator. Similarly, we hope that a correct form of the Abelian operator for a double-winding Wilson loop can be found in the similar way. When we change the line integral to the surface integral, our considerations of the diagrams which give the leading contribution to the strong coupling expansion seems to be useful to construct the NAST for a double-winding Wilson loop. These results will be discussed in a forthcoming paper. In order to perform the strong coupling expansion in the lattice gauge theory, we must calculate the following integrations for the polynomials of group matrix elements over each links: where U ij (i, j = 1, 2, · · · , N ) denotes a matrix element of a matrix U ∈ SU (N ) belonging to the SU (N ) group with the property U −1 = U † , and dU is an invariant measure (Haar measure) on the compact group which is left-invariant and right-invariant We can normalize the measure such that By using properties of the invariant measure, Creutz has shown that eq.(A.1) can be evaluated by the following formula [29,30]: where J is a source variable and is an arbitrary N × N matrix, |J| = det(J), ∂ ji ≡ ∂/∂J ji , and cof(∂) is a cofactor of ∂, respectively. We list some of explicit results from the above formula as The last eq.(A.11) consist for N > 2. For N = 2, Following relation can be shown by using property of invariant measure, From this relation, we also obtain, The following more practical formulae are useful to calculate the expectation value of double-winding Wilson loop by using strong coupling expansion. Let X, Y, A, B be elements of SU (N ) group. From eq.(A.9), we find From eq.(A.8), we find From eq.(A.11), we find for N > 2, In this section, we show explicitly how eq.(16) is obtained. From eq.(8) and eq.(9), a contribution to a coplanar double-winding Wilson loop average W (C 1 × C 2 ) from the bottom panel of Fig.7 is expressed as B.1: Diagrammatic representation of the integration ruleW2 eq.(B.4) for the product of two double-plaquettes with the same clockwise orientation: Integration is performed over the link variables U on the link which is common to two double-plaquettes with the same clockwise orientation. By decomposing the path-ordered product of the link variables along the loop, the plaquette variables for the single plaquette p1 and p2 to the left and right of U is respectively represented by tr(U † p 1 ) := tr(U † X) and tr(U † p 2 ) := tr(Y U ). Here X and Y represent the products of the link variables along staple-shaped paths with the same orientations.
where U † pj and U † p k denote respectively plaquette variables on (S 1 − S 2 ) and S 2 areas. Here note that U † p represents the plaquette variable for the plaquette p with the clockwise orientation.
First, integration with respect to the link variables {U } on the (S 1 − S 2 ) area can be performed with the same technique of the strong coupling expansion as that for the fundamental Wilson loop to obtain where we have defined Next, we perform the integration in eq.(B.3) over the link variables {U } inside of the S 2 area, which excludes the links on the loop C 2 = ∂S 2 (the boundary of S 2 ). As shown in Fig.B.1, performing the integration with respect to the link variables U on the link which is common to two double-plaquettes with the same clockwise orientation using eq.(A.19), we obtaiñ .
Here  From the above consideration, definingW n by the result of connecting n adjacent double-plaquettes one after another by integrating over the link variables inside the S 2 area, we can conclude thatW n is written as This statement is proved by the mathematical induction. Indeed, by applying the same procedures as those given in eq.(B.4) and eq.(B.7) to eq.(B.8), we find the relationship Therefore, we have obtained the recurrence relation which holds for the coefficients α n and β n for n ≥ 1: Solving this recurrence relation with the initial condition eq.(B.5), we obtain the explicit form for the coefficients α n and β n : Because the expansion coefficient is applied to n double-plaquettes.
Appendix C: Explicit calculation of the coefficient p3 In this section, we show explicitly how eq.(22) is obtained. From eq.(8) and eq.(9), a contribution to a coplanar double-winding Wilson loop average W (C 1 × C 2 ) from the top panel of Fig.7 is expressed as dU W (C 1 × C 2 ) · pj ∈(S1−S2) 1 g 2 tr(U † pj ) · p k ∈S2 1 g 2 tr(U p k ) , (C.1) where U † pj and U p k stand respectively for plaquette variables on the (S 1 − S 2 ) and S 2 areas. Here note that U † p and U p respectively represent the plaquette variables for the plaquette p with clockwise and counterclockwise orientations. In this section, we focus on the N = 3 case.
First, the integration with respect to the link variables {U } on the (S 1 − S 2 ) area can be performed with the same technique of the strong coupling expansion as that for the fundamental Wilson loop to obtain where we have defined Next, we perform the integration in eq.(C.3) over link variables {U } inside of the S 2 area, which excludes the links on the loop C 2 as the boundary of the S 2 area. As shown in Fig.C.1, performing the integration over the link variable U using eq.(A.18) for two plaquettes that have a common link U , we obtain dU tr(XU ) · tr(U † Y ) = 1 N tr(XY ). (C.4) From this observation, we conclude that one factor of 1/N appears if two plaquettes are connected after common links are integrated. When S 2 plaquettes are connected one after another by using eq.(C.4), a factor of (1/N ) S2−1 is FIG. C.1: Diagrammatic representation of the integration rule eq.(C.4) for the product of two plaquettes with the same counterclockwise orientation: Integration is performed over the link variable U on the link which is common to two plaquettes with the same counterclockwise orientation. The plaquette variables for the plaquette p1 and p2 to the left and right of U is respectively represented by tr(Up 1 ) := tr(XU ) and tr(Up 2 ) := tr(U † Y ). Here X and Y represent the products of the link variables along staple-shaped paths with the same orientations.
applied, and after that only the path ordered product of the link variables on the loop C 2 as the boundary of S 2 is left unintegrated. Therefore, eq.(C.2) becomes where the integral is only for the link variable on the loop C 2 . As shown in Fig.B.3, by using the decomposition W (C 2 ) := tr(AX) and W (C 2 × C 2 ) := tr(AXAX), and by repeatedly using eq.(A.9), we obtain where we have used the cyclicity of the trace in the second equality. Note that this result is meaningful only when N = 3, because we have used eq.(A.9) in the above calculation, eq.(A.10) holds for M = 0 (mod3). For N = 3, thus, we obtain W (C 1 × C 2 ) p3 = −3 1 3g 2 S1 , (C.7) which indeed yields p 3 = −3.