*Article* **A Possible Resolution to Troubles of SU(2) Center Vortex Detection in Smooth Lattice Configurations**

**Rudolf Golubich \* and Manfried Faber**

**\*** Correspondence: rudolf.golubich@gmail.com

**Abstract:** The *center vortex model* of quantum-chromodynamics can explain confinement and chiral symmetry breaking. We present a possible resolution for problems of the vortex detection in smooth configurations and discuss improvements for the detection of center vortices.

**Keywords:** quantum chromodynamics; confinement; center vortex model; vacuum structure; cooling

**PACS:** 11.15.Ha; 12.38.Gc

#### **1. Introduction**

The center vortex model assumes that the relevant excitations of the QCD vacuum are *center vortices*, closed color magnetic flux lines evolving in time. It can explain Confinement [1] and chiral symmetry breaking [2–4]. In four-dimensional space–time, the flux lines form closed surfaces in dual space, see Figure 1. In the low-temperature phase, they percolate space–time in all dimensions.

**Figure 1.** The geometric relation between piercings, the flux line and the vortex surface is schematically shown. Left: A flux line can be traced by following non-trivial plaquettes (depicted in orange with a "−1") after transformation to maximal center gauge and projection to the center degrees of freedom. Middle: Each non-trivial plaquette belongs to four elementary cubes, where the flux enters and has to leave through another plaquette. The depicted grey rectangles correspond to the same plaquette. For each cube, the three involved coordinates are indicated. Right: Due to the evolution in time, the flux line (depicted as orange line) forms a closed two-dimensional surface in four-dimensional spacetime.

Within lattice simulations, the center vortices are detected in *maximal center gauge* after projection to the center degrees of freedom. The procedure is described in more detail in Section 2. As long as the detected vortices reproduce the relevant physics, we speak of a valid *vortex finding property*. During the analysis of the color structure of vortices in smooth configurations [5] one is confronted with a loss of the vortex-finding property. Problems in

**Citation:** Golubich, R.; Faber, M. A Possible Resolution to Troubles of SU(2) Center Vortex Detection in Smooth Lattice Configurations. *Universe* **2021**, *7*, 122. https://doi.org/10.3390/ universe7050122

Academic Editor: Dmitri Antonov

Received: 24 March 2021 Accepted: 25 April 2021 Published: 29 April 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

Atominstitut, Technische Universität Wien, 1040 Wien, Austria; faber@kph.tuwien.ac.at

detecting center vortices due to ambiguities in the gauge-fixing procedure were already found by Kovacs and Tomboulis [6]. They also point out that the thickness of vortices is of importance for the extraction of properties related to confinement. We found that this thickness can cause troubles in the vortex detection, resulting in a loss of the string tension. In search for improvements in the vortex detection, the cause of this loss is analyzed and a possible resolution discussed. We model the influence of cooling on the vortex thickness and the corresponding loss of the vortex density. An upper limit for the lattice spacing and a lower limit for the lattice size is presented. These limits are derived from measurements of the vortex density and estimates of the cross-section of flux tubes.

#### **2. Materials and Methods**

Our lattice simulations of the SU(2) Wilson action cover an interval of inverse coupling *β* ∈ [2.1, 3.6] in steps of 0.05. We start with low *β* values to identify discretization effects and to detect the onset of finite size effects. To check how far the compatibility of our model reaches, we expand the calculations to relatively large values of *β*. The lattice spacing *a* corresponding to the respective values of *β* is determined by assuming a physical string tension of (440 MeV)<sup>2</sup> via a cubic interpolation of the literature values given in Table 1. This is complemented by an extrapolation according to the asymptotic renormalization group equation for *β* > 2.576

$$a(\beta) = \Lambda^{-1} e^{-\frac{\beta}{8\beta\_0}} \text{ with } \beta\_0 = \frac{11}{24\pi^2} \text{ and } \Lambda = 0.015(2) \,\text{fm}^{-1},\tag{1}$$

with Λ obtained by fitting this equation to the values of *a* for *β* ≥ 2.6 in Table 1.

**Table 1.** The indicted dependence of the lattice spacing *a* in fm and the string tension *σ* in lattice units on the inverse coupling *β* is taken from references [7–11].


The analysis is performed on lattices of size 8<sup>4</sup> and 104 with 0, 1, 2, 3, 5 and 10 Pisa-Cooling [12] steps with a cooling parameter of 0.05. We have chosen these small lattice sizes because, in bigger lattices, the finite-size effects are expected at higher values of *β* and, as we will show, the detection of center vortices becomes increasingly difficult with rising values of *β*.

A central part of our analysis consists of identifying non-trivial center regions, regions whose is perimeter evaluated as close to non-trivial center elements, using the algorithms presented in references [13–15]. In the gauge-fixing procedure, we look for gauge matrices Ω that maximize the functional

$$\mathcal{R}^2 = \sum\_{\mathbf{x}} \sum\_{\mu} |\operatorname{Tr}[\acute{l}\acute{l}\_{\mu}(\mathbf{x})]|^2 \quad \text{with} \quad \acute{l}\acute{l}\_{\mu}(\mathbf{x}) = \Omega(\mathbf{x} + \boldsymbol{\epsilon}\_{\mu}) \amalg \chi(\mathbf{x}) \Omega^\dagger(\mathbf{x}). \tag{2}$$

The non-trivial center regions are used to guide this procedure to prevent the problems found by Bornyakov et al. [16]. The detection of such non-trivial center regions is based on enlarging regions until they get as close as possible to non-trivial center elements. This quite calculation-intensive procedure is depicted in Figure 2.

As not all resulting regions are evaluated sufficiently near to a non-trivial center element, we take only those into account which are sufficiently near to non-trivial center elements.

During the gauge-fixing, only such gauge matrices Ω are allowed that preserve the sign of the non-trivial center regions. As this causes the rejection of some gauge matrices, the number of required simulated annealing steps until convergence of the gauge functional might increase.

**Figure 2.** The non-trivial center regions, used for the gauge-fixing procedure, are detected by repeating the depicted procedure until every plaquette either belongs to an identified region or has already been used as the seed to grow a region. The direction of enlargement of the respective regions is marked by an arrow. Plaquettes that belong to a region are colored; plaquettes already used as seed are shaded. In the final determination of the non-trivial center regions enclosing a thick vortex, no collision handling is performed.

After gauge fixing and projection, plaquettes are identified that evaluate non-trivial center elements. These are dubbed *P-plaquettes* and considered to be pierced by a P-vortex.

If the number of P-plaquettes is smaller than the number of non-trivial center regions used to guide the gauge-fixing procedure, this is a clear indication of a failing vortex detection. For each value of *β*, the proportion of configurations where this is the case is determined. This allows to quantify the loss of the vortex-finding property besides quantifying it directly via the string tension of the center-projected configurations.

The further analysis is performed in the full SU(2) configurations. For each P-plaquette, a non-trivial center region that encloses the P-plaquette is identified. This center region is considered to be pierced by the thick vortex detected by the P-vortex. Figure 3 depicts the relation between P-vortices, thick vortices and the non-trivial center regions.

**Figure 3.** Two-dimensional slices through a four-dimensional lattice are depicted. The vortex detection starts as a best-fit procedure of P-Vortices to thick vortices, indicated by the first arrow. Then, starting from the detected P-plaquettes, non-trivial center regions are identified to reconstruct the thick vortex. These non-trivial center regions are, in general, not rectangular.

The cross-section of the flux building thethick vortex, *A*vort, is measured by counting the plaquettes that build up the non-trivial center regions enclosing the corresponding thick vortex. In each configuration, we determine minimal, average and maximal cross-sections.

The string tension *σ* is determined via Creutz ratios *χ* calculated in the center projected configurations

$$\sigma \approx \chi(\mathbb{R}, T) = -\ln \frac{\langle \mathcal{W}(\mathbb{R} + 1, T + 1) \rangle}{\langle \mathcal{W}(\mathbb{R}, T + 1) \rangle} \frac{\langle \mathcal{W}(\mathbb{R}, T) \rangle}{\langle \mathcal{W}(\mathbb{R} + 1, T) \rangle} \,\tag{3}$$

with *R* × *T* Wilson loops *W*(*R*, *T*). As the Coloumb-part of the potential is strongly suppressed after projecting to the center degrees of freedom, the linear part corresponding to a non-vanishing string tension is already reproduced with small loop sizes, as we saw in references [14,15,17]. Symmetric Creutz ratios are used and the average of *χ*(1, 1) and *χ*(2, 2) is taken to determine the string tension. Our study is based on the data generated in reference [5], where we did not save a sufficiently wide range of Wilson loop data.

Assuming independence of vortex piercings, the string tension can also be related to the vortex density *vort*, the number of P-plaquettes per unit volume, via

$$
\sigma \approx -\ln(1 - 2 \times \varrho\_{\text{vort}}).\tag{4}
$$

The requirement of uncorrelated piercing is only fulfilled if the vortex surface is strongly smoothed, otherwise this simple equation overestimates the string tension.

The working hypothesis is that the loss of the vortex-finding property, observed via a loss of the string tension, when cooling is applied, can be related to a thickening of the vortices. We will try to model the loss of the vortex density based on an analysis of the geometric structure of center vortices.

#### **3. Results**

The different measurements are performed for a lot of different values of *β* and several cooling steps. So as not to overload the visualizations only a part of the intermediate results is depicted, showing only specific numbers of cooling steps and restricting to a smaller interval of *β*-values. Those parts of the data that are dominated by finite size effects are identified and excluded from the further analysis.

Starting with the quantification of the vortex-finding property presented in Figure 4, some troubles are brought to light. The proportion of configurations where fewer Pplaquettes have been identified than non-trivial center regions exist, rises rapidly when passing a specific value of *β*. This specific value depends on the lattice size and the number of cooling steps.

When reducing the lattice size or increasing the number of cooling steps, the loss of the vortex-finding property occurs at lowered values of *β*. The proportion depicted seems to saturate at about 30%, except for 10 cooling steps at a lattice of size 84, where it reaches higher values. The fact that some non-trivial center regions have no corresponding P-plaquettes after gauge fixing and projection to the center degrees of freedom hints at a possible explanation for part of the lost string tension. The gauge functional given in Equation (2) is local in the sense that each gauge matrix Ω is solely based on the eight gluonic links connected to the specific lattice point. Farther distances than a single lattice spacing are not directly taken into account. In contrast, the detection of the non-trivial center regions is, in a sense, more physical as it is based solely on gauge-independent quantities, that is, the evaluation of arbitrary big Wilson loops. When detecting P-vortices in smooth configurations and high lattice resolutions, the center flux can be distributed over many link variables. Each of these links can evaluate arbitrarily close to the trivial center element, although a Wilson loop build by the links can evaluate arbitrarily near to the non-trivial center element. In such a scenario, a gauge-fixing procedure, only taking the vicinity of lattice points into account, will likely fail and result in an underestimated string tension.

**Figure 4.** The proportion of configurations is depicted where less non-trivial plaquettes have been identified than non-trivial regions exist. The datapoints are joined to guide the eye. Due to the logarithmic scaling of axes, only non-vanishing values are depicted: all lines start with 0% at lower values of *β*. The interruption of the green line corresponding to the lattice of size 10<sup>4</sup> at 10 cooling steps at *β* = 2.55 results from a vanishing percentage at the respective *β*-value. Observe that the curves rise at different values of *β* for different number of cooling steps and different lattice sizes.

Looking at the Creutz ratios depicted in Figure 5 two possibly intertwined effects can be observed.

**Figure 5.** The string tension *σ* is estimated via an average of the Creutz ratios *χ*(1, 1) and *χ*(2, 2) calculated in center-projected configurations for different numbers of cooling steps and lattice sizes. The datapoints are joined by lines to guide the eye. The literature values correspond to those listed in Table 1; the asymptotic line is given by Equation (1). Observe that, in the low *β*-regime, an underestimation of the string tension correlates to the number of cooling steps. This underestimation is independent of the lattice size. At higher values of *β*, finite size effects set in.

At sufficiently low values of *β*, the string tension is independent of the lattice size, but decreases with an increasing number of cooling steps. Of interest is that, for sufficiently small values of *β*, the deviation from the asymptotic prediction decreases with a rising value of *β*—for example, the 104-lattice starts at 10 cooling steps with an underestimation of the asymptotic string tension of 50% ± 1% at *β* = 2.1, improving to 40.8% ± 0.4% at *β* = 2.25. At higher values of *β*, the independence from the lattice size no longer holds. For different lattice sizes, a sudden decrease in the string tension occurs at different values of *β*. The respective *β*-values are compatible for different numbers of cooling steps. The dependency on the lattice size and the independence on the number of cooling steps

hint at finite size effects, but finite size effects do not give a direct explanation of the reduction in the string tension at lower values of *β*: We do not observe a dependency on the lattice size in the low *β*-regime. Based on the deviations of the string tensions for different lattice sizes, we expect finite size effects to occur at length scales around 1.3 fm, independent of cooling: observing that the lattice of size 8<sup>4</sup> deviates from the 104-lattice at *<sup>β</sup>* <sup>≈</sup> 2.3, corresponding to a lattice spacing of *a* ≈ 0.165, we acquire a physical lattice extend around 1.32 fm for the smaller lattice. The finite size effects on the bigger lattice set in at *β* between approximately 2.35 and 2.4, resulting in a length scale between approximately 1.2 fm and 1.4 fm. This length scales are compatible with the findings of Kovacs and Tomboulis [18]. In Ref. [5] we also found color-homogeneous regions embedded in the vortex surface with roughly the same diameters. Similar distances can also be found between neighbouring piercings of a Wilson loop, extracted from the vortex density, as will be seen in Table 4.

A relation to the thickness *A*vort of center vortices is suspected and points towards possible further analysis. The possibility of a thick vortex expanding due to a spreading of the center flux was already suggested by Kovacs and Tomboulis in [19]. Assuming a circular cross-section of the flux tube, its diameter can be calculated as

$$d\_{\text{flux}} = 2 \times \underbrace{\sqrt{\frac{A\_{\text{vort}}}{\pi}}}\_{r\_{\text{flux}}}.\tag{5}$$

with *A*vort being the area of the flux cross-section. That flux lines are closed requires that within each two-dimensional slice through the lattice at least two vortex piercings can find place. This give a criteria on the lattice extent *L*

$$L \gg \mathbf{2} \* d\_{\text{flux}}.\tag{6}$$

If *A*vort measured by a plaquette count exceeds 19 for a lattice of size 104, or 12 for a lattice of size 84, we can expect finite size effects to step in. These thresholds are of relevance for the average, minimal and maximal flux tube cross-section depicted in Figures 6–8. The mean flux tube cross-section presented in Figure 6 shows that we have to restrict to relative low values of *β* to stay away from finite size effects.

**Figure 6.** The average cross-sections of the flux tubes, measured by counting plaquettes, increases when cooling is applied. It reaches a threshold at which finite size effects are expected to become problematic, shown as a dashed line for the two lattice sizes. Measurements performed on lattices of different size have good compatibility.

