Different faces of confinement

In this review, we provide a short outlook of some of the currently most popular pictures and promising approaches to non-perturbative physics and confinement in gauge theories. A qualitative and by no means exhaustive discussion presented here covers such key topics as the phases of QCD matter, the order parameters for confinement, the central vortex and monopole pictures of the QCD vacuum structure, fundamental properties of the string tension, confinement realisations in gauge-Higgs and Yang-Mills theories, magnetic order/disorder phase transition among others.


Introduction
Quantum Chromodynamics (QCD) based upon SU(3) c gauge theory of color represents a real-world example of a fundamental Yang-Mills (YM) theory applied to description of strong interactions and is an organic part of the Standard Model (SM) of particle physics. This theory is extremely successful in predicting various measurable phenomena at particle colliders. The class of phenomena that originate from (or driven by) strong interactions is extremely wide and covers such areas as nuclear physics, hadron physics, physics of quark-gluon plasma, high-temperature and high-density QCD, high-energy particle production and hadronisation. Depending on characteristic length scales, QCD behaves very differently. At short space-time separations, e.g. once we zoom into distances much shorter than the proton radius, QCD appears as a weakly-coupled theory enabling a precise Perturbation Theory (PT) analysis. Much of its success has been achieved in this asymptotic freedom or ultraviolet (UV) regime where the quark-gluon interaction strength recedes. Thus success highlights the QCD theory as the correct theory of strong interactions at the fundamental level precisely matching all the existing observations up to very high momentum transfers reached by the Large Hadron Collider (LHC) so far. However, on the opposite side of length-scales in the infrared (IR) limit, QCD enters entirely different, strongly-coupled domain, rendering the PT inapplicable and creating substantial problems for making reliable predictions at intermediate and low momentum transfers, i.e. at large distances. While it is conventionally believed that QCD should remain the correct theory of strong interactions also at large distances, in the so-called confined regime, deriving reliable predictions remains a big theoretical challenge. For one of the broadest and comprehensive overviews of many phenomenological and theoretical aspects of QCD and QCD-like gauge theories spanning from IR to UV, from dilute to dense regimes, see Ref. [1].
The problem of confinement concerns the strongly-coupled sector of QCD composed of interacting colored partons (quarks and gluons). In virtue of color confinement, the colored particles appear to be always trapped (confined) inside colorless composites. The latter emerge as asymptotic states rendering the long-distance regime of hadron physics described by Effective Field Theory (EFT) approaches such as the chiral PT, as well as a variety of non-perturbative techniques realised in numerical simulations on the lattice. Much of the discussion in the current review is devoted to highlighting main ideas and possible existing ways to address the confinement problem that is known as the main unsolved problem in the SM framework. Despite the major efforts of the research community and tremendous progress made over last few decades, it does not appear to be fully and consistently resolved yet. There are several important subtleties in formulation of this problem to be discussed in what follows. One of the standard ways of formulating the problem is that there is no complete understanding of why these fundamental degrees of freedom (DoFs) of QCD (or, generically, of any strongly-coupled YM theory) do not emerge in the physical spectrum of asymptotic states and how the composite hadrons are dynamically produced starting from the fundamental DoFs in the initial state. In a phenomenological sense, there is a fundamental mismatch between the underlined DoFs of QCD in its short-and long-distance regimes manifest in experimental measurements, and there is not a single consistent theoretical framework that goes beyond the framework of PT and treats both weakly-and strongly-coupled regimes on the same footing. For practical purposes, various phenomenological approaches have been proposed that characterise the long-distance effects of QCD absorbing them into universal elements of a given scattering process, such as non-perturbative matrix elements, fragmentation functions or parton distributions. As a commonly adopted picture, a color-electric flux tube (also known as a color string) is stretched among the partons produced in a high-energy collision. A string-like picture emerges in the limit of large number of colors already in D = 1 + 1 dimensions as has been advocated by t'Hooft back in early 70'es -see e.g. Ref. [2]. As produced partons move away from each other at large enough distances, those flux tubes fragment into composite particles such as mesons and baryons where initial (anti)quarks and gluons get necessarily combined with newly emerged ones from the vacuum into color-neutral configurations. In a nutshell, the basic problem concerns a first-principle derivation of long-distance hadron spectrum and dynamics from an underlined strongly-coupled gauge theory. More specifically, a successful model of confinement is expected to provide a first-principle dynamical description of the string formation, its basic characteristics and string-breaking effects, also connecting those unambiguously to dynamics of the fundamental DoFs of the underlined gauge theory and deducing the phase structure of the theory at various densities and temperatures. While there are no compelling solutions yet available, there are several distinct approaches to confinement treatment being actively developed in the literature. Not only a large variety of treatments of confinement has hit the literature in past decades, but also a proper definition of confinement, what we actually mean by this word, posses a notorious difficulty as was thoroughly discussed in Refs. [3,4]. In this review, we will try to summarise some of the existing attractive treatments of confinement and ideas and why confinement occurs in the way it does in a conceptual and qualitative manner, without pretending to provide an exhaustive overview of all relevant details and corresponding references.
The review is organised as follows. In Sect. 2 we discuss the basic ingredients of the QCD phase diagram at different temperatures and values of the baryon chemical potential. In Sect. 3 we provide a brief description of magnetic order/disorder phases and introduce the basic notions of the lattice gauge theory that would be used in follow-up discussions. In Sect. 4 we overview basic concepts and ideas that lead to different asymptotic behaviors of Wilson loop VEV as an order parameter for confining phase. Such distinct properties of QCD scattering amplitudes as the Regge trajectories and the associated picture of a color string have been outlined in Sect. 5. In Sect. 6, we provide a detailed outlook on the complementarity between the Higgs and confining phases and describe such a common feature for both phases as color confinement. In Sect. 7, a brief description of the string hadronisation picture realised in the Lund model is given. Sect. 8 elaborates on why confinement criteria based upon gauge symmetry remnants (un)breaking may be spoiled by gauge-fixing artefacts, highlighting the need for a gauge-invariant description of confinement. Sect. 9 introduces the basics of the center symmetry based confinement criterion and its implications. Sect. 10 gives a brief outlook on another order parameter of confinement, the Polyakov loop, particularly suitable for confinement description at finite temperatures. In Sect. 11, yet another important order parameter of confinement probing the vortex structure of the QCD vacuum, the t'Hooft loop, has been introduced and the basic features of the center vortices have been described. Sect. 12 elaborates on the most important characteristics of the string tension as the probes for a confining phase. The foundations and implications of the center vortex mechanism of confinement, with its basic tests performed in the literature, have been discussed in Sect. 13. Sect. 14 connects the chiral symmetry breaking and the topological charge to the existence of vortex configurations. In Sect. 15, we briefly describe the Gribov-Zwanziger scenario of confinement relating it to the non-perturbative behavior of propagators and describing how a color string could emerge in this scenario by considering constituent gluons in the gluon chain model. An renown dual superconductivity picture of confinement and the fundamental role of magnetic monopoles have been briefly described in Sect. 16. A novel generalisation of the confinement criterion applicable in gauge theories with matter in the fundamental representation has been briefly discussed in Sect. 17. Sect. 18 highlights an important recent development in understanding the confining property of the gauge-field vacuum and Higgs-confinement transitions via a novel non-local order parameter. A summary and concluding remarks are given in Sect. 19.

Phase structure of QCD matter
Following the discovery of asymptotic freedom in QCD [5,6], it has been realised that phase transitions in the hot and dense QCD matter between the hadronic (confined) and quark-gluon (deconfined) phases are crucial for understanding the cosmological evolution as well as the state of matter and dynamics of neutron stars [7][8][9][10][11][12][13][14]. Besides, the idea of experimental measurements through heavy-ion collisions has been offered as a tantalising opportunity for explorations of this interesting physics. In those early times, a hypothetical state of QCD matter at characteristic temperatures of around 100 MeV has been envisaged as existing in two possible states of "hadronic plasma" [9] and "quark-gluon plasma" (QGP) [10], with an energy density of order 1 GeV/fm 3 . Later on, it has been understood that the QCD phase diagram has a much richer structure, particularly, at high baryon number densities, with a lot of important implications for understanding, for instance, neutron star physics as well as heavy-ion collisions at particle colliders.
Strongly-interacting QGP was first discovered at RHIC collider in 2005 [15][16][17][18] and later has been confirmed at much higher energies at the CERN LHC (for a detailed review, see e.g. Refs. [19,20] and references therein). In the QGP phase, as the name suggests, the strong interactions between constituents of the plasma, "dressed" light quarks and gluons being its collective excitations, is driven by their SU(3) c color charges. For a comprehensive review of early developments and key ideas in analysis of strongly-coupled QCD phenomena and QGP in particular, see e.g. Ref. [21], while an overview of more recent theoretical and experimental studies can be found in Refs. [19,20,22].
In a weakly-interacting QCD gas at very high T, the microscopic quark-gluons interactions are relatively weak and should obey the predictions of asymptotic freedom. The leading-order perturbative QCD coupling that determines the strength of QCD interactions at asymptotically short distances, is given in terms of the QCD energy scale Λ QCD ≈ O(1 GeV), momentum transfer Q Λ QCD and the active quark flavors' number N f . In a perturbative domain of QCD, when going towards shorter distances l Λ −1 QCD , the color charge is being diluted compared to the "soft" and non-perturbative domain of QCD at larger distances l ∼ Λ −1 QCD where the charge is being built-up effectively, due to the phenomenon called the color charge "anti-screening" [5,6]. This is quite an opposite effect to what happens in QED. This behavior of the coupling is demonstrated in Fig. 1 (left panel), together with experimentally measured values. As soon as α s (Q) hits large values entering the strongly-coupled (confined) regime at lower T, the PT ceases to work such that effective and non-perturbative methods  next-to-leading order; NNLO: next-to-next-to-leading order; NNLO+res.: NNLO matched to a resummed calculation; N 3 LO: next-to-NNLO).