Taking a look at the maximal flux tube cross-section depicted in Figure 7, we can expect finite size effects at even lower values of *β*: None of the data with 10 cooling steps can be expected to be free of finite size effects.

**Figure 7.** The maximal cross-sections of the flux tubes hint at finite size effects. Within our *β*-interval, only the lattice of size 10<sup>4</sup> stays below the threshold when cooling is applied. With cooling, the different lattice sizes become more and more incompatible.

The lattice of size 8<sup>4</sup> could be too small even without any cooling applied. With cooling and increasing *β* the different lattice sizes become more and more incompatible. This may be caused by finite size effects and insufficient statistics. Still, the overall behaviour with cooling is qualitatively reproduced and allows the gain of another estimate on the growth rate.

Looking at the minimal tube size depicted in Figure 8 an even more sudden rise in the cross-section can be observed.

We expect that the minimal flux tube cross-sections starts to grow with a certain *β*, where the high action density of non-trivial plaquettes leads to a suppression within the path integral. This causes a dependency on the lattice size due to the reduced statistics. For sufficiently low values of *β* and sufficiently low numbers of cooling steps, the minimal flux tube cross-section is given by exactly one plaquette, independent of *β* and the number of cooling steps. We restrict further analysis to 104-lattices with *<sup>β</sup>* <sup>≤</sup> 2.3 and, at most, five cooling steps. Nevertheless, we depict the full data in all relevant figures to allow for a check of the plausibility of our model by looking at the specific deviations of the data from our prediction.

Assuming an exponential growth in the flux tubes' cross-section with an increase in the number of cooling steps, a model of the form

$$A\_{\rm vort}(N\_{\rm cool}) \;= A\_{\rm vort}(0) \; e^{N\_{\rm cord} \; (\not\!\!g\_{\rm cord} + \not\!\!g\_{\rm cord} \; \text{at} \;)} \tag{7}$$

is fit to the data, with *N*cool being the number of cooling steps and *a* the lattice spacing. The fit-parameter *g*cool corresponds to the exponential growth in the flux tube with cooling. As the tube size is measured by counting plaquettes, we have to account for discretization effects. This is done by adding another fit-parameter *g*discret in the exponent, related to the lattice spacing and the number of cooling steps. The two parameters are not necessarily constant as they can depend on the specific structure of interest. We restrain from carrying along another index: In the following, the values of these two parameters are to be considered only with respect to the specific context. They differ for the average cross-sections and the maximal cross-sections of flux tubes. The fit of this model to the average flux tube sizes is shown in Figure 9 in physical units. The fit is done for small *β* and cooling steps indicated by black points.

**Figure 9.** The measured data of the average flux tube cross-section for various numbers of cooling steps and several *β* are shown by black and orange points. The dashed lines depict the fits according to Equation (7), where only the black datapoints were used. The corresponding fit parameters are given in Table 2. Deviations of the data from the fits can be related to finite size effects.

The fit, dashed lines reproduce the data well until the expected onset of finite size effects for cross-sections, increasing with the number of cooling steps and *β*. This onset is compatible to the estimates in Equation (6) and will be discussed later. At present, we concentrate on the growth in the flux tube cross-section described by the fit parameters given in Table 2. The suspected exponential growth of *A*vort is confirmed by the good quality of the fit for positive *g*cool, even for larger values of *β* and cooling steps.

**Table 2.** The parameters of the model described by Equation (7) and depicted in Figure 9 for average cross-sections are shown.


The negative value of *g*discret reflects the decreasing slope of the dashed lines with increasing *β*, indicating an influence of the lattice resolution: a coarser lattice reduces the growth of *A*vort. The overall behaviour of *A*vort is qualitatively reproduced by the maximal cross-sections, as depicted in Figure 10.

**Figure 10.** The measured data of the maximal flux tube cross-section for various numbers of cooling steps and several *β* are shown by black and orange points. The dashed lines depict the fits according to Equation (7), where only the black datapoints were used. The corresponding fit parameters are given in Table 3. Deviations of the data from the fits can be related to finite size effects.

Only the growth has slowed down, as can be seen in the values given in Table 3.

**Table 3.** The parameters of the model described by Equation (7) and depicted in Figure 10 for maximal cross-sections are shown.


This implies that the growthin *A*vort with increased cooling is limited.

A further influence of cooling is a smoothing of the vortex surface. We will now model this smoothing and show that the vortex flux tubes can be thickened without pushing each other apart. The vortex density vort allows to gain information about the distance of the vortex centers. Here, we have to take into account that some of the P-plaquettes belong to correlated piercings and can be attributed to short-range fluctuations. We define the quantity *A*max as the non-overlapping area around vortex centers.

The vortex density vort is usually calculated by dividing the number P-plaquettes by the total plaquette number. Given enough statistics, it can be determined by counting the number of piercings *N*vort within a sufficiently large Wilson loop of Area *A*loop build by *N*loop plaquettes

$$\rho\_{\text{vort}} = \frac{N\_{\text{vort}}}{N\_{\text{loop}}} = \frac{N\_{\text{vort}}}{A\_{\text{loop}} \ast a^{-2}} = \frac{N\_{\text{vort}}}{(A\_{\text{free}} + N\_{\text{vort}} \ast A\_{\text{max}}) \ast a^{-2}}.\tag{8}$$

In the last identity, we have split the area of the loop into two non-overlapping parts: each piercing is enclosed by circular area given by *A*max and *A*free covers the remaining part of the loop. When cooling is applied, we have to take into account that *A*max grows.

$$\varphi\_{\text{vort}}(N\_{\text{cool}}) = \frac{N\_{\text{vort}}}{(A\_{\text{free}} + N\_{\text{vort}} \* (A\_{\text{max}}(0) + \delta A\_{\text{max}}(N\_{\text{cool}}))) \* a^{-2}}.\tag{9}$$

Using *A*loop = *A*free + *N*vort ∗ *A*max(0) and a model of the form given in Equation (7) for *<sup>A</sup>*max(*N*cool) we attain *<sup>δ</sup>A*max <sup>=</sup> *<sup>A</sup>*max(0)(*eN*cool (*g*cool+*g*discrete *<sup>a</sup>*) <sup>−</sup> <sup>1</sup>). It follows

$$\varrho\_{\text{vort}}(N\_{\text{cool}}) = \frac{\varrho\_{\text{vort}}(0)}{1 + \varrho\_{\text{vort}}(0) \ A\_{\text{max}}(0) a^{-2} \left( e^{N\_{\text{cool}} \left( \frac{\eta}{\text{cool}} + \frac{\eta}{\text{discreto}} \right. \right. \right)}} \tag{10}$$

We fit *g*cool, *g*discrete and *A*max(0) to the measurements of vort. The measured data and the fit are shown in Figure 11.

**Figure 11.** The vortex density is depicted for different values of *β* and different numbers of cooling steps. For the model prediction, shown as dashed lines, only the black datapoints were used. That the datapoints fall below the model prediction at specific numbers of cooling steps for different values of *β* can be explained by finite size effects. The corresponding parameters of the model are given in Table 4.

The respective fit parameters are listed in Table 4.

**Table 4.** The parameters of the model described by Equation (10) for the loss of the vortex density during cooling.


The value of *A*max(0) is larger than the flux tube cross-sections depicted in Figure 9. This, and the fact that the value of *g*cool, for the vortex density is smaller than those of the vortex flux tube cross-sections indicate that the majority of piercings remain separated from

one another even when cooling is applied. Assuming circular geometry, we can calculate the minimal possible distance between vortex centers

$$d\_{\text{center}}(N\_{\text{cool}}) = 2\sqrt{\frac{A\_{\text{max}}(N\_{\text{cool}})}{\pi}}.\tag{11}$$

To determine how many cooling steps are possible, we need to know how much the vortices can grow by cooling without getting into conflict. We estimate the minimal available separation by

$$s\_{\text{flux}}(N\_{\text{cool}}) = \underbrace{2\sqrt{\frac{A\_{\text{max}}(0)}{\pi}}}\_{d\_{\text{other}}(0)} - \underbrace{2\sqrt{\frac{A\_{\text{vort}}(N\_{\text{cool}})}{\pi}}}\_{d\_{\text{flux}}(N\_{\text{cool}})}.\tag{12}$$

We use *d*center(0), the average distance between piercings when no loss of the vortex density occurred, and subtract the average diameter of the flux tubes *d*flux(*N*cool) with cooling applied. If *s*flux(*N*cool) becomes smaller than one lattice spacing, our methods of center vortex detection are likely to fail: we can no longer find two non-overlapping non-trivial center regions enclosing the thick vortex flux tubes. This allows for a limit for the lattice spacing to be derived, *a*, given in Equation (13) together with a limit on *L* based on Equation (6)

$$a < s\_{\text{flux}} \quad \text{and} \quad L > \text{Max}(2\overline{d\_{\text{flux}}}, \text{Max}(d\_{\text{flux}})). \tag{13}$$

The requirement for the lattice extent *L* is based on the fact that two vortex piercings have to fit in every two-dimensional slicing through the lattice. Assuming a vanishing minimal flux tube size, the limit is given either by two times the average diameter *d*flux or one times the maximal diameter Max(*d*flux)—whatever is bigger. The assumption of a vanishing minimal flux tube size is an approximation: on the lattice, the minimal size is given by exactly one plaquette, which is normally negligible in comparison to the lattice extent.

Using what we learned so far, we can evaluate these inequalities and find numerical values for the upper limit of *a* and the lower limit of *L*. These are depicted in Figure 12 and will now be discussed. Discretization effects are neglected by setting *g*discret = 0. Fitting the average flux tube cross-sections for configurations without cooling for 2.1 ≤ *β* ≤ 2.3, see Figure 9, by a polynom up to quadratic order with respect to the lattice spacing *a* gives

$$\overline{A\_{\text{vort}}}(0) \approx 3.367(38) \ a^2 + 0.200(9) \text{ fm } a,\tag{14}$$

compatible with the values we found in [20]. A fit to the maximal cross-sections without cooling for 2.1 ≤ *β* ≤ 2.3, see Figure 10, results in higher fit parameters

$$\text{Max}(A\_{\text{vort}}(0)) \approx 11.3(2) \ a^2 + 0.224(37) \text{ fm } a. \tag{15}$$

Using this fit and Equation (12) with *A*max(0) from Table 4 we obtain an upper limit for the lattice spacing that depends on the number of cooling steps and *g*cool. With this limit, we can determine a lower limit for the required lattice extent. Both limits are shown in Figure 12 for the two different values of *g*cool resulting from average and maximal flux tube sizes from Tables 2 and 3. Let us remember how these limits were derived. Closed flux lines require sufficient room for two piercings within each two dimensional slice through the lattice—a lower limit for the lattice extent arises.

Taking the stronger limits with *g*cool = 0.14, we determine the corresponding limits of *β* for given lattice size and number of cooling steps. In Table 5 some numerical values are shown.

**Figure 12.** Based on the growth in the flux tubes and the reduction in the vortex density in dependency of the number of cooling steps, an upper limit for the lattice spacing (left) and a lower limit for the lattice extent (right) can be derived, as given in Equation (13). The weaker limit depicted in red is based on the slower growth in the maximal sized flux tubes with *g*cool = 0.0999 (see Table 3), the stronger limit, depicted in orange, is based on the faster growth in average-sized flux tubes with *g*cool = 0.14 (see Table 2).


**Table 5.** For different numbers of cooling steps and different lattice extents, the table gives a lower and an upper limit for *β*. "None" indicates that the limits exclude one another.

We now look at the meaning of these limits for the string tension. In Figure 5, we observe that the deviation from the asymptotic prediction decreases with increasing *β* in the low *β*-regime. We believe that this behaviour holds within the *β*-intervals of Table 5. The upper limit of *β* can be extended by increasing the lattice size. It would be interesting to see if this alone suffices to restore full compatibility with the asymptotic string

tension with modest cooling, but the required computational power might exceed our present capabilities.

#### **4. Discussion**

Using non-trivial center regions we analyzed how Pisa-cooling influences the crosssections of thick center vortices. We found an exponential growth that slows down with increasing cross-sections. By geometric arguments, we derived an upper limit for the lattice spacing above which discretization effects trouble the vortex detection and a lower limit for the lattice extent where finite size effects set in. This window gets smaller with cooling and decreasing lattice extent. Cooling results in deviations from the asymptotic behaviour: an underestimation of the string tension occurs. Within the window, increasing *β* leads to better agreement with the asymptotic behaviour. It would be interesting to see whether the string tension calculated on the projected lattice is, in fact, fully restored with sufficiently large *β* or if only a partial restoration occurs.

By improving the method of center vortex detection, it might be possible to soften the aforementioned limits. The method of vortex detection used in this work was based on the direct maximal center gauge guided by non-trivial center regions [13–15]: we identify regions whose perimeter evaluates to the non-trivial center element and preserve their evaluation during gauge-fixing and center projection. This approach comes with three possibilities of improvement.

The growth in the flux tube due to cooling results in the non-trivial center factors within the evaluation of Wilson loops being spread over more and more links. In the original direct maximal center gauge the contribution to the gauge functional at a given site x is determined by its attached links only. By taking farther links into account, the troubles arising from the spread of the center flux may be counteracted.

When two thick vortices are not separated by at least one lattice spacing, the identification of non-trivial center regions enclosing the single piercings might fail. The original method used for the detection of non-trivial center regions is based on enlarging the perimeter of Wilson loops while preventing overlaps of the resulting regions: if overlaps occurred, the region that evaluates to a higher trace is deleted. By allowing overlaps, an improvement might be possible: more non-trivial center regions are kept to guide the further gauge-fixing procedure.

With rising number of cooling steps, more non-trivial center regions than P-plaquettes were found: The direct maximal center gauge failed to preserve some of the non-trivial center regions. This could be counteracted by inserting non-trivial factors before starting the simulated annealing procedure used to maximize the gauge functional. These non-trivial factors should guarantee that each non-trivial center region evaluates to the non-trivial center element when evaluated in the center projected configuration.

**Author Contributions:** Conceptualization, R.G.; methodology, R.G.; software, R.G. and M.F.; writing original draft preparation, R.G.; writing–-review and editing, R.G. and M.F.; supervision, M.F.; project administration, M.F.; All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** The data presented in this study are available on request from the corresponding author.

**Acknowledgments:** We thank the company *Huemer-Group* (www.huemer-group.com, accessed on 3 March 2021) and Dominik Theuerkauf for providing the computational resources speeding up our calculations.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


## *Review* **From Center-Vortex Ensembles to the Confining Flux Tube**

**David R. Junior 1,2,\*,†, Luis E. Oxman 1,† and Gustavo M. Simões 1,†**


**\*** Correspondence: davidjunior@id.uff.br

† These authors contributed equally to this work.

**Abstract:** In this review, we discuss the present status of the description of confining flux tubes in SU(N) pure Yang–Mills theory in terms of ensembles of percolating center vortices. This is based on three main pillars: modeling in the continuum the ensemble components detected in the lattice, the derivation of effective field representations, and contrasting the associated properties with Monte Carlo lattice results. The integration of the present knowledge about these points is essential to get closer to a unified physical picture for confinement. Here, we shall emphasize the last advances, which point to the importance of including the non-oriented center-vortex component and non-Abelian degrees of freedom when modeling the center-vortex ensemble measure. These inputs are responsible for the emergence of topological solitons and the possibility of accommodating the asymptotic scaling properties of the confining string tension.

**Keywords:** confinement; ensembles and effective fields; topological solitons

#### **1. Introduction**

Our knowledge about the elementary particles, as well as three of the four known fundamental interactions, is successfully described by the standard model of particle physics. In particular, the quantitative behavior of the electromagnetic, weak, and strong interactions is encoded in the common language of gauge theories. In the strong sector, an important and intriguing phenomenon regarding the possible asymptotic particle states takes place. When quarks and gluons are created in a collision, they cannot move apart. Instead, they give rise to jets of colorless particles (hadrons) formed by confined quark and gluon degrees of freedom. Although confinement is key for the existence of protons and neutrons, a first-principles understanding of the mechanism underlying this phenomenon is still lacking. At high energies, the detailed scattering properties between quarks and gluons are successfully reproduced by QCD perturbative calculations in the continuum, which are possible thanks to asymptotic freedom. This is in contrast with the status at low-energies, where the validity of quantum chromodynamics (QCD) is well-established from computer simulations of the hadron spectrum which successfully make contact with the observed masses. This review focuses on this type of non-perturbative problem in pure *SU*(*N*) Yang–Mills (YM) theory, which is a challenging open problem in contemporary physics. Here again, Monte Carlo simulations provide a direct way to deal with the large quantum fluctuations and compute averages of observables such as the Wilson loop, which is an order parameter for confinement in pure YM theories. As usual, the lattice calculations, as well as the center-vortex ensembles we shall discuss, consider an Euclidean (3d or 4d) spacetime. Unless explicitly stated, this is the metric that will be used throughout this work. For heavy quark probes in an irreducible representation D, the Wilson loop is given by:

$$\mathcal{W}\_{\rm D}(\mathcal{E}\_{\rm c}) = \frac{1}{\mathcal{\widetilde{\mathcal{Y}}}} \operatorname{tr} \mathbf{D} \Big( P \Big\{ e^{i \int\_{\mathcal{E}\_{\rm g}} d\mathbf{x}\_{\mu} A\_{\mu}(\mathbf{x})} \Big\} \Big) \,, \tag{1}$$

**Citation:** Junior, D.R.; Oxman, L.E.; Simões, G.M. From Center-Vortex Ensembles to the Confining Flux Tube. *Universe* **2021**, *7*, 253. https:// doi.org/10.3390/universe7080253

Academic Editor: Dmitri Antonov

Received: 7 June 2021 Accepted: 15 July 2021 Published: 21 July 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

where D is the dimensionality of D. The closed path C<sup>e</sup> can be thought of as associated to the creation, propagation, and annihilation of a pair of quark/antiquark probes. From a rectangular path with sides *T* and *R*, information about the static interquark potential was obtained from the large *T* behavior WD(Ce) <sup>∼</sup> *<sup>e</sup>*−*T V*D(*R*). An area law, given by the propagation time *T* multiplied by the interquark distance *R*, corresponds to a linear confining potential [1] (for a review, see [2]).

There are many model-independent facts that point to the importance of the center of the group *SU*(*N*) to describe the confining properties of YM theory. In this regard, the first ideas relating the possible phases to the *Z*(*N*) properties of the vacuum were developed in [3]. There, disorder vortex field and string field operators were introduced in (2 + 1)d and (3 + 1)d Minkowski spacetime, respectively. At equal time, they satisfy

$$
\hat{\mathcal{W}}\_{\rm F}(\mathcal{C}\_{\rm e})\,\hat{\mathcal{V}}(\mathbf{x}) = e^{i2\pi \,\mathrm{L}\left(\mathbf{x}, \mathcal{C}\_{\rm e}\right)/N} \,\hat{\mathcal{V}}(\mathbf{x}) \,\hat{\mathcal{W}}\_{\rm F}(\mathcal{C}\_{\rm e})\,,\quad\text{in } (2+1)\mathrm{d},\tag{2}
$$

$$
\hat{\mathcal{W}}\_{\rm F}(\mathcal{C}\_{\rm c}) \, \hat{\mathcal{V}}(\mathcal{C}) = e^{i2\pi \, \mathrm{L}(\mathcal{C}, \mathcal{C}\_{\rm c}) / N} \, \hat{\mathcal{V}}(\mathcal{C}) \, \hat{\mathcal{W}}\_{\rm F}(\mathcal{C}\_{\rm c}) \, \, \text{ in (3+1)} \mathbf{d}, \tag{3}
$$

where the subindex F denotes the fundamental representation, **<sup>x</sup>** <sup>∈</sup> <sup>R</sup><sup>2</sup> (C ∈ <sup>R</sup>3) is a point (curve) in real space where a thin pointlike (looplike) thin center vortex is created in three (four) dimensional spacetime. *L*(**x**, Ce) and *L*(C, Ce) are the corresponding linking numbers. An explicit realization of *<sup>V</sup>*<sup>ˆ</sup> was given by the action *<sup>V</sup>*<sup>ˆ</sup> <sup>|</sup>*<sup>A</sup>* <sup>=</sup> <sup>|</sup>*A<sup>S</sup>* , where |*A* are quantum states with well-defined shape *A*<sup>0</sup> = 0, *Ai* (*i* = 1, 2, 3) at a given time. The field *A<sup>S</sup> <sup>μ</sup>* has the form of a gauge transformation, but performed with a singular phase *S* ∈ *SU*(*N*). To define the operator *<sup>V</sup>*ˆ(**x**) (respectively *<sup>V</sup>*ˆ(C)), *<sup>S</sup>* must change by a center element when going around any spatial closed loop that links **x** (respectively C). Spurious singularities may be eliminated by using the adjoint representation Ad(*S*), which leaves a physical effect only at the point **x**, or closed path C, where Ad(*S*) is multivalued. Arguments in favor of characterizing confinement as a *magnetic Z*(*N*) spontaneous symmetry breaking phase (center-vortex condensate),

$$
\langle \hat{\mathcal{V}}(\mathbf{x}) \rangle \neq 0, \qquad \langle \hat{\mathcal{V}}(\mathcal{C}) \rangle \sim e^{-\mu \text{Perimeter}(\mathcal{C})}, \tag{4}
$$

were also given in that work.

The lattice also provides direct information about the role played by the center of *SU*(*N*) in the confinement/deconfinement phase transition. This is observed in the properties of the Polyakov loops *P***x**(A ), which are given by Equation (1) computed on a straight path located at a spatial coordinate **x** and extending along the Euclidean time-direction. Due to the finite-temperature periodicity conditions, these segments can be thought of as circles. By considering the fundamental representation, *P***x**(A ) was analyzed in the lattice [4]. When changing from higher to lower temperatures, the distribution of the phase factors of *P***x**(A ), for typical Monte Carlo configurations, shows a phase transition. At higher temperatures, for most **x**, the phase factors are close to one of the center elements *<sup>e</sup>i*2*πk*/*N*, *<sup>k</sup>* <sup>=</sup> 0, ... , *<sup>N</sup>* <sup>−</sup> 1. On the other hand, below the transition, they are equally distributed on *Z*(*N*), as a function of the spatial site **x**. As a result, the Monte Carlo calculation gives a transition from a non-vanishing to a vanishing gauge-field average *P***<sup>x</sup>** , which is in fact **x**-independent, where the *electric Z*(*N*) symmetry is not broken. This corresponds to a transition from a deconfined phase at higher *T*, where the quark free energy is finite, to a confined phase below *Tc*, where the free energy diverges.

In the full Monte Carlo simulations, the relevance of *Z*(*N*) is also manifested in general Wilson loops at asymptotic distances. In this regime, the string tension only depends on the *N*-ality *k* of D, which determines how the center *Z*(*N*) of *SU*(*N*) is realized in the given quark representation [5],

$$\mathcal{D}\left(e^{i\frac{2\pi}{N}}I\right) = \left(e^{i\frac{2\pi}{N}}\right)^{k}I\_{\mathcal{G}^{\dagger}}.\tag{5}$$

Regarding the confinement mechanism, lattice calculations aimed at determining the relevant degrees of freedom have been performed for many years. In particular, procedures have been constructed to analyze Monte Carlo *Uμ*(*x*) ∈ *SU*(*N*) link-configurations and extract center projected configurations *Zμ*(*x*) ∈ *Z*(*N*) [6–9] (for recent techniques to improve the detection of center vortices, see [10]). A given plaquette is then said to be pierced by a thin center vortex if the product of these center elements along the corresponding links is non-trivial. Observables may then be evaluated by considering vortex-removed and vortex-only configurations. The confining properties are only well described in the latter case [6,7,11–18]. In the lattice, the analysis and visualization of center-vortex configurations [19] led to important insights regarding the origin of the topological charge density in the YM vacuum. In 3d (4d), thin center vortices are localized on worldlines (worldsheets) *ω*. In this case, the Wilson loop in Equation (1) yields a center element

$$\mathcal{W}\_{\rm D}(\mathcal{E}\_{\rm e}) = \mathcal{Z}\_{\rm D}(\mathcal{E}\_{\rm e}) = \frac{1}{\mathcal{\mathcal{Y}}} \text{tr} \left[ \mathbf{D} \left( e^{i \frac{2\pi}{N}} I \right) \right]^{\mathcal{L}(\omega, \mathcal{E}\_{\rm e})},\tag{6}$$

where L(*ω*, Ce) is the total linking number between *ω* and Ce. This result also applies to thick center vortices, when their cores are completely linked by Ce. In this case, *ω* refers to the thick center vortex guiding centers. In the scaling limit, where the lattice calculations make contact with the continuum, the density of thin center vortices detected at low temperatures is finite [7,20]. Furthermore, center vortices percolate and have positive stiffness [21,22], while the fundamental Wilson loop average over *Zμ*(*x*) displays an area law. This is in accordance with center-vortex condensation and the Wilson loop confinement criteria. For *SU*(2), a model based on the projected thin center-vortex ensemble captures 97.7% of the fundamental string tension. On the other hand, the percentage drops to ∼62% for *SU*(3) [23]. One of the most important features of the center-vortex scenario is that it naturally explains asymptotic *N*-ality: the center element contribution in Equation (6) only depends on the *N*-ality of D. For these reasons, it is believed that the confinement mechanism should involve these degrees of freedom. For a recent discussion about this area of research, see [24].

When it comes to accommodating the model-independent full Monte Carlo calculations, some questions arise. In 3d, the full asymptotic string tension dependence on D is very well fitted by the Casimir law [25]

$$
\sigma\_k^{(3)} = \frac{k(N-k)}{N-1} \,\, \, \, \, \, \, \, \tag{7}
$$

which is proportional to the lowest quadratic Casimir among those representations with the same *N*-ality *k* of D, which corresponds to the antisymmetric representation. In addition, it is precisely at asymptotic interquark distances where a model based on an ensemble of thin objects should be more reliable. This is different at intermediate distances, where finite-size effects allowed for an explanation of the observed scaling with the Casimir of D [26,27]. Then, one question is: how to capture the asymptotic law in Equation (7) from an average over percolating thin center-vortices? In 4d, where the available data cannot tell between a Casimir or a Sine law [28]

$$
\sigma\_k^{(4)} = \frac{k(N-k)}{N-1} \quad \text{vs.} \quad \sigma\_k^{(4)} = \frac{\sin k\pi/N}{\sin \pi/N} \;/ \tag{8}
$$

is there any ensemble based on center-vortices that could reproduce one of these behaviors? More importantly, how can one explain this together with the formation of the confining flux tube observed in the lattice? This means reproducing the Lüscher term [29–31] and the observed transverse field distributions (see [32–34], and references therein). Here, we shall review some developments aimed at providing a possible answer to these questions.

In Section 2, we shall discuss the simplest Abelian center-vortex ensembles. In Section 3, we summarize, from different points of view, additional non-Abelian information and correlations that could be natural ingredients to be taken into account. In Section 4, we review ensembles of percolating oriented and non-oriented center vortices in 3d and 4d, their effective field description, as well as the possibility to accommodate the asymptotic properties of the confining string. Finally, in Section 5, we discuss recent lattice results in the light of our effective description, and present some perspectives.

#### **2. Center-Vortex Ensembles**

The idea that center vortices are the dominant degrees of freedom in the infrared regime means, in practice, that the Wilson loop average at asymptotic distances may well be captured by modeling the average of the center-elements in Equation (6). This line of research was mainly explored in the lattice [35] by considering an ensemble of fluctuating worldlines (in 3d) or worldsurfaces (in 4d) with tension and stiffness (see also the discussion at the beginning of Section 3.2). For example, in 4d, a theory of fluctuating center-vortex worldsurfaces in four dimensions was introduced by considering the lattice action [35]

$$S\_{\rm latt}(\omega) = \mu \mathcal{A}(\omega) + cN\_{\rm p} \, , \tag{9}$$

where A(*ω*) is the area of the vortex closed worldsurface *ω*, formed by a set of plaquettes, and *Np* is the number of pairs of neighboring plaquettes of the surface lying on different planes. The latter term, as well as the lattice regularization, contribute to the stiffness of the vortices. This model, initially introduced for *SU*(2), and then generalized for *SU*(3) [36], is able to describe important features, such as the confining string tension for fundamental quarks and the order of the deconfinement transition. This type of model can be also formulated in the continuum. The objective is the same, that is, looking for natural ensemble measures to compute center-element averages and compare them with the asymptotic information extracted from the full Monte Carlo average WD(C*e*) . A successful comparison is expected to give important clues about the underlying mechanism of confinement. When computing center-element averages in the continuum, the simplest model has the form:

$$
\langle \mathcal{Z}\_{\rm D}(\mathcal{E}\_{\rm t}) \rangle = \mathcal{N} \sum\_{\omega} e^{-S(\omega)} \frac{1}{\mathcal{J}} \operatorname{tr} \left[ \mathbf{D} \left( e^{i \frac{2\pi}{N}} I \right) \right]^{\mathbf{L}(\omega; \mathcal{E}\_{\rm v})},\tag{10}
$$

where ∑*<sup>ω</sup>* represents the sum over different configurations in a diluted gas of closed worldlines (in 3d) or worldsurfaces (in 4d). The weight factor *e*−*S*(*ω*) implements the effect of center-vortex tension (*μ*) and stiffness (1/*κ*) observed in the lattice [21,22]. More precisely, *S*(*ω*) contains a term proportional to the length or area of *ω*, and another one proportional to a power of the absolute value of the curvature of *ω*. See Equation (A3) for an explicit formula in 3 dimensions. *S*(*ω*) could also contain interactions with a scalar field *ψ* that, when integrated with a corresponding weight *W*(*ψ*), generates interactions among the variables *ω*.

Extended models can also be introduced where the defining elements are not only given by *ω* but also by additional labels. At the level of the gauge field variables *Aμ*, the center-vortex sectors can be characterized by different mappings *S*<sup>0</sup> ∈ *SU*(*N*) containing defects (see Section 4.2). A center vortex with guiding center *ω* and magnetic weight *β* is characterized by *<sup>S</sup>*<sup>0</sup> <sup>=</sup> *<sup>e</sup>*−*iχβ*·*T*, *<sup>β</sup>* · *<sup>T</sup>* <sup>≡</sup> *<sup>β</sup>*|*qTq*, where *<sup>χ</sup>* is a multivalued angle that changes by 2*π* when going around *ω*, and *Tq*, *q* = 1, ... , *N* − 1 are the Cartan generators. As they carry a single weight, these vortices are known as oriented (in the Cartan subalgebra). For elementary center vortices, the tuple *β* is one of the magnetic weights *β<sup>i</sup>* (*i* = 1, ... , *N*) of the fundamental representation. In the region outside the vortex cores, *Aμ* is locally a pure gauge configuration constructed with *S*0. Then, for fundamental quarks, the contribution to a large loop contained in that region is *i*-independent and given by the elementary center-element (1/*N*)tr\$ *e*−*i*2*πβi*·*T*% <sup>=</sup> *<sup>e</sup>i*2*π*/*<sup>N</sup>* to the power L(*ω*, <sup>C</sup>e). Different elementary fluxes may join to form more complex configurations, provided this is done in a way that conserves the flux. For example, *N* center-vortex guiding centers associated with different

magnetic weights *β<sup>i</sup>* can be matched. For simplicity, let us consider the *SU*(3) case in three dimensions and a configuration characterized by *S*<sup>0</sup> = *eiχ*1*β*1·*Tei<sup>χ</sup>*2*β*2·*T*, where *χ*<sup>1</sup> and *χ*<sup>2</sup> are multivalued when going around the closed worldlines *ω*<sup>1</sup> and *ω*2, respectively. These worldlines could meet at a point, then follow a common open line *γ*, and again bifurcate to close the corresponding loops. In this case, we would have a pair of fluxes entering the initial point, carrying the fundamental weights *β*1, *β*2, and a flux leaving along *γ*, carrying the weight *β*<sup>1</sup> + *β*2. In *SU*(3), this sum is an antifundamental weight −*β*3. In other words, there are three fluxes entering the initial point, which carry the three different fundamental weights *β*1, *β*2, *β*3. This can be readily generalized to *SU*(*N*), where *N* fluxes carrying the different fundamental weights can meet at a point, as these weights satisfy ∑*<sup>i</sup> β<sup>i</sup>* = 0. Vortices may also be non-oriented [37], in the sense that they may not be described by a single weight. In this case, the center-vortex components with different fundamental weights are interpolated by instantons in 3d and monopole worldlines in 4d. These lower dimensional junctions, which carry a flux of the form *β<sup>i</sup>* − *βj*, should be weighted with additional phenomenological terms in *S*(*ω*). Furthermore, in the 4d case, three monopole worldlines carrying fluxes *β<sup>i</sup>* − *βj*, *β<sup>j</sup>* − *βk*, *β<sup>k</sup>* − *β<sup>i</sup>* can be matched at a spacetime point. Similar higher-order matching rules are also possible. In what follows, we shall discuss the different ensembles, starting with the simplest possibilities in 3d and 4d.

#### **3. Abelian Effective Description of Center Vortices**

In this section, we shall briefly discuss center-vortex ensembles formed by diluted closed worldlines in 3d (Section 3.1) or worldsurfaces in 4d (Section 3.2), characterized by no other properties than tension, stiffness, and vortex–vortex interactions. No additional degrees of freedom, matching rules or correlations with lower dimensional objects will be considered here.

#### *3.1. Three Dimensions*

In a planar system, thin center vortices are localized on points, so they are created or annihilated by a field operator *V*ˆ(*x*). The emergence of this order parameter can be clearly seen by applying polymer techniques to center-vortex worldlines [38]. In [39], the center-element average for fundamental quarks, over all possible diluted loops, was initially represented in the form

$$
\langle \mathcal{Z}\_{\mathbb{F}}(\mathcal{C}\_{\mathfrak{c}}) \rangle = \mathcal{N} \int [D\psi] \, e^{-\mathcal{W}[\psi]} \, e^{\int\_0^\infty \frac{d\mathcal{L}}{L} \int dx \int du \, Q(\mathbf{x}, \boldsymbol{\mu}, \mathbf{x}, \boldsymbol{\mu}, \mathbf{L})} \,, \tag{11}
$$

where *Q*(*x*, *u*, *x*0, *u*0, *L*) is the integral over all paths with length *L*, starting (ending) at *x*<sup>0</sup> (*x*) with unit tangent vector *u*<sup>0</sup> (*u*), in the presence of scalar and vector sources *ψ* and <sup>2</sup>*<sup>π</sup> <sup>N</sup> sμ*, and weighted by tension and stiffness. The factor *W*[*ψ*] = *<sup>ζ</sup>* 2 *d*3*x ψ*2(*x*) generates, upon integration of the auxiliary scalar field *ψ*, repulsive contact interactions between the loops with strength given by the parameter <sup>1</sup> *<sup>ζ</sup>* . Indeed, as in the exponential we have *x* = *x*0, *u* = *u*0, its expansion generates the diluted loop ensemble. As usual, the factor 1/*L* is to avoid loop overcounting when choosing *x*<sup>0</sup> on a given loop. The external source *s<sup>μ</sup>* is localized on a surface *S*(Ce) whose border is the Wilson loop. As a consequence, it generates the intersection numbers between the loop-variables in *Q* and *S*(Ce), which coincide with the different linking-numbers. Using the large-distance behavior of *Q*(*x*, *u*, *x*0, *u*0, *L*), which satisfies a Fokker–Planck diffusion equation (given by Equations (A1) and (A7), with *b<sup>μ</sup>* Abelian, and *D*(Γ*γ*[*bμ*]) being the complex number Γ*γ*[*bμ*]) we then showed that the ensemble average of center elements becomes represented by a complex scalar field *V*(*x*),

$$\langle \mathcal{Q}\_{\mathbf{F}}(\mathcal{L}\_{\mathbf{t}}) \rangle \approx \mathcal{N} \int [\mathcal{D}V][\mathcal{D}\bar{V}] \, e^{-\int d^3 x \left[\frac{1}{\hbar^2} \sum \nabla D\_{\mu} V + \frac{1}{2\xi} (\nabla V - v^2)^2\right]} \, \_\mu$$

$$v^2 \propto -\mu \kappa > 0, \qquad D\_{\mu} = \partial\_{\mu} - i \frac{2\pi}{N} s\_{\mu} \, . \tag{12}$$

This was obtained for small (positive) stiffness 1/*κ* and repulsive contact interactions. The scalar field *V* is originated due to the approximate behavior of *Q*(*x*, *u*, *x*0, *u*0, *L*) in Equation (A12), which turns the exponential in Equation (11) into a functional determinant. The squared mass parameter of this field is proportional to *κμ*, where *μ* is the center-vortex tension. For percolating objects (*μ* < 0), the *U*(1) symmetry of the effective field theory is spontaneously broken (*κμ* < 0). Among the consequences, we have:


$$S\_{\rm latt}^{(3)} = \beta \sum\_{\mathbf{x}, \mu} \text{Re} \left[ 1 - e^{i\gamma(\mathbf{x} + \boldsymbol{\mu})} e^{-i\gamma(\mathbf{x})} e^{-i\mathbf{x}\_{\boldsymbol{\mu}}(\mathbf{x})} \right] . \tag{13}$$

The external source in Equation (12) translates into the frustration *eiαμ*(**x**) = *e<sup>i</sup>* <sup>2</sup>*<sup>π</sup> <sup>N</sup>* if *S*(Ce) is crossed by the link and is trivial otherwise;


This point of view will be useful to propose other ensemble measures relying on lattice models, as in the case where the derivation of the effective description is not known, see for example Sections 3.2 and 5.3. It is also interesting to see that the initial ensemble properties encoded in Equation (11) are recovered close to the 3d XY model critical point, as expected. Indeed, using the same techniques reviewed in [40] for the case without frustration, the partition function may be formulated in terms of integer-valued divergenceless currents, originated after using the Fourier decomposition

$$e^{\beta \cos \gamma} = \sum\_{b=-\infty}^{\infty} I\_b(\beta) e^{ib\gamma} \, , \tag{14}$$

at every lattice link. The resulting expression turns out to be equivalent to a grand canonical ensemble of non-backtracking closed loops formed by currents of strength |*bμ*| = 1. In the model without frustration, close to the critical point *β*<sup>c</sup> ≈ 0.454 (continuum limit), the relevant configurations are known to be formed by large loops rather than by multiple small loops, and multiple occupation of links is disfavored, thus making contact with the initial properties parametrized in the ensemble (see Table 1 below).

**Table 1.** The correspondence between the effective field and 3d XY model representations of the Abelian center−vortex ensemble.


**Figure 1.** The Wilson loop and the frustration are represented in red and green, respectively. Configurations of type (**a**), which involve sites joined by open lines, do not contribute to the partition function. Only site configurations joined by loops, like the one in (**b**), contribute (with a center-element).

#### *3.2. Four Dimensions*

Regarding the effective description of 4d ensembles based on random surfaces, as in the 3 + 1 dimensional world center vortices are one-dimensional objects spanning closed worldsurfaces, the emergent order parameter would be a *string field*. However, unlike the 3d case, a derivation starting from the ensemble of closed worldsurfaces with stiffness is still lacking. Such generalization should initially describe a growth process where a surface is generated, and then derive a Fokker–Planck equation for the lattice loop-to-loop probability. Similarly to what happens with end-to-end probabilities for polymers, where stiffness is essential to get a continuum limit when the monomer size goes to zero [41,42], curvature effects are expected to be essential for the continuum limit of triangulated random surfaces. Indeed, ensembles of surfaces which consider only the Polyakov (or Nambu-Goto) action leads to a phase of branched polymers [43,44]. On the other hand, in [45], the phase fluctuations of an Abelian string field with frozen modulus were approximated by a lattice *field* theory: the *U*(1) gauge-invariant Abelian Wilson action. In other words, the Goldstone modes for a condensate of one-dimensional objects are gauge fields. Motivated by this enormous simplification and by an analogy with the 3d case, in [46] we proposed a Wilson action with frustration as a starting point to define a measure for percolating center vortices in four dimensions. This proposal will be discussed in Section 5.3. For the time being, we summarize the main initial steps, which are analogous to items 1–4 in Section 3.1:


Thus, the main simplification in 4d is that, in a condensate, the effective description can be captured by a local field. Similarly to 3d, where the soft modes can be read in the phase of the vortex field *<sup>V</sup>*(*x*) <sup>∼</sup> *v eiγ*(*x*), the natural soft modes in 4d are given by a compact gauge field,

$$V(\mathbb{C}) \sim v \, e^{i\gamma\_{\Lambda}(\mathbb{C})}, \qquad \gamma\_{\Lambda}(\mathbb{C}) = \oint\_{\mathbb{C}} d\mathbf{x}\_{\mu} \, \Lambda\_{\mu} \,. \tag{15}$$

**Figure 2.** (**a**) Configurations formed by link variables distributed on plaquettes organized on an open surface do not contribute, as the *V<sup>μ</sup>* link−variables at the surface edges cannot form singlets; (**b**) when they are organized on closed surfaces, singlets can be formed and the group−integral is non−trivial.

#### **4. Center-Vortex Gauge Fields, Matching Rules, and Correlations**

The simplest center-vortex ensembles discussed in Section 3 could provide an important basis to understand the confinement mechanism at asymptotic distances. However, they do not contain enough ingredients to reproduce more intricate properties. In this section, we shall discuss the center-vortex gauge fields and typically non-Abelian elements that could characterize the associated ensembles.

#### *4.1. Thick Center Vortices and Intermediate Casimir Scaling*

Before discussing generalized center-vortex ensembles with matching rules and nonoriented components, let us recall how the consideration of center-vortex thickness and the natural non-Abelian orientations in the gauge group can account for the observed Casimir scaling at intermediate distances. Some ideas along this line were initially pursued in [47]. In [26,27] (see also [48]), a simple model was put forward in the lattice, where the contribution to a planar Wilson loop along a curve C<sup>e</sup> was modeled. The starting point is to postulate an ensemble of thick center vortices whose total flux, as measured by a fundamental holonomy, have different possibilities *<sup>z</sup><sup>j</sup>* <sup>=</sup> *<sup>e</sup>i*2*πj*/*N*, *<sup>j</sup>* <sup>=</sup> 1, ... , *<sup>N</sup>* <sup>−</sup> 1. When a thick center vortex is partially linked, the contribution to the Wilson loop is given by the insertion of a group element *Gj*(*x*, *S*) that depends on the location (*x*) of the center-vortex midpoint (or guiding center) with respect to Ce. It also depends on a group orientation *S*,

$$G\_{\vec{\jmath}}(\mathbf{x}, \mathbf{S}) = \mathcal{S} \mathcal{G}\_{\vec{\jmath}}(\mathbf{x}) \mathbf{S}^{\dagger} \,, \tag{16}$$

where <sup>G</sup>*<sup>j</sup>* <sup>=</sup> exp *i α<sup>j</sup>* · *T* is in the Cartan subgroup and the tuples *α<sup>j</sup>* are formed by model-dependent scalar profiles. These profiles implement the natural condition that *Gj*(*x*, *I*) = *z<sup>j</sup> IN*, if the thick center vortex is fully enclosed by Ce, it is *IN* if it is not enclosed at all, and it gives an interpolating value otherwise. After averaging over random group orientations in [26,27], they arrived at

$$\sigma\_{\mathcal{C}\_{\mathbf{v}}}(\mathbf{D}) \equiv -\sum\_{\mathbf{x}} \frac{1}{A} \ln(1 - \sum\_{j=0}^{N-1} f\_j(1 - \frac{1}{\mathcal{G}} \text{Tr} \, \mathbf{D}\left(\mathcal{G}\_j\right)))\,,\tag{17}$$

where *fj* is the probability that a given plaquette of the planar surface enclosed by C<sup>e</sup> be pierced by the midpoint of a center-vortex of type *<sup>j</sup>*, *<sup>σ</sup>*C*<sup>e</sup>* (D) is the string tension in representation D, and *A* is the minimal area of Ce. At intermediate distances, after some natural approximations, an appropriate choice of profiles, and using the key formula

$$\text{Tr}\left(\mathcal{D}(T\_q)\mathcal{D}(T\_p)\right) = \mathcal{O}\,\delta\_{qp}\frac{\mathcal{C}\_2(\mathcal{D})}{N^2 - 1}\,'\,\tag{18}$$

the Casimir Scaling

$$\frac{\sigma\_{\rm I}(\rm D)}{\sigma\_{\rm I}(\rm F)} = \frac{\rm C\_{2}(\rm D)}{\rm C\_{2}(\rm F)}\tag{19}$$

was obtained. In [26,27], based on a specific choice of probabilities and profiles, it was also possible to reproduce different asymptotic behaviors, such as the Casimir and the Sine law. In Section 5, we shall review a different line based on oriented and non-oriented center vortices, which naturally lead to an asymptotic Casimir law. As these models are generated from weighted center-element averages, they are expected to be applicable in the asymptotic region.

#### *4.2. Center-Vortex Sectors in Continuum YM Theory*

Center vortex correlations were considered for the first time in [3]. In (2 + 1)d Minkowski spacetime, the order–disorder algebra in Equation (2) says that the action of *<sup>V</sup>*ˆ(**x**) on <sup>|</sup>*<sup>A</sup>* gives

$$\mathcal{W}\_{\rm F}(\mathcal{C}\_{\rm e}) \left( \mathcal{V}(\mathbf{x}) | A \right) \rangle = e^{i \frac{2\pi}{N}} \mathcal{W}\_{\rm F}(\mathcal{C}\_{\rm e}) \left( \mathcal{V}(\mathbf{x}) | A \right) \,, \tag{20}$$

if **x** is encircled by Ce, and it leaves the state |*A* unaltered otherwise. Here, |*A* is a state with well-defined shape in the Weyl gauge *<sup>A</sup>*<sup>0</sup> <sup>=</sup> 0, That is, *<sup>V</sup>*ˆ(**x**)|*<sup>A</sup>* is a state where a thin center-vortex is created on top of *Ai*. In particular, the action of *V*ˆ *<sup>N</sup>*(**x**) is trivial. Then, the possible phases were effectively described by a model with magnetic *Z*(*N*) symmetry

$$\mathcal{L} = \partial^{\mu}\mathcal{V}\,\partial\_{\mu}V + m^{2}\,\mathcal{V}V + \frac{\lambda}{2}\left(\mathcal{V}V\right)^{2} + \xi\left(V^{N} + \mathcal{V}^{N}\right). \tag{21}$$

This includes quadratic and quartic correlations, as well as the *N*-th order terms that capture the possibility that *N* vortices may annihilate. The case *m*<sup>2</sup> > 0 would correspond to a Higgs phase where center vortices are in the spectrum of asymptotic states. The case *m*<sup>2</sup> < 0 corresponds to a center-vortex condensate, with *N* degenerate classical vacua, so that *Z*(*N*) is spontaneously broken. For a detailed analysis of this effective description, see [49,50]. In [3], based on the center-vortex operator definition *<sup>V</sup>*ˆ(**x**)|*<sup>A</sup>* <sup>=</sup> <sup>|</sup>*A<sup>S</sup>* , discussed in Section 1, 3d Euclidean vortex Green's functions *<sup>V</sup>*¯(*y*)*V*(*x*) were defined. This was done by considering the YM path-integral over configurations *Aμ* with boundary conditions around the pair of points *<sup>x</sup>*, *<sup>y</sup>* <sup>∈</sup> <sup>R</sup>3, such that a vortex is created at *x*, it is then propagated, and finally annihilated at *y*. When |*x* − *y*| → ∞, an exponential decay would correspond to a Higgs phase and *V* = 0, because of the clustering property. This agrees with the discussion above, where the Higgs phase *m*<sup>2</sup> > 0 is characterized by a *Z*(*N*) symmetric vacuum. On the other hand, a condensate would correspond to a Green's function that tends to a constant.

Now, from the definition of the operator *V*ˆ(**x**), it is clear that it introduces singularities in the gauge fields. If *A* is smooth, the configuration *A<sup>S</sup>* is singular, with a field strength containing a delta-singularity at the center vortex location **x**. As pointed out by 't Hooft, the operator's definition could be made more precise by smearing the singularities over an infinitesimal region around **x**. Otherwise, we would be working with singular infinite action gauge fields. Although this direction was not pursued in that work, the smeared Green's functions could depend on the choice of boundary conditions, for the mapping *<sup>S</sup>* <sup>∈</sup> *SU*(*N*), around *<sup>x</sup>* and *<sup>y</sup>*. In other words, the vortex field *<sup>V</sup>*<sup>ˆ</sup> could hide non-Abelian degrees of freedom which are not evidenced by the algebra in Equation (2), which only depends on properties with respect to the Wilson loop.

In [51], we proposed a partition of the full configuration space of *smooth* gauge fields {*Aμ*} into sectors V(*S*0) ⊂ {*Aμ*} characterized by topological labels *S*0. For this objective, we introduced *N*<sup>f</sup> auxiliary adjoint scalar fields *ψ<sup>I</sup>* by means of an identity in the YM path integral, which constrain them to be a solution to a classical equation of motion for the minimization of an auxiliary action *S*aux(*ψI*, *A*). Imposing regularity and boundary conditions, the solution *ψI*(*A*) is unique, and can be decomposed by means of a gener-

alized polar decomposition *<sup>ψ</sup>I*(*A*) = *SqIS*−1, where *<sup>S</sup>*(*x*) <sup>∈</sup> *SU*(*N*) and (*q*1, ... , *qN*<sup>f</sup> ) is a "modulus" tuple. The phase defects cannot be eliminated by regular gauge transformations *U*, which act on the left *S* → *US*. A gauge field is then said to belong to a given sector V(*S*0) if *S*(*A*) is equivalent to a class representative *S*0. The continuum of possible labels *S*<sup>0</sup> are characterized by the location of oriented and non-oriented center-vortex guiding centers, with all possible matching rules (see the discussion in Section 2). Although a possible label for an oriented center-vortex would be *S*<sup>0</sup> = *eiχβ*·*T*, a typical non-oriented configuration is characterized by *S*<sup>0</sup> = *eiχβ*·*TW*. In 3d, close to some points (instantons) on the center-vortex worldline generated by *eiχβ*·*T*, the mapping *W* behaves as a Weyl transformation that changes the fundamental weight *β* to *β* . Similarly, in 4d, the change occurs at some monopole worldlines on the center-vortex worldsurfaces generated by *eiχβ*·*<sup>T</sup>* (see [46]). The full YM partition function and averages of observables were then represented by a sum over partial contributions,

$$Z\_{\rm YM} = \sum\_{\rm S\_0} Z\_{(\rm S\_0)^\star} \quad \text{ (O)}\\ \chi\_{\rm YM} = \frac{1}{Z\_{\rm YM}} \sum\_{\rm S\_0} \int\_{\mathcal{V}(\mathcal{S}\_0)} [D\mathcal{A}\_\mu] \, O \, e^{-S\_{\rm YM}} \,. \tag{22}$$

Here, ∑*S*<sup>0</sup> is a short-hand notation for the contribution originated from the continuum of labels *S*0. These ideas provided a glimpse of a path connecting first principles Yang–Mills theory to an ensemble containing all possible center-vortex configurations. In addition to addressing this important conceptual issue, the partition into sectors may circumvent the well-known Gribov problem when fixing the gauge in non-Abelian gauge theories, as Singer's no go theorem [52] only applies to global gauges in configuration space (see [53] for a detailed discussion). In [51], the gauge was locally fixed by a regular gauge transformation that rotates *S* to the reference *S*0, which is imposed by a sector dependent condition *fS*<sup>0</sup> (*ψ*) = 0. Furthermore, the theory was shown to be renormalizable in the vortex-free sector [54]. The extension of the renormalization proof to sectors labeled by center vortices is under way, and will be presented elsewhere. An interesting consequence of this construction is that a new label may be generated by the right multiplication, *<sup>S</sup>*<sup>0</sup> <sup>→</sup> *<sup>S</sup>*0*U*˜ <sup>−</sup>1, with regular *U*˜ , which is not necessarily connected to *S*<sup>0</sup> by a regular gauge transformation. That is, given a center-vortex sector, there is a continuum of physically inequivalent sectors characterized by non-Abelian d.o.f. where the defects are located at the same spacetime points. In the context of effective Yang–Mills–Higgs models, which describe the confining string as a smooth topological classical vortex solution, the presence of similar internal d.o.f. was previously noted in a large class of color-flavor symmetric theories [55–64].

#### **5. Mixed Ensembles of Oriented and Non-Oriented Center Vortices**

The general properties of center vortices discussed so far motivate the search for a natural ensemble that captures all the asymptotic properties of confinement. Among them, the formation of a confining flux tube is the most elusive one in this scenario. The formation of this object would also explain the Lüscher term, which has not been observed in projected center-vortex ensembles. Furthermore, the asymptotic Casimir law (cf. Equation (7)) should be reproduced in 3d, while in 4d we would like to understand the coexistence of *N*-ality with the Abelian-like flux tube profiles [32–34]. It is clear that a confining flux tube requires an ensemble whose effective description contains topological solitons, namely, a confining domain wall in (2 + 1)d and a vortex in (3 + 1)d. However, the simple models of oriented and uncorrelated center vortices discussed in Section 3 do not have the conditions to support these topological objects1. In what follows, we shall review how the inclusion of the center-vortex matching rules and correlations with lower dimensional defects (see Sections 2 and 4.2) could fill the gap between center-vortex ensembles and the formation of a flux tube. In [65,66], lattice studies showed that the 4d Abelian-projected lattice is not represented by a monopole Coulomb gas, but rather by collimated fluxes attached to the monopoles. In the continuum, these configurations correspond to the previously discussed non-oriented center vortices. While in 4d the lower dimensional defects on center-vortex worldsurfaces are monopole worldlines, in 3d they

are instantons. The relevance of non-oriented center vortices to generate a non-vanishing Pontryagin index was shown in [37]. Now, although oriented and non-oriented center vortices, located at the same place, would contribute to a large Wilson loop with the same center-element, it is natural to weight them with different effective actions. In the second case, the measure should also depend on the location of the lower-dimensional defects.

#### *5.1. 3d Ensemble with Asymptotic Casimir Law*

In this section, we review the mixed ensembles formed by oriented and non-oriented center-vortices with *N*-line matching rules introduced in [67]. In that reference, to prepare the formalism so as to include the different correlations, we initially wrote the contribution to the Wilson loop of a thin center-vortex loop *l* as

$$\mathcal{W}\_{\rm D}(\mathcal{E}\_{\rm e})|\_{\rm loop} = \frac{1}{N} \operatorname{Tr} \Gamma\_l[b\_{\mu}^{\mathcal{E}\_{\rm e}}], \quad \Gamma\_{\gamma}[b\_{\mu}] = P\{e^{i\int\_{\gamma} dx\_{\mu} b\_{\mu}}\} \,, \tag{23}$$

where *b*C<sup>e</sup> *<sup>μ</sup>* = 2*πβ<sup>e</sup>* · *T s*C<sup>e</sup> *<sup>μ</sup>* , *β<sup>e</sup>* is the highest magnetic weight of D, and *s* Ce *<sup>μ</sup>* is a source localized on Ce. Here, we use the notation *β<sup>e</sup>* · *T* = *βe*|*qTq*, with *Tq*, *q* = 1 ... , *N* − 1 being the Cartan generators of *SU*(*N*). Then, after weighting each loop with a phenomenological factor *e*−*S*(*l*) accounting for tension and stiffness (cf. Equation (10)), and summing over all possible diluted loops, we obtained the center-element average

$$
\langle \mathcal{Z}\_{\rm D}(\mathcal{C}\_{\rm e}) \rangle = e^{\int\_0^{\infty} \frac{d L}{\mathcal{L}} \int d\mathbf{x} \int d\mathbf{u} \,\mathrm{tr} \, Q(\mathbf{x}, \mu, \mathbf{x}, \mu, \mathbf{L})} \,, \tag{24}
$$

where *Q*(*x*, *u*, *x*0, *u*0, *L*) is the integral over all the paths with length *L* that begin at *x*<sup>0</sup> with unit tangent vector *u*0, and end at *x* with orientation *u*. This is given by Equation (A1), using as D the fundamental representation. This object satisfies a non-Abelian diffusion equation whose large *κ*-limit (small stiffness) solution (cf. Equation (A12)) led to approximate Equation (24) by

$$\langle \mathcal{Z}\_{\rm D}(\mathcal{C}\_{\rm e}) \rangle \approx Z\_{\rm loops} = \mathcal{N} \int [d\phi] \, e^{-\int d^3 x \, \phi^\dagger \mathcal{O} \phi} , \quad \mathcal{O} = -\frac{1}{3\kappa} (I\_N \partial\_\mu - i b\_\mu^{\mathcal{C}\_\mathbf{e}})^2 + \mu I\_N \, . \tag{25}$$

where *φ* is an emergent complex scalar field in the fundamental representation.

One basic defining property of center vortices is that *N* such objects can be virtually created out of the vacuum at *x*<sup>0</sup> and then annihilated at *x*. At the level of the gauge fields, this is related to the possibility of matching *N* guiding centers each one carrying a different fundamental magnetic weight *βi*, *i* = 1, ... , *N*, which satisfy *β*<sup>1</sup> + ... + *β<sup>N</sup>* = 0. Then, to incorporate all possible oriented center-vortex line matchings (see Section 5.1.1), we expanded the loop ensemble in Equation (25) considering the *N* types of weights, each one represented by a fundamental field *φi*, *i* = 1, ... , *N*. At this point, the center-element average over loops was generated from the partition function

$$Z\_{\rm loops}^{N} = \int [D\Phi^{\dagger}][D\Phi] \, e^{-\int d^{3}x \left[\frac{1}{N\pi} \text{Tr}\left((D\_{\mu}\Phi)^{\dagger}D\_{\mu}\Phi\right) + \mu \text{Tr}(\Phi^{\dagger}\Phi)\right]}\,\tag{26}$$

where Φ is a complex *N* × *N* matrix with components Φ*ij* = *φj*|*i*.

#### 5.1.1. Including *N*-Vortex Matching

The contribution to the Wilson loop of *N* center-vortex worldlines starting at *x*<sup>0</sup> and ending at *x*, and carrying different weights, was rewritten as

$$\mathcal{W}\_{\rm D}(\mathcal{E}\_{\rm e})|\_{N-\text{lines}} = \frac{1}{N!} \boldsymbol{\varepsilon}\_{\dot{i}\_1 \ldots \dot{i}\_N} \boldsymbol{\varepsilon}\_{\dot{i}\_1' \ldots \dot{i}\_N'} \boldsymbol{\Gamma}\_{\gamma 1} [\boldsymbol{b}\_{\mu}^{\mathcal{E}\_{\rm e}}]|\_{\dot{i}\_1 \dot{i}\_1'} \ldots \boldsymbol{\Gamma}\_{\gamma N} [\boldsymbol{b}\_{\mu}^{\mathcal{E}\_{\rm e}}]|\_{\dot{i}\_N \dot{i}\_N'} \,. \tag{27}$$

By weighting each line in Equation (27) with the factor *e*−*S*(*γi*), and integrating over paths with fixed endpoints and over all the lengths *Li* (cf. Equation (A13)), we obtained

$$\mathcal{L}\_N \propto \int d^3 \mathbf{x} d^3 \mathbf{x}\_0 \epsilon\_{i\_1 \dots i\_N} \epsilon\_{j\_1 \dots j\_N} G(\mathbf{x}, \mathbf{x}\_0)\_{i\_1 j\_1} \dots G(\mathbf{x}, \mathbf{x}\_0)\_{i\_N j\_N} \tag{28}$$

where *G*(*x*, *x*0) is the Green's function of the operator O. In this manner, the *N*-line contribution in Equation (28) and similar processes were generated by adding a term ∝ (det Φ + det Φ†). The effective description thus obtained is separately invariant under local and global *SU*(*N*) symmetries *Sc*(*x*), *Sf* ∈ *SU*(*N*)

$$\begin{aligned} \Phi &\rightarrow S\_{\mathfrak{c}}(\mathbf{x})\Phi, \quad b\_{\mu} \rightarrow S\_{\mathfrak{c}}(\mathbf{x})b\_{\mu}S\_{\mathfrak{c}}^{-1}(\mathbf{x}) + iS\_{\mathfrak{c}}(\mathbf{x})\partial\_{\mu}S\_{\mathfrak{c}}^{-1}(\mathbf{x}),\\ \Phi &\rightarrow \Phi S\_f \,. \end{aligned} \tag{29}$$

In the effective description, other natural terms compatible with these symmetries, like the vortex–vortex interaction Tr(Φ†Φ)2, should also be included, thus leading to the center-element average ZD(Ce) = *Z*v[*b*C<sup>e</sup> *<sup>μ</sup>* ]/*Z*v[0],

$$Z\_{\rm V}[b\_{\mu}^{\mathcal{L}\_{\Phi}}] = \int [D\Phi^{\dagger}][D\Phi]e^{-\int d^{3}x \left[\frac{1}{2}\text{Tr}\left((D\_{\mu}\Phi)^{\dagger}D\_{\mu}\Phi\right) + \mu\text{Tr}(\Phi^{\sigma}\Phi) + \frac{\lambda\_{0}}{2}\text{Tr}(\Phi^{\sigma}\Phi)^{2} - \xi\_{0}(\text{det}\Phi + \text{det}\Phi^{\sigma})\right]}\,. \tag{30}$$

This effective description has some similarities with the 't Hooft model (cf. Equation (21)). More specifically, they coincide for configurations of the type Φ = *V IN*. However, there is no reason for the path-integral to favor this type of restricted configuration. Up to this point, in the percolating phase (*μ* < 0), the quadratic and quartic terms tend to produce a manifold of classical vacua labeled by *U*(*N*), while the addition of the det Φ-interaction reduces this manifold to *SU*(*N*). Then, unlike the 't Hooft model, in the SSB phase this effective description has a continuum set of classical vacua which precludes the formation of the stable domain wall. It is interesting to formulate the Goldstone modes *V*(*x*) ∈ *SU*(*N*) in the lattice, which leads to

$$S\_{\rm latt}^{(3)}(b\_{\mu}^{\mathcal{L}\_{\rm r}}) = \vec{\beta} \sum\_{\mathbf{x}, \mu} \text{Re} \left[ \mathbb{I} - \vec{\mathcal{U}}\_{\mu} V(\mathbf{x} + \boldsymbol{\hat{\mu}}) V^{\dagger}(\mathbf{x}) \right],\tag{31}$$

where *<sup>U</sup>μ*(*x*) = *<sup>e</sup>i*2*πβ*e·*<sup>T</sup>* <sup>∈</sup> *<sup>Z</sup>*(*N*), if the link *<sup>x</sup>*, *<sup>μ</sup>* crosses *<sup>S</sup>*(Ce), and it is the identity otherwise. As expected, in the expansion of the partition function, besides the contribution of sites distributed on links that form loops, there is also one originated from *N* lines that start or end at a common site *<sup>x</sup>*. In the former case, the singlets are included in *<sup>N</sup>* <sup>⊗</sup> *<sup>N</sup>*¯ , while in the latter they are in the products of *N V*(*x*) or *V*†(*x*) (compare with the Abelian case in Section 3.1). In this way, the rules originating Equation (30) can be recovered in the lattice. This type of cross-checking is useful to better understand proposals of lattice ensemble measures in situations where it is harder to derive the effective field description, like in 4d spacetime.

#### 5.1.2. Including Non-Oriented Center Vortices in 3d

In terms of Gilmore–Perelemov group coherent-states (see [68,69] for a complete discussion or [46] for a summary of the main ideas) |*g*, *ω* = *g*|*ω* , *g* ∈ *SU*(*N*), Equations (23) and (27) became

$$\begin{split} \mathcal{W}\mathsf{W}\big(\mathcal{C}\_{\mathsf{e}}\big)|\_{\mathsf{loop}} &\propto \int d\mu(\boldsymbol{g}) \left< \boldsymbol{g}, \omega \boldsymbol{\nu} \middle| \Gamma\_{\mathbb{I}}[\boldsymbol{b}\_{\mu}^{\mathcal{C}\_{\mathsf{e}}}] \big| \boldsymbol{g}, \boldsymbol{\omega} \right> \,, \\ \mathcal{W}\big(\mathcal{W}\_{\mathsf{e}}\big)|\_{\mathsf{N}-\mathsf{linear}} &\propto \int d\mu(\boldsymbol{g}) d\mu(\boldsymbol{g}\_{0}) \left< \boldsymbol{g}, \omega \boldsymbol{\nu} \middle| \Gamma\_{\mathbb{I}\mathbb{I}}[\boldsymbol{b}\_{\mu}^{\mathcal{C}\_{\mathsf{e}}}] \big| \boldsymbol{g}\_{0}, \omega\_{1} \right> \,, \dots \left< \boldsymbol{g}, \omega\_{\mathsf{N}} \middle| \Gamma\_{\mathbb{I}\mathbb{N}}[\boldsymbol{b}\_{\mu}^{\mathcal{C}\_{\mathsf{e}}}] \big| \boldsymbol{g}\_{0}, \omega\_{\mathsf{N}} \right> \,. \end{split} \tag{32}$$

The first contribution can be thought of as associated to the creation of a center-vortex with initial fundamental weight *ω* and group orientation *g*, which is propagated along the closed worldline *l*, and is then annihilated. The second corresponds to *N* vortices with different magnetic weights *β<sup>i</sup>* = 2*N ωi*, *i* = 1, ... , *N*, created out of the vacuum at a spacetime point *x*0, that follow separate worldlines *γ<sup>i</sup>* and then annihilate at *xf* . Following a similar interpretation, and recalling that the center-vortex weights change at the instantons, we introduced non-oriented center vortices. When a closed object is formed by n parts *γ*<sup>1</sup> ◦ *γ*<sup>2</sup> ◦ ... ◦ *γ<sup>n</sup>* with *n* instantons at points *x*<sup>1</sup> ... *xn*, we considered the contribution

$$\mathbb{C}\_{\mathfrak{n}} = \int d\mu(\mathfrak{g}\_1) \dots \int d\mu(\mathfrak{g}\_{\mathfrak{n}}) \langle \mathfrak{g}\_{\mathfrak{l}\mathfrak{l}} \omega | \mathfrak{g}\_{\mathfrak{2}\mathfrak{l}} \omega' \rangle \langle \mathfrak{g}\_{\mathfrak{2}\mathfrak{l}} \omega | \mathfrak{g}\_{\mathfrak{3}\mathfrak{l}} \omega' \rangle \dots \langle \mathfrak{g}\_{\mathfrak{n}\mathfrak{l}} \omega | \mathfrak{g}\_{\mathfrak{l}\mathfrak{l}} \omega' \rangle \times \dots \quad \text{(33)}$$

$$\times \langle \mathfrak{g}\_{\mathfrak{l}\mathfrak{l}} \omega' | \Gamma\_{\gamma\_{\mathfrak{n}}} [b\_{\mu}^{\mathfrak{C}\_{\mathfrak{q}}}] | \mathfrak{g}\_{\mathfrak{n}\mathfrak{l}} \omega \rangle \dots \langle \mathfrak{g}\_{\mathfrak{3}\mathfrak{l}} \omega' | \Gamma\_{\gamma\_{\mathfrak{2}}} [b\_{\mu}^{\mathfrak{C}\_{\mathfrak{q}}}] | \mathfrak{g}\_{\mathfrak{2}\mathfrak{l}} \omega \rangle \langle \mathfrak{g}\_{\mathfrak{2}\mathfrak{l}} \omega' | \Gamma\_{\gamma\_{\mathfrak{1}}} [b\_{\mu}^{\mathfrak{C}\_{\mathfrak{q}}}] | \mathfrak{g}\_{\mathfrak{1}\mathfrak{l}} \omega \rangle \dots$$

Here, a center vortex is propagated along *γ*<sup>1</sup> from *x*1, with orientation *g*<sup>1</sup> and weight *ω*, up to *x*2, with orientation *g*<sup>2</sup> and weight *ω* . At *x*2, keeping the orientation *g*2, the weight changes to *ω* , and then *γ*<sup>2</sup> is followed, etc. This precisely characterizes a non-oriented center vortex, where the flux orientation along the Cartan subalgebra changes. Additionally, notice that |*ω ω*| is the root vector *Eα*, which is in line with the presence of pointlike defects carrying adjoint charge. Moreover, when the chain configuration links the Wilson loop Ce, one of the holonomies Γ*γ*<sup>1</sup> , ... , Γ*γ<sup>n</sup>* gives a center element, while all the others are trivial, thus leading to the expected center-element for a chain, up to a positive and real weight factor. Performing the integrals on the group, we arrived at an additional vertex and the final formula for the ensemble average of WD(Ce), incorporating all the configurations discussed so far,

$$\langle \mathcal{Z}\_{\rm D}(\mathcal{L}\_{\rm e}) \rangle = \frac{Z[b\_{\mu}^{\downarrow\_{\rm e}}]}{Z[0]}, \quad Z[b\_{\mu}] = \int [D\Phi] \, e^{-S\_{\rm eff}(\Phi, b\_{\mu})} \,, \tag{34}$$

$$S\_{\rm eff}(\Phi, b\_{\mu}) = \int d^3 \mathbf{x} \left( \text{Tr} (D\_{\mu} \Phi)^{\dagger} D\_{\mu} \Phi + V(\Phi, \Phi^{\dagger}) \right), \quad D\_{\mu} = \partial\_{\mu} - i b\_{\mu} \,, \tag{35}$$

$$V(\Phi, \Phi^{\dagger}) = \frac{3}{2}\lambda\_0 \text{tr}(\Phi^{\dagger}\Phi + \frac{\mu}{\lambda\_0}I\_N)^2 - \xi\_0(3\kappa)^{\frac{N}{2}}(\det \Phi + \det \Phi^{\dagger}) - 3\theta\_0 \text{tr}\left(\Phi^{\dagger}T\_A \Phi T\_A\right), \tag{36}$$

where *<sup>λ</sup>*0, *<sup>ξ</sup>*0, *<sup>ϑ</sup>*<sup>0</sup> <sup>&</sup>gt; 0, and we have made the redefinition <sup>Φ</sup> <sup>→</sup> <sup>√</sup>3*κ*<sup>Φ</sup> of the field. When vortices with positive stiffness percolate (1/*κ* > 0, *μ* < 0), a condensate is formed. In the parameter region *λ*<sup>0</sup> , *ξ*<sup>0</sup> >> *ϑ*0, the most relevant fluctuations will be parametrized by Φ ∝ *S*, *S* ∈ *SU*(*N*). It is interesting to check in the lattice how the different configuration types are recovered. The additional non-oriented component in the discretized theory is generated from the product of an adjoint variable arising from the new term

$$\text{Tr}\Big(\Phi^\dagger T\_A \Phi T\_A\Big) \sim \text{const.} \text{Tr}(\text{Ad}(S))\Big),\tag{37}$$

at a lattice site *<sup>x</sup>*, with the adjoint contribution in *<sup>N</sup>* <sup>⊗</sup> *<sup>N</sup>*¯ associated with *<sup>V</sup>*(*x*) and *<sup>V</sup>*†(*x*).

#### *5.2. Saddle-Point Analysis in 3d*

For non-trivial *ϑ*, the *SU*(*N*) classical vacua degeneracy is lifted, and the possible global minima become discrete:

$$\Phi = \upsilon \mathcal{Z}\_N, \quad \mathcal{Z}\_N = \left\{ \varepsilon^{i\frac{2\mu\mu}{N}} ; n = 0, 1, \dots, N - 1 \right\} \; , \tag{38}$$

$$\left(\mathfrak{H}\lambda\_0 \mathfrak{x} N\left(v^2 + \frac{\mu}{\lambda\_0}\right) - 2\mathfrak{f}\_0(3\kappa)^{\frac{M}{2}} N v^{N-2} - 3\kappa \theta\_0(N^2 - 1) = 0\right). \tag{39}$$

Thus, the presence of instantons opens the possibility of stable domain walls that interpolate the different vacua. In this case, the calculation may be approximated by a saddle-point expansion. Considering a large circular Wilson Loop C<sup>e</sup> centered at the origin of the *x*<sup>2</sup> − *x*<sup>3</sup> plane, the effect of the source is simply to impose the boundary conditions

$$\lim\_{\mathbf{x}\_1 \to -\infty} \Phi(\mathbf{x}\_1, \mathbf{x}\_2, \mathbf{x}\_3) = \upsilon I\_{N\prime} \qquad \lim\_{\mathbf{x}\_1 \to \infty} \Phi(\mathbf{x}\_1, \mathbf{x}\_2, \mathbf{x}\_3) = \upsilon e^{2\pi \boldsymbol{\beta}\_{\boldsymbol{\theta}} \cdot \mathbf{T}}, \qquad (0, \mathbf{x}\_2, \mathbf{x}\_3) \in \mathcal{S}(\mathcal{C}\_{\boldsymbol{\theta}}) \,. \tag{40}$$

In [67], we showed that the Ansatz

$$\Phi = (\eta I\_N + \eta\_0 \beta \cdot T)e^{i\theta \beta \cdot T}e^{i\alpha} \tag{41}$$

closes the equations of motion, yielding scalar equations for the profiles *η*, *η*0, *θ*, *α*. Due to the relation *ei*2*πβe*·*<sup>T</sup>* = *e*−*<sup>i</sup>* <sup>2</sup>*k<sup>π</sup> <sup>N</sup>* , the boundary conditions (40) may be imposed either by a solution where *α* varies with *θ* constant, or vice versa. The first possibility is closely related to the 't Hooft model (cf. Equation (21)). In the second case, the *θ* variation is governed by the Sine-Gordon equation

*∂*2 *<sup>x</sup>*<sup>1</sup> *<sup>θ</sup>* <sup>=</sup> <sup>3</sup>*κϑ*<sup>0</sup> <sup>2</sup> sin(*θ*) . (42)

In this manner, for quarks with *N*-ality *k*, we obtained the asymptotic Casimir Law

$$
\epsilon\_k = \frac{k(N-k)}{N-1} \epsilon\_1 \,\,\,\,\,\tag{43}
$$

where <sup>1</sup> is proportional to the Sine-Gordon parameter 3*κϑ*0.

#### *5.3. A 4d Ensemble with Asymptotic Casimir Law*

Here, we review the ensembles of oriented and non-oriented center vortices in four dimensions as proposed in [46]. In that study, instead of deriving the effective description of center-vortex ensembles with negative tension and positive stiffness, we started the discussion from the natural Goldstone modes defined on the lattice (see also Section 3.2). The missing steps are expected to be implemented by deriving diffusion loop equations including the effect of stiffness. The lattice description of an Abelian ensemble of worldsurfaces coupled to an external Kalb–Ramond field in the form

$$\int d\sigma\_1 d\sigma\_2 \, B\_{\mu\nu}(X(\sigma\_1, \sigma\_2)) \Sigma^{\mu\nu}(X(\sigma\_1, \sigma\_2)), \qquad \Sigma^{\mu\nu} = \frac{\partial X^\mu}{\partial \sigma\_1} \frac{\partial X^\nu}{\partial \sigma\_2} - \frac{\partial X^\nu}{\partial \sigma\_1} \frac{\partial X^\mu}{\partial \sigma\_2} \,, \tag{44}$$

where *Xμ*(*σ*1, *σ*2) is a parametrization of the worldsurface, was obtained in [45]. This was done in terms of a complex-valued string field *V*(*C*), where *C* is a closed loop formed by a set of lattice links. The associated action is

$$S\_V = -\sum\_{\mathbb{C}} \sum\_{p \in \eta(\mathbb{C})} \left[ \mathcal{V}(\mathbb{C} + p) \mathcal{U}\_p V(\mathbb{C}) + \mathcal{V}(\mathbb{C} - p) \mathcal{U}\_p V(\mathbb{C}) \right] + \sum\_{\mathbb{C}} m^2 \mathcal{V}(\mathbb{C}) V(\mathbb{C}) \,. \tag{45}$$

*η*(*C*) is the set of plaquettes that share at least one common link with *C*, while *C* + *p* is the path that follows *C* until the initial site of the common link, then detours through the other three links of *p*, and continues along the remaining part of *C*. In addition, the coupling (44) originates the plaquette field *Up* = *<sup>e</sup>ia*2*Bμν*(*p*) . Then, the following polar decomposition was considered

$$V(\mathbb{C}) = w(\mathbb{C}) \prod\_{l \in \mathbb{C}} V\_{l\prime} \qquad V\_l \in \mathcal{U}(1) \; , \tag{46}$$

with a phase factor that has a "local" character, as it was written in terms of the holonomy along *C* of gauge field link-variables *Vl*. Finally, when a condensate is formed (*m*<sup>2</sup> < 0), it was argued that the modulus is practically frozen2, so that *<sup>w</sup>*(*C*) <sup>≈</sup> *<sup>w</sup>* <sup>&</sup>gt; 0. By using this fact in Equation (45), the only links whose contribution do not cancel are those belonging to *p*:

$$\mathcal{V}(\mathbb{C} + p)! I\_p V(\mathbb{C}) = w^2 \prod\_{l \in \mathbb{C} + p} \prod\_{l' \in \mathbb{C}} \mathcal{V}\_l! I\_p V\_{l'} = w^2 \mathcal{U}\_p \prod\_{l \in p} \mathcal{V}\_l \,. \tag{47}$$

Thus,

$$\mathcal{S}\_{\text{latt}}^{(4)}(a\_p) = \vec{\beta} \sum\_p \text{Re} \left[ \mathbb{I} - \mathcal{U}\_p \prod\_{l \in p} V\_l \right]. \tag{48}$$

where the sum is over all plaquettes *p* and a constant was added such that the action vanishes for a trivial plaquette. Then, the description of a loop condensate, where loops are expected to percolate, is much simpler than that associated with a general phase. The string field parameter gives place to simpler gauge field Goldstone variables *V<sup>μ</sup>* = *ei*Λ*μ*(*l*) , governed by a Wilson action with frustration *Up*. This was the starting input used in [46]. An external Kalb–Ramond field that generates the center elements when the simplest center-vortex worldsurface link <sup>C</sup><sup>e</sup> is obtained by replacing *<sup>B</sup>μν* <sup>→</sup> <sup>2</sup>*π<sup>k</sup> <sup>N</sup> sμν*, where *k* is the *N*-ality of the quark representation D and

$$s\_{\mu\nu} = \int\_{S(\mathbb{C}\_{\varepsilon})} d^2 \vec{\sigma}\_{\mu\nu} \delta^{(4)}(\mathbf{x} - X(\sigma\_1, \sigma\_2)) \, , \tag{49}$$

$$d^2 \overline{\sigma}\_{\mu\nu} = \frac{1}{2} \epsilon\_{\mu\nu a\beta} \left( \frac{\partial X^a}{\partial \sigma\_1} \frac{\partial X^\beta}{\partial \sigma\_2} - \frac{\partial X^\beta}{\partial \sigma\_1} \frac{\partial X^a}{\partial \sigma\_2} \right) d\sigma\_1 d\sigma\_2 \tag{50}$$

is localized on *S*(Ce). In the lattice, this localized source corresponds to a frustration *Up* <sup>=</sup> *<sup>e</sup>iα<sup>p</sup>* , where *<sup>α</sup><sup>p</sup>* <sup>=</sup> <sup>−</sup>2*πk*/*<sup>N</sup>* if *<sup>p</sup>* intersects *<sup>S</sup>*(Ce) and it is trivial otherwise. Similarly to the 3d case, we can check a posteriori that the lattice expansion involves an average of center elements over closed worldsurfaces (see Section 3.2). This is a consequence of the properties of *U*(1) group integrals. This also applies to the non-Abelian extension *V<sup>μ</sup>* ∈ *SU*(*N*), governed by

$$S\_V^{\text{latt}}(\mathfrak{a}\_{\mu\nu}) = \tilde{\beta} \sum\_{\mathbf{x}, \boldsymbol{\mu} \preccurlyeq \boldsymbol{\nu}} \text{Re } \text{tr} \left[ \boldsymbol{I} - \boldsymbol{\bar{U}}\_{\mu\nu} \boldsymbol{V}\_{\boldsymbol{\mu}}(\mathbf{x}) \boldsymbol{V}\_{\boldsymbol{\nu}}(\mathbf{x} + \boldsymbol{\hat{\mu}}) \boldsymbol{V}\_{\boldsymbol{\mu}}^{\boldsymbol{\dagger}}(\mathbf{x} + \boldsymbol{\hat{\nu}}) \boldsymbol{V}\_{\boldsymbol{\nu}}^{\boldsymbol{\dagger}}(\mathbf{x}) \right],$$

where plaquettes are denoted as usual. The closed surfaces are generated because *<sup>N</sup>* <sup>⊗</sup> *<sup>N</sup>*¯ contain a singlet. Interestingly, the *SU*(*N*) version has additional configurations where *N* open worldsurfaces meet at a loop formed by a set of links. This is due to the presence of a singlet in the product of *N* link variables. Therefore, the associated normalized partition function

$$\frac{Z\_{\rm v}^{\rm lat}[\boldsymbol{a}\_{\mu\nu}]}{Z\_{\rm v}^{\rm lat}[0]}, \quad Z\_{\rm v}^{\rm lat}[\boldsymbol{a}\_{\mu\nu}] = \int [\mathcal{D}V\_{\mu}] \, e^{-S\_{\rm V}^{\rm lat}(\boldsymbol{a}\_{\mu\nu})} \tag{51}$$

is an average of the center elements generated when a Wilson loop in representation D is linked by an ensemble of oriented center-vortex worldsurfaces with matching rules.

#### *5.4. Including Non-Oriented Center Vortices in 4d*

Although thin oriented or non-oriented center vortices contribute with the same centerelement to the Wilson loop, they are distinct gauge field configurations, with different Yang–Mills action densities. It is then important to underline that the ensemble measure could depend on the monopole component. In order to attach center vortices to monopoles, we included dual adjoint holonomies defined on a "gas" of monopole loops and fused worldlines. In this case, because of the integration properties in the group there are additional relevant configurations like those of Figure 3a,b. The use of adjoint holonomies is in line with the fact that monopoles carry weights of the adjoint representation (the difference of fundamental weights), see [46,51].

**Figure 3.** Non−oriented center vortices containing monopole worldlines. We show a configuration that contributes to the lowest order in *β*˜ (**a**), and one that becomes more important as *β*˜ is increased (**b**). A non-oriented center vortex with three matched monopole worldlines is shown in (**c**).

Then, partial contributions with *n*-loops were generated by

$$\left. \mathcal{Z}\_{\text{mix}}^{\text{lat}} \left[ \boldsymbol{a}\_{\mu\nu} \right] \right|\_{\text{P}} \propto \int \left[ \mathcal{D}V\_{\mu} \right] e^{-S\_{\text{V}}^{\text{lat}}(\boldsymbol{a}\_{\mu\nu})} \left. \mathcal{W}\_{\text{Ad}}^{(1)} \dots \mathcal{W}\_{\text{Ad}}^{(n)} \right. $$

$$\mathcal{W}\_{\text{Ad}}^{(k)} = \frac{1}{N^2 - 1} \text{tr} \left( \prod\_{(\mathbf{x}, \boldsymbol{\mu}) \in \mathcal{C}\_{k}^{\text{lat}}} \text{Ad} \left( V\_{\mu}(\mathbf{x}) \right) \right) . \tag{52}$$

In addition to the matching rules of *N* worldsurfaces, which in the continuum occur as *N* different fundamental magnetic weights add to zero, monopole worldlines carrying different adjoint weights (roots) can also be fused. For example, when *N* ≥ 3, three worldlines carrying different roots that add up to zero can be created at a point. For this reason, we also considered partial contributions to the ensemble like

$$\left. Z\_{\rm mix}^{\rm latter}[a\_{\mu\nu}] \right|\_{\rm P} \propto \int \left[ \mathcal{D}V\_{\mu} \right] \left. e^{-S\_{\rm V}^{\rm latter}(a\_{\mu\nu})} D\_{\rm 3}^{\rm latter} \right. \tag{53}$$

where *D*latt <sup>3</sup> is formed by combining three adjoint holonomies Ad(Γlatt *<sup>j</sup>* ) (see Figure 3c). Other natural rules involve the matching of four worldlines. Then, weighting the monopole holonomies with the simplest geometrical properties (tension and stiffness), the lattice mixed ensemble of oriented and non-oriented center vortices with matching rules can be pictorially represented as

$$Z\_{\rm mix}^{\rm latt}[a\_{\mu\nu}] = \int [\mathcal{D}V\_{\mu}] \, e^{-S\_{\rm \xi}^{\rm att}(a\_{\mu\nu})} \times \dots \, \tag{54}$$

where the dots represent possible combinations of holonomies as illustrated in Figure 4.

Then, noting that *ei*2*πk*/*<sup>N</sup>* = *e*−*<sup>i</sup>* <sup>2</sup>*π β*·*w*<sup>e</sup> , where *β* is a fundamental magnetic weight and *w*<sup>e</sup> is a weight of the quark representation D, we considered the naive continuum limit, *Vμ*(*x*) = *eia*Λ*μ*(*x*) , Λ*<sup>μ</sup>* ∈ su(*N*),

$$Z\_{\rm mix}[s\_{\mu\nu}] = \int [\mathcal{D}\Lambda\_{\mu}] \, e^{-\int d^4 \mathbf{x} \, \frac{1}{4\hbar^2} \left(F\_{\mu\nu}(\Lambda) - 2\pi s\_{\mu\nu} \mathbf{\hat{g}}\_{\mathbf{v}} \cdot \mathbf{T}\right)^2} \times \dots \,\tag{55}$$

The dots represent all possible monopole configurations to be attached to center-vortex worldsurfaces (see Figure 5). Each contribution was obtained using the methods in the Appendix A. The first factor in Figure 5 (monopole loops) generates emergent adjoint fields coupled to the effective field Λ*μ*.

**Figure 4.** Natural combinations of holonomies that can be used to model the mixed ensemble of oriented and non-oriented center vortices. Each contribution is weighted with tension and stiffness.

**Figure 5.** Continuum limit of the monopole sector. The worldline contributions are obtained from the solution to a Fokker–Planck diffusion equation.

For example, a diluted ensemble of a given species of monopoles, with tension *μ*˜ and stiffness <sup>1</sup> *<sup>κ</sup>*˜ , is generated by

$$\frac{1}{2} \int\_{0}^{\infty} \frac{d\mu}{\text{L}} \quad \int d^4 \mathbf{x} \text{ } d\mu \text{ tr } \mathbb{Q}(\mathbf{x}, \mu, \mathbf{x}; \mu, \mathbf{L}) \Big|\_{\mathbf{x}} \tag{56}$$

where *Q* is given by Equation (A1) and D corresponds to the adjoint representation. In the small-stiffness approximation, the non-Abelian diffusion equation for *Q* is solved by Equation (A12), with

$$\mathcal{O} = -\frac{\pi}{12\overline{\kappa}} \left( \partial\_{\mu} - i \operatorname{Ad} \left( \Lambda\_{\mu} \right) \right)^{2} + \overline{\mu} I\_{\overline{\gg}\_{\operatorname{Ad}}} \,. \tag{57}$$

Therefore, the factor in Equation (56) was approximated by

$$e^{-\text{Tr}\,\ln O} = \int [\mathcal{D}\zeta][\mathcal{D}\zeta^{\dagger}] \, e^{-\int d^{4}\mathbf{x}\left((D\_{\mu}\zeta^{\dagger}, D\_{\mu}\zeta) + m^{2}(\zeta^{\dagger}, \zeta)\right)}$$

$$\pi^{2} = (12/\pi)\,\text{j}\mathbb{R}\mathbf{\tilde{\tau}}, \quad D\_{\mu}(\Lambda)\,\zeta = \partial\_{\mu}\zeta - i\left[\Lambda\_{\mu\nu}\zeta\right] \, \, \, \,\tag{58}$$

where *ζ* is an emergent complex adjoint field, and we have introduced the Killing product between two Lie algebra elements *X*,*Y* as (*X*,*Y*) ≡ tr(Ad(*X*)Ad(*Y*)). In the continuum, the path-integral of Ad(Γ[Λ]) over shapes and lengths led to the Green's function for the operator *O*, so that fusion rules like the one in Equation (53) became effective Feynman diagrams. Indeed, to differentiate the monopole lines that can be fused, the monopole loop ensemble was extended to include different species. At the end, a set of real adjoint fields

*ψ<sup>I</sup>* ∈ su(*N*) emerged (*I* is a flavor index). This, together with the non-Abelian Goldstone modes (gauge fields), led to a class of effective Yang–Mills–Higgs (YMH) models,

$$Z\_{\rm mix}[s\_{\mu\nu}] = \int [\mathcal{D}\Lambda\_{\mu}][\mathcal{D}\Psi] \, e^{-\int d^4x \left[\frac{1}{4\pi^2} \left(F\_{\mu\nu}(\Lambda) - 2\pi s\_{\mu\nu}\mathfrak{g}\_{\bullet} \cdot T\right)^2 + \frac{1}{2}(D\_{\mu}\mathfrak{g}\_{I}D\_{\mu}\mathfrak{g}\_{I}) + V\_{\rm It}(\mathfrak{g})\right]} \,. \tag{59}$$

The vertex couplings weight the abundance of each fusion type. Percolating monopole worldlines (positive stiffness and negative tension) favor a spontaneous symmetry breaking phase that can easily correspond to S*U*(*N*) → *Z*(*N*) SSB. This pattern has been extensively studied in the literature (see [55–60,70–72] and references therein).

#### *5.5. Analysis of the Saddle Point in 4d*

In [73–75], we investigated a possible model containing *<sup>N</sup>*<sup>2</sup> <sup>−</sup> 1 real adjoint scalar fields *ψ<sup>I</sup>* and Ad(*SU*(*N*)) flavor symmetry,

$$V\_{\rm H}(\psi) = c + \frac{\mu^2}{2} (\psi\_{A\prime} \psi\_A) + \frac{\kappa}{3} f\_{ABC}(\psi\_A \wedge \psi\_B, \psi\_C) + \frac{\lambda}{4} (\psi\_A \wedge \psi\_B)^2 \,, \tag{60}$$

where *X* ∧*Y* ≡ −*i*[*X*,*Y*]. This model includes some of the correlations previously discussed. The case *μ*˜ = 0 is specially interesting. At this point, the classical vacua are

$$
\Lambda\_{\mu} = \frac{i}{\mathcal{S}} S \partial\_{\mu} S^{-1}, \quad \psi\_A = v S T\_A S^{-1}. \tag{61a}
$$

Then, the Higgs vacua manifold is Ad(*SU*(*N*)) and the system undergoes *SU*(*N*) → *Z*(*N*) SSB, which leads to stable confining center strings. Interestingly, at *μ*˜ = 0, we were able to find a set of BPS equations that provide vortex solutions whose energy is

$$
\epsilon = 2\pi \mathfrak{g} v^2 \mathfrak{\beta} \cdot 2\delta \text{ ,} \tag{62}
$$

where *δ* is the sum of all positive roots of the Lie algebra of *SU*(*N*). Using an inductive proof based on the Young tableau properties, we showed that the smallest *β* · 2*δ* factor is given by the *k*-A weight, the highest weight of the totally antisymmetric representation with *N*-ality *k*. Then, for a general representation D(·) with *N*-ality *k*, the asymptotic string tension satisfies *<sup>σ</sup>*(D)

$$\frac{\sigma(\mathbf{D})}{\sigma(\mathbf{F})} = \frac{\mathbb{C}\_2(k \cdot \mathbf{A})}{\mathbb{C}\_2(\mathbf{F})} = \frac{k(N - k)}{N - 1} \; , \tag{63}$$

which is one of the possible behaviors observed in lattice simulations. Furthermore, the radial energy distribution transverse to the string is *k*(*N* − *k*) times the distribution for a Nielsen– Olesen vortex. For *k* = 1, this agrees with the YM energy distribution of the fundamental confining string, recently obtained from lattice Monte Carlo simulations [32–34].

#### **6. Discussion**

We reviewed ensembles formed by oriented and non-oriented center vortices in 3d and 4d Euclidean spacetime that could capture the confinement properties of *SU*(*N*) pure Yang–Mills theory. Different measures to compute center-element averages were discussed. In 3d and 4d, they include percolating oriented center-vortex worldlines and worldsurfaces that generate emergent Goldstone modes, which correspond to compact scalar and gauge fields, respectively. The models also have the natural matching rules of *N* center vortices, as well as the non-oriented component where center-vortex worldlines (worldsurfaces) are attached to lower-dimensional defects, i.e., instantons (monopole worldlines) in 3d (4d). In addition to the weighting center vortices with tension and stiffness, it is also natural to include additional weights for the lower dimensional defects. In 4d, monopole matching rules are also included. The corresponding effective field content and the SSB pattern may lead to the formation of a confining center string, represented by a domain wall (vortex) in two-dimensional (three-dimensional) real space. The Lüscher term is originated as usual, from the string-like transverse fluctuations of the flux tube. An asymptotic Casimir law can also be accommodated. This asymptotic behavior was observed in 3d, while in 4d it is among the possibilities.

More recently, the transverse distribution of the 4d YM energy-momentum tensor *Tμν* and the field profiles have been analyzed at intermediate and nearly asymptotic distances [32,34,76]. In [34], it was numerically shown that the *Tμν* tensor of the Abelian Nielsen–Olesen (ANO) model cannot fit the *SU*(3) data at the vortex guiding center for *L* = 0.46 fm (intermediate distance) and *L* = 0.92 fm (near asymptotic distance) at the same time. In fact, in [34], it was shown that the components of the energy-momentum tensor at the origin may not be accommodated for *L* = 0.46 fm. Then, on this basis, an ANO effective model to describe the fundamental string was discarded. However, while it is clear that an effective model for the confining flux tube should work at asymptotic distances, it is not that obvious that the same model could be extrapolated to intermediate distances. By intermediate distances we mean those where the string tension scales with the quadratic Casimir of the quark representation. In particular, this is the region where adjoint quarks are still confined by a linear potential, before the breaking of the adjoint string. On the other hand, in the asymptotic region, gluonic excitations around external quarks in a given irreducible representation *D*(·) may be created, so as to produce an asymptotic scaling law that only depends on the *N*-ality of *D*(·). As discussed in this review, the effective field descriptions were derived by considering the (weighted) average of center elements over oriented and non-oriented center vortices, which is expected to be applicable at asymptotic distances. In other words, we wonder if it is meaningful to discard possible effective models on the basis of the lack of adjustment to lattice data on a wide range that includes the intermediate region, where these models are not expected to fully capture the physics. Additionally, note that the known mechanism to explain intermediate Casimir scaling is based on including center-vortex thickness. In turn, these finite-size effects are not included in the ensemble definition that leads to our effective model. Interestingly, while the lattice data rule out the ANO model at intermediate distances *L* = 0.46 fm, such profiles are still among the possibilities at the nearly asymptotic distance *L* = 0.92 fm. Accordingly, the 4d *SU*(*N*) → *Z*(*N*) models we discussed in this review have a point in parameter space where the infinite flux tube profiles Abelianize, while keeping all the required *N*-ality properties. Additionally, the ideas presented in this review imply that not only an asymptotic Casimir law should be observed, but also that the transverse confining flux tube profiles for quarks in different representations should be the same, up to the asymptotic scaling law. This is true for both 3d and 4d, with the profiles being of the Sine-Gordon type in 3d. It would be interesting to test these predictions with lattice simulations.

**Funding:** This research was funded by Conselho Nacional de Desenvolvimento Científico e Tecnológico (CNPq), the Coordenação de Aperfeiçoamento de Pessoal de Nível Superior (CAPES), and the Deutscher Akademischer Austauschdienst (DAAD).

**Conflicts of Interest:** The authors declare no conflicts of interest.

#### **Appendix A. Non-Abelian Diffusion**

Center vortices in 3 dimensions and monopoles in 4 dimensions are propagated along worldlines in Euclidean spacetime. Then, the corresponding ensembles will naturally involve the building block *Q* associated to a worldline with length *L* that starts at *x*<sup>0</sup> with orientation *u*<sup>0</sup> and ends at *x* with final orientation *u*. This is given by

$$Q(\mathbf{x}, \boldsymbol{\mu}, \mathbf{x}\_0, \boldsymbol{\mu}\_0, \boldsymbol{L}) = \int [d\mathbf{x}(\mathbf{s})]\_{\mathbf{x}\_0, \boldsymbol{\mu}\_0}^{\mathbf{x}, \boldsymbol{\mu}} e^{-S(\boldsymbol{\gamma})} \mathbf{D}(\bar{\Gamma}\_{\boldsymbol{\gamma}}[b\_{\boldsymbol{\mu}}]) \, \, \, \, \, \tag{A1}$$

$$\Gamma\_{\gamma}[b\_{\mu}] = P\{e^{i\int\_{\gamma} d\mathbf{x}\_{\mu}b\_{\mu}}\}\,,\tag{A2}$$

where *S*(*γ*) is a vortex effective action, and an interaction with a general non-Abelian gauge field *bμ* was considered. We are interested in the specific form

$$S(\gamma) = \int\_0^L ds \left(\frac{1}{2\kappa} \dot{u}\_\mu \dot{u}\_\mu + \mu\right), \quad u\_\mu(s) = \frac{dx\_\mu}{ds},\tag{A3}$$

which corresponds to tension *μ* and stiffness 1/*κ*. These objects were extensively studied in [46,77]. In what follows, we review the results obtained.

For the simplest center-vortex worldlines in 3d, D is the defining *SU*(*N*) representation, while for monopole worldlines in 4d, D corresponds to the adjoint. To derive a diffusion equation for this object, the paths were discretized into *M* segments of length Δ*L* = *L*/*M*. In this case, the path ordering was obtained from

$$P\{\boldsymbol{\varepsilon}^{-}\}^{\mathrm{l.}\mathrm{d}sH(\mathbf{x}(s),\boldsymbol{\mu}(s))}\} = \boldsymbol{\varepsilon}^{-H(\mathbf{x}\_{M}\boldsymbol{\mu}\_{M})\Delta L} \dots \boldsymbol{\varepsilon}^{-H(\mathbf{x}\_{1}\boldsymbol{\mu}\_{1})\Delta L},\tag{A4}$$