Deep-inelastic scattering and global PDF fits:
Studies of DIS final states have led to a number of precise determinations of -s: a combination [501] of precision measurements at HERA, based on NLO fits to inclusive jet cross sections in neutral current DIS at high Q 2 , provides combined values of -s at di erent energy scales Q, as shown in Fig [23]. On the left panel, the QCD evolution of the characteristic parton (quark and gluon) density and length-scale with respect to rapidity Y = ln(1/x) and ln Q 2 . The figure is taken from Ref. [24].
are applied, being however often vastly disconnected from the microscopic QCD theory. One could perform a consistent matching of the fundamental QCD to the effective Lagrangian of chiral PT at the "soft" scale Q 4π f π 1 GeV where both descriptions are expected to be valid and overlap. Such a matching provides a clue about the IR behavior of α S (Q) that tends to get "frozen" at the value of α S IR 0.56 [25].
Besides the weakly-coupled short wavelength modes of partonic DoFs with Q = 2πT dominating thermodynamic evolution at very high T, the QGP features also long wavelength (non-perturbative) modes, with length scales of l > T −1 . The latter modes dominate the evolution at not-so-high T forming a liquid effectively turning QGP into ideal fluid [22,26,27]. The latter fundamental property of QGP has been discovered first at RHIC [15][16][17][18] and then confirmed at the LHC. Other effects of such strongly-interacting QGP are manifested through a collective flow phenomenon [27] as well as in an effective suppression of high-energy partons transiting through a hot and dense deconfined medium [28,29] (for a review, see Ref. [20] and references therein).
Taking the ratio of interaction-to-kinetic energy of the QGP constituents and assuming equal contributions from chromo-electric and chromo-magnetic interactions, one introduces the so-called plasma parameter [30] expressed in terms of the fundamental (quark) and adjoint (gluon) Casimir invariants of SU(3) c , C q and C g , respectively, and T-dependent average distance between the partons a satisfying The latter evolves in T only through N f (T). Weakly-interacting (ideal) plasmas have a very low Γ < 10 −3 , while a strongly-interacting plasma typically has a much larger Γ 1. Taking a nearly-ideal (weakly-coupled) massless QCD gas, for instance, one obtains Γ ∼ α S d 1/3 F serving as a lower estimate for the plasma parameter as it ignores the partonic interactions in the ideal gas approximation. In a realistic case of QGP created in heavy-ion collisions at RHIC, one finds T ≈ 200 MeV and α S = Universe 2021, 7, 9 5 of 57 0.3 − 0.5 with only two relevant active flavours, N F = 2, leading to a value of Γ 1.5 − 6, indeed being deeply inside the strongly-coupled plasma regime.
The QCD evolution of partonic matter in terms of basic kinematic parameters of resolved partons in the medium is illustrated in Fig. 1 (right panel). For instance, developing the partonic cascades in typical momentum transfer Q one resolves the partons with a transverse area 1/Q 2 , such that at larger Q and T ∼ Q one observes a dilution of the parton density controlled by the DGLAP evolution equations (see for instance Refs. [31,32]). One may also observe how the parton density evolves with energy or, more conveniently, with a fraction of light cone momentum taken by a given radiated parton out of a parent particle, x = k + /P + . One may visualise the partonic cascade off the initial particle effectively as Brownian-like motion in the transverse plane that can be considered as the Gribov diffusion process in the evolution "time" Y = ln(1/x). The latter parameter is simply a rapidity difference between the radiated and parent partons, while the diffusion constant is D ∼ α S . Such an evolution is controlled by BFKL equations (for more details, see e.g. Refs. [31,32] and references therein).
The partonic cascade is dominated essentially by soft gluons at high energies or at very small fractions x 1, and they are of the same size at a fixed scale Q. As soon as the parton scattering cross section ∼ α S /Q 2 multiplied by the probability to find a parton at a given Q with a fraction x, xG A (x, Q 2 ), becomes of the order of the geometrical cross-section of an area A occupied by the gluons, ∼ πR 2 A , the gluons start to overlap effectively. Due to a repulsive interaction between gluons, however, their occupation number saturates at f g ∼ 1/α S . In particular, this occurs for gluons with transverse momenta below a certain emergent scale Q s (x), k ⊥ ≤ Q s (x), known as a saturation or "close packing" scale [33] (see also Refs. [34,35]), thus, representing a fixed point in the parton x-evolution. Such a saturation phenomenon is rather generic as an analogical scaling of the density ∼ α −1 characterises various Bose-Einstein condensation phenomena, in particular, those in the Higgs mechanism and in superconductivity [36]. Such a highly coherent gluonic state of matter has properties of a classical field [34] and is known in the literature as the Color Glass Condensate (CGC) [24,35,37], or glasma [38]. Indeed, in the path integral formulation of the SU(N) gauge theory, for instance, one sums over all gauge-field configurations weighted with exp(−iS g /h), where the action can be written as such that g 2 s multipliesh in the exponent. Here, f abc , (a, b, c) ∈ {1, . . . , N 2 − 1} are the SU(N) structure constants. The path integral would be dominated by the classical configurations forh → 0 (classical limit) which is therefore equivalent to taking the weak coupling limit of the theory g 2 s → 0 where the action is large, S g h, and so is the number of quanta in these configurations, f g ∼ S g /h [34]. There are certain reasons to believe that such classical-field configurations should describe the state of cold nuclear matter in initial stages of ultra-relativistic heavy-ion collisions [24,37].
Needless to mention, strongly-interacting QCD exhibits a variety of emergent collective effects and phenomena other then those of QGP that are very difficult to understand and to predict starting from the first-principle microscopic theory of QCD. Observable predictions of the hot/dense QCD theory depend on the equation of state (EoS) of compressed nuclear matter but the latter has not been fully understood yet. This situation is analogical to emergent phenomena in atomic and condensed-matter physics driven by the QED interaction theory at the microscopic level. Notably Universe 2021, 7, 9 6 of 57 enough, besides the hadronic and QGP phases, QCD matter features also other distinct phases predicted in various approaches [39,40]. Among important examples of various realisations of confining non-abelian gauge-field dynamics in cosmology are the relaxation phenomena in real-time cosmological evolution of the QCD vacuum [41] and a possibility of phase transitions in a "dark" strongly-coupled SU(N) gauge sectors [42], both potentially testable via detection of stochastic primordial gravitational-wave spectra in future measurements. The homogeneous gluon condensates in the effective SU(N) theory (such QCD gluodynamics) have also been found to play an important role in generation of the observable cosmological constant [20,41,[43][44][45][46]. For a recent review of implications of the quantum YM vacuum for the Dark Energy problem, see Ref. [47].
Systematic explorations of QCD matter at high densities and temperatures including the search for critical end point (CEP) in the middle of the phase diagram at µ B ∼ 0.4 GeV shown in Fig. 2 has started only about ten years ago. The CEP is located at the end of the first-order phase transition boundary between the hadronic phase and QGP, where a second-order phase transition is predicted to occur. One expects a number of new phenomena in a vicinity of that point [48][49][50][51] that have been searched for by the RHIC Beam Energy Scan program.
Currently, a number of different studies of QCD phases in various parts of the (T, µ B )-diagram are being deployed, both experimentally and theoretically, and a high complexity has started to emerge. Particularly intense are explorations of low µ B 0 [52][53][54][55] and high µ B 100 − 600 MeV [40,49,51,55] domains, with possible transitions in between, also indicated in Fig. 2. Other CEPs may also be expected to emerge such as those for chiral (crossover at low-T, not shown in the figure) and nuclear liquid-gas (in the nuclear matter ground-state at nearly-zero T and µ B = 0.93 GeV) transitions.
More specifically, looking at the QCD phase diagram in Fig. 2 along the direction of increasing baryon chemical potential µ B , we notice that at energies close to the binding energy of bulk nuclear matter the so-called cold nuclear matter phase is found. Interactions between nucleons (quark bound Universe 2021, 7, 9 7 of 57 states) may lead to pairing and di-baryon condensation that spontaneously break the U(1) B baryon number symmetry (see e.g. Refs. [56,57]). Physically, this also means that the system is in a superfluid (confining) phase. This system has a close analogy, for instance, with liquid helium where one also finds the Bose-Einstein condensation and Goldstone modes, both associated with the superfluidity property. The same physics emerges in ordinary nuclear matter based upon nuclear many-body theory which is applicable at not too large densities.
There is still a substantial lack of knowledge on a transition between the cold nuclear matter and high-density QCD phases, particularly relevant for physics of neutron stars. Since QCD is asymptotically free, one can go to very high µ B in the quark-matter phase and employ weak-coupling techniques [58]. In this quark-matter phase of dense QCD, such calculations predict a nearly Fermi-liquid with residual interactions that lead to pairing among quarks in a gauge-dependent way. This is described by means of a gauge-dependent di-quark condensate qq playing a role of an order parameter in the dynamical Higgs mechanism such that we deal with a Higgs phase. Indeed, such a di-quark condensate emerges due to long-range attractive forces between the quarks through a Cooper-like pairs' condensation [59,60]. Such a high-density (baryon) superfluid phase where the SU(3) c gluon field is fully "Higgsed" is known in the literature as a color superconductor 1 (CSC) (for a comprehensive review on key aspects of dense QCD, see Refs. [61,62]). Formation of such Cooper pairs of quarks can be seen in QCD with three massless u, d, s flavours at large baryon number densities featuring the following color and flavor symmetries' reduction [40,63] down to a diagonal subgroup SU(3) c+L+R . The corresponding symmetry transformations involve a simultaneous "rotation" of color and flavor group representations known as the color-flavor locking (CFL). Such a CFL phase is known not to be topologically ordered [64]. Then, in the CFL quark-matter phase one could also find an order parameter for U(1) B symmetry breaking (down to Z 2 ) in analogy to the di-baryon condensate in the nuclear-matter phase -it can be viewed as a cubic power of the di-quark condensate thus being associated with a superfluid flow. In fact, both quark matter and nuclear matter phases were found to be relevant for the EoS of neutron stars (see e.g. Refs. [62,65,66]), and the signatures of possible phase transitions might show up in mass-radii relations for neutron stars and gravitational-wave spectra from neutron star collisions. As at high temperatures no baryon number symmetry breaking occurs, one supposedly crosses the line where U(1) B gets restored when the system heats up. As we noticed above, at low temperatures, both low-and high-density phases have the same order parameter w.r.t. U(1) B breaking, and one of the fundamental open questions is whether a boundary between the quark-matter (Higgs) and nuclear-matter (confinement) phases actually exists. Following Refs. [67][68][69], one could consider a simplified picture of pure QCD and include three massless flavors in a maximally symmetric realisation, such that there is no distinction in symmetry realisations between the hadronic phase and asymptotically high-density phase. The latter means there may be no phase transition that is consistent with identical global symmetry realisations in both regimes, t'Hooft anomalies' matching and with smoothly connecting low-lying excitations (see e.g. Refs. [67,70,71]). Such an assumption has become a working one for many phenomenological studies modelling the EoS for neutron star physics (see e.g. Ref. [72] and references therein). Below, following the recent results of Ref. [58] one may conclude, however, that the Schäfer-Wilczek conjecture about quark-hadron continuity at large µ B may be largely oversimplified. The reality may be even more complex than what emerges in existing theoretical approaches. The basic problem is that there are no well-justified theoretical 1 Real superconductors have observable phenomena such as persistent currents distinguishing the superconducting phase from a normal one. The name "color superconductor" come out rather misleading in a sense that there are no observable persistent color-charge currents associated to this phase [58] Universe 2021, 7, 9 8 of 57 methods available for treatment of the strong-coupling regime of QCD, with a non-zero chemical potential, where lattice simulations may not be very reliable.
Finally, yet another QCD phase which is believed to be located somewhere between the chirally restored and confined phases is known as quarkyonic matter [73] that may also have some relevance for neutron star physics [74]. In the limit of large number of color charges N c , the gluons' contribution scales as ∼ N 2 c compared to that of quarks ∼ N c such that this phase is assumed to have energy densities well beyond Λ 4 QCD . Since gluons are bounded in glueballs, one ends up with N c DoFs in this phase.
Let us now turn to a discussion of methods of the lattice gauge theory that became the main tool for explorations of non-perturbative physics in gauge theories and, in particular, QCD in the strongly-coupled regime and the associated dynamics of confinement, at least, at not too large chemical potentials.

Ising model and lattice gauge theory
To what extent one can expect to derive precision results for low-energy observables from the first-principle QCD theory? A default answer to this question is that we should not expect that, at least, analytically. The collective phenomena that are manifest in the strongly-coupled regime of a gauge theory are so complex that none of the existing analytic approaches captures all the relevant dynamics and yields satisfactory results. At the same time, a theory may remain to be correct even if methods of extracting observable information from it are not perfect or suitable. Often though, we start with a simplified model that hopefully captures the same physics as a realistic one but where we have a better control, and then we abstract the lessons that we learn from such a model back to more complicated theories such as QCD.
Luckily, a precise and reliable analysis is possible, but only numerically. The best available framework so far is the lattice gauge theory providing a first-principle numerical approach for strongly-coupled theories like QCD. In fact, this framework is often considered as a "numerical experiment" and may be regarded as a black-box whose results need to fit a certain theoretical picture of real underlined physical phenomena and objects providing means to understand those phenomena qualitatively. Whether or not the lattice results fit a particular picture of confinement is an ongoing long-standing debate in the literature. For relatively recent detailed reviews on non-perturbative physics and the confinement problem, see e.g. Refs. [3,4,[75][76][77] and references therein. Here and below, we follow the notation adopted in Ref. [3] unless noted otherwise, acknowledging that the latter reference represents one of the most complete, pedagogical and sophisticated reviews available in the literature on what the confinement problem actually is from various perspectives and approaches.
In order to built a consistent picture of confinement, we need to elaborate on such important notions as ordered and disordered systems. One of the simplest examples of lattice field theory follows the basic principles of statistical mechanics where the most relevant properties of these systems are readily seen in the Ising model of ferromagnetism. For illustration, consider a simple system -a square (D = 2), cubic (D = 3) or hypercubic (D > 3) array (or lattice) of atoms, each with two spin states -in external magnetic field h. This system is described by the Hamiltonian, Universe 2021, 7, 9 9 of 57 where s(x) = +1 and −1 would correspond to an atom at a point x with spin up and down, respectively, and we denote here the total number of spins as N. Probability for a specific configuration of spins, {s(x)}, at a given temperature T can be written as In the case of zero external field, h = 0, the system apparently possesses a global Z 2 symmetry w.r.t. transformations such that the mean magnetisation (average spin) vanishes. This is a system in a so-called disordered state. Assume that the spins in the initial state are aligned. The exact Z 2 symmetry means that at any given temperature any finite system would end up in a disordered state provided that one waits for long enough for that to occur. This leads to non-existence of permanent magnets as any alignment of the spins would be destroyed by thermal fluctuations. However, for large N, i.e. for macroscopic magnets, the time between sizable fluctuations that could flip a lot of spins would grow exponentially and eventually exceeds the lifetime of the Universe. For non-zero h, however, the Z 2 symmetry appears to be explicitly broken enabling s = 0 at any temperature. In this case, the system appears to be in an ordered state where a large amount of spins point in the same direction. Now, consider magnetisation of a large system in the limit of vanishing h. One could show that, in general, this quantity is non-vanishing yielding the so-called spontaneous symmetry breaking (SSB) of the global Z 2 symmetry which occurs particularly at low temperatures (ordered state). A global symmetry is said to be broken spontaneously when the Hamiltonian and the corresponding equations of motion are symmetric but the solutions for physical observables (such as the magnetisation introduced above) are not. At high T above a certain critical temperature (Curie temperature), the averaged spin vanishes, and the spin system appears again in a symmetric (disordered) state. Considering the vacuum expectation value (VEV) of a product of two spins we notice, G(r) ≡ s(0)s(r) ∼ exp(−r/l), i.e. it falls off exponentially with the distance between atoms r in a disordered state, where l is the correlation length. There is a phase transition between the ordered and disordered phases of the system at the Curie temperature for any D > 1, while for D = 1 the system is in a disordered phase at any T. The existence of such phase transitions associated with a global symmetry breaking is a generic property of many different systems and is also manifest in strongly-coupled gauge theories as will be discussed below.
Let us further promote the global Z 2 symmetry to a local one whose transformation parameter depends on position of the associated DoFs, ξ(x) = ±1, and can be chosen independently at each site (gauge transformations). For this purpose, let us consider the links of the lattice s µ (x) along each dimension µ = 1 . . . D as dynamical DoFs subjected to the gauge transformation and write down the Hamiltonian of the gauge-invariant Ising model Thereby we arrive at the simplest example of the Z 2 lattice gauge theory. In order to describe such systems, one considers observables that are invariant under gauge transformations. A particularly important class of observables can be obtained by taking the VEV of the so-called Wilson loop -a product of links on the lattice around a given closed contour C [78], The Hamiltonian (13) is given by the simplest Wilson loop given by a plaquette, the minimal closed loop on the lattice. In analogy to the gauged Ising model, in a generic lattice gauge theory described by a certain (discrete or continuous) gauge group G, one starts with the Euclidean action where the link variables are the elements of the gauge group. For instance, in the case of a non-abelian group G ≡ SU(2) the group elements in discretized spacetime are in terms of the lattice spacing a, the gauge coupling g, the Pauli spin matrices σ a , a = 1, 2, 3, and the SU(2) gauge field A a µ (x). By convention, the link variable U µ (x) is associated with a line running from site x on the lattice to a neighbor site x +μ in the positive direction µ. The probability distribution of lattice configurations of the gauge field is found in full analogy to that of the Ising model, namely, where the Euclidean action, also known as the Wilson action, is invariant under local gauge transformations We used the fact that the trace of any SU(2) group element is real. A straightforward extension to SU(N) gauge theory leads to with suitably generalised group elements U µ (x). Expanding the latter in powers of A a µ (x), taking β = 2N/g 2 and turning to the continuum limit of vanishing lattice spacing a → 0, one arrives at the standard expressions for the action and gauge transformations in Euclidean spacetime in terms of the field strength tensor F µν and a gauge group element G(x). Here, the repeated indices are summed over as usual.
The formulation of the lattice gauge theory in Euclidean spacetime has quickly become the cornerstone and the main reference for numerical analysis of basic characteristics of the corresponding quantum field theory (QFT) in Minkowski spacetime (such as its low lying spectrum and the static potential). This is due to a single most important fact that the Euclidean formulation of the field theory is conveniently considered as a statistical (not quantum) system whose analysis can be performed using the power of the lattice Monte Carlo methods. For a detailed description of these methods, see e.g. Ref. [79].
The Euclidean formulation is particularly designed for studies of QFT at finite temperatures in equilibrium and works in Euclidean space with periodic time direction for bosonic fields while fermion fields fulfill antiperiodic boundary conditions in the time direction (for a recent review, see e.g. Refs. [80,81]). A finite T theory is then constructed from its zero-temperature counterpart by replacing bosonic and fermionic four-momenta k µ in Euclidean integrals by 2πnT and (2n + 1)πT, respectively, and then switching from k µ integration to summation over n. In a hot medium, an average momentum transfer is the given in terms of temperature, Q = 2πT. The study of thermodynamics and phase transitions is performed in the Hamiltonian formalism starting from the thermal partition function, and the "time" is Euclidean in the path integral formalism from the beginning, at any temperature. The order of deconfinement phase transition in Euclidean SU(3) lattice gauge theory has been studied in this approach by Monte Carlo methods in Ref. [82].
In the continuum limit, in order to obtain the Minkowski action of the corresponding QFT starting from the thermal theory action in Euclidean spacetime one conventionally adopts the Wick rotation t → −it and A 0 → iA 0 relying on analyticity property of the vector-potential. Then, an assumption that a numerical simulation successfully set up in Euclidean spacetime yields relevant results to the corresponding QFT in Minkowski spacetime would be justified only for smooth transitions between short-distance to long-distance physics enabling analytic (in physical time and in A 0 ) continuations of amplitudes from Minkowski to Euclidean spacetime, and backwards. Indeed, such an assumption is violated in the most general case as stated by the so-called Maiani-Testa no-go theorem [83] related to the "failure" of the Wick rotation mentioned above. Indeed, when going out from thermodynamics approaching the study of bound states, the Wick rotation is applicable only to compute static characteristics of the QCD medium such as vacuum condensates as well as masses of stable particles which are the minority of the QCD spectrum. Resonances such as the majority of mesons, charmed and stranged baryons, tetraquarks, pentaquarks, and hadron molecules are accessible in Euclidean space only indirectly and only under restrictive assumptions. For more details on the associated problems in treatment of two-particle systems, see Ref. [84], while a review on the status of three-particle systems can be found e.g. in Ref. [85].
A manifestation of non-analytic structures (domain walls) in the YM vacuum in physical time has also been discussed recently in the context of the non-stationary background of expanding Universe in Ref. [46]. Such structures were found as attractor cosmological solutions at sufficiently large physical times asymptotically matching the YM dynamics on the Minkowski background. In the essence of Maiani-Testa theorem, such non-analytic (domain-wall) solutions found in (nearly) Minkowski background would in general not match the corresponding lattice simulations in Euclidean spacetime, so their implications for confinement are yet unclear and should be studied separately. As long as such solutions are concerned, one may conjecture that the Euclidean YM field theory predictions match those in Minkowski spacetime only in regions sufficiently far away from the non-analytic phase boundaries. This conjecture however requires further in-depth studies of implications of these novel solutions for confinement dynamics.
Another crucial limitation of Monte-Carlo lattice simulations concerns the thermal gauge theory with non-vanishing chemical potential. Indeed, the action becomes complex if the temperature T and the chemical potential µ are both non-zero, meaning that standard Monte-Carlo methods fail in this case (for a thorough review on this issue, see e.g. Ref. [86]). In particular, due to the sign problem the lattice simulations of QCD at µ B > 0 exhibit difficulties in reproducing the quark-gluon plasma as observed in heavy-ion collisions, even under an assumption of thermal equilibrium. The situation becomes even worse when considering the nuclear matter in neutron stars or collapsing black holes at very large density in the curved spacetime. The way to proceed is to expand the pressure in µ B /T and calculate the physical observables as Taylor expansions in this quantity, see e.g. Ref. [87]. In practice, this requires calculating operators of high order, which are noisy and require very large statistics [88]. Recently an alternative summation scheme for the equation of state of QCD at finite real chemical potential was proposed in [89], designed to overcome those shortcomings. Using simulations at zero and imaginary chemical potentials the extracted LO and NLO parameters describing the chemical potential dependence of the baryon density were extrapolated to large real chemical potentials. Proposed expansion scheme converges faster than the Taylor series at finite density, thus, leading to an unprecedented coverage up to µ B /T ≤ 3.5 and to more precise results for the thermodynamic observables.

Asymptotic behavior of large Wilson loop VEVs
Different phases of a gauge theory are classified based on the behavior of Wilson loop VEVs at large Euclidean times compared to spacial separations, i.e. T E R. Computing those in Euclidean spacetime provides a direct access to the interaction energy between the static field sources in Minkowski QFT when the mass of the sources (and hence the fundamental energy scale of a confining gauge theory) is taken to infinity. Introducing a massive scalar field (a "scalar quark") in an arbitrary representation r to the gauge theory on D-dimensional lattice, the corresponding action is invariant under gauge transformation of the scalar field: µ (x), and the gauge-field holonomy is U(p) for a given plaquette p.
Consider an operator that creates a particle-antiparticle pair in a color-singlet state at a given time T E and separation R, that also creates a color-electric flux tube (or string) stretched between the charges. In the limit of heavy static color-charged sources, m 1 in lattice units, the second term in Eq. (22) may be considered as a small perturbation, so the string-breaking effect can be neglected to a first approximation. Indeed, as matter fields are very heavy in this limit, it would take an infinite energy to pull them out of the vacuum and to place them on mass shell in order for them to bind to the sources and hence to screen their charge. This means that one would stretch the flux tube to an infinite length before it can ever break apart, which is of course an unrealistic but still useful picture to test the confinement property of the quantum vacuum.
Thus, by integrating out φ in the functional integral one finds for the VEV to the leading order in 1/m 2 expansion, where is the VEV of the Wilson loop written in terms of the time-like holonomy U(R, T E ) of the pure gauge theory. Here, the link variables run counterclockwise on a time-like rectangular contour C = R × T E , Universe 2021, 7, 9 13 of 57 the group character is χ r and the sum runs over states with two static charges. In the continuum limit, the corresponding holonomy is given by the path-ordered exponential So the Wilson loop (holonomy) operator in this case represents a rectangular time-like loop describing the creation, propagation and, finally, destruction of two static quark and antiquark placed at certain fixed spacial points. The time-like links in a given Wilson loop, hence, can be considered as the worldlines of static heavy charges.
On the other hand, in the operator formalism one deduces that [3] where ∆E n is the energy of the nth excited state above the vacuum, and in the last part of this relation only the dominant contribution (at large T E ) from the minimum-energy eigenstate has been taken into account. In this case, ∆E min = V r (R) corresponds to the energy difference between two static charges, being in other words the interaction (static) potential between them V r (R). Hence, the VEV of the rectangular Wilson loop is characterised by the potential V(R), which can be inverted as Now consider, for instance, a planar non self-intersecting Wilson loop in U(1) gauge theory, and using the Stokes law, it can be written as where the areal integration represents the magnetic flux and proceeds through the minimal area of the large Wilson loop. Thus, due to additive nature of the flux, such a planar Wilson loop can be arbitrarily split into a product of smaller loops whose areas add up to the one of the large loop Here, the orientations of the smaller loops are chosen in such a way that neighboring contours run in opposite directions to each other. In the case of magnetic disorder, the magnetic fluxes through smaller loops C i (e.g. plaquette variables, in the case of smallest loops) are completely uncorrelated, such that the VEV factorises as where A and A are the larger and smaller Wilson loop areas, respectively. Assuming the absence of light matter fields that could, in principle, screen the color charge of the massive sources, and considering a rectangular Wilson loop with C = R × T E , the magnetically disordered state is characterised by the linear growth of the interaction potential with distance R between the static charges asymptotically, hence, represents a potential of a linear string. Here, V 0 is interpreted as a self-energy contribution, and σ r has the meaning of the string tension in a given group representation r that does not depend on the subloop area A . For an illustration of the total potential interpolating small-R (Coulomb) and large-R (confining) regimes, see is then reproduced for T E R as expected, or for a generic contour enclosing a large minimal area including also a dependence on the perimeter of the contour P(C). Note, the gluon propagator is singular in the UV regime in the continuum limit which generically induces a singular term that is interpreted as a divergent self-energy V 0 of the charged particles and antiparticles propagating in the loop. The latter produces a perimeter-law contribution to the large Wilson-loop VEV in the above expression. So, usually a kind of smearing of the loop via a superposition of nearby loops is required to regularise the Wilson loop in the continuum limit (see e.g. Ref. [90]), while on the lattice such a short-distance regularisation is always implicit. It is straightforward to show that for any gauge group and D = 2, only magnetically disordered phase is realised, reproducing the area-law falloff due to the absence of a Bianchi constraint on the components of the field strength tensor [91]. It is, however, a much harder problem to prove the area-law falloff of large Wilson loop VEVs in a generic YM theory with a non-trivial center symmetry which represents the basic confinement problem (for more details, see below). A remarkable property of a Wilson loop is that it characterises vacuum fluctuations of the gauge field, i.e. without the presence of any external sources, with a spacelike loop C, in terms of the ground-state Ψ 0 of the Hamiltonian of the pure gauge theory. As the spacelike and timelike loops are related by a Lorentz transformation, one deduces that the potential energy of interaction between static charges is directly connected to the gauge-field vacuum fluctuations in the absence of color-charged sources.
In D > 2 lattice, the Bianchi constraint emerges that correlates the field strength values at neighbor sites so that those no longer fluctuate independently from one point to another [92]. The absence of those correlations among the smallest Wilson loops, the plaquette variables, is the single most important requirement that provides the area-law relation for Wilson loops of arbitrary sizes. For D > 2 such correlations disappear and the area-law is established in the strong-coupling limit only, i.e. in the leading order in β 1. In the weakly-coupled regime β 1 in D = 3 + 1 electrodynamics, this property does not hold and one recovers the massless phase instead, with the potential [93] corresponding to a perimeter-law falloff of the Wilson loop VEV, where P(C) = 2T E for a rectangular loop C (with T E R), while the coupling g(R) is a slow function of R that approaches a constant in the Coulomb phase. In non-abelian theories the magnetically disorder phase has been established for sufficiently large Wilson loops using the non-abelian Stokes law (see e.g. Refs. [94][95][96][97][98][99][100][101]), and also employing a finite-range behavior of field strength correlators [102,103]. Let us now briefly discuss one of the most distinctive features of long-range dynamics of QCD associated with Regge trajectories.

Regge trajectories and QCD strings
We have seen that the magnetic disorder phase manifests itself through linear dependence of the static potential, and this behavior is inherent to that of a string. What is the nature of such a "color string", and how is it formed? Which phenomenological implications such strings may have?
In hadronic scattering processes, the t-channel exchanges of QCD resonances are considered to be important at high energies. As suggested by quantum mechanics, a given scattering amplitude can be represented as a series expansion in partial waves, in terms of the Legendre polynomials of the first kind and of order l, P l (cos θ), the scattering angle θ and the partial wave amplitudes a l . For a 2 → 2 process and particles of equal mass, for instance, Considering an exchange of a single resonance only, with spin l 0 and at large s → ∞, the amplitude behaves as A(s, t) ∝ s l 0 , such that by means of the optical theorem, the corresponding total cross section, σ tot ∝ s l 0 −1 . This result does not work very well against the experimental data for an integer value of l 0 . The way out is to adopt that there are several resonances being exchanged in the t-channel that should all be taken into account. This is consistently done in the formalism of the Regge theory operating with an analytical continuation of partial amplitudes a l to the complex angular momentum plane (for a thorough discussion of Regge theory principles and applications, see e.g. Ref. [104]). The poles in this plane are traced out by straight lines known as Regge trajectories, l = α(t), and are associated with particles. The squared mass of an exchanged resonance with spin l corresponds to those t at which l is an integer. As a result of the Regge theory, the asymptotic energy dependence of the scattering amplitude reads As a striking feature of QCD that has not been observed e.g. in the electroweak (EW) theory, the Regge trajectories appear to be almost linear functions, and one of the big questions is which dynamics could provide such a simple behavior confirmed experimentally. Namely, hadrons of a given flavour quantum number appear to lie at almost parallel Regge trajectories. It is clear that such a behaviour must be specific to confining dynamics of QCD. Apparently, the potential that binds the quark and anti-quark together into a meson and rises with the interquark separation linearly should be responsible for such a behaviour. One adopts the physical picture of a string stretched between q andq as a narrow color-electric flux tube, which carries the energy E = σr, so that one can neglect the quark masses. Considering for simplicity the leading Regge trajectory that maximises l at a given t, the flux tube of length r rotates about its center such that its end points move with the speed of light, and √ t = in terms of the string tension σ, and the transverse velocity v ⊥ = 2x/r. Analogously, the angular momentum of such a system providing us finally with the Regge slope l/t = 1/2πσ ≡ α = const. The latter can be extracted by fitting to the experimental data α 0.9 GeV −2 yielding the string tension value of σ 0.18 GeV 2 = 0.91 GeV/fm.
The fundamental question is how non-local string-like objects emerge from the local microscopic parton (quark and gluon) dynamics in QCD. For some peculiar reasons, the gluon field between a static quark and anti-quark gets "squeezed" into a narrow cylindrical domain whose transverse area is nearly independent on the interquark distance -the main effect of the magnetic disorder phase. In a color electric flux tube picture, the energy stored in such a QCD string is proportional to the string tension σ that can be found in terms of the color electric field E a i ≡ F a 0k as an integral over the transverse area of the flux tube as [3] Such a string then wildly fluctuates in transverse directions, and the energy of such fluctuations tends to grow with the distance between the static sources. At some critical distance, the strong fluctuations destabilise the flux tube making the longer strings less energetically favourable than the shorter ones. So, instead of indefinitely (and linearly) rising energy stored in a flux tube with its length, one encounters a string breaking effect realised due to the presence of quarks in QCD or, in a general YM theory, matter fields in fundamental representation of the gauge group. Let us elaborate on this point in some more details in what follows.

Color confinement and Higgs-confinement complementarity
A traditional and rather generic question one may ask here is what we actually mean by confinement in a gauge theory with and without matter fields that transform in the fundamental representation of the gauge group. As was discussed above, in pure non-abelian gauge theories without dynamical matter fields, the existing attempts to prove confinement consist in demonstrating the area-law dependence of W(C), or equivalently, in showing linear dependence of the static quark potential at large separations 2 . As we will elaborate in more formal details below, confinement in 2 One should make a side remark here: considering static charges in fundamental representation, with a non-zero coupling to the gauge field in the action, automatically implies that the theory is not a pure non-abelian gauge theory. Obviously, a pure gauge theory feature neither "static" nor "dynamic" quark fields, and moreover as such the latter are not distinguished by a pure YM theory is associated with an unbroken center symmetry. Thus, the non-perturbative vacuum of QCD or, in general, a non-abelian gauge theory in the range of length-scales where the static potential satisfies a linearly-rising behavior is considered to be in a confined phase.
In the presence of dynamical quarks in the theory, there would actually be no a linear static potential between heavy test quarks at asymptotically large R. Indeed, if one attempts to pull them apart one eventually observes pair creation (out of the vacuum), thus, ending up with formation of mesons at very large distances. In this picture, such a dynamical quark-antiquark pair creation occurs at the ends of the two shorter strings at the breaking point of the larger one such that the color charge of the static charges gets effectively screened off. Such a string breaking or fragmentation phenomenon in QCD causes flattening out of the static quark potential at large distances in consistency with the Regge trajectories of QCD and with vast phenomenology of particle physics processes with hadronic final states. Such a picture has become the cornerstone of the hadronisation modelling when long strings loose their stability and decay into shorter strings yielding the spree of hadrons measurable by experiments at long distances. As we will discuss more later on, no exact center symmetry can be found in such a theory since it generically gets broken by the presence of matter in fundamental gauge-group representation. There are reasons to expect a finite range in intermediate distances where the potential could be seen as approximately linear and hence string-like. So even as confinement is unquestionably useful way of thinking about long-range physics of QCD, it is by far a more complex phenomenon than an assumption about an asymptotically linear static potential associated with unbroken center symmetry.
The phenomenological reality is that coloured quarks and anti-quarks at long distances are always bind together into composite states -mesons and baryons -and do not exist as isolated color charges. This is realised in an effective string-based hadronisation picture that is proven to work very well phenomenologically in a variety of high-energy scattering processes with hadron final states (see below). The corresponding dynamics has been studied in lattice gauge theory simulations in the strong-coupling regime when matter fields are present in the action [105][106][107]. The resulting hadrons are automatically color-neutral and are the true asymptotic states of QCD, not the colored quarks and gluons. Hence, sometimes QCD confinement is naively identified with color confinement (known also as C-confinement) due to the color charge being effectively screened away at large distances by dynamical matter fields such that the colored partons may only propagate at short distances. However, one must be a little more careful with such an identification. If color confinement were the only property of the confining phase, than typical Higgs theories (such as the weak interactions' theory in the SM) should also be considered as confining [3] although they do not feature such phenomena as flux tube formation and Regge trajectories [108,109]. This is why "true confinement" appears to be a more complex phenomenon and, in addition to C-confinement, it should also be connected to other distinct properties of the quantum ground state such as magnetic disorder associated with an unbroken global symmetry [3]. It does appear indeed rather obvious that C-confinement always accompanies the magnetic disorder phase while the opposite may not necessarily be always true [110].
Indeed, consider an even simpler SU(2)-invariant gauge-Higgs theory [111], with Yukawa-type interaction term that can be straightforwardly deduced from Eq. (22). Here, the confinement regime is reproduced for small β, γ 1 characterised by the linear rise of the static potential followed by its flattening at large separations due to string breaking. So, this regime is very similar to the long-range the action unless the static ones are made very much heavier than the dynamic ones. So any statements about the linear static potential in pure non-abelian gauge theories should be taken with reservations and does make sense only when taking a limit of heavy (static) matter fields that can be effectively integrated out in the corresponding path integral of the theory. However, the latter procedure, formally, eliminates such heavy charges from asymptotic states of the resulting EFT entirely making it impossible to use them as probes for vacuum dynamics and hence confinement in pure gauge theories. So, in practice, one does not eliminate them from the asymptotic states of a gauge theory but rather retains them as heavy sources but with a finite mass. dynamics of real QCD. However, at large values of β, γ 1 one enters the Higgs regime characterised by the presence of massive vector bosons, analogues to those in the EW theory. This is the so-called massive phase characterised by a Yukawa-type potential for T E R corresponding to a perimeter-law for a generic large planar loop C, W(C) ∼ exp[−V 0 P(C)] with R 1/m. In fact, in both confinement and Higgs (massive) regimes the color field is not detectable far from its source. Indeed, while in the confinement regime there are only color-singlets in the physical spectrum of this theory, in the Higgs regime the gauge forces are the short-range ones, such that one charge screening mechanism transforms into another as the couplings change. This is due to the fact the gauge-invariant operators in SU(2) theory that create color-singlet states in the confinement domain are also responsible for creation of massive vector bosons in the Higgs domain (for a first discussion on role of the EW theory operators for generation of particle spectra, see e.g. Ref. [108]), and those states evolve into each other with varying the model parameters. Whether this happens continuously or via a first-order phase transition is a subject of ongoing research in the literature, to be discussed below.
Referring to the EW theory as a particularly important example one should be also very careful about what one actually means by the Higgs phase and the associated Higgs mechanism. Conventionally, the Higgs phase is described in terms of a Mexican-hat shape potential emerging due to formation of classical scalar fields' (Higgs) condensates in a weakly coupled regime and, as a cause, leading to spontaneous breaking of a given symmetry. While the gauge symmetry is manifest at the Lagrangian level, due to its spontaneous breakdown by means of the Higgs condensate, it is not a symmetry of solutions of the corresponding equations of motion. Note, however, that it is meaningless to talk about spontaneous breaking of a gauge symmetry without specifying a certain gauge-fixing condition. Indeed, the Higgs vacuum VEV depends on the gauge choice that we make in practical calculations and can be fixed to any value by an appropriate choice of the gauge while the actual physical observables and physical states must be gauge-invariant and do not depend on this choice. The gauge symmetry SSB phase cannot be regarded as a true physical system provided that the gauge symmetries are redundancies of description and cannot actually break spontaneously. The latter is the statement of the so-called Elitzur's theorem [112]. Indeed, according to this theorem, a local gauge symmetry, in variance to less powerful global symmetries, can not break spontaneously such that VEVs of any gauge-noninvariant observables must be zero.
In general, in a gauge theory with fundamental-representation matter fields such as a gauge-Higgs theory, for instance, one typically does not expect to physically identify a local order parameter which would distinguish between the Higgs and confinement phases as qualitative descriptions of the corresponding field configurations. If there is no gauge-invariant way to distinguish between these regimes than it would be justified to attribute them to a single phase as mentioned earlier. A discussion of this issue known as Higgs-confinement complementarity goes back to as early as late 70'es and early 80'es. In Refs. [109,113,114], by varying parameters in relatively simple lattice gauge-Higgs theories with a global symmetry, analyticity over a set of observables has been rigorously proven when going from a confining regime in the phase diagram to a regime characteristic for the Higgs phase. Although at certain large values of β such a phase boundary emerges (see e.g. Ref. [115]), one can find an analyticity line continuously connecting any two points in the parameter space except γ = 0 3 . In other words, in those models where this is true there would indeed be no thermodynamical phase transitions (or phase boundaries) along this path that separate the two regimes suggesting a possible existence of a single, massive phase all along the phase diagram (see Ref. [3] for a more elaborate discussion). Can this statement be applied only for some specific models, or it is always true?
This important result, first obtained in specific models, was then conjectured by some of the authors into a kind of "folk theorem" (known also as the Fradkin-Shenker-Banks-Rabinovici theorem) stating that the corresponding conclusion is expected to be always correct. Namely, if there is no local order parameter distinguishing different symmetry realisations one should probably expect continuity of phases. There are many examples where such a continuity has indeed been confirmed in simulations such as in transition from low-to high-temperature QCD when turning from physics of dilute gas of hadronic resonances to the physics of quark-gluon plasma (at low µ B ). Indeed, in Euclidean description of real QCD there are certain reasons to believe that there is no thermodynamic phase transition that separates these two regimes. However, as will be discussed below, the analyticity conjecture may not actually be always true. As was argued in Ref. [58] considering a discontinuity in a non-local order parameter, the Fradkin-Shenker-Banks-Rabinovici theorem does not apply to models where a global symmetry is broken in the same way in both the Higgs and confinement regimes, i.e. where the Higgs fields are charged under global symmetries.
In fact, already in the string-breaking picture of hadronisation, by construction, the gluon vector-potential cannot retain its analyticity and is inherently discontinuous in the effective string-length (or string-time) scale as the string breaks apart and no gluon field is expected to retain between the daughter strings. Whether or not the observables still remain analytic upon such a string breaking is one of the big questions for confinement models. One interesting example of the analyticity breakdown is associated with the notion of "dense QCD" or QCD at large baryon chemical potentials in the phase with broken U(1) B . We will elaborate on this aspect in the end of this review.

String hadronisation and the Lund model
One of the existing successful realisations of the string hadronisation picture is the so-called the Lund string fragmentation model [116] implemented in Monte-Carlo event generators widely-used in phenomenology of particle physics such as Pythia [117,118]. It realises the basic picture of linear confinement described above, where a flux tube is stretched between the color-charged endpoints of the back-to-back qq system that is characterised by the string tension σ 1 GeV/fm and the transverse size close to that of the proton, r p 0.7 fm. In the simplest formulation of the hadronisation model, the quarks at the endpoints are assumed to be massless and to have zero transverse momenta. As the energy transfers between the endpoint quarks and the flux tube, they move along the light cone experiencing the "yo-yo"-type oscillations. As the quarks move apart and pair-creation of dynamical qq pairs is enabled, there is non-zeroth probability for the initial "quark-string-antiquark" system to break up into smaller strings. For a simple illustration of this phenomenon, see Fig. 4. Ordering the newly produced pairs as q iqi , with i = 1, . . . , n − 1, into a chain along the string, depending on the initial energy of q andq one eventually ends up with production of a set of n mesons, {qq 1 , q 1q2 , . . . , q n−2qn−1 , q n−1q } moving along x axis of the initial string. The q iqi production vertices with coordinates (t i , x i ) have a spacelike separation, with no unique time-ordering, satisfying the constraint that the produced ith meson must be on its mass shell, In a more elaborate formulation, quarks have mass m q while the color string wildly fluctuates not only in longitudinal but also in transverse directions, and the amplitude of those fluctuations tend to grow with the string length and may eventually destabilise the system causing the string to break-up. The transverse momenta p ⊥ of the (anti)quarks are then naturally incorporated by giving q andq opposite kicks in the transverse plane, with the mean square p 2 ⊥ = σ/π ≡ κ 2 (0.25 GeV) 2 , such that the produced meson receives p 2 ⊥had = 2κ 2 . The virtual (anti)quarks tunnel over a distance m ⊥ /σ, with m ⊥ = m 2 q + p 2 ⊥ the transverse quark mass, before they become on-shell, dE dz = dp z dz = dE dt = dp z dt = κ so energy-momentum quantities can be read off from space-time ones Motion of quarks and antiquarks in a qq system: z t q q gives simple but powerful picture of hadron production ⊥ /σ). In the framework of Lund model, a consistent selection of the produced DoFs is performed according to the probability distribution [116], implying an equilibrium distribution of the production vertices on the string where Γ = σ 2 (t 2 − x 2 ), a, b are free parameters, and z is the light-cone momentum fraction carried away by a produced meson. The remaining (1 − z)-part of the momentum is kept by the string and is then redistributed among other mesons in its subsequent fragmentation. Even though the hadron masses do not enter this approach directly, a good description of the produced particle spectra can be reached with only a few free parameters. More complicated qqgg . . . topologies can be introduced considering a gluon as a state with separate color and anticolor indices, well justified in the large-N c limit [119]. The string gets then stretched between q andq as usual while each of the gluons attach at intermediate points along the string respecting the color flow that goes in and out of each gluon. Notably, the fragmentation procedure of such a string does not require any extra free parameters [120]. The fact that there is no string that connects q andq directly in this case leads to asymmetries in the produced particle spectra in consistency with experimental observations [121]. At last, baryon production can be conceptually tackled by enabling a diquark-antidiquark breaking e.g. via sequential qq production stages (for more details on this mechanism, see e.g. Refs. [122,123]).

Gauge symmetry remnants and confinement criteria
Due to the Elitzur's theorem [112] described above the phases of a gauge theory cannot be distinguished by means of the breaking of any local gauge symmetry. Thus, there must be an additional, global symmetry whose breaking enables us to identify those phases, at least, when a local order parameter is concerned. In the Ising model, the role of such a global symmetry is played by the Z 2 symmetry as we have noticed earlier. Fixing a covariant gauge, in general, does not eliminate the gauge freedom entirely, but leaves certain remnant (both dependent and independent on spacetime coordinates) symmetries that can in principle get spontaneously broken since the Elitzur's theorem does not apply to those.
One of the examples of a possible confinement criterion known as the Kugo-Ojima condition [124,125] states that the full residual gauge symmetry in Landau gauge ∂ µ A a µ = 0 must remain unbroken in order to ensure that the expectation value of color charge operator ψ|Q a |ψ vanishes in any physical state ψ. The spacetime-dependent (but global) part of such a full residual gauge symmetry w.r.t. gauge transformations A µ → GA µ G † in Landau gauge is known to take the following form [126,127] where in terms of a finite number of arbitrary parameters a µ , and the SU(2) gauge coupling constant g. Besides, for confinement to hold yet another, spacetime-independent, part of the full residual gauge symmetry is required to be unbroken in addition to that in Eq. (48). An analogical criterion of confinement has also been formulated in Coulomb gauge [128,129].  Thus, according to the Kugo-Ojima and Coulomb confinement criteria the phase boundary between the confining and de-confining regimes of a gauge theory is associated with the boundary between the unbroken and broken full (x-dependent and independent) remnants of the gauge symmetry in Landau and Coulomb gauges, respectively. However, a problem highlighted by lattice simulations and demonstrated in Fig. 5 is that these criteria predict transitions between confinement and deconfinement phases where actually no such transitions appear in the exact numerical analysis [131]. Different remnant symmetries emergent in different gauges break at different values of the couplings, so the resulting phase boundary is in fact gauge-dependent and might indeed emerge even when there is no actual change in the physical state of the system [131]. In order to distinguish a confining state from a non-confining one, one should instead come up with a gauge-invariant criterion whose violation would indicate a true boundary between the magnetic order and disorder states that is the same in any gauge. As was discussed earlier, there is no such criterion in a gauge-Higgs theory. This might indicate that there is no such gauge-invariant separation of phases that can be attributed to a spontaneous breaking of a given symmetry and the system is in the massive phase characterised by the string-breaking effects and a perimeter-law behaviour of Wilson loop VEVs [3].
There are compelling reasons to believe that the same picture is realised in QCD at very large separations, supported also by lattice simulations. In the gauge-Higgs theory, only in the limit of Higgs decoupling, γ = 0, the state of magnetic disorder emerges indicated by the area-law falloff of large Wilson loops at arbitrary large spacetime separations. The same occurs in the infinite quark mass limit in QCD such that it takes an infinite amount of energy in order to put an infinitely heavy quark-antiquark pair on its mass-shell from the vacuum such that the area-law persists to arbitrarily large string lengths.

Center symmetry
So, when one talks about the true (gauge-invariant) separation of phases, one implies a strong first-order (non-analytic) phase transition between the magnetic order (massive) and disorder states that exists at a well-defined (unique!) combination of model parameters in any gauge. Such a non-analytic behaviour is associated with a spontaneous breaking of a certain symmetry and, to comply with the Elitzur's theorem [112], such a symmetry must be global. This type of a symmetry exists and is called the center symmetry -a specific subgroup of a given gauge symmetry group which is defined as a subset of the gauge group elements that commutes with all the elements of the gauge group. For instance, the center of the SU(N) gauge symmetry group is its Z N subgroup {exp 2πin/N1 N }, with n = 0, . . . , N − 1.
Each of an infinite number of SU(N) representations can be separated into N possible subsets or N-alities depending on the corresponding representation of Z N (there are only N of those). Hence, each SU(N) representation is characterised by N-ality k that is found as the number of boxes in the associated Young tableau mod N. In other words, N-ality reflects how a given representation transforms under the center symmetry subgroup of the gauge group. For instance, if for a matrix representation M[g] of an SU(N) group element g, M[zg] = z k M[g] for a center Z N element z, one says that g belongs to a representation of N-ality k (for a more detailed pedagogical discussion, see e.g. Ref. [3]). In the lattice formulation, one could show that the action (19) of a pure gauge theory is invariant under the timelike link transformation on a fixed time slice t = t 0 . This transformation is a particular case of the singular gauge transformation defined on a time-periodic lattice with a period L t as where G( x, t) is a periodic function up to a center symmetry transformation, i.e.
G( x, L t + 1) = z * G( x, 1) , that also leaves Wilson loops invariant on the lattice. Such a transformation corresponds to an "almost" gauge transformation in the continuum limit, where the second term is dropped for t = L t and for µ = 0 when it turns into delta-function. Matter fields in the fundamental representation of the gauge group SU(N), or any other fields with N-ality k = 0, break the center symmetry Z N explicitly if they are not decoupled from the theory -like the Higgs field for a non-zero coupling γ in the example discussed above or the quark sector of real QCD (with k = 1), with finite quark masses. Such a breaking, which is also a necessary ingredient of the string hadronisation model (see above), causes the static potential to flatten out instead of growing linearly at asymptotically large distances as the matter fields are, in fact, responsible for the string breaking phenomenon. Gluons or other particles in the adjoint representation having N-ality k = 0 do not break the center symmetry so they cannot screen the color charge of a static source if the latter has a non-zero N-ality. A well-known exception is the G 2 gauge symmetry which has a trivial center subgroup, with a single unit element only, such that the gluons can bind to any source producing a color-singlet state.
An important criterion of confinement is thus associated with unbroken center symmetry in a pure YM theory implying an asymptotically and infinitely rising static quark potential and signalling the area-law falloff of large Wilson line VEVs and hence the presence of the magnetic disorder state. The center symmetry can also be spontaneously broken by thermal effects i.e. at high temperatures, in pure YM theories causing the same effect of flattening out the static potential asymptotically as that of the matter fields. Other possible sources of the center symmetry breaking should also be considered, in order to reconstruct a full picture of phases in the underlined gauge theory.

Polyakov loop
Consider a finite (in space) lattice which is periodic in time. Such a lattice is used, in particular, in quantum statistical mechanics at finite temperatures T where the partition function reads In the continuum limit of a field theory and in Euclidean time T E , the latter generalises to a path integral where the periodic boundary condition in time φ( x, 0) = φ( x, β T ) is imposed through an implicit delta-function. Upon lattice regularisation, the temperature is related to the lattice period in time L t as T = 1/(L t a), with a being the lattice spacing as usual, and hence β T = L t a is the total time extension of the lattice. While neither the gauge-field action nor Wilson loops are affected by the Z N center symmetry transformation (50), the trace of the following holonomy winding in time around the periodic-time lattice, known as the Polyakov loop [132], is Z N non-invariant, i.e. it transforms as P → zP. In the continuum limit, one can represent the Polyakov loop holonomy as follows in terms of an SU(N) matrix S(x) that diagonalises P ( x).
One can show that in the case of P ( x) being a center element, all µ j are equal and such a holonomy determines the finite-temperature classical instanton solutions known from Refs. [133,134].
In fact, in center-projected configurations that will be discussed below the Polyakov loop holonomies P ( x) are the only center elements. In general, the Polyakov loop is non-trivially charged under Z N meaning that its expectation value plays a role of an order parameter for spontaneous breaking of the center symmetry. Hence, the Polyakov loop is yet another important characteristics of the confined (magnetic disorder) phase of the gauge theory, and vacuum fluctuations of the gauge field are responsible for formation of this phase and in some ways are associated with the center symmetry.
Indeed, a difference between free energies of two states, one containing a single isolated (heavy) static charge q and the other one defined in a pure gauge theory is as follows which is obtained by integrating out the massive quark field (in the m → ∞ limit) in the path integral for Z q over the period of the lattice 0 ≤ T E < β T . Indeed, F q for a single quark q would be infinite if P( x) → 0, i.e. in the case of unbroken center symmetry. At high temperatures (small β T ), the center symmetry is in general spontaneously broken such that the isolated charges are described by finite-energy states (deconfining phase). A magnetic disorder-to-order phase transition associated with thermal breaking of the center symmetry is expected to occur at a critical temperature [135].

t'Hooft loop and center vortices
The singular gauge transformation in the continuum limit (53), unlike the ordinary center symmetry transformation, leaves the action non-invariant. As a result of such a transform, a singular loop of magnetic flux, the so-called thin center vortex, is being created. For instance, as was mentioned earlier the holonomy for a closed spacelike loop C in U(1) gauge theory is given in terms of the magnetic flux Φ B through the loop. For a loop winding around a solenoid oriented along the z-axis, it is possible that Φ B = 0 even for a zeroth magnetic field along the closed loop which can be obtained as a result of a singular gauge transformation applied to A µ = 0 with a discontinuous G(x). If in cylindrical coordinates {r, θ, z, t}, the corresponding transformation function G has a discontinuity in θ for r > 0, then where exp(±ieΦ B ) is an element of U(1) group, the sign ± depends on the orientation of the loop C, such that a singular line of magnetic flux (thin vortex) is produced along the z-axis. Instead of z-axis, one could introduce yet another closed contour C topologically linked to C such that the singular gauge transformation operator G that creates a magnetic flux along C would satisfy Switching over to SU(N) YM theory, the U(1) group element that multiplies a transformation operator in Eq. (61) should be replaced by a center-group Z N element in order for such a transform to create a thin vortex (and hence to affect the action) on the (D − 2)-dimensional hypersurface only, and not on the Dirac D − 1 region that it envelops. Above, the spacelike Wilson loop C is topologically linked to the (D − 2)-dimensional thin vortex, with the corresponding linking number l. Upon quantisation of the non-abelian magnetic flux, its quanta are known in the literature as the thin center vortices, while a regularisation of the singular color-magnetic field by smearing it out in the transverse directions to the (D − 2) hypersurface leads to a vortex with finite thickness, or a thick center vortex. For a more detailed description of the vortex configurations and properties, see e.g. Ref. [3] and references therein. Consider an operator B(C) that creates a thin center vortex at a fixed time t 0 along a given loop C in a D = 3 + 1 gauge theory [136]. If C and another closed loop C are topologically linked (with l = 1) in a three-dimensional surface, then is valid. In this case, the operator B(C) is known as the t'Hooft loop. As was demonstrated in Ref. [136], the VEV of a Wilson loop W(C) ≡ U(C) and a t'Hooft loop B(C) may satisfy either a perimeter-law or an area-law falloffs, but not simultaneously. Indeed, the confined (magnetic disorder) phase corresponding to an unbroken center symmetry is realised when while the opposite case, implies a spontaneously broken center symmetry (magnetically-ordered phase). Indeed, the Wilson and t'Hooft loop operators can be considered dual to each other as the first one creates a closed loop of color-electric flux, while the second one creates a closed loop of color-magnetic flux (thin center vortex), at a fixed time t in both cases. One could introduce a vortex on a finite lattice in D = 4 by replacing U(p ) → ξU(p ) for a given plaquette p in the SU(N) gauge-field action [137] which can be viewed as a change in periodic boundary conditions, also referred to as twisted boundary conditions. Such a change creates a thick center vortex on the lattice parallel to (z − t)-plane, satisfying the ordinary periodic boundary conditions in z, t coordinates.
In the simplest case of SU(2) gauge symmetry, the (magnetic) free energy of a Z 2 -center vortex F m can be found as in terms of the partition functions with ordinary and twisted boundary conditions, Z + and Z − , respectively, while the free energy of closed color-electric flux F e is It was shown in Ref. [138] that the VEV of a rectangular Wilson loop C with area A(C) is bounded from above as A sufficient condition for the existence of a magnetic-disorder phase, and hence confinement, in terms of the behavior of the magnetic vortex free energy then reads i.e. it falls off exponentially at large L x L y area, such exp(−F e ) F m 1. Indeed, the latter limit, together with Eq. (70), implies an area-law upper bound for a large Wilson loop and, hence, the asymptotic string tension. In Ref. [139] it has been pointed out that quark confinement emerges from a vortex condensate supported by the mass gap.

Fundamental properties of the string tension
One of the fundamental characteristics of confinement is an non-vanishing asymptotic string tension or, equivalently, the asymptotic linearity of the static potential [3,4]. As was proven in Ref. [140], the potential is always convex and is saturated by a straight line from above. At not too large distances, the string tension for a quark in a given representation r of the gauge group interacting with an antiquark can be approximated as This is the property known as the Casimir scaling which is strictly valid in the large-N limit. Here, σ F is the string tension for the defining (fundamental) representation. Such a scaling can be proven in a two-dimensional theory and then to a good precision can be found also in 4D by means of the dimensional reduction [141], supported also by numerical simulations [142]. For a more recent analysis of the Casimir scaling in D = 2 + 1 SU(N) theory in the vortex picture, see Ref. [143]. Asymptotically at very large distances, the Casimir scaling does not hold (apart from N = 2 and large-N cases), and can be effective at intermediate distances only. The dimensional reduction is a specific (approximate) property of quantum state of the theory Ψ 0 [A] emergent at large length scales. According to this property a calculation of the VEV of a large Wilson loop W(R, T) in fundamental representation in a D = 4 gauge theory can be sequentially reduced to that in D = 3 theory [144,145] and then down to D = 2 case [91]. In this case, where in last relation corresponds to the fact that in D = 2 the Wilson loop VEVs obey an area-law falloff. For this property to hold in the strong coupling limit, the vacuum functional should take the same form in D = 2, 3, 4 at large length-scales: Note, this form can not be correct at short distances in PT, so it should be regarded as approximate and generically valid in the non-perturbative regime only. It is also not correct for Wilson loops in the adjoint representation which follow a perimeter-law due to the color screening effect. An elaborate form for the vacuum functional that matches both the dimensional reduction form and the correct free-field limit has been proposed in Refs. [146] predicting the glueball mass spectrum in D = 2 + 1 in consistency with the lattice calculations. For other proposals, see e.g. Refs. [147][148][149][150][151].
Another fundamental property of the string tension, presumably closely related to confinement, is the observation that the string tension depends only on N-ality of the gauge group representation. For static quarks in the adjoint representation, for instance, gluons screen their charges at large distances causing the string to break at separations R satisfying 2E < σ A R, where E is the gluonic energy of the produced "gluelump" state, σ A = C A /C F σ F is the string tension in the adjoint representation, valid at intermediate distances. For numerical studies of the adjoint string tensions, see e.g. Ref. [152]. While the precise form of the N-ality dependence is not known there are several models widely used in the literature. Among them, for instance, the "Casimir scaling" proposal assumes that the string tension for the lowest dimensional representation (k-string tension), behaves asymptotically for the SU(N) gauge theory as If the true confinement phenomenon implies the formation of an electric flux tube in the form of a quantum Nambu-like string, typical predictions of the string model such as subleading deviations from linearity of the potential as well as the spectrum of string excitations should find their evidence in a first principle analysis of confining gauge theories. In particular, one such prediction is a subleading 1/R correction term to the static quark potential emerging due to transverse fluctuations of the string known as the Lüscher term [153,154] such that the VEV of a large rectangular Wilson loop can be generically parameterised as where the second term in the exponent is a self-energy contribution that diverges in the continuum limit as was mentioned above. On the lattice, one may extract the asymptotic string tension σ as the following ratio computed at large loop areas known as the Creutz ratio.
Another property of the Nambu string is that the cross section area of the string grows logarithmically with the quark separation, the effect knows as roughening [155,156]. An agreement with Nambu string model predictions was found earlier in the analysis of closed string excitations in the D = 2 + 1 SU(N) gauge theory in Ref. [157].

Center vortex mechanism of confinement
The center vortex mechanism of confinement is strongly supported by the fact that the static potential slope depends only on N-ality, while N-ality zero (or adjoint) string tensions vanish at asymptotically large distances. Also, when adopting a picture of pair creation of particles out of the vacuum at a certain distance causing the string to break, one implies a microscopic perturbative language of particle states in a particular configuration. While an extrapolation of perturbative particle states towards large distances may not necessarily work out well in confining theories, an effective particle picture of string breaking is still considered to adequately reflect the reality, at least qualitatively. In proper path integral computations, one sums over all possible field configurations that should provide the same result for the gauge-invariant observables (such as Wilson loop VEVs) as the phenomenologically successful effective particle picture of the string breaking. Ultimately, one would like to find out how the vacuum field fluctuations induce N-ality dependence of the asymptotic string tension and describe colour screening of the static sources [3].
While instantons [158] are saddle points of the classical gauge-field action, vortices are interpreted to be saddle points of the effective one-loop action [159,160] that incorporates the vacuum polarisation effects, and hence have a pronounced fundamental meaning (see also Ref. [4]). Fluctuations of center vortices that can be identified as solitonic objects in typical field configurations are known to give rise to an area law of Wilson loops. A remarkable property is that Wilson loops in different representations but with the same N-ality get the same contributions from center vortices, while loops of N-ality zero are not affected. This follows from the simple fact that the creation of a vortex linked to the loop C affects the loop holonomy of a given N-ality k as U(C) → zU(C) and its VEV as W r (C) → z k W r (C) for z from the center group Z N of SU(N).
In a more generic case, consider a set of vortices linked to a given loop C, with linking numbers l 1,2,3,... , having the center elements z 1,2,3,... . Then, creation of this set modifies the Wilson loop VEV as W r (C) → Z k (C)W r (C), where Z(C) = z l 1 1 z l 2 2 z l 3 3 . . . . In the vortex picture of confinement [136,159,161,162] (see also Ref. [3] and references therein) the gauge-field vacuum configuration is considered to be a set of vortices superimposed on a non-confining configuration. Then random fluctuations in number of vortices in the system as well as in their linking numbers to a given Wilson loop C induces the area law dependence of the corresponding Wilson loop VEV. The loop holonomy can be represented in a factorised form U(C) = Z(C)u(C), where u(C) is a contribution from a non-confining background, and Z(C) ⊂ Z N is a center-valued holonomy. Then, the vortex mechanism implies factorisation of the Wilson loop VEV A detailed proof relies on a weak correlation between Z(C) and U(C), as well as between Z(C 1 ) and Z(C 2 ) for any large loops C, C 1,2 , and can be found for instance in Ref. [3]. It manifestly demonstrates that the string tension computed for smaller loops is the same as that for the larger ones provided that the above assumptions hold and Z(C i ) experience independent fluctuations. Numerical estimates [163] suggest that the thickness of the vortex is close to one fermi so, in principle, the Wilson loops with an extension below this scale may get affected. As was demonstrated in Refs. [164,165], such a vortex thickness plays an important role for generating the Casimir scaling at intermediate distances. At large distances dominated by large Wilson loops the N-ality dependence of the linear static potential is reproduced as expected. From this point of view, vortices are non-local objects that represent specific field configurations that lead to an asymptotic string tension as a function of N-ality.
The link configurations U µ (x) = g(x)z µ (x)g −1 (x +μ) that produce Z(C) holonomies can be transformed into the link configurations z µ (x) of Z N lattice gauge theory responsible for confinement by means of a specific SU(N) gauge transformation g(x). The thin vortices then have a meaning of excitations of the center-group Z N lattice gauge theory. The original link variables U µ (x) get separated into a product of center elements z µ (x) and the link variables of the non-confining background V µ (x) by the g(x) transform The main aim of the vortex mechanism of confinement is to find a specific g(x) for a given non-confining background V µ (x) typically assumed to be a small fluctuation about the unity. Locations of center vortices then can be extracted from z µ (x) after the above factorisation U → zV has been performed [166]. One such g(x) transforms the DoFs into a specific gauge known as the direct maximal center gauge where the deviation of the links in the adjoint representation from the identity matrix is minimal, or where the quantity with the adjoint link U A µ (x), is maximal. Locations of center vortices then can be extracted in a dedicated Monte-Carlo procedure from identified center elements z µ (x) once the center mapping (projection) U µ (x) → z µ (x) has been performed. If the product Z(p) of z µ (x) on the projected Z N lattice around a plaquette p satisfies Z(p) = 1, a thin vortex (or P-vortex) is then located on that plaquette. The vortex picture of confinement then reduces to a consideration of P-vortices as random surfaces percolating through the spacetime volume. Uncorrelated piercings by the P-vortices on a given planar surface correspond to uncorrelated large center-projected loops. The numerical procedures, however, may fix the projected lattice to only one out of a large amount of local maxima of the gauge-fixing functional K known as the Gribov copies [167], not straight to its global maximum, which is considered to be a problem in several widely used center-gauge fixing approaches.
The problem of Gribov copies is one of the main obstacles for a consistent treatment of the confinement problem. Considering a set of gauge-equivalent configurations of the gauge field known as a gauge orbit and imposing a gauge-fixing condition as a certain hypersurface in a space of gauge field configurations, the Gribov copies can be visualised as many possible intersections of the gauge orbit with the gauge-fixing hypersurface. Summing over all contributions from Gribov copies in a path integral, the latter may actually vanish since those contributions come with opposite signs and may mutually eliminate each other in a given observable. This is the statement of the Neuberger's theorem [168] rendering BRST quantization not well defined in the non-perturbative regime of a gauge theory (see a detailed discussion in Ref. [3]). One possibility is to restrict the functional integral to a subspace of gauge configurations with positive Faddeev-Popov determinant, the so-called Gribov region, and its boundary containing also the lowest non-trivial eigenmode with zeroth eigenvalue is called the Gribov horizon. An instructive example of such a hypersurface restricted to the Gribov region is the Landau gauge fixing condition that minimises the functional such that the corresponding Gribov region consists of all possible minima of R = R[A] for a given gauge orbit. However, various gauge orbits might cross the Gribov region a different number of times leading to different weights assigned to different gauge orbits. A proposal to consider only unique global minima of R[A] functional for each gauge orbit [169] may be very difficult to realise in practical calculations. Also, there is no any reason to believe a particular Gribov copy with the global minimum for R[A] is more physical than other local minima. Lattice procedures, in general, assume that a particular choice of a Gribov copy would not make a big difference on the numerical results.
In order to establish a direct connection between the existence of P-vortices and magnetically disordered phase, following the reasoning of Ref. [3], let us first consider whether the center-projected Z µ (x) link variables (extracted, for instance, in a maximal center gauge) are responsible for confinement. For this purpose, it is instructive to consider the VEV of the rectangular R × T E Wilson loop, W(R, T E ), defined in Eq. (77). If such a loop is constructed from Z µ (x) links on a center-projected lattice, the corresponding Creutz ratio (78) appears to converge much faster to σ than for the unprojected Wilson loop VEV. Already at R = 2 the static potential becomes linear -the property of the so-called precocious linearity. The fact that the asymptotic string tensions extracted from center-projected and unprojected Wilson loop VEVs at large R are the same is known as the center dominance. There is also an excellent agreement of the Creutz ratios on the center-projected lattice with the well known predictions of the asymptotic freedom for large β (small gauge couplings).
A slow convergence of the Creutz ratio (78) to the string tension at large R in the unprojected case means that here we deal with thick vortices linked to large Wilson loops. The center projection effectively shrinks the thickness of the vortices down to a single lattice spacing, so the linking appears to be relevant already for small center-projected loops. Indeed, as was deduced earlier, P-vortex piercings are totally uncorrelated on a planar surface causing the linearity of the potential already at small distances. One naturally wishes to establish that each thin P-vortex in the projected configurations matches a thick center vortex in unprojected lattice in order to prove that the P-vortices do not carry artefacts of the gauge fixing procedure and indeed are responsible for the underlined physics of magnetic disorder (and hence confinement).
As thoroughly described in Ref. [3], one way of proving the relevant correlation of P-vortices with gauge-invariant observables (unprojected Wilson loops) is to compute a so-called vortex-limited Wilson loop VEV defined as an expectation value of an ordinary unprojected loop holonomy W n (C) but taken in the ensemble of configurations where the minimal surface area of the loop is pierced by n P-vortices. Then, considering for simplicity the SU(2) theory, if the ratios asymptotically behave as W n (C)/W 0 (C) → (−1) n provided that Z(C) = (−1) n (−1 per each vortex piercing) then the procedure of finding thin P-vortices on the center-projected lattice effectively locates thick center vortices on the unprojected lattice. This indeed has been confirmed by lattice simulations, see e.g. Ref. [170].
Another test proposed in Refs. [171] suggests to insert a thin vortex found by the center projection operation into a thick vertex on the unprojected lattice and then to check if their disordering effects, due to center dominance, cancel out asymptotically at large distances. Indeed, an explicit calculation shows that this procedure eliminates the string tension and hence the disorder effect. It was also checked in Ref. [172] that P-vortex density is independent on the lattice spacing in the continuum limit as expected for physical objects. An additional observation of Ref. [173] revealed that the continuum action density appears to be singular at the location of P-vortices which, together with their constant density, signals an intricate cancellation between action and entropy at a surface of infinite action associated with a vortex.
As was discussed above, at finite temperatures T in a time-periodic lattice the Polyakov loop VEVs determine the quark free energy F q . In SU(2) gauge theory, at T > T c = 220 MeV a deconfinement transition occurs when F q becomes finite and the static quark potential goes flat. However, even at large T > T c space-like Wilson loops retain their area-law falloff such that vacuum fluctuations inherit some of the key properties of the confined phase.
This observation fits well with the center-vortex mechanism of confinement [3]. At low T, due to uncorrelated piercings of the minimal loop areas, one finds P( x) = 0 and an exponential falloff of the Polyakov loop correlators for large interquark separation P( x)P( x + R) ∼ exp[−σ(T)L t R], with σ(T) the T-dependent string tension of a flux tube stretched between q andq. Since the vortices running in spacelike directions have a finite diameter, as temperature rises, they get squeezed by the reduced finite lattice extension in time L t until they effectively stop percolating eliminating the exponential falloff of Polyakov loop correlator and hence P( x) is no longer zero [163,174]. The asymptotic behaviour of the space-like Wilson loop, however, is determined by the piercings of center vortices oriented in periodic time (i.e. running in timelike directions), and their cross-section is not limited by a small extension in the time direction at large T. Thus, the corresponding P-vortices keep percolating on a time slice in the spacial directions such that the exponential falloff of spacelike Wilson loops remains unaffected in the deconfined regime [175,176].
As we already discussed above, the center symmetry turns out to be explicitly broken by the dynamical fields in the fundamental representation. The center dominance in the confinement region in SU(2) gauge-Higgs theory has been tested in Ref. [177]. In a region where the screening effects by the matter fields become important the center vortices do not disappear but somehow rearrange themselves in order to allow for asymptotically vanishing string tension but still generating a linear slope in the potential at intermediate distances (no signature of linearity has been found in the Higgs region at any scale). In the presence of matter fields, Dirac volume shrinks and the vortex piercings of the Wilson loop minimal area are expected to become correlated at large distances, but to the best of our knowledge there is no full consensus on exactly how this occurs.

Chiral symmetry breaking and topological charge
The global chiral symmetry of QCD light u, d quark sector SU(N f ) R × SU(N f ) L (with the number of flavors, say, N f = 2) is broken spontaneously by the order parameter known as the quark (or chiral) condensate qq = 0. In addition, it is also broken explicitly by light current quark mass turning the Goldstone bosons, the pions, into massive pseudo-Goldstone states. Another less known mechanism based upon the linear sigma model of effective quark-meson interactions introduces yet another source of the global chiral symmetry breaking, through a linear term in σ-field proportional to the quark condensate. Such a breaking is also explicit and as such it provides an additional finite contribution to the pion mass. A symmetry breaking due to quark condensation phenomenon is often referred to as dynamical symmetry breaking and is considered a baseline for Technicolor models of EW symmetry breaking [178,179] (for a detailed review of existing concepts, see e.g. Ref. [180]).
As was discussed earlier in the case of Ising model, in order to get a nontrivial value of the order parameter one should perform two limits in a certain order -first, take volume to infinity, and then set the quark masses to zero. This procedure leads to the well-known Banks-Casher relation [181] between the chiral condensate as the trace of the quark propagator and the value of the density of the close-to-zero eigenvalues of the Dirac operator characterised by vacuum field configurations. The latter density receives no perturbative contributions, and hence the dynamical chiral symmetry breaking is an intrinsically non-perturbative phenomenon.
Provided that in real QCD with light quarks the string tension vanishes asymptotically due to colour-screening and string breaking, the chiral condensate by itself is not tied to the area-law falloff of large Wilson loop VEVs and does not even require the presence of gauge fields, in analogy to the effective Nambu-Jona-Lasinio model [182]. Naively, one might think that these observations indicate no immediate connection between the chiral symmetry breaking mechanism and the confinement phenomenon. As was emphasised in Ref. [183], the low-lying Dirac eigenmodes which are crucial for chiral symmetry breaking provide vanishingly small contributions to the string tension and to the Polyakov loop in both confined and deconfined phases. These observations provided no indication of an immediate correspondence between chiral symmetry breaking and confinement.
Interestingly enough though, the critical temperatures of chiral and deconfinement phase transitions appear to be the same or close to each other as suggested by lattice simulations, motivating a further search for possible hidden connections between the two transitions. In particular, a connection between the Polyakov loop, center symmetry, and the chiral condensate may be due to the fact that, after integrating out fermions, the chiral condensate is basically a complex expectation value of many Wilson loops, including those wrapping around compact dimensions. As was elaborated in detail in Ref. [184], the spectral properties of the Dirac operator are affected by confinement, in particular, causing the correlators of Dirac eigenvector densities to decay exponentially instead of a power law in the deconfined phase. Ultimately, one would need to establish a link between the spectral properties of the Dirac operator in the infrared regime presumably, responsible for chiral symmetry breaking with those in the ultraviolet regime tightly connected to confinement.
Remarkably, in Ref. [171,185] it was shown that the chiral condensate vanishes as soon as vortices are removed from the underlined field configurations, while the chiral condensate values are notably larger in center-projected configurations than those on the unmodified lattice. This observation shows that the center vortices are responsible not only for magnetic disorder but also determine the chiral symmetry breaking -thus, both phenomena are tightly connected [186].
It is well known that the axial symmetry U(1) A of the classical QCD action is broken by chiral anomaly at quantum level. The topological charge given by the integral of the divergence of the axial current, receives contributions from finite action configurations known as instantons [158]. Due to the Atiyah-Singer Index theorem, integer Q value has a meaning of a difference of numbers of zero modes of the Dirac operator with positive and negative chiralities. The η meson, which would have been a (pseudo-)Goldstone boson of U(1) A breaking, appears to be way too heavy phenomenologically (above 1 GeV). Its mass is found to be proportional to the topological susceptibility found in pure gauge theory in the chiral and large-N c limits, i.e.
-the relation known as the Veneziano-Witten formula [187,188]. Here, V → ∞ is a large volume, and f π is the pion decay constant. For lattice calculations of the topological susceptibility and tests of the Veneziano-Witten formula, see e.g. Refs. [189,190]. The topological susceptibility χ is characterised by the vacuum quantum-field fluctuations in a pure gauge theory, without any quark fields. Like the chiral condensate, the density of topological charge may not seem to immediately connect to the IR property of confinement, and naively one would guess that it may be determined by non-confining configurations such as instantons in the standard picture. However, as was shown in Ref. [191] a P-vortex acquires a fractional topological charge at "writhing" points, and it is possible to get a correct topological susceptability in certain vortex models [192]. Moreover, the results of Ref. [171] actually demonstrate that the topological charge tends to vanish upon vortex removal, while in Ref. [193] it was shown that χ computed from P-vortices appears to be consistent with the measurements. So the initial naive guess do appear to be wrong, and confinement plays a crucial role here as well.
Yet another, more recent, test of the vortex mechanism considering the effective quark propagator in Landau gauge in the following IR form has been performed in Refs. [186,195] (see also Refs. [3,4] for a pedagogical discussion). With an appropriate smoothing ("cooling") procedure in the SU(3) gauge theory that eliminates short-distance fluctuations, the effective mass M(k) and renormalisation Z(k) functions have been computed for the full, vortex-only and vortex-removed configurations and compared to each other. Removing the vortices causes the mass function to plummet dramatically -see Fig. 6, while the full and vortex-only results have appeared to be essentially the same, hence demonstrating a critical role of the vortices in dynamical mass generation and chiral symmetry breaking. The maxima of the action for vortex-only configurations appear to resemble those of instantons, while the number density of those objects is notably similar for the full and vortex-only configurations and by far much larger than that for the vortex-removed case. It seems likely that vortices and instantons are indeed connected in some very non-trivial way.
Remarkably, the center vortices thus appear to describe a number of fundamentally important IR phenomena in non-abelian gauge theories in a gauge-invariant way. Nevertheless, there are also weak points in the vortex mechanism of confinement that require further clarification and, in a perspective, a more complete understanding of, for instance, the Gribov copies problem and a lack of natural  Smoothing is necessary for the center-projected configurations, bec isn't suited to rough Z N configurations. Smoothing is carried out in Z N subgroup, so it is more appropriate to call the result "vortex-only" A comparison of M(p), Z(P) for the full and vortex removed co Obviously vortex removal causes a severe reduction of the mass func zero. By contrast, M(p), Z(p) for the full and vortex-only configurati 5), and M(p) is of course an effect due to chiral symmetry breaking.   [194]. explanation of the Lüscher term. A further, more complete theory of vortices should address these issues, hopefully, on a first-principle basis.
A lack of a perfect consistency of the vortex scenario with full numerical results for the SU(3) gauge theory has also emerged in the literature. For instance, center projection in SU(3) case yields 2/3rd of the asymptotic string action computed on the full lattice [196]. However, a consistency has been substantially improved by means of a certain gauge-field smoothing procedure [195], so this may not be regarded as a critical problem.
In order to make the next step in our understanding of the vortex dynamics, it may enlightening to suggest an EFT of vortices as non-local dynamical objects -fluctuating surfaces -in D dimensions where all the IR phenomena described above would emerge naturally among its key predictions. Such a theory known as the random surface model that resembles a string theory on lattice has been proposed and elaborated e.g. in Refs. [191,192,[197][198][199].
In order to build the simplest D = 4 action density of vortices in this framework, one considers an extrinsic curvature of the vortex worldsheet multiplied essentially by a single coupling while the additional Nambu-like string term proportional to the area of the vortex worldsheet appears to be redundant and can be omitted. In the SU(2) version of this model, one assigns (−1) n to the Wilson loop holonomy for the number of vortex piercings n per minimal loop area and then one averages it over an ensemble of center vortex configurations. The latter can be generated by Monte-Carlo methods for a lattice action density given by the number of cases when a single link is shared by two adjacent orthogonal vortex plaquettes.
In order to compute the topological charge density in this model, for instance, one employs a weighted stochastic procedure of introducing the monopole lines to the surface of each vortex plaquette (see Sect. 16 below for a brief description of the magnetic monopoles' scenario of confinement). The topological susceptibility appears to be insensitive to the monopole lines' density --a sign of strong predictive power of the model. Besides, the model correctly predicts the emergence of the chiral condensate at T < T c and restoration of the chiral symmetry at T > T c , with a critical temperature of the transition T c . A variation in the lattice time extension can provide a temperature dependence, and the second-order deconfinement phase transition has been found. The single dimensionless coupling and the lattice spacing a determine a wide range of long-distance non-perturbative phenomena and were fixed through a matching to the physical T c / √ σ and σ/a 2 = (440MeV) 2 , in terms of the string tension σ. Upon such a matching, the temperature-dependent values of σ, the chiral condensate and χ are shown to be in agreement with the full theory. Remarkably, in the case of SU(3) gauge theory the random surface model predicts the electric flux tubes in a form of Y-shaped string junctions for baryons (three-quark systems) [200], also in agreement with numerical results of Refs. [201,202].
An alternative EFT approach to dynamics of vortices was suggested in Ref. [162] that is based upon a gauge theory with an adjoint matter field and its gauge-invariant mass term which provides a mass for the gauge field via the Higgs mechanism. Besides the vortex solutions it also naturally reveals another type of solutions with magnetic monopoles running along the vortex sheets that are necessary to generate a topological charge.
For a more thorough discussion on the existing vortex-based scenarios, we refer the reader to Ref. [3]. Now, we turn to alternative scenarios of confinement, yet trying to connect them with the existence of vortices whenever possible.

Gribov-Zwanziger scenario, non-perturbative propagators and gluon chains
Starting from the Coulomb gauge, in Refs. [167,169,203] it has been suggested that very small eigenvalues of the Faddeev-Popov operator that are located close to the Gribov horizon contribute the most to the Coulomb potential V C (R) and could in principle enhance it to a linear form (see also Ref. [204]). This is the so-called Gribov-Zwanziger scenario of confinement. As it should be for a confined phase, numerical analysis on the lattice demonstrates the linear rise of Coulomb potential V C (R) which is basically a separation-dependent part of the interaction energy of the physical qq state defined as with the ground-state of the theory Ψ 0 , the vacuum energy E 0 and with self-energy contribution neglected at large R. However, the slope of the extracted Coulomb potential V C (R) is significantly (for a factor of 2-3, depending on the gauge coupling) larger than that of the static quark potential V(R) lim T→∞ V(R, T) obtained by gauge-invariant methods [129]. Although the latter is in agreement with Zwanziger inequality [205], V(R) ≤ V C (R), the potential is overconfining prompting discussions in the literature on whether the Coulomb potential in this formulation actually is the full story of confinement or some crucial ingredients are still missing. It is worth mentioning however that the asymptotic string tension of the Coulomb potential appears to vanish as soon as vortices are removed from the underlined gauge field configurations rendering the importance of the vortices for understanding the confinement phenomenon in the Coulomb gauge [206]. Such configurations without vortices in fact behave as perturbations of the free gauge theory, in consistency with expectations.
In the confined phase, the Coulomb self-energy of an isolated static charge E is expected to be infinite, and the main condition for that reads Universe 2021, 7, 9 35 of 57 where the first relation relies on the continuum limit of small eigenvalues λ → 0 of the Faddeev-Popov operator, with the corresponding eigenstates φ λ and density of the eigenvalue distribution ρ(λ). By using the lattice methods it was found that [204] ρ(λ) ∼ λ 0. 25 , yielding a divergent E → ∞ and hence satisfying the confinement criterion (87). An enhancement of ρ(λ) and F(λ) close to the Gribov horizon λ → 0 seems to be associated with the role of a center vortex ensemble. However, as was advocated in Ref. [129] the Coulomb force appears to be confining also at temperatures above the deconfinement phase transition temperature which contradicts to the fact that a confining potential must be associated with a phase of magnetic disorder. The linear confining Coulomb potential in the Gribov-Zwanziger scenario can be associated with the instantaneous part of the two-gluon correlator. So, confinement could be effectively considered as an emergent property due to a gluon exchange with a non-perturbative (dressed) gluon propagator. A naive calculation shows that a linear potential may arise if the propagator of the gluon exchange scales with momentum transfer as ∼ 1/k 4 at k → 0, at least, in one of the possible gauges [207]. One typically attempts to analyse the IR behavior of the effective gluon and ghost propagators and vertices using the formalism of the Dyson-Schwinger equations following from the disappearance of the functional integral of a total derivative, with subsequent differentiation over the sources {j k }. For a review on phenomenological implications of the Dyson-Schwinger approach, see e.g. Ref. [208] and references therein. The full gluon and ghost propagators in Euclidean spacetime are conventionally represented in terms of form factors as respectively, such that their IR behavior, as the virtuality of the exchange vanishes k 2 → 0, is controlled by where κ gh and κ gl are the so-called IR critical exponents (or anomalous dimensions), to be determined in the calculations. A necessary condition for the Kugo-Ojima confinement criterion is that the ghost propagator features an enhanced (stronger than 1/k 2 ) IR singularity, i.e. lim k→0 [J(k 2 )] −1 = 0, known as the horizon condition [209]. The second condition is the vanishing gluon propagator, lim k→0 [Z(k 2 )/k 2 ] = 0. This is the exactly case for the so-called scaling solution [210][211][212] that implies a specific relation between κ gl and κ gh in D-dimensions [209,210,212,213] For D = 4 case, the values are found to be κ gh 0.595 and κ gl −1.19, such that the gluon propagator indeed tends to vanish at k → 0. In order to explain confinement, it was argued in Ref. [214] that the quark-gluon vertex should be sufficiently singular in the long-distance limit, such that its combination with a non-singular gluon propagator gives rise to the confining potential. The scaling solution has been confirmed by a lattice analysis of Ref. [215] in SU(2) gauge theory in Landau gauge and only in D = 2 dimensions, but it was not observed for D > 2 [216,217].
Another well-known solution, the so-called decoupling solution, with has been proposed e.g. in Refs. [218][219][220]. This solution corresponds to a saturated form of the IR gluon propagator tending to a constant and, hence, effectively decouples from the dynamics in analogy to a massive particle. It is worth noticing here that the non-perturbative gluon propagator does not behave as a propagator for a massive state. Indeed, from numerical simulations one observes indications for violation of positivity, in consistency with the fact that no colored gluons exist in the asymptotic spectrum of a gauge theory that is traditionally connected to gluon confinement [221,222]. Besides, the decoupling solution implies a simple 1/k 2 pole for the ghost propagator. This solution appears to be favoured by known lattice simulations for D > 2 which also indicate a disagreement with the Kugo-Ojima criterion. A more generic criterion for quark confinement applicable in arbitrary gauges relying on the IR behaviour of ghost and gluon propagators has been proposed in Ref. [223].
One would remark here that the primary probe for the magnetic disorder phase is, of course, the area-law falloff of gauge-invariant observables, Wilson loop VEVs, and not the gluon propagator itself which is not a gauge-invariant object. So one should be extra careful in interpreting the IR behavior of the propagator in order to avoid spurious results. For recent comprehensive effort to obtain a linear static potential in the framework of Dyson-Schwinger formalism in Coulomb gauge, see Ref. [224]. For a thorough analysis of the Polyakov line VEVs and effective potential based upon the formalism of Functional Renormalisation Group [225] has been performed in Ref. [226,227], and an agreement with lattice results has been found. However, the search for the area-law dependence of large Wilson loops' VEVs with these methods has not been successful so far.
The picture of strongly collimated color-electric flux tubes stretched between the color-charged static sources does not seem to apply to the distribution of color-electric field in the Coulomb gauge [3]. Indeed, there is a significant long-range dipole contribution to the Coulomb electric field that would cause rather strong van der Waals type forces between hadrons at large distances. This would immediately contradict to the mass gap existence [139] that requires only short-range forces between composite color-neutral states. This problem generically emerges in any confinement scenario such as the Dyson-Schwinger type approaches where a confining force is associated with a single (dressed) gluon exchange at large distances. While providing a linear potential, such one-gluon exchange scenarios (including the Coulomb confinement one) imply a spread out of the electric field towards large distances, possibly with flux collimation to some extent [228].
A possible development that may eventually address the shortcomings of the Coulomb confinement scenario discussed earlier is to notice that the qq state defined in Eq. (86) is not necessarily a minimum-energy state of a system containing a single qq pair, and lower energy states could in principle be constructed using operators Q j i -functionals of the lattice links -that effectively create "constituent" coupled gluons asΨ where schematically, The resulting state effectively represents a chain of gluons bound by attractive forces, with a q and q at the end of the chain, at large R, that could, in principle, provide a necessary suppression of the long-range dipole fields. Hence, such a gluon chain may be viewed as a color-electric flux tube itself [229,230]. Indeed, as q andq get separated, more and more constituent gluons get pulled out of the vacuum to minimise the energy of the system [231,232]. This picture rather naturally emerges by expanding the Wilson line stretched between q andq in powers of the gluon field and actually implies the absence of dipole fields at large R. In the limit of large number of colors N in the SU(N) theory, such a chain of gluons on a given time slice is dominated by a high-order planar Feynman amplitude that can be, in principle, tackled by analytic methods. Among remarkable features of the gluon chain model are the Casimir scaling in the leading order of 1/N expansion and a subleading 1/N 2 string breaking effect at some critical length-scale leading to a correct N-ality dependence of the string tension asymptotically. In the case of heavy (static) charges in the adjoint representation of SU(N), for instance, in the limit N → ∞ two gluon chains instead of one are formed between the charges, leading to twice larger adjoint string tension compared to the one in the fundamental representation, i.e. σ A = 2σ F . The latter is defined only at intermediate distances, but must disappear at asymptotic distances due to color screening by N-ality zero gluons in the vacuum. Although gluons do not break the center symmetry as such, they take part in the color screening on the same footing as light quarks in QCD such that both quarks and gluons are absent in the asymptotic spectrum in the virtue of C-confinement and the string hadronisation model. This suggests a non-trivial but less explored and speculative possibility that the non-perturbative gauge-field vacuum somehow rearranges itself at large distances in such a way that the center symmetry might get broken somehow even without the presence of matter fields in the fundamental representation 4 .
Indeed, pulling the two gluons (or adjoint matter states) apart from each other, eventually the virtual gluons from the QCD vacuum are prompted to bind to the octet-charged sources yielding color-singlet states -gluelumps -at asymptotically large distances. Such a gluon color-screening mechanism is very similar to that driven by dynamical virtual quarks being brought on mass-shell to screen the charge of heavy static quarks as the latter move apart, and the energy accumulated in the string is partially spent for that purpose. So, the color-screening and hence the string-breaking phenomenon is not particularly sensitive to N-ality but rather to the color charge itself being the necessary prerequisite for C-confinement. While formation of a flux tube between the two gluons at intermediate distances applies for the confining phase in the strongly-coupled regime, C-confinement as an asymptotic phenomenon occurs also in the Higgs phase but without formation an intermediate flux tube.
Note, an adjoint string breaks via 1/N 2 suppressed but very important (at large R) interaction between the gluon chains enabling them to transform into a pair of gluelumps as described above (see also Ref. [230]). This correctly generalises for sources in arbitrary gauge group representation giving rise to N-ality dependence of the asymptotic string tension. Such an important string-like property of the gluon chain as the Lüscher term appears due to fluctuations in the gluons' positions on the chain [230]. As was demonstrated in Ref. [233], introducing up two gluons in a chain preserves the linearity of the Coulomb potential but that is already enough to bring its slope much closer to the true static potential (i.e. obtained by gauge-invariant methods). In this calculation, it was shown that multi-gluon configurations in the chain appear to be increasingly important at large R, also strongly reducing the sensitivity of the results to the lattice volume. This means that the long-range dipole field becomes strongly suppressed indicating a possible formation of a localised color flux tube (for more discussion on this aspect, see Ref. [3]).

Dual superconductivity and magnetic monopoles
As was proposed a long ago in Refs. [132,[234][235][236][237], the QCD vacuum could be viewed as a "dual" superconductor an analog of type-II superconductor where the electric and magnetic fields are interchanged. These studies have pioneered the developments of a beautiful theory of what 4 By construction, a pure gauge theory does not contain any fundamental-representation charges. So instead of heavy quarks, the use of "constituent gluons" as static color charges to probe the formation and properties of the flux tubes of finite lengths, color screening, string breaking mechanism and the phases of the theory would be the most natural approach to study confinement in pure non-abelian theories.
Universe 2021, 7, 9 38 of 57 is sometimes called the dual superconductor picture of confinement. In a usual superconductor one deals with a condensate of electric charges (in fact, bosonic Cooper pairs) and, due to repulsion (or confinement) the magnetic fields get squeezed into magnetic flux tubes with a constant energy density (Abrikosov vortices). In a "dual" superconductor one instead works with a condensate of magnetic charges known as magnetic monopoles where electric field of static charges would be squeezed (confined) into electric flux tubes. The latter realisation is what we often regard as ordinary confinement in QCD. Both, the static potential of magnetic monopoles in type-II superconductor and the static potential of color-electric charges in "dual" superconductor would rise linearly with the charge separation.
This effect gives rise to a very simple picture of confinement essentially based upon a suitable generalization of the Landau-Ginzburg superconductivity theory. Indeed, starting from relativistic abelian Higgs model on recovers the magnetic flux-tube Abrikosov-like solutions dubbed as the Nielsen-Olesen vortices [238]. Attributing a non-trivial winding number n to the Higgs complex phase, a Nielsen-Olesen vortex carries the magnetic flux 2πn/e. In the dual version, such vortex carries an electric flux that confine the electric charges. A particular model, where the dual abelian Higgs model with confinement is realised, is the N = 2 supersymmetric YM theory known as the Seiberg-Witten model [239,240] having several distinct types of electric flux tubes. In this model, a continuous set of distinct vacua is spanned by the "moduli" space of certain scalar field operators. Soft supersymmetry breaking then reduces the theory down to an effective N = 1 theory where confinement of electric charge is realised due to condensation of the monopole field and electric flux tube formation. This happens in full analogy to confinement of magnetic charge due to magnetic flux tube formation in usual type II superconductors and in the ordinary abelian Higgs model. The duality transformation in the Seiberg-Witten model inverts a certain combination of the effective coupling constant and the θ angle enabling one to obtain the effective action of light fields at any value of the gauge coupling from the detailed knowledge about the weak-coupling regime of the theory and its infrared singularities (for a detailed review of the underlined concepts and formalism, see e.g. Refs. [241][242][243]). Such a duality is due to an exact symmetry of the abelian effective theory manifest at low energies, and not of the original SU(2) theory. In fact, this duality is a proper generalisation of the famous electric-magnetic duality of the Dirac formulation of Maxwell electrodynamics (with magnetic monopoles) exchanging the electric charge q e and its magnetic counterpart q m = 2π/q e . Hence, by means of such duality transformation one hopes to learn about strong-coupling (or long-distance) dynamics of a given from the weak-coupling regime of its dual formulation. In Ref. [244] it was shown that k-string tensions in SU(N) version of the Seiberg-Witten model obey the Sine law being numerically not very different to that of the Casimir scaling (c.f. Eq. (75)).
In itself, the superconductivity picture of confinement is an abelian mechanism which has been explored originally by Polyakov [245] in the context of confinement of electric charges in D = 2 + 1 compact U(1) gauge theory. This theory turns out to be an important starting point to approach QCD confinement. While in D = 2 + 1 case, the compact QED features monopoles (topological excitations), in D = 3 + 1 those monopoles are point-like defects in spacetime, i.e. they are also instantons. Effectively integrating out all the DoFs except monopoles in D = 2 + 1 compact QED, it was shown in Refs. [132,237,245] that the action of the monopole gas interacting by means of Coulomb force on the lattice reads with i, j = 1 . . . N for N monopoles, and the lattice Coulomb propagator G at large distances behaves as G ∼ 1/4π|r i − r j |. A Wilson loop in this approach can be expressed in terms of the monopole density and appears as a current loop that generates its magnetic field being effectively screened away by the (anti)monopoles from the background. Such an effect causes the area-law falloff for the Wilson loop VEVs. Polyakov has explicitly demonstrated that even arbitrarily low density of these monopoles is sufficient to produce confinement and the mass gap of the theory. This happens in a regime when the entropy related to the size and shape of large Wilson loops wins over the cost in the monopole action for a large loop. The latter effect occurs at any coupling for D = 3 QED, but only for large enough couplings in the D = 4 case.
In the case of YM theories one needs to extract an abelian subgroup from the gauge group e.g. by means of an adjoint Higgs field. An important realisation in the case of SU(2) gauge theory is the Georgi-Glashow model where in the minimum of the Higgs potential and in unitary gauge there is a residual U(1) local gauge symmetry. Due to this symmetry, the model exhibits magnetic ('t Hooft-Polyakov) monopoles [246,247] as instanton solutions of the classical equations of motion in D = 3, or as static solutions (solitons) in D = 3 + 1. The Higgs field that is used to fix the unitary gauge necessarily vanishes at the center of each 't Hooft-Polyakov monopole making the unitary gauge fixing ambiguous at those sites. The Wilson loop VEVs are then computed in a similar way as was done in compact D = 3 QED resulting in a finite string tension σ ∼ exp(−S m ) [237]. In D = 4, the Georgi-Glashow theory has both confining and non-confining phases, however, stable monopole solutions exist only in the non-confining phase where they do not form a Coulomb plasma.
An important caveat in D = 3 theory is that one cannot simply neglect of effects of W bosons at large distances (and hence in the analysis of confinement) in the long-range effective action. Indeed, the string tensions cannot acquire a correct N-ality dependence without W bosons. The Coulomb monopole gas approximation can be justified in a certain intermediate range below a string-breaking length-scale where a W bosons carrying two units of electric charge are pair-produced and screen the charges of the static sources, also possessing two units of electric charge. Analogically, the dual abelian Higgs model that ignores the effect of W bosons predicts a wrong N-ality dependence of the Wilson loop VEVs. Thus, it is unable to consistently describe long-range physics of vacuum fluctuations at characteristic distances exceeding the color screening length-scale. Non-abelian supersymmetric versions of the dual Higgs model have been proposed in a number of existing works yielding specific non-abelian vortex solutions; see for a detailed review on these aspects e.g. Refs. [75,248] and references therein.
Dynamical "abelization" of SU(N) gauge fields can be achieved even without an adjoint Higgs field. Instead of using an adjoint Higgs field, another way to extract a Cartan (abelian) subgroup U(1) N−1 of SU(N) suggested in Ref. [249] is the so-called abelian projection, using a composite operator which transforms like a matter field in adjoint representation and fixing a gauge in which this operator is diagonal. The same effect emerges also with adjoint fermions fields [250,251], or by adding a trace deformation term to the action [252], both methods have been successfully explored by lattice simulations (see e.g. Refs. [253][254][255]).
The gluons from the coset of abelian projection are charged under U(1) N−1 while the monopole condensation would describe their confinement in a way similar to the dual abelian Higgs model. The basic idea then is to look for a specific gauge in which the quantum fluctuations of the U(1) N−1 -charged gluons are strongly suppressed compared to the fluctuations of "photons" from the Cartan subgroup U(1) N−1 . In such a gauge called maximal abelian gauge [256] the link variables would be close to a diagonal form. For instance, in the SU(2) gauge theory this is achieved by means of requiring ∑ Tr[U µ (x)σ 3 U † µ (x)σ 3 ] to be maximal while leaving the residual U(1) symmetry w.r.t. gauge transformations This enables one to decompose U µ (x) = C µ (x)u µ (x), where C µ (x) matrix is expressed in terms of a "matter" field c µ (x) with two units of U(1) charge, while the diagonal u µ (x) = diag(exp(iθ µ (x)), exp(−iθ µ (x))) is given in terms of abelian U(1) gauge field θ µ (x), the "photon", coupled to the "matter" field c µ (x). One therefore obtains the abelian-projected lattice by means of U µ (x) → exp(iθ µ (x)) projection. Note in the case of SU(3) theory the maximal abelian projection is not unambiguously defined as has been discussed for instance in Ref. [257].
In the monopole dominance approximation [258,259], one then replaces the link variables by the monopole links constructed from the Dirac string variables and the Coulomb propagator and then one computes the VEVs of the Wilson loops over an ensemble of such monopoles. This procedure leads to (almost) the same values for the asymptotic string tensions of the single-charged Wilson loops in the SU(2) lattice gauge theory as in the gauge-invariant approach. Also, the single-charged Polyakov loops computed in the abelian-projected configurations and in the monopole dominance approximation agree with each other and both vanish below the critical temperature of the deconfinement transition, in consistency with expectations. However, these results do not agree for double-charged Polyakov loops. Vanishing VEVs of the latter, and hence the confining disorder, are found in the monopole dominance approximation which is inconsistent with the charge screening effect that must be in place for double-charged static sources, and for that matter -with the N-ality requirement. This means that in the case of magnetic disorder dominated by abelian gauge field configurations the abelian flux can not be distributed according to the Coulomb monopole-gas approximation. 138 9 Monopoles, Calorons, and Dual Superconductivity  insensitive to linking with such vortices, and this will ensure that the string tensions σ n will satisfy (9.54) as required. The picture outlined above can be tested by numerical simulation. We work in the indirect maximal center gauge, which uses maximal abelian gauge as an intermediate step, and identify the locations of both monopoles (by abelian projection) and vortices (by center projection). We also measure the excess action The latter problem is not present in the abelian-projected configurations yielding a correct asymptotic behavior of large Wilson loop VEVs in the fundamental representation. Fixing an abelian projection gauge in the SU(2) gauge theory arranges the monopoles and antimonopoles coupled to each other into a chain with the total monopole flux of ±2π. At a certain fixed time, such a flux can be squeezed into center vortex structures on the abelian-projected lattice [260] -for an illustration of this effect, see Fig. 7. Indeed, the numerical analysis that locates both the (anti)monopoles through abelian projection and the center vortices through center projection showed that almost all (anti)monopoles are located on the vortex sheets arranging themselves into alternating order in a chain (for an inspiring discussion, see Refs. [3,4]). The fact that the double-charged (Wilson and Polyakov) loops do not get contributions from a linking with such vortices on the abelian-projected lattice is reflected in a vanishing asymptotic string tension in this case, in agreement with the charge screening effect.
It is instructive to introduce a specific order parameter, the VEV of a monopole creation operator µ( x) , that would signal an emergence of the dual superconducting phase in a non-abelian gauge theory [261,262]. The operator µ( x) effectively inserts a monopole configuration at a certain position into the system such that it does not commute with the total magnetic charge operator, and hence its VEV would break the corresponding dual U(1) gauge symmetry, a remnant of the gauge symmetry. According to the monopole condensation mechanism of confinement, the system is in a confining phase if and only if µ( x) = 0 while a transition to a non-confining configuration occurs when µ( x) → 0 which indeed coincides with a more generic numerical analysis in the full theory (also at finite temperatures). There are, however, severe ambiguities such that µ( x) may vanish also in the absence of any thermodynamic transition to a deconfined phase [263]. Indeed, as was already briefly discussed earlier the breaking of a gauge symmetry remnant cannot be utilised as a correct signature of the magnetic disorder phase.
In a pure non-abelian gauge theory in D = 4, classical instanton solutions can not be responsible for magnetic disorder of the vacuum field configurations since their field strength falls off too fast at large distances. However, at finite temperature the instanton solutions as saddle points of the Euclidean gauge fields' action called calorons can be relevant for confinement. The latter solutions were found in Refs. [264][265][266] and are known in the literature as KvBLL solutions. They may contain monopole constituents sourcing both electric and magnetic fields, also known as dyons or Bogolmolny-Prasad-Sommerfield (BPS) monopoles [267,268], which can be widely separated. The thermal approach to pure 4D YM theories based upon nonperturbative results on a thermal ground state in the deconfined phase derived from (anti)caloron ensemble has been thoroughly discussed in Ref. [269] and in references therein. Among important corollaries to this approach is, for instance, the derivation of the 3D critical exponent of the Ising model for correlation length cirticality.
The early work of Ref. [134] made important contributions to understanding the trivial-holonomy calorons in SU(2) Euclidean gauge theory based on Ref. [270], while nontrivial holonomy solutions have been studied e.g. in Refs. [264,266,[271][272][273]. Considering for instance the maximally non-trivial Polyakov loop holonomy P( x) introduced in SU(N) in Eq. (57), where µ j are ordered and spaced with a maximal distance from each other the probability density of calorons in the vacuum would be peaked at TrP( x) = 0. As was discussed earlier, in the center vortex mechanism of confinement the vanishing Polyakov loop expectation value computed on an ensemble, where positive and negative fluctuations in vortex configurations cancel out, is a signature of unbroken center symmetry and, hence, that of the confinement property.
Notably, in the caloron configurations the maximally non-trivial Polyakov loop holonomy vanishes by itself before any averaging as the basic property of such configurations. Due to this property, a system of widely separated dyons, the dyon gas, whose free energy is minimal for Tr P( x) = 0, has been considered as the basis for description of the magnetic disorder in YM theories [274]. Indeed, it was shown that the k-string tensions extracted from space-like Wilson loops are in agreement with those that determine the asymptotic behavior of Polyakov loop correlators, and follow the Sine law (97). In the high temperature regime, a phase transition to the deconfinement phase occurs with T c / √ σ values being in perfect agreement with numerical lattice results. Despite of such tremendous success, the path integral measure for the multi-dyon configurations appears to be not positively definite thus violating the basic property of the exact measure [275]. However, numerical simulations of Refs. [276] with a suitable parameterisation of the integration measure confirmed that the confining static potential indeed emerges in the dyon gas approximation. Thus, one can conclude that the monopole mechanism based upon the caloron classical solutions is one of the most promising scenarios of confinement in non-abelian gauge theories.
There are some critical points to be made regarding the dyon gas picture of confinement neatly summarised in Ref. [3]. The same question as for the monopole Coulomb gas applies also for the dyon gas regarding the asymptotic string tension of double-charged Wilson loops that should disappear due to a screening by gluons. Another question concerns the probability distribution of Polyakov loop holonomies which is peaked for a vanishing maximally-nontrivial holonomy in the fundamental representation i.e. Tr P( x) = 0. If this is indeed true, it implies a negative expectation value of the Polyakov loop in the adjoint representation. However, if the latter is positive, the probability distribution would be peaked at the center-element holonomy as suggested by the center vortex scenario of confinement. Remarkably, the expectation value of the adjoint Polyakov line in the phase of magnetic disorder has been found to be positive in the SU(3) theory in Ref. [277].
Besides, considering the asymptotic behavior of double-winding Wilson loops VEVs' (i.e. Wilson loops winding around closed co-planar loops C 1 and C 2 ), one reveals a dramatic difference in predictions of the monopole and vortex mechanisms of confinement [278] (see also [3,4]). The vortex scenario provides the "difference-of-areas" law behavior for such loops, ∼ exp[−σ|A(C 1 ) − A(C 2 )|, where the "−" sign is due to a vortex linking to the largest loop, correctly reproducing the full lattice results. In this case, the monopole scenario predicts the "sum-of-areas" falloff as ∼ exp[−σ(A(C 1 ) + A(C 2 ))], which is disfavoured by numerical simulations. The latter observation indicates that the vacuum cannot be in a dual-superconducting state of monopole/dyon plasma. In fact, the "difference-of-areas" law is restored by heavy W bosons that are present in the full YM theory, but not in an abelain part of it. As was suggested in Ref. [4], upon integrating out the W bosons' states, one expects the monopole-antimonopole lines to get collimated into Z 2 vortices, as is illustrated in Fig. 7, i.e. in a similar fashion to what has been seen on an abelian-projected lattice. This effectively turns the monopole ensemble into a configuration of Z 2 vortices, offering an intricate connection between the two pictures of confinement.
Another observation of Ref. [279] in the case of G 2 gauge theory has suggested that the Polyakov loop expectation value exactly vanishes in the dyon gas picture. The color screening in G 2 , however, requires to bind a static source in the fundamental representation to a minimum of three gluons likely leading to a very small, but non-zero, Polyakov loop expectation value that would be difficult to identify numerically [3]. If the color screening mechanism is a valid approach, the center symmetry breaking at large distances would be manifest in G 2 , whereas the dyon gas approximation, like the Coulomb gas approximation discussed earlier, might be lacking something relevant in the asymptotic regime.

Separation-of-charge confinement criterion
A clear symmetry-based distinction between the confining and Higgs phases in a gauge theory with fundamental representation matter fields has been recently proposed in Refs. [110,280,281]. An important generalised criterion of confinement valid in both pure YM theories and YM theories with matter in fundamental representation states that with Ψ V being the qq state connected by a Wilson line, for any choice of the gauge bi-covariant non-local operator V ab ( x, y; A). The latter depends only on the gauge field, thus, eliminating any possibility for a string breaking by means of dynamical matter fields. This criterion is a necessary and sufficient condition for separation-of-charge confinement (or S c -confinement, for short), which is meaningful only in gauge theories with a non-trivial center symmetry. In Eq. (101), H is the Hamiltonian, E vac is the vacuum energy, E 0 (R) ∼ σR at R → ∞ is an asymptotically linear function which has the meaning of the ground-state energy of the qq in a pure SU(N) gauge theory (but not in the one with matter fields) where the above criterion is equivalent to the area-law falloff of Wilson loop VEVs.
In the SU(2) gauge-Higgs theory, the confining phase is found for γ β 1 and γ 1/10 where the S c -confinement condition (101) is satisfied. But deeply in the Higgs phase, for other couplings' ranges, this criterion is not fulfilled, hence we deal with only a weaker C-confinement situation there. In Refs. [110,280,281] it has been shown that a transition between the Cand S c -confinement phases must take place in the gauge-Higgs theory, and the unbroken custodial symmetry has been found to separate the S c -confining (if not massless) phase from the Higgs phase corresponding to a C-confined spin glass state where the custodial symmetry is actually broken. It would be very interesting to see how such a new concept of S c -confinement can be applied for more realistic theories such as QCD.

Separating the Higgs and confinement phases: vortex holonomy phase
As was discussed earlier, the results of Refs. [109,113,114] state that it is possible to identify a continuous path between the Higgs and confinement regimes where no first-order phase transitions occur like what is believed to happen in physics of high-to-low temperature QCD (smooth crossover) transition and in several other specific models. Indeed, one would naively expect that such a continuity always takes place unless the phases are separated by different realisations of global symmetries.
An important counter-example to this statement has been recently explored in a whole class of models in Refs. [58,64]. Namely, it has been demonstrated that even in the case of a spontaneous global U(1) symmetry breaking in both Higgs and confinement regimes (an analog to the baryon symmetry breaking in dense QCD at low T) it is still possible to identify a novel non-local order parameter (a vortex holonomy phase) that separates the two phases leading to a thermodynamical phase transition. This proof provides an important argument against the quark-hadron continuity (Schäfer-Wilczek) conjecture in dense QCD mentioned above in Sect. 2.
For this purpose, the authors of Ref. [58] started with the Polyakov's D = 4 compact U(1) gauge theory discussed in the previous section, then imposed a single global U(1) G symmetry (an analog to U(1) B of QCD) and added three complex scalar fields that effectively mimic dynamical quark fields in QCD, with the following assignments under U(1) × U(1) G group Then, the product φ + φ − appears to be an analog to the gauge invariant baryon operator, φ 0 can couple to φ + φ − product and can be considered as a baryon interpolating field or a source for baryons. A VEV in φ 0 would bring the theory into a superfluid phase of spontaneously broken U(1) G , such that the global symmetries' realisation is the same in the confinement and Higgs phases.
In the monopole-driven confinement picture, the monopoles in this theory would have a finite action that can be UV-completed through an SU(2) symmetry just as in the original Polyakov's description. Finally, additional discrete Z 2 (charge conjugation and flavor-flip) symmetries have also been introduced, and can be considered as analogous to the flavor symmetry in QCD. The effect of monopoles in this model induces an additional term in the Lagrangian [237] V m (σ) ∝ e −S I cos σ , which is a potential for the "dual photon" -a periodic scalar field σ → σ + 2π related to the field strength by the abelian duality relation As one varies the adjustable mass parameters three different regimes of the theory emerge. One of them corresponds to compact 3D U(1) theory with heavy scalar quarks and confinement where no symmetries are spontaneously broken ("gapped confined" regime). The second regime features the Higgs mechanism with a non-zero VEV φ + φ − = 0, such that a cubic φ + φ − φ 0 term in the potential drives condensation of φ 0 , φ 0 = 0 and hence spontaneous breaking of U(1) G (Higgs phase with a single massless Goldstone boson). The third regime has monopole-driven confinement, while φ 0 = 0 spontaneously breaks U(1) G symmetry and the heavy charged scalars are very heavy and may be disregarded. A key question here is that whether the two phases with spontaneously broken U(1) G ("Higgs" and "confining") with no distinguishing local order parameter are really two distinct phases or might be continuously connected.
The main claim of Ref. [58] is that these phases are distinct and can be distinguished only by a new non-local order parameter that is connected with topological excitations. Physically, both phases with spontaneously broken U(1) G should be considered as superfluids that have vortices. As a consequence of U(1) G breaking, the theory possesses a gapless Nambu-Goldstone mode which is a phase of φ 0 condensate, as well as topologically stable vortex excitations when the phase of the φ 0 condensate winds around the unit circle when one goes around a given loop. In ordinary superfluids in three spacial dimensions this provides vortex loops, but in D = 2 + 1 vortices act like point particles, the point around which the condensate phase winds. Winding of that phase, in the language of superfluids, is exactly what is called quantized minimal "circulation". The winding number can then be found in terms of a contour integral of the gradient of the phase, or simply as The charge particles in the superfluid phase would then interact with (minimal-energy) vortices through acquiring an Aharonov-Bohm phase as one sends a charged particle into a loop that links with the worldline of the vortex. This phase is measured by the Wilson loop holonomy whose expectation value Ω(C) ∼ e −mP(C) .
This means that short-range quantum fluctuations in the phase (with dynamical fundamental representation charges) automatically lead to a perimeter-law decay of the VEV of a large Wilson loop. In fact, this represents the same physics as that of string breaking that turns the area-law behavior of a pure YM theory into a perimeter-law behavior present in real QCD and in the Higgs phase as was discussed in detail in previous sections. Let us consider a Wilson loop expectation value in the presence of a minimal-energy vortex which can be thought in terms of a constrained functional integral where one integrates over all the field configurations in the theory but with a constraint forcing the presence of a vortex along some large worldline C in D = 3 spacetime. The same short-range physics guarantees it is going to have the same perimeter-law behavior but it can have an additional phase factor Ω(C) w=1 ∼ e iΦ e −mP(C) . (109) In order to extract the phase, one defines the ratio [58] -the vortex order parameter, taking the size of the Wilson loop and its separation from the vortex r v arbitrarily large simultaneously. The symmetries of the theory guarantee that Wilson loop VEV Ω(C) is real, and it can be made positive. While the charge conjugation and reflection symmetries are broken, the flavor-flip symmetry is preserved in the presence of a vortex. The latter also flips the sign of the gauge field and hence conjugates the holonomy guaranteeing that the vortex-constrained expectation value of the Wilson loop is also real, but can be either positive or negative, i.e. O Ω = ±1. This means that the vortex order parameter cannot vary smoothly under variations of model parameters when moving between the two phases. In Ref. [58] it was demonstrated by means of a semi-classical analysis and through minimization of the vortex energy that in the Higgs phase the vortex order parameter must be equal to −1.
Integrating out the heavy charged scalars in the confining phase, one can conclude that the vortex-constrained Wilson loop expectation value hardly knows about the presence of the vortex, thus O Ω = 1. It was argued in Ref. [58] that if one varies parameters of the theory and at some point the magnetic flux carried by vortices suddenly jumps, which is what the vortex order parameter is really probing, that surely is going to change the core energy density of the vortex. In this case, it does change the probability of having vortex excitations in the ground state wave function affecting the ground state energy density. In other words, a sudden change in the vortex properties really should be reflected in genuine thermodynamic phase transition.
Given the close analogies of the considered model with QCD, by construction of the model above, this conclusion may be straightforwardly generalised to a D = 4 non-abelian theory with fundamental-representation matter charged under a given global symmetry. This is the case of dense QCD with broken U(1) B , where the SU(3) vortex order parameter can be shown to take two distinct values in two phases [58] O Ω ≡ lim TrΩ(C) = e 2πi/3 CFL/Higgs phase +1 nuclear/confining phase , (111) such that the latter O Ω = +1 is understood as the characteristic signature of confined QCD phase with spontaneously broken baryon symmetry. This indeed illustrates the main point suggesting that the quark-hadron continuity between the nuclear and quark-matter phases may not hold. Ref. [70] has argued that despite the noted discontinuity in the vortex order parameter the continuity of phases may still be intact due to a continuous connection between the vortex in the CSC/CFL phase and the corresponding one in the nuclear phase. In response to this claim, Ref. [58] has explicitly proven the existence of thermodynamical phase transition at the interface between the two phases connected to the manifest discontinuity in the non-local vortex order parameter. While the debate about this important issue will likely continue in the literature, it once again reveals the surprising underlined complexity of the non-perturbative QCD vacuum and the associated approaches to the confinement problem may still not be in their final form. It would be very instructive to find possible connections between the vortex holonomy phase and its discontinuity with the S c -confinement criterion briefly discussed in the previous section, both pursuing the same goal of sharply separating the Higgs and confining phases.

Summary
To summarise, the mass gap and color confinement that are realised already in a gauge-Higgs theory may not be connected to an asymptotically rising static potential and Regge trajectories and, hence, they do not necessarily represent an emergence of the magnetic disorder state. On the other hand, the magnetic disorder and the associated area-law behaviour of Wilson line VEVs imply color confinement and the mass gap automatically. In this sense, color confinement and the mass gap represent only a small part of a bigger picture of confinement and should be considered as a consequence of the confined magnetic disorder state and flux tubes formation corresponding to a phase with unbroken non-trivial center symmetry. If there is no a non-analytic boundary between the massive and magnetic disorder phases at finite values of coupling constants and at some critical length-scale, i.e. a first-order phase transition, one should talk about a single massive phase at all scales, like for instance in the gauge-Higgs theories. The flux tubes formation is only an approximate picture in this case, roughly consistent with reality at some intermediate distances, but it does not necessarily represent an emergence of a new phase.
One of the big questions for real QCD though, i.e. with physical quarks and gluons, which would distinguish it from EW theory is then whether a magnetic disorder phase really exists within some finite interval of characteristic length-scales that would abruptly transit to a massive phase at asymptotically large length-scales (due to string-breaking), or not. If not, then real QCD would always be considered on the same footing with a gauge-Higgs theory as existing in the massive phase only which is one of the basic options actively discussed in the literature. One thing, however, that distinguishes real QCD from the EW theory is the existence of experimentally observed Regge trajectories in QCD with light quarks that, in fact, may indicate the presence of a non-analytic phase boundary at moderately large distances in QCD, in variance to EW theory, while the both would be in the massive phase asymptotically. Numerical values of the coupling constant here should play decisive role here, and for weak couplings the magnetic disorder may not emerge at all.
In fact, light "sea" quarks, i.e. with masses way below the confinement energy scale of QCD, emerge due to gluon splitting G a → qq such that correlated qq pairs could be viewed as effective gluons as long as the resolution length-scale is above the wave-length of such a pair. If at such length-scales the strong-coupling constant is large enough, one can view physics in such a regime as that of an effective pure YM theory in a magnetically disordered phase with unbroken center symmetry. By pulling q andq apart from each other at length-scales larger than the resolution scale, the center symmetry gets effectively broken and the theory enters the massive phase. In the infinite quark mass limit, however, the string-breaking length-scale grows indefinitely making the magnetic disorder phase valid for asymptotically large distances.
Depending on the values of the gauge coupling constant, the same can occur in a gauge-Higgs theory in a strongly coupled regime which is supported by lattice simulations. Thus, we arrive to radically different phase structure of a gauge theory depending on whether it is in a strongly-coupled or in a weakly-coupled regime. But even without matter fields involved, the major problem of confinement remains, namely, to understand why pure YM theories with a non-trivial center symmetry in D ≤ 4 dimensions can only exist in a state of magnetic disorder. Once this key problem is solved, it will become clearer under what conditions realistic theories like a gauge-Higgs theory or real QCD may exist the magnetically disordered phase, and a first-order phase transition towards the massive phase may occur, if at all.
As was elaborated in this review, the phase structure and properties of the quantum QCD vacuum is still under intense explorations, both experimentally and theoretically, numerically and analytically, and is far from its complete and satisfactory description. However, a tremendous progress has been made and some basic contours of the fundamental picture of confinement have started to emerge. We do understand a confining phase as an asymptotic magnetic disorder phase with unbroken non-trivial center symmetry that manifests itself through the area-law behavior of large Wilson loop VEVs and, hence, a linear rise of the corresponding static (string) potential. However, real QCD features such a phase only pre-asymptotically where color-electric flux tubes exist at not-too-large distances, while they break apart at length-scales beyond an inverse to the lightest meson mass (pion) scale yielding a massive phase asymptotically. Quanta of magnetic flux, vortices, have proven to play a crucial role at all stages, from formation of a flux tube to its breaking. Given the overwhelming qualitative and quantitative evidence collected in vast amounts of studies in the literature, the center vortex mechanism remains among the most favored scenarios of confinement so far. Other ideas such as the monopole scenario highlight the underlined complexity of confining phase and phase transitions and offer different perspectives, but in one way or another connect to the vortex picture. Various order parameters briefly described in this review probe the confining phase and are capable of separating it from a non-confining (Higgs) phase, with the latter remaining under a continuous debate in the literature.