where *H*(*x*, *u*) = −*i*D(*uμbμ*(*x*)). The relation between the building block *QM* associated to a discretized path containing *M* segments of length Δ*L* and that associated with a path of length *L* − Δ*L* is given by:

$$Q\_M(\mathbf{x}, \boldsymbol{\mu}, \mathbf{x}\_0, \boldsymbol{\mu}\_0, L) = \int d^n \mathbf{x}' d^{n-1} \boldsymbol{\mu}' e^{-\mu \Delta L} \boldsymbol{\psi}(\boldsymbol{\mu} - \boldsymbol{\mu}') \times$$

$$e^{-\mu \Delta L} e^{-H(\mathbf{x}, \boldsymbol{\mu}) \Delta L} \delta(\mathbf{x} - \mathbf{x}' - \boldsymbol{\mu} \Delta L) \, Q\_{M-1}(\mathbf{x}', \mathbf{x}\_0, \boldsymbol{\mu}', \boldsymbol{\mu}\_0) \,. \tag{A5}$$

with

$$\psi(\mathfrak{u} - \mathfrak{u}') = \mathcal{N}\mathfrak{e}^{-\frac{1}{2\tilde{\kappa}}\Delta L\left(\frac{\mathfrak{u} - \mathfrak{u}'}{\Delta L}\right)^2} \tag{A6}$$

arising from the discretization of the stiffness term. It acts like an angular distribution in velocity space, which tends to bring *u* close to *u*. Expanding Equation (A5) to first order in Δ*L*, and taking the limit Δ*L* → 0, the diffusion equation

$$\left(\partial\_L - \frac{\kappa \sigma}{2} \hat{L}\_{\mu}^2 + \mu + u\_{\mu} (\partial\_{\mu} - i \mathcal{D}(b\_{\mu})) \right) \mathcal{Q}(\mathbf{x}, \mu, \mathbf{x}\_0, \mathbf{u}\_0, \mathbf{L}) = \mathbf{0} \,, \tag{A7}$$

was obtained, to be solved with the initial condition

$$Q(\mathbf{x}, \boldsymbol{\mu}, \mathbf{x}\_0, \boldsymbol{\mu}\_0, \mathbf{0}) = \delta(\mathbf{x} - \mathbf{x}\_0)\delta(\boldsymbol{\mu} - \mathbf{u}\_0)I\_{\mathbb{R}^d} \,. \tag{A8}$$

D is the dimension of the quark representation D and *L*ˆ <sup>2</sup> *<sup>u</sup>* is the Laplacian on the sphere *Sn*−1. The constant *σ* is given, in *n* spacetime dimensions, by

$$\sigma = \frac{\sqrt{\pi}}{2^{n-3}} \frac{\Gamma\left(\frac{n-2}{2}\right) \Gamma\left(\frac{n+1}{2}\right)}{\Gamma^2\left(\frac{n-1}{2}\right) \Gamma\left(\frac{n-3}{2}\right)} \left(\frac{4\Gamma(n-3)}{\Gamma\left(\frac{n-3}{2}\right)} - \frac{\Gamma(n-1)}{\Gamma\left(\frac{n+1}{2}\right)}\right) \,. \tag{A9}$$

For the cases considered in this review (*n* = 3, 4), *σ* = 1, 2/*π*, respectively. In the limit of small stiffness, there is practically no correlation between *u* and *u*0, which allowed for a consistent solution of these equations with only the lowest angular momenta components:

$$Q(\mathbf{x}, \boldsymbol{\mu}, \mathbf{x}\_0, \boldsymbol{\mu}\_0, L) \approx Q\_0(\mathbf{x}, \mathbf{x}\_0, L), \quad \partial\_L Q\_0(\mathbf{x}, \mathbf{x}\_0, L) = -O Q\_0(\mathbf{x}, \mathbf{x}\_0, L) \,, \tag{A10}$$

$$O = -\frac{2}{(n-1)\sigma\kappa n}(\partial\_{\mu} - iD(b\_{\mu}))^2 + \mu\_{\prime} \quad Q\_0(\mathbf{x}, \mathbf{x}\_0, \mathbf{0}) = \frac{1}{\Omega\_{n-1}}\delta(\mathbf{x} - \mathbf{x}\_0) \tag{A11}$$

<sup>Ω</sup>*n*−<sup>1</sup> being the solid angle of *<sup>S</sup>n*−1. This implies,

$$Q(\mathbf{x}, \boldsymbol{\mu}, \mathbf{x}\_0, \boldsymbol{\mu}\_0, L) \approx \left< \mathbf{x} | \boldsymbol{\varepsilon}^{-L\mathcal{O}} | \mathbf{x}\_0 \right> . \tag{A12}$$

Then, in this limit, we also have

$$\int\_{0}^{\infty} dL \, du \, du\_{0} \int [D\mathbf{x}]\_{\mathbf{x}\_{0}, \mathbf{u}\_{0}}^{\mathbf{x}, \mathbf{u}} e^{-S(\gamma)} \, \mathcal{D}(\Gamma[b]) = \int\_{0}^{\infty} dL \, du \, du\_{0} \, \mathcal{Q}(\mathbf{x}, \mathbf{u}, \mathbf{x}\_{0}, \mathbf{u}\_{0}, L)$$

$$\approx \langle \mathbf{x} | \mathcal{O}^{-1} | \mathbf{x}\_{0} \rangle, \quad \mathcal{O} \, \mathcal{G}(\mathbf{x}, \mathbf{x}\_{0}) = \delta(\mathbf{x} - \mathbf{x}\_{0}) \, I\_{\mathcal{O}} \,. \tag{A13}$$

#### **Notes**


#### **References**

