Confinement on ℝ3 × 𝕊1 and double-string collapse

We study confining strings in N\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ \mathcal{N} $$\end{document} = 1 supersymmetric SU(Nc) Yang-Mills theory in the semiclassical regime on ℝ1,2× 𝕊1. Static quarks are expected to be confined by double strings composed of two domain walls — which are lines in ℝ2 — rather than by a single flux tube. Each domain wall carries part of the quarks’ chromoelectric flux. We numerically study this mechanism and find that double-string confinement holds for strings of all N-alities, except for those between fundamental quarks. We show that, for Nc≥ 5, the two domain walls confining unit N-ality quarks attract and form non-BPS bound states, collapsing to a single flux line. We determine the N-ality dependence of the string tensions for 2 ≤ Nc≤ 10. Compared to known scaling laws, we find a weaker, almost flat N-ality dependence, which is qualitatively explained by the properties of BPS domain walls. We also quantitatively study the behavior of confining strings upon increasing the 𝕊1 size by including the effect of virtual “W-bosons” and show that the qualitative features of double-string confinement persist.


Introduction, summary and outlook
Quark confinement is an old and difficult problem in quantum field theory. 1 A direct analytical attack in non-supersymmetric gauge theories on R 1,3 is beyond current capabilities. The study of "toy" theories approximating different aspects of the physical world is thus a worthwhile endeavour. Supersymmetric Yang-Mills (SYM) theory stands out in this regard. It has been the subject of many studies due to its tractability and its similarity to non-supersymmetric pure Yang-Mills (YM) theory, to which it is connected by decoupling the gaugino.
Of special interest to us in this paper is the theory formulated on R 3 × S 1 , with supersymmetric boundary conditions on S 1 [2,3]. Here, the dynamics is drastically simplified when LN Λ 1; L is the S 1 circumference, N is the number of colors, and Λ is the strong coupling scale of the SU(N ) gauge group [4,5]. Remarkably, in this regime, center stability, confinement, and discrete chiral symmetry breaking become intertwined and are all due to the proliferation of various topological molecules in the Yang-Mills vacuum: magnetic [6,7] and neutral "bions" [8][9][10][11][12], composite objects made of various monopole-instantons [13,14]. The magnetic bions, in particular, provide a fascinating and non-trivial locally four-dimensional generalization [6,7] of the Polyakov mechanism [15] of confinement. Despite looking three-dimensional (3D) at long distances, the theory remembers much about its four-dimensional (4D) origin. In many cases, it is known or believed that the small-L theory is connected to the R 4 theory "adiabatically," i.e. without a phase transition. Lattice studies [16][17][18] have already provided evidence for the adiabaticity.
The R 3 × S 1 setup described above provides a rare example of analytically tractable nonperturbative phenomena in four dimensional gauge theories and many of its aspects have been the subject of previous investigations, see [19] for a review. Here, we continue the study of the properties of confining strings in this calculable framework. Our goal is to study the detailed structure of confining strings in SYM in the semiclassical regime: we want to determine the string tensions, their N -ality dependence, and get at least a glimpse at how the confining string configurations evolve upon increase of L towards the R 4 limit.
In ref. [20], it was shown that heavy fundamental quarks are confined by particular "double-string" configurations: a quark-antiquark pair has its chromoelectric flux split into two -rather than a single -flux tube. These flux tubes are actually domain walls (DW) connecting the N different vacua of SYM. Static DWs are lines in R 2 of finite energy per unit JHEP01(2021)044 length. Along their length, they carry a fraction of the chromoelectric flux of quarks. These fractional fluxes are precisely such that a pair of DWs form a "double string" stretched between a quark and antiquark, leading to quark confinement with a linearly rising potential. 2 The double-string nature of confining strings is closely tied to the recently discovered mixed 0-form chiral/1-form center anomaly [21][22][23][24][25][26] in SYM. 3 The anomaly-inflow aspect and the quark deconfinement on DWs was studied in detail in ref. [28]. There, the chromoelectric fluxes carried by Bogomolnyi-Prasad-Sommerfield (BPS) DWs connecting the N vacua of SYM were determined. In our study of confining strings, we shall make extensive use of these results.

Summary of results and outline
In this paper, we study in detail the double-string confinement mechanism, in the framework of the abelian EFT valid at LN Λ 1, also accounting for the leading effect of virtual W -boson loops. This EFT is described in some detail in section 2.
For the reader familiar with the R 3 Polyakov model, but less familiar with confinement in circle-compactified theories at small L, we review, in section 3, the main differences between Polyakov confinement and confinement in SYM on R 3 × S 1 , for gauge group SU (2). This discussion does not suffer from too much Lie-algebraic notation needed to study SU(N ), but is sufficient to stress the main peculiarities of small-L confinement in SYM theory, which persist to arbitrary values of N .
For general SU(N ) gauge groups, the double-string confinement mechanism is explained in sections 4.1, 4.2. The taxonomy of the double-string configurations of sections 4.2.1-4.2.4 captures most qualitative aspects of the physics for any number of colors N and for quark sources of any N -ality k. It is found, via numerical simulations, to fail only for N ≥ 5, k = 1, when DWs exhibit attraction and the double string collapses.
Our numerical simulations were performed for SU(N ) gauge theories with N ≤ 10. The numerical computations of the two-dimensional classical equations of motion with quark sources (see section 5) can become quite demanding, especially when sources of all Nalities for up to N = 10 are considered. Including the effect of virtual W -bosons demands even more resources due to the broadening of the double strings. Our computations for the confining string configurations were primarily performed on the Niagara cluster at the SciNet HPC Consortium [29,30]. Each quark source configuration was run on a single hyperthreaded processor core, allowing us to run 80 quark sources simultaneously on a single compute node. In total, our simulations utilized approximately 15 core-years of compute time on Niagara. 2 A typical example of a double string is shown on figure 5 in the bulk of this paper. 3 Historically, within quantum field theory, the double-string picture was the first used to argue that quarks become deconfined on DWs -the quark's chromoelectric flux is absorbed by the DWs, of equal tension, to the left and the right. The equal tensions of different walls in SYM is due to supersymmetry, as many of the walls are BPS objects (in non-supersymmetric generalizations the unbroken 0-form center symmetry also plays a crucial role). The anomaly inflow connection, understood later, gives an additional general argument valid beyond semiclassics [21][22][23][24][25][26]. Earlier arguments within string theory originated in [27]. Recall that the sine law approximately holds in softly-broken Seiberg-Witten theory and the MQCD embedding of SYM. The square-root law was found in the MIT Bag Model, as well as, approximately, in deformed Yang-Mills theory on R 3 ×S 1 . Our results show that the N -ality dependence of f in SYM on R 3 × S 1 at small L is flatter than any of the other known scaling laws. A qualitative explanation follows from the BPS DW properties and is given in the text. The flat N -ality dependence shows that in the abelian large-N limit, confining strings remain interacting.

JHEP01(2021)044
Our main results are as follows: 1. We first argue, in section 4.2, that for static quark sources of N -ality k, the lowest tension confining strings are double-string configurations composed of two of the BPS 1-walls of SYM theory. These BPS 1-walls carry electric fluxes that add up to match the flux of the quark source w k , the highest weight of the k-index antisymmetric representation.
We also show that quarks of N -ality k with charges different from w k are confined by double-string configurations made of two BPS p-walls, with 4 0 ≤ p ≤ k. These metastable configurations are expected to decay to the lowest tension ones by the pair production of W -bosons, an effect that cannot be seen by the small-L abelian EFT.
2. Numerically (see section 6) the above qualitative expectations are found to hold for N ≤ 4 for all k. However, for N ≥ 5 double strings form only for k > 1, i.e. for all JHEP01(2021)044 but the fundamental representation quarks. Fundamental quark sources for N ≥ 5 were found to be confined by a single-string configuration, formed upon a collapse of the two BPS 1-walls into a non-BPS bound state. See section 6.1.
The collapse of the double string for k = 1 is the most unexpected result of our findings and, at present, we have no simple qualitative explanation for it. To make sure it is not a numerical artifact, we perform many checks on this collapse. One includes simulating the largest possible distance between the sources, making sure that the collapse is not a finite separation effect (see figure 9). Most importantly, we show that the fact that this collapse occurs only for k = 1, N ≥ 5 can already be seen by studying the interactions of parallel BPS DWs with appropriate fluxes. These onedimensional configurations are studied in great detail in section 7, confirming that the collapse occurs only for BPS 1-walls whose fluxes add to that of a fundamental quark.
3. The results of the two items above imply that the N -ality dependence of the string tensions is rather weak. Qualitatively, this is because the k-strings of lowest tension with k > 1 are all double strings made of BPS 1-walls carrying appropriate fluxes. Hence, one expects that the k-string tension is approximately twice the tension of the BPS 1-wall, for all k > 1. This expectation is borne out in our numerical simulations, see figure 1 for SU (10), as well as section 6.2. The observed slight increase of the string tension with k seen in the data is consistent with the increase of the length of the strings due to the larger transverse separation for larger N -ality. On the other hand, the fact that for k = 1 the tension is significantly lower is due to the double-string collapse and reflects the binding energy of the non-BPS bound state.
The N -ality dependence found here does not conform to the usual large-N argument that T (k, N ) should approach kT (1, N ). The difference between the abelian and nonabelian large-N has been emphasized many times [8,[31][32][33]: the confining strings discussed here do not become free in the N → ∞, LN -fixed abelian large-N limit. 4. We find that, whenever confinement is due to a double string, the separation between the strings increases as a logarithm of the distance between the sources. This qualitatively matches with a simple model of exponential DW repulsion proposed earlier for an SU(2) theory [20]. The fit of the parameters to the logarithmic growth is in section 6.3.

5.
We probe larger values of L by including the effect of virtual W -boson (and superpartner) loops on double-string confinement for N ≤ 9 in section 6.4. These effects become more important as the S 1 size increases. The main lesson is that, while still in the semiclassical regime, these effects do not change the qualitative picture of doublestring confinement for all k and N ≤ 9: no separation of collapsed k = 1 strings for N > 4 occurs, nor do the N ≤ 4 k = 1 strings collapse. The quantitative changes due to increasing L are the larger transverse separation of the double string and the change of the k-string tension ratio, see appendix C.3 and C.4. Upon increasing L within the abelianized regime, the range of variation of f (N, k) is seen to be even smaller than that shown in figure 1 (see figure 15).

JHEP01(2021)044
The expectation for the string-tension ratio on R 4 , based on large-N arguments and the M-theory embedding of SYM, is that f (N, k) is smilar to the sine-law and Casimirscaling curves on figure 1. The behavior we find indicates that the evolution of the string-tension ratio towards R 4 may not be monotonic. It is also possible that at the finite values of the coupling (see eqs. (6.8), (6.11)) we chose to study, the unknown higher-order corrections to the Kähler potential qualitatively affect the string-tension ratios.

Outlook
In this paper, we studied the details of double-string confinement in the semiclassical regime on R 3 × S 1 , in the framework of SYM. While originally found independently, double-string confinement in semiclassical circle-compactified theories is, in fact, required by anomaly inflow of the mixed anomalies involving the spontaneously broken discrete 0-form (e.g. chiral or CP) symmetry and the 1-form center symmetry [28]. This confinement mechanism transcends SYM and also holds in other theories that have similar higher-form anomalies. Two examples are deformed Yang-Mills theory and QCD with adjoint fermions, and more can be found in [34,35]. It would be of interest to study double-string confinement in more detail in those theories as well. We note that QCD with adjoint fermions stands out among the many theories, since it has gapless modes in the bulk.
The question of the behavior of confining strings as L is increased beyond the ΛN L < 1 regime is, of course, of great interest. Lattice studies of adiabatic continuity in SYM, as in [17], and in other theories should be able to shed light on the transition between the abelian confinement studied in this paper and the full-fledged non-abelian confinement in the R 4 theory.
There are several other features of circle-compactified theories, their domain walls, and their relation to confinement that would be nice to understand better. All our studies here, and in virtually all papers on the subject, are done in the framework of the 3D EFT of the lightest fields, the Kaluza-Klein zero modes of the 4D fields. The small-L theory, however, is weakly coupled at all scales and it should be possible -although it may be complicated, see [36] for related work and the recent remarks in [37] -to study more interesting dynamical questions that necessitate the incorportation of the Kaluza-Klein modes in the description. One such question would be to understand dynamically how quarks "become anyons" on the domain walls. This aspect of anomaly inflow on domain walls, see e.g. [26], is not possible to see in our 3D EFT. Further, the pair production of "Wbosons," which we allude to many times and which makes the string tensions N -ality-only dependent, is also outside of the realm of validity of the 3D EFT. Naturally, incorporating the Kaluza-Klein modes in the dynamics might shed further light on the large-L behavior.
We hope to return to these interesting questions in the future.

JHEP01(2021)044
2 SYM on R 3 × S 1 : effective field theory, scales and symmetries We shall not describe in detail the microscopic dynamics leading to the long-distance theory 5 described below. Our starting point here is the result that the infrared dynamics of SYM on R 3 ×S 1 is described by an effective field theory (EFT) of chiral superfields. These superfields are the duals of the Cartan subalgebra photons and their superpartners. The simplicity of the EFT is due to the fact that the theory abelianizes, breaking spontaneously the gauge group SU(N ) → U(1) N −1 at a scale of order the lightest W -boson mass, m W = 2π LN . Weak coupling is then owed to the small-L condition m W Λ. The lowest components of the chiral superfields of the EFT are the Cartan-subalgebra valued bosonic fields, the dual photons and holonomy scalars combined into dimensionless complex scalar fields: Here φ a are real scalar fields describing the deviation of the S 1 holonomy eigenvalues from its center-symmetric value 6 and σ a are the duals of the Cartan-subalgebra photons. The more precise relation to 4D fields is as follows: 7 Here F a µ3 denotes the mixed R 3 -S 1 component of the Cartan field strength tensor of the 4D theory, F a µν is the field strength along R 3 , all taken independent of the S 1 -coordinate x 3 , and g 2 is the 4D SYM gauge coupling at a scale of order 8 1 L . Physically, eq. (2.2) implies that spatial derivatives of φ a are 2D duals of the 4D magnetic field's components along R 2 ; likewise, spatial derivatives of σ a are 2D duals of the 4D electric field's components along R 2 . 9 As indicated in (2.1), the target space of the dual photons fields is the unit cell of the SU(N ) weight lattice spanned by w k = (w 1 k , . . . , w N −1 k ), k = 1 . . . N − 1, the fundamental 5 In historical order, see [2,3] and the instanton calculation of [4,5], completed recently in [9,38] by the calculations of the non-cancelling one-loop determinants in monopole-instanton backgrounds. 6 Explicitly, the N eigenvalues of the S 1 -holonomy in the fundamental representation, labeled by where ρ and νA denote the Weyl vector and the weights of the fundamental representation (see our notations and conventions below). We give this explicit formula in order to stress that O(1) excursions of φ away from the origin do not lead to a breakdown of semiclassics. Such a breakdown requires that two eigenvalues of the holonomy coincide, causing the appearance of massless W -bosons, but this requires O(1/g 2 ) deviations of φ away from the origin. However, all our numerical results for DWs and strings do not involve such large excursions -instead, the range of variation of φ in the classical solutions of our EFT is O(1). 7 The spacetime metric is (+, −, −) and 012 = 012 = 1. 8 See the discussion around eq. (2.8) for more detail on the renormalization scale. 9 For example, labeling the spacetime coordinates as (0, 1, 2)=(t, y, z), as we do later in this paper (reserving x for the complex field (2.1)), the spatial gradients g 2 4πL ∂zσ a = −F a 0y = −E a y and g 2 4πL ∂yσ a = F a 0z = E a z are 2D duals to (or a 90 o rotation of) the electric field.

JHEP01(2021)044
weights. 10 A simple way to understand the σ-field periodicity is that it allows for non vanishing monodromies corresponding to the insertion of probe electric charges (quarks) of any N -ality, as appropriate in an SU(N ) theory. 11

EFT
At small LN Λ 1, the bosonic part of the long distance theory is described by the weakly-coupled R 1,2 Lagrangian: Here, W (x) is the holomorphic superpotential: α A is the affine (or lowest) root of the SU(N ) algebra. 12 The Kähler metric and its inverse are denoted by K ab and K ab in (2.3). 13 As opposed to the superpotential W , whose form can be determined by holomorphy and symmetries, the Kähler potential is under control only at small L. In fact, only one term in its weak-coupling expansion at ΛN L 1 has been computed. The expression for the inverse Kähler metric is where the one-loop correction is k ab . Physically, the matrix k ab determining the inverse Kähler metric represents the mixing of the electric Cartan photons due to the non-Cartan "W -boson" and superpartner loops. The form of the inverse Kähler metric K ab (2.5) is justified in the semiclassical LN Λ 1 limit and the leading correction k ab was computed in various ways in [9,38,39]. The Lagrangian (2.3) in terms of the dual photon fields is obtained after a linear-chiral superfield duality transformation, which takes into account the one-loop mixing. The duality is performed in an expansion in powers of g 2 and is 10 Vectors in the Cartan subalgebra will be sometimes denoted by bold face: σ = (σ 1 , . . . , σ N −1 ), φ = (φ 1 , . . . , φ N −1 ) and the complex field x = (x 1 , . . . , x N −1 ) and similar for their complex conjugatesx. The dot product used throughout is the usual Euclidean one. We shall interchangeably also use a notation σ instead of σ (and likewise for group-lattice vectors) to denote these Cartan subalgebra vectors. 11 Explicitly, from (2.2), the monodromy of the σ field around a spatial loop C ∈ R 2 (the y-z plane), taken in the counterclockwise direction, is, within our conventions, dσ a = − 4πL where n is an outward normal to C and the arrows denote vectors in R 2 . Thus, the monodromy dσ measures the total electric charge bounded by C. It does, in fact, equal 2πλ, where λ is the vector of U(1) N −1 charges inside C. See eq. (4.4) and figure 2 below. 12 Roots are normalized to have square length 2; roots and coroots for SU(N ) are identified, and αA ·wB =

JHEP01(2021)044
described in detail in [38]. 14 The one-loop correction to the (inverse) Kähler metric is currently known at the center symmetric point: Here ψ is the logarithmic derivative of the gamma function and β AB denote the positive roots of SU(N ). These are N −1 dimensional vectors in the root lattice (whose components are denoted ( β AB ) a , a = 1, . . . , N − 1); recall that the simple roots α A = β A,A+1 , A = 1, . . . N − 1 are a subset of the positive roots. We note that (2.7) is a purely group theoretic factor; its implications will be discussed after we introduce the relevant scales.

Scales
The scales appearing in the long-distance theory (2.3) are determined by the microscopic dynamics of the underlying 4D SYM theory on S 1 of size L. We begin with the strong coupling scale Λ of the SU(N ) theory ΛL 4π Here, and everywhere in this paper, g denotes the 4D SU(N ) gauge coupling taken at the scale 4π L (of order the mass of the heaviest W -boson). The mass of the lightest W -boson is m W = 2π LN and the condition for validity of our EFT, usually stated as ΛLN 1, can more precisely be phrased as Λ m W = 2π LN . All scales appearing in the long-distance theory (2.3) can be expressed in terms of Λ and L (or Λ and m W ). The scale M appearing in the kinetic term is while the nonperturbative scale m M , which determines the strength of the instanton effects generating the superpotential W (see [12] for the SU(N ) normalizations) is: In both (2.9) and (2.10), we used (2.8) to express g 2 N through the strong coupling scale Λ and L, to leading accuracy at g 2 N 1 (as indicated by the dots). The one-loop mixing matrix k ab can be diagonalized using a discrete Fourier transform (as done in unpublished work with A. Cherman), but for our purposes a numerical evaluation will suffice. As already noted, k ab from (2.7) is purely group-theoretic in nature and it 14 For use below, we only note that when the Cartan-photon mixing is accounted for, the photon-dual photon duality relation (2.2) is replaced by

JHEP01(2021)044
is easy to see that it has only negative eigenvalues, the smallest of which is about −6, for the range of N we study. This means that, upon increasing g 2 N , the eigenvalues of K ab decrease away from unity (their vanishing would signal a strong coupling singularity where our EFT breaks down). As (2.8) shows, at fixed strong coupling scale Λ, an increase of g 2 N is tantamount to increasing L. Thus, the decrease of the K ab eigenvalues shows that a strong coupling regime is approached, as expected, in the large-L limit. Conversely, a tiny value of g 2 N corresponds to small ΛLN and implies that the corrections k ab are negligible.
We shall begin our study of confining strings with values of g 2 N such that the effect of k ab is negligible, i.e. we shall first take K ab = δ ab . Later on, we shall also use the one-loop correction (2.7) as a probe of the L-dependence of the confining strings, while keeping the theory still well in (or at least on the "safe" side of) the calculable regime.
Following the usual terminology, we often call the scale m the "dual photon mass." We should, however, keep in mind that m is really the mass of the heaviest of the N − 1 dual photons, whose mass spectrum is given by m k ∼ m sin 2 πk N , k = 1, . . . , N − 1. 15 The widths of DWs are generally determined by the lightest dual photon mass.

Symmetries and vacua
Of special interest to us here are the discrete chiral and center global symmetries of SYM: 2N . This is the usual discrete R-symmetry of SYM that acts on the fermionic superpartners of (2.1) (and on the superpotential), by a phase rotation. However, of most relevance to us is that it also shifts the dual photons: 16 where we also show the action of Z 2N on the terms appearing in the superpotential (2.4). Here, ρ denotes the Weyl vector, ρ = w 1 + . . . + w N −1 . As we shall see, this symmetry is spontaneously broken to Z 2 (fermion number), leading to the existence of N ground states and DWs between them. 15 We cannot help but mention that the EFT (2.3), in addition to providing a semiclassically calculable framework to study confinement, exhibits many other, perhaps more unusual, phenomena. For example, in the abelian large-N limit [32], the 3D theory (2.3) becomes effectively four dimensional by generating an emergent latticized dimension, with many features similar to T -duality of string theory. Further, as shown in [40,41], the theory generates doubly exponentially small nonperturbative scales, ∼ e −e 1/g 2 , which are manifested in the existence of exponentially weakly bound states of the dual photons of (2.3) -the 3D remnant of 4D glueball supermultiplets. 16 Notice that Z2N acts as ZN on the dual photon and the superpotential.

JHEP01(2021)044
2. A "0-form" center symmetry Z (1),S 1 L N . The notation is chosen to emphasize that this is the dimensional reduction of the S 1 L -component of the 4D 1-form center symmetry. As explained in detail in [32,40,42], the 0-form center symmetry acts on the dual photons and their superpartners: x → Px, (2.12) e αa·x → e α a+1(modN ) ·x , (2.13) and is an exact unbroken global symmetry of the theory.
In a basis independent way, the operation denoted by P in (2.12) is the product of Weyl reflections with respect to all simple roots [42], while in the standard Ndimensional basis for the weight vectors it is a Z N cyclic permutation of the vector's components [32,40]. This clockwise action is evident from the way Z The 0-form and 1-form center symmetries discussed above are parts of the Z (1),R 4 N 1form center of the R 4 theory. As recently realized [21], SYM on R 4 has a mixed 't Hooft anomaly between the 0-form Z 2 , as in (2.14) below. Anomaly inflow on the resulting DWs implies that the DW worldvolume is nontrivial, supporting the phenomena of quark deconfinement (and braiding of Wilson lines, for DWs with three-dimensional worldvolume). Aspects of this inflow has been studied in many works [22,23,[43][44][45][46]. In the semiclassical small-L setup of this paper, anomaly inflow and the deconfinement of quarks on DWs were studied in [28]. For a recent discussion of these symmetries see also [47].
In the k-th vacuum, the dual photon field σ has a nonzero expectation value σ k , while φ is fixed at the center symmetric point: We also showed the expectation value W k of the superpotential in the k-th ground state. The N vacua (2.14) are interchanged by the action of the spontaneously broken Z symmetry is also unbroken in the bulk of SYM, corresponding to the confinement of quarks. 17 This follows from P 2πρ In words, the action of the 0-form center Z (1),S 1 L N on the vacua (2.14) is a weight-lattice shift of σ k , which is an identification, as per (2.1).

JHEP01(2021)044
There are DWs between the various N vacua. A k-wall is one that connects vacua k units apart, i.e. σ changes from σ 0 to σ k . The BPS k-walls are the ones of lowest tension among all DW connecting the same two vacua. The tension of a BPS k-wall depends only on the asymptotics of the superpotential on the two sides of the DW, W q and W q+k(modN ) , and is independent on the details of the Kähler metric, see ref. [48]. The BPS DW tensions in SYM are These BPS k-walls and their tensions play an important role in the constructions of candidate minimum-tension double-string confining configurations.

Relation to (and differences from) confinement in the R Polyakov model
Before we continue our detailed study of confinement in SU(N ) SYM, we pause to discuss the similarities and differences between confinement in SYM on R 3 ×S 1 and the well-known 3D Polyakov model. This discussion is intended for the reader who is familiar with the Polyakov model, but has not delved into the intricacies of confinement in abelianized gauge theories on R 3 × S 1 .
In this section, we outline the difference between confinement in the SU(2) Polyakov model on R 3 and SU(2) SYM on R 3 × S 1 . The discussion of SU(2) suffices to show the necessity to consider double-string configurations even for fundamental quarks in the simplest setting, without much of the Lie-theoretic technology necessary to study SU(N ) SYM.

SU(2)
Polyakov model on R 3 . The 3D Polyakov model is a nonsupersymmetric SU(2) gauge theory with an adjoint Higgs field, breaking the gauge group to U(1). Like the SYM theory discussed earlier, the long-distance physics is described by a single dual photon (in SYM, there are additional fields, the dual photon's superpartners). The dual photon is a dimensionless field, whose periodicity is in the weight lattice, σ ≡ σ + 2πw 1 , where w 1 = 1/ √ 2 is the fundamental weight of SU(2) and the simple root of SU(2) is α 1 = √ 2 (in order to have a unified notation with the rest this paper, we keep these unusual normalizations for SU (2)). Thus the U(1) charge of a fundamental quark is 18 Thus, in our normalization of the dual photons from Footnote 11, the monodromy (denoted by ∆σ nF ) of the dual photon σ around a charge equal to n times the fundamental charge (n = ±1, ±2, . . .) is: 18 Usually one takes the charge of W-bosons to be ±1 and the charges of fundamental quarks to be ±1/2; clearly the different normalization does not affect the physics. We stick with using roots and weights for SU (2) in order to emphasize the fact that our comparison of the SU(2) Polyakov model and SU(2) SYM generalizes to SU(N ).

JHEP01(2021)044
With this preamble and normalization, the long-distance effective Lagrangian, first written by Polyakov [15], is: 19 The Polyakov model (3.2) has a single vacuum at α 1 σ = 0. The periodic nature of the dual photon, σ ≡ σ + 2πw 1 , and the relation w 1 α 1 = 1 imply that the vacuum at α 1 σ = 2π is equivalent to the one at σ = 0 (as are all periodic repetitions thereof). Domain wall-like configurations interpolating between these two values of σ have a jump ∆σ = 2π √ 2 upon crossing the wall. Thus, from (3.1), we conclude that the jump of σ across the "domain wall" means that the "DW" configuration carries flux appropriate to a flux tube confining a single fundamental quark/anti-quark pair.
If one considers, within the semiclassical regime, higher charges in the Polyakov model, obtaining confining flux configurations with n > 1, as per (3.1), clearly requires double-as well as multiple-string configurations. For example, confining a massive W -boson would require two n = 1 "DWs" (as in the "gluon chain" model [49]); naturally, if stretched too long these adjoint strings will necessarily decay via W -boson pair production as discussed already in [50], see also our related discussion of the breakdown of the abelianized EFT in section 7.5.
Our main point, as we discuss below for SU (2), and in the rest of the paper for N > 2, is that in SYM on R 3 × S 1 , one has to consider double-string configurations even for fundamental quarks. This is, in fact, necessitated by the new mixed discrete chiral/center symmetry anomaly, as discussed already in [28]. SU(2) SYM on R 3 × S 1 . Let us now contrast Polyakov's effective Lagrangian (3.2) with the long distance Lagrangian of SYM on R 3 × S 1 . The latter, setting the φ field to its vanishing vev (recall (2.14)), has the form There is, however, one striking difference between (3.3) and (3.2): the factor of 2 inside the cosine potential. This innocent-looking factor is due to the nature of the confining objects -the "magnetic bions", composite instantons made of BPS monopole-instantons 19 We have shifted the vacuum energy of L Polyakov to zero and have written the lagrangian up to a normalization constant in the kinetic term (a number of order unity, not essential for our discussion); further, a power-law dependence on v/g 2 3 of the potential term is absorbed into the coefficient c. Here g3 is the 3D gauge coupling of mass dimension 1/2 and v g 2 3 is the scale of breaking SU(2) → U(1). The potential terms are due to a dilute gas of instantons (∼ e iα 1 σ ) and anti-instantons (∼ e −iα 1 σ ) of action ∼ v/g 2 3 1, as shown on the first line in (3.2).
Before we continue with the description of confining strings via (3.3), let us note that if we take L → 0, we arrive to the 3D version of SYM, rather than to the Polyakov model (3.2). The difference between the two is elucidated in, e.g. [3]: the SYM on R 3 has a "runaway vacuum", a peculiarity occuring in supersymmetric theories. To see this, one has to study the behaviour of the φ-field, the bosonic superpartner of σ, which has been set to its vev in (3.3). The expectation value of φ, which is not a compact field in the R 3 limit, as it is not related to the S 1 holonomy, is then seen to run off to infinity.
The Lagrangian (3.3) implies that SYM has two vacua in the fundamental domain for σ (σ ≡ σ + 2πw 1 ): at α 1 σ = 0 and α 1 σ = π. Notice that these two vacua are not related by a 2πw 1 shift of σ, thus they are not equivalent (however, the vacuum at α 1 σ = 2π is equivalent to the one at σ = 0 by the weight-lattice periodicity of σ). Rather, these two vacua correspond to the spontaneous breaking of the Z 4 → Z 2 chiral symmetry of SYM (2.11). The domain wall between the vacua at σ = 0 and σ = π/ √ 2 thus has a jump ∆σ = π/ √ 2. On the other hand, the domain wall between σ = 2π/ √ 2 (equivalent to the σ = 0 vacuum) and σ = π/ √ 2 also carries flux |∆σ| = π/ √ 2. (Notice that these are distinct DWs between the same vacua and the fluxes they carry are in accordance of our general results for BPS k-wall fluxes given in eq. (4.6) of the next section.) Thus, the domain walls in the SU(2) SYM theory carry electric flux ∆σ = π/ √ 2exactly half the electric flux, eq. (3.1) with n = 1, needed to confine fundamental quarks. A single DW does not have sufficient flux to confine a single fundamental quark and we arrive at the necessity to study "double string" configurations made of two domain walls, even for the confinement of fundamental quarks, in contrast with the Polyakov model. 20 It was numerically shown, already in [20], that it is double string configurations that confine fundamental quarks in SU(2) SYM on R 3 × S 1 (i.e. they do not collapse to a single flux tube). The novelty of this paper is to make use of the analytical results of [28] classifying the BPS DW fluxes (and the further insight into their properties) in SU(N > 2) SYM to study the double-string confinement mechanism for SU(N ) SYM.
Apart from finding the weakest N -ality dependence, among any known confining theory, the unexpected feature we found is the collapse of the double string confining fundamental quarks for N ≥ 5, but not for N = 2, 3, 4. This has to do with the interactions of the appropriate DWs, which changes from repulsion to attraction upon increase of N . As we discuss in section 6.1, this is not analytically understood, but has been subjected to many tests (see section 7).
We now continue with describing the details of our study for general SU(N ) SYM thoeries.

Setup
To study confinement, we add a static source term to the action for quarks of weights 21 ± λ at fixed positions r 1,2 ∈ R 2 : We denote, recalling footnote 9, the R 2 coordinates as (z, y), and take r 2 =(z=0, y=0) and r 1 =(z=0, y = R). Then, using the duality relation (2.6), we find We then add the action of the EFT (2.3), S EFT = dtdzdyL EFT , to the source action S source and introduce dimensionless coordinatest,z,ỹ, via t =t m , z =z m , y =ỹ m , to obtain, after dropping the tilde from the dimensionless coordinates 22 We numerically solve the static equation of motion 23 following from (4.3) Distance scales in the equations above are all dimensionless numbers, which are converted to physical units upon multiplication by 1/m, i.e. the (heaviest, at K ab = δ ab ) dual photon Compton wavelength, with m defined in (2.10). Likewise, the dimensionless static energy, 21 We stress that since the long-distance dynamics is abelian, Wilson loops can be classified by the weights λ. These are SU(N ) weight lattice elements, representing the charges of the quarks under the unbroken U(1) N −1 . These abelian sources can be given a fully gauge invariant description using Wilson loops representing static sources in SU(N ) representations R, extending in R 1,2 , but "decorated" by the insertion of appropriate powers of the S 1 holonomies. In the abelianized phase where the holonomy has a Z (1),S 1 L N -preserving expectation value, individual weights can be selected by a discrete Fourier transform, see [33] for details. 22 Here, ∇ denotes the gradient operator in R 2 and not a vector in the group lattice. 23 This equation implies that σ exhibits a 2πλ jump upon crossing the horizontal (y) axis in the positive z direction, for y between 0 and R. See also the original derivation in [15,51], also given in [33].
where λ 1 and λ 2 are weight-lattice vectors. Assuming that the DWs repel, a double-string configuration like the one shown will then be created. Among all DWs connecting vacua k units apart, the ones of least tension are the BPS ones. Ignoring the possibility of attractive interactions between the DWs, in the text we describe the candidate lowest tension double-string configurations confining quarks in representations of any N -ality. We argue that quarks with weights in the k-index N -ality k representation, with k ≤ N 2 , are confined by double strings made of BPS p-walls with 0 ≤ p ≤ k (we use p = 0 to denote appropriate non-BPS walls). Most of these N k double strings are metastable and are expected to decay into double strings of lowest tension, made of BPS 1-walls. The dotted line inside the double string denotes a weight-lattice discontinuity of the dual photon field, which is a gauge identification. The relation of this picture to deconfinement on DWs can be seen by noting that when the two k-walls are BPS, they have equal tensions. Thus, upon "opening up" the picture, we find that quarks can move along the DWs with no energy cost, giving a microscopic realization of deconfinement of quarks on DWs between the 0-th and k-th vacua, as argued for by anomaly inflow. is converted to the physical energy E = MẼ. A dimensionless string tension (or a DW tension (2.15))T computed from the simulation, is converted to physical units as T = M mT . When we study the effect of changing L at fixed Λ (i.e. changing g 2 N in (2.5) according to (2.8)), being mindful of these scalings and the definitions of M and m in terms of L and Λ will be important.

BPS DW fluxes and double-string configurations
The double-string confinement mechanism is illustrated on figure 2 and is explained in more detail in the caption. As already alluded to, the confining string consists of two DWs between the N vacua of SYM. Its essence is that due to confinement, the chromoelectric flux of color sources can only propagate along one dimensional lines, the DWs of SYM theory. As explained at length below, the fluxes carried by the DWs are such that one needs two DWs in order to absorb the quark fluxes.
We shall consider quark sources in different representations of SU(N ) and discuss what kind of double-string configurations are expected. But first, we remind the reader of the results of [28], which are essential for our arguments. In that paper, the electric fluxes carried by BPS k-walls in SYM theory were determined. It has long been known, see [48], that there are N k BPS walls between vacua k units apart. Ref. [28] showed that the BPS k-walls fall into two classes, labeled by I and II below, which have monodromies of σ, or electric fluxes Φ k,I i 1 ,...,i k and Φ k,II i 1 ,...,i k−1 , given by: The total number of BPS k-walls is, as already stated, We stress that this is because the weights ( w i 1 , . . . , w i k ) are to be taken all different, with indices ranging from 1 to N − 1; likewise all weights in the second set, ( w i 1 , . . . , w i k−1 ), are to be taken different. 24 The above spectrum of BPS k-wall fluxes is invariant under k → N − k up to reversal of the overall sign of the electric flux (a parity transformation, as in [44]).

Quarks in the highest weight of the k-index antisymmetric representation
We shall now use the results for the BPS wall tensions to propose double-string configurations confining quarks of various weights. We begin by first concentrating on quark sources of N -ality k and of weights w k -the highest weight of the k-index antisymmetric representation. One result (ignoring DW interactions) that we can immediately obtain is that quarks of weights w k can be confined by double-string configurations made out of BPS 1-walls only, for all k. This is because the fluxes of BPS 1-walls are 2π as per (4.6). As shown on the figure, the total monodromy is Φ 1,II − Φ 1,I k = 2π w k , as required for a string confining quarks of weight w k . Since the BPS 1-walls are the ones of lowest tension (recall (2.15)), we expect that the minimum energy confining configuration for quarks in an N -ality k representation will be due to the double-string configuration made of two BPS 1-walls. To

JHEP01(2021)044
determine the actual N -ality dependence of the string tensions, one requires a numerical evaluation of the energy of the double-string configuration at large quark separation, which will incorporate the domain wall interactions.
For N -ality k quarks of weights different from w k , we shall argue below that it is not possible (for every weight) to engineer a double-string configuration made out of BPS 1-walls. In fact, as we shall see below it is not possible to construct a double-string configuration made of BPS p-walls for all weights of a given N -ality. Thus, non-BPS DWs become important for some weights. However, since any N -ality k weight is related to the highest weight w k by the addition of root vectors, it is clear that at sufficiently large quark separation, any such double-string configuration made of non-BPS walls (or of BPS p-walls with p > 1) will decay to the BPS 1-wall double-string configuration by the creation of W -boson and superpartner pairs. The screening of any N -ality k weight down to w k is possible because W -boson charges take values in the root lattice.
The fact that the "long distance" EFT (2.3) describing an abelian confining theory breaks down when distances between probe quarks are taken sufficiently large has been noted before: see [1] and [52] for a more recent discussion and references. The screening of static sources by the pair creation of massive "W -bosons" is responsible for "restoring justice" and ensuring that asymptotic string tensions are only N -ality dependent. The importance of this screening, and the associated EFT breakdown, is not new and holds for other models with abelian confinement, as discussed for Seiberg-Witten (SW) theory [31] and deformed Yang-Mills (dYM) theory [33].
Below, we consider in some detail how quarks of N -ality k, but of general weights, i.e. weights different from the fundamental weight w k are confined.

Quarks of general weights of N-ality k = 1
Fundamental quarks of any weight ν A , A = 1 . . . , N , of the fundamental representation (with ν 1 = w 1 ) are all confined by BPS 1-walls only and have the same tension. As we stress in (4.8) below, this is due to the unbroken 0-form center symmetry (which holds also in dYM but not in SW theory). It can also be seen from the fact that ν A = w A − w A−1 , with w 0 = w N = 0, that the difference between two appropriately chosen BPS 1-wall fluxes from (4.6) is always a weight of the fundamental representation. Explicitly, with Φ 1,I A taken from (4.6), we have that showing that for all weights of N -ality 1, quarks are confined by double strings of lowest tension made of two BPS 1-walls (when A = 1 or N , one of the weights must be taken Φ 1,II from (4.6)). Note that we could have also used BPS 2-walls of fluxes Φ 2,II A (from the second line of from (4.6)) instead of Φ 1,I A to engineer a double-string configuration confining fundamental quarks; however, these are expected to have higher tension due to the higher tension of the BPS 2-walls.
Finally, we stress that the fact that all string tensions of N -ality 1 are the same is because the unbroken Z : and relates the strings confining the different weights of fundamental quarks to each other. This Z N symmetry action (4.8) will be important in our discussions of higher N -alities below.

Quarks of general weights of N-ality k = 2
We next consider the two-index representation obtained by taking two fundamental quarks close to each other. The weights of the resulting two-index representations are ν A + ν B .
There are a total of N 2 such weights split between the two-index symmetric and antisymmetric representations (which share many weights). Let us arrange the N 2 weights into an N × N matrix with entries ν A,B ≡ ν A + ν B (for use below, we take the indices to cyclically "wrap around", i.e. we take ν N +1,l ≡ ν 1,l and ν l,N +1 ≡ ν l,1 ). Further, because of the Z N action (4.8), it is sufficient to restrict to A = 1, since all strings confining quarks of weights ν A,A+p(modN ) , p = 1, . . . N , are related by a symmetry and so should share the same properties, including equal tensions. Thus, we shall say that the weights ν A,A+p , with A = 1, . . . , N , form a Z N orbit and shall focus on its representative ν 1,1+p . Clearly, as p = 1, . . . N , there are N such orbits for the k = 2 two-index representations. Consider first the p = 1 nearest off-diagonal Z N orbit of ν 1,1+1 . Quarks with these weights are all confined by double strings made of BPS 1-walls; one can see this using an expression essentially identical to (4.7). This is because ν 1,2 = w 2 − w 0 = w 2 (recall we defined w 0 = 0) is proportional to the difference between two BPS 1-wall fluxes from (4.6), Φ 1,II − Φ 1,I 2 . Further, we take p > 1 and consider the p-nearest off-diagonal Z N orbit of ν 1,1+p with p = N . Now, we argue that the lowest tension double strings for p > 1, p = N , can be constructed of two BPS 2-walls. To see this, we note that the difference between two BPS 2-wall fluxes from (4.6). The only two-index N -ality 2 case left to consider is the p = N diagonal weights in the Z N orbit of ν 1,1 = 2 w 1 . Because the weights determining the fluxes of BPS walls on the right-hand side of (4.6) are all to be taken different, one cannot arrange this weight to be a difference of two BPS wall fluxes (our numerical simulations show that quarks of the diagonal weights are confined by non-BPS double-string-like configurations). The N -ality k = 2 strings lying in different Z N orbits are studied in more detail in section 7.4. As already stated, we expect that for 1 < p ≤ N all N -ality two strings of higher tension are metastable and decay to the lowest tension double strings made of BPS 1-walls, the ones in the Z N orbit of w 2 , by the production of W -boson pairs. 25 25 We also recall that adjoint quarks are also confined in the abelian regime. Consider a W -boson whose weight is a simple root αa = νa − νa+1 = 2 wa − wa+1 − wa−1. This cannot be written as the difference of two BPS k-wall fluxes, for any k, hence one expects that other configurations will be responsible for confinement of non-Cartan gluons in the abelian regime. See section 7.5, where such a triple-string configuration confining a W -boson of weight α1 is found.

Quarks of general weights of N-ality 2 < k ≤ N 2
We can now similarly consider k-index representations of N -ality k made by putting together k fundamental quarks. 26 Clearly, now there are N k weights. Similarly to our discussion of k = 2 above, we label them as ν a 1 ,a 2 ,...,a k = ν a 1 + . . . + ν a k , all a i ∈ {1, . . . N }. (4.10) The highest weight of the k-index antisymmetric representation is w k and it lies in the Z N orbit of ν a,a+1,a+2,...,a+k−1 , where the a = 1 element is w k = ν 1,2,3,...,k , as implied from the above definitions.
The N k weights (4.10) of the k-index N -ality k representations split into Z N orbits, labeled by k − 1 integers n ≡ (n 1 , . . . , n k−1 ) all from 1 to N , a+n 1 ,a+n 2 ,...,a+n k−1 , a = 1, . . . , N, (4.11) where a = 1 and n = (1, 2, 3, . . . , k−1) corresponds to w k . Again, without loss of generality, we can focus on the a = 1 element of a given Z N orbit, for which we can write, from the definition (4.10) and recalling (4.7): For generic 27 weights, the above equation implies that meaning that for such weights of an N -ality k k-index representation, quarks are confined by double strings made of two BPS k-walls whose fluxes Φ k,II n 1 ,...,n k−1 and Φ k,I 1,1+n 1 ,...,1+n k−1 , as per (4.6), combine to provide the monodromy required by (4.12).
An extreme non-generic case was already seen: for n = (1, 2, 3, . . . , k − 1), quarks are confined by double strings made of BPS 1-walls, instead of BPS k-walls. More generally, if some of the numbers n i are neighboring, there are cancellations between the two groups of weights on the right-hand side of (4.12) reducing the "rank" of the BPS walls required to construct a double string of the right monodromy. Like in the k = 2 case, there are also quarks (e.g. n = (N, N, . . . , N )) that cannot be confined by double-string configurations made of BPS-walls.

Numerical simulation of double-string configurations
To study double-string confinement, we numerically solve the static equation of motion (4.4), reproduced here with the minor modification of symmetrizing the locations of the two quarks: While focusing on the lowest string tension of a given N -ality, it suffices to consider N -alities k ≤ N 2 . Strings of N -ality N 2 < k < N have tensions identical to those of strings of N -ality N − k, since the BPS wall tensions (2.15) obey T k = T N −k . 27 Such that the numbers n are all distinct, non-neighboring, and not 1, N − 1, or N .

JHEP01(2021)044
It is often of interest to also solve for the associated one-dimensional domain walls, whose static eqation of motion can be thought of as discarding the source term ( λ → 0) and eliminating one of the spatial dimensions (∇ 2 → d 2 dz 2 ) in (5.1): Note that both equations should also satisfy the boundary condition that x is equal to one of the vacua at each connected piece of the boundary at infinity.
In what follows, we will describe the numerical methods used to solve these two equations 28 efficiently. The methods are mostly identical to those described in sections 4 and 5 of [28], with some extensions and modifications. Whenever there is overlap, we summarize the main points needed to reproduce the work here, but refer the interested reader to the aforementioned paper for detailed description.

Gauss-Seidel finite-difference relaxation
The core method used is the Gauss-Seidel finite-difference method [53]. For the 2D EOM, we first discretize our space into a rectangular grid 29 with "pixels" being squares of size h × h. In all cases, we have used h = 0.1 in units of the dual photon mass m.
Next, the Laplacian is converted to a finite difference using the centered-difference approximation: (5.4) We solve this large system of coupled linear equations using the Gauss-Seidel method: start from some initial field configuration, which can be any random 30 guess, with the only constraint being that the boundary of the grid, which we approximate to be spatial infinity, must equal one of the vacua (2.14) of the theory. Then, we iteratively update the value of the field at each point of the grid (with the exception of the boundary points) by the average of its four neighboring points, with a Laplacian term subtracted: (5.5) If the process converges, then x a new ≈ x a old , and so (5.5) becomes the same as (5.4). Since the boundary points were not touched, the boundary value problem is solved. Once the 28 From this point on, we will refer to (5.1) and (5.2) as the 2D EOM and 1D EOM, respectively. 29 The dimension of the grid in the direction of the quark separation must obviously be larger than R. In addition to this constraint, care should be taken to ensure that the grid size is large enough to accommodate the non-constant part of the solution, usually by trial-and-error. However, the edge effect is in general not a big concern in a confining theory such as ours, in contrast with a non-confining theory; see appendix A. 30 Literally, as we will soon see.

JHEP01(2021)044
solution is obtained, we can easily numerically compute various quantities of interest, such as numerically integrating 31 to get the static dimensionless energy (4.5). Convergence can be measured by an error function, 32 which we have chosen to be The notation here refers to the maximum value among all points and all vector components. We have experimentally found an error tolerance of e < 10 −9 to be a good convergence criterion, which was used in producing all the results in this paper. However, it is computationally costly for the standard Gauss-Seidel method to comply with such a small tolerance, especially when the computation must be repeated many times as in the present case. Two methods for speeding up the Gauss-Seidel process that are specific to the present problem are provided in appendix B. Finally, to implement the right-hand side of (5.1), one needs to use the numerical derivative of Dirac delta function, given by For a derivation, see section 5 of [28].
As for the 1D EOM, the exact same method suitably adapted to one dimension also works. We will only make two comments about the 1D case and refer the interested reader to section 4 of [28] for more details. First, unlike the 2D case, the 1D EOM has two disconnected pieces of boundary (namely, z = ±∞), and the nontrivial solutions are the DWs interpolating between two different vacua. Secondly, upon obtaining a solution to the 1D EOM, one can check whether it is BPS (i.e. whether it saturates the lower bound on energy) by either comparing its numerical energy to the analytic expression (2.15) or checking that a plot of its first derivative matches the BPS equation [48]:

Initial configurations
Although the initial configuration can be arbitrary, a good choice that is close to the real solution can speed up the convergence process. So we use the "BPS initial configuration," which we now illustrate with the most frequently-used example. Recall from the discussion in section 4 that quarks of weights w k are expected to be confined by double-string configurations made of two BPS 1-walls, for any k. We will use the same fluxes as the previous discussion to achieve the monodromy, except that for ease of numerical implementation we shift the overall field by i 2π N ρ. Concretely, we impose the boundary condition x| ∞ = i 2π N ρ, and let the central region between the two quarks have field values w N = 0 and w k , separated by a field discontinuity 31 While integrating, care must be taken to avoid including the discontinuous jump inside the double string, introduced by the monodromy (see figure 5). For example, one can use one-sided derivatives. 32 Note that this is different from the error function used in [28]. In between the quarks are many copies of pairs of BPS DWs, interpolating from the boundary value to 0 and from i2πw 2 back to the boundary value. The field discontinuity has a value that matches the monodromy and its location and shape match the Dirac delta function in the 2D EOM (5.1). We have used Kähler metric without quantum correction. Using the Gauss-Seidel iterations, this initial configuration relaxes to figure 5. along the y-axis, 33 consistent with the location of the Dirac delta term in (5.1). We then solve for the pair of BPS DWs 34 which interpolate between i 2π N ρ → 0 and i2πw k → i 2π N ρ. To construct the initial configuration, we cover the region between the two quarks by many pairs of BPS DWs, while leaving the external region at the constant boundary value. We show a typical example in figure 3.
There is a valid concern regarding the choice of initial configuration that needs to be addressed. In principle, there could be more than one solution that minimizes the action. Then the numerical solution might depend on the choice of initial configuration. We argue that this is not the case by testing against the most arbitrary initial condition possible: a completely random field, as shown in figure 4. When the Gauss-Seidel method is applied with a tolerance of 10 −9 , both this random field 35 and the BPS initial configuration (figure 3) relax to figure 5, with the exact same numerical energy.

A typical example of a double string
Finally, we illustrate what a typical double-string solution looks like, before diving into the various interesting results we learned from mass-producing such computations. Note that 33 In all diagrams in this paper, the y-axis is drawn to be horizontal, while the z-axis is vertical. 34 Note that there is a freedom to choose the positions of the two kinks, which corresponds to a guess of the string separation. There is no way to know this value apriori. 35 Readers who wish to reproduce this random initial configuration study should note that it takes expo-  although the double-string confinement prediction [20] was found to hold for most quarks, the fundamental quarks break this picture, as we shall see in the next section. We present the result for the gauge group N = 5, with a quark-antiquark pair of weight ±λ = ±w 2 , separated by a distance R = 15, all while ignoring the one-loop correction in the Kähler metric (2.5). Figure 5a and 5b show the real and imaginary component of the third field, respectively. Compare this with the schematic prediction of figure 2. We emphasize again that the field discontinuity at the y-axis is completely consistent with the location of the Dirac delta term in (5.1) and matches the monodromy of the quarks. Figure 5c and 5d show the potential and gradient energy density, which are the first and second term in the integrand of (4.5), respectively. The total dimensionless energy is just the integral of the sum of these two pictures. As is clear from the scale, the total energy density is dominated by the gradient term, but a plot of the potential energy density most clearly reflects the structure of a double string, while that of the gradient energy is concentrated at the two quarks.
At the center of the double string, the two quarks are relatively far away, and so recalling the discussion of the relation between (5.1) and (5.2) at the beginning of section 5, one would expect that a cross-section corresponds to the solution of the 1D EOM. As a consistency check, we take a vertical cross section of the solution at y = 0 and compare it to a plot of the corresponding pair of BPS solutions in figure 6. To a first order approximation, they do match quite well.

Properties of confining strings
Upon mass production of the above computations, various interesting properties of confining strings were found. We will first present all results without the quantum correction, i.e. ignoring the second term in the Kähler metric (2.5), and only afterward describe the effect JHEP01(2021)044 of the quantum correction. Additionally, we are mainly interested in quarks of weights w k . As explained in section 4.2.1, quarks of N -ality k are screened down to w k by the creation of W -bosons and superpartner pairs. Hence, we will often concentrate on these quarks without further notice.

Collapse of the double string for fundamental quarks
The picture of double-string confinement proposed in [20] was found to hold in most cases. However, it breaks down for fundamental quarks. To describe this, it helps to separate theories with N ≥ 5 and N < 5, because it turns out the latter suffer from low-rank peculiarities.
For N ≥ 5, quarks of weights w k with k = 1 have double-string solutions 36 just like that of the "typical example" in figure 5. But for w 1 , the two strings attract and form a 36 As far as we could tell from simulations up to N = 10. bound state. In other words, the quarks are confined by a single string, as shown in figure 7 for N = 8. The cross section appears to be some nonlinear merger of the corresponding pair of BPS solutions, as shown in figure 8. This is in contrast with the double-string case for which the cross section (near the midpoint between the quarks) is almost identical to a pair of BPS solutions, as shown in figure 6.

JHEP01(2021)044
We checked that this collapse is not just a small R effect by going to as far as R = 200, see figure 9. Stronger evidence comes from the fact that the attractive interaction between JHEP01(2021)044 two DWs of total flux w 1 can even be seen in 1D solutions. A comprehensive study of DW interactions in 1D is performed in section 7.
For N < 5, this collapse does not occur and w 1 quarks are confined by double strings. We emphasize here that the above general assertions of whether the strings collapse only pertain to large quark-separation behavior. For instance, for N = 4 and small quark separation R, the DWs of w 1 quarks do feel the attraction that would eventually become insurmountable for their higher N cousins; but for large R, the repulsion wins over and a double-string is formed. The gradual switch starts around R = 17, and the progression is shown in figure 10. Similar gradual repulsion as R increases also occurs in the high N -ality cases for N > 8.
We currently do not have an analytical understanding of the reason for the collapse or of why unit N -ality is special starting at N = 5. However, we noticed that this non-generic behavior of DWs in low dimensions and a turning point at N = 5 has a precedent. It was found in [28] that BPS DWs can have certain components of their dual photons possessing a "double-dipping" profile; σ 1 for both of the BPS DWs in figure 6 is a good case in point. This physically means that the corresponding double string has most of its electric flux flowing from a quark to its antiquark, but along a small portion of the DW cross section, the electric flux flows backwards from the antiquark to the quark (of course the total flux is fixed by the quark charge). Crucially, this is a generic phenomenon for N ≥ 5 but does not occur for N < 5. We suspect there are some connections or common causes between the field-reversal phenomenon and the double-string collapse. Perhaps these hints can help future investigations along this direction.

N-ality dependence of string tensions
It is expected that for large quark separation R, the static energyẼ (4.5) grows linearly as a function of R. Furthermore, after the effect of screening is taken into account, the string JHEP01(2021)044 tensionT which is the slope of anẼ vs R graph, should depend only on the rank of the gauge group N and N-ality k of the quarks: Since N -ality k quarks are screened down to weight w k , we can numerically test (6.1) and compute the string tensions using simulations of the confinement of w k quarks. We indeed found the asymptotic linear relation to hold for all cases tested: for N up to 10 and for quarks of any N -ality. We also numerically extracted the tensionsT (N, k) using a least-squares fit. Typical plots for k = 1 and k = 1 are shown in figure 11. To a first order approximation, the general features of the values of the tensions can be explained using what we already know. For quarks of weight w k with k = 1, we know that they are confined by double strings made of two BPS 1-walls. Of course, what we mean by this statement is that far away from the quarks, such as at the center of the double string, the one-dimensional cross section matches the shape of two BPS 1-walls (figure 6), whose combined monodromy produces the weight w k .
In fact, from figure 5 and our other numerical simulations, it appears that the strings look like (rotated/deformed) 1-walls for most of their length, aside from the area immediately adjacent to the quarks. If we make the approximation that the string tension per unit length at every point on the string is the same as that at the center and ignore the curvature of the string altogether, then we would expect the overall dimensionless string tension to be roughly twice the BPS 1-walls tension, 2T 1-wall = 4N sin π N , as in (2.15), for the double-string configurations that do not collapse. On the other hand, for fundamental strings, w 1 , we have seen that the double string collapses to form a bound state, so we expect its string tension to be less than 2T 1-wall . Indeed, just like the case for SU(10) in figure 11, all k = 1 tensions are less than 2T 1-wall , reflecting the nonzero binding energy, and all k = 1 tensions are slightly above or below 2T 1-wall .

JHEP01(2021)044
In addition to the k-string tensions, of particular interest to us is their unitless ratio This quantity is considered a probe of the confinement mechanism and has been compared among different toy models of confinement as well as to lattice simulations. Some wellknown scaling laws from other models are the Casimir scaling (due to one-gluon exchange, or in the strong coupling limit, see e.g. [1]), the square root of Casimir scaling (first found in the MIT Bag Model [54]; it was found to also approximately hold in deformed Yang-Mills theory on R 3 × S 1 , see [33]), and the Douglas-Shenker sine law (which holds in softly broken Seiberg-Witten theory [31], MQCD [55], and in some two-dimensional models [56,57]), 37 .
With these previous results in mind, we present our numerical results ofT (N, k), f (N, k), and the uncertainties arising from their least squares fit in table C.1 in appendix C. Note thatT (N, k) =T (N, N −k) due to charge conjugation symmetry, and so we only need to compute N -ality up to N 2 . A plot of f (N, k) is shown in figure 12, and a comparison between the present theory and other known scaling laws for SU(10) is shown in figure 1. Among the known scaling laws, the values of f (N, k) are closest to the square root of Casimir scaling, but the curve is much flatter. In fact, the curve is relatively flat apart from a large jump between the k = 1 and k = 2 points. This is because the tensions for k = 1 deviate only slightly from 2T 1-wall as compared to k = 1, which is a collapsed string. Because of the lack of string collapse, the f (4, k) plot (not shown) is the only special case 38 where the curve almost stays constant at 1.

Growth of the string separation
The string separation d was predicted to grow logarithmically as a function of R in [20]. The argument goes as follows. Approximate the double string to be composed of mostly parallel DWs separated by a distance d, and take the double-string configuration to be a rectangle of dimensions R and d. The energy must have a term proportional to the string 37 Ref. [58] found yet another different scaling law, similar to (6.5), a sine-squared one which holds for the tension of strings wrapped around the S 1 , inferred from the correlator of two Polyakov loops, at small gaugino mass m in softly broken SYM on R 3 × S 1 . This theory is conjectured to be continuously connected to pure Yang-Mills theory [9]. 38 N -ality dependence only makes sense starting at N = 4. tension per unit length (taken to be M m, recall section 4.1) times the length of the double string, which is ∼ 2M m(R + d). The one-dimensional DW repulsion 39 per unit length of the DW can be modelled by an exponential term ∼ M me −md , and the corresponding term in the energy must be proportional to M mRe −md . Thus, a larger d increases the energy due to the increased length of the string, but is compensated by a smaller repulsion energy. Up to numerical factors, the energy is

JHEP01(2021)044
The minimum occurs at a separation d that scales as From our simulations, all of the non-collapsing double strings indeed have a logarithmic growth in d. We define d to be the distance between the two maximum points in a y = 0 cross-section of the potential energy density. A typical fit to a logarithm is given in figure 13. A list of all the parameters fitted to a logarithmic description of d is attached in table C.2 in appendix C. Despite the success of the logarithmic growth prediction, this simple model does not capture the N -ality dependence of string separation, which turns out to be significant. The separation generally gets larger as N -ality increases, an effect that becomes obvious for large N , and most strikingly for N = 10, as shown in figure 14.  (5), k = 2). The step-like nature of the data is entirely a numerical side-effect: the pixels of the discretization have size h = 0.1, while R is increased in intervals of 1. When log(R) grows slower than 0.1 upon an increase of ∆R = 1, the width will have a step-like growth.
(b) k = 5. Figure 14. N -ality dependence of string separation. In SU(10), the string separation is significantly larger for k = 5 than k = 2 (compared at the same R, of course).

Towards larger L: the effect of W-boson loops
As discussed in section 2, in the leading semiclassical approximation, up to now we neglected the one-loop correction to the Kähler metric (2.5). In this approximation, all our results for the string tensions are expressed in dimensionless form, as given in (2.15), in terms of the parameters m and M defined in (2.10) and (2.9). In this section, we shall describe the results of simulations that take into account the one-loop correction to the Kähler metric. This introduces a dependence of the simulation results on the dimensionless parameter: should recall that M and m also depend on . To take this into account, we use (6.8), the definitions of M and m, and the definition of Λ, (2.8), now rewritten as to rewrite these dimensionful parameters and the dimensionful string tension as functions of , Λ, and N . Our goal is to keep Λ fixed for theories with different N and at different circle sizes L. It is clear from (6.9) that increasing is equivalent to increasing L at fixed Λ. Thus, we find 40 T ( ). (6.10) We study two values of : 1 = .09 and 2 = .12. These are chosen to be small enough that the Kähler metric does not have a strong coupling singularity for the values of N we study, yet large enough to produce visible deviations from the simulations with K ab = δ ab . From (6.9), these values of correspond to S 1 -sizes L 1 and L 2 such that L 1 Λ 4π 0.0086, for 1 = .09, and L 2 Λ 4π 0.0314, for 2 = .12, (6.11) corresponding to a roughly fourfold increase of L, L 2 3.6L 1 . We also notice that for up to N ≤ 10, the semiclassical expansion parameters ΛL 1 N 2π and ΛL 2 N 2π are smaller than unity, though hardly infinitesimal -the latter is about .6 for N = 10. Thus the validity of the semiclassical approximation is at best qualitative. As will be clear from our subsequent discussion, the string tensions are seen to not vary wildly as a function of ; this small variation can be taken to qualitatively support the use of semiclassics.
The results for the string tensions and the fitted string separations for these two values of are given in tables. (C.3) and (C.4), respectively. Before discussing them, let us make some comments on the simulations with = 0. First, we note that the equation of motion (5.1) and the behavior of the Kähler metric (2.5) discussed earlier imply that, after including W -boson loops, the BPS walls are expected to become wider. This is because the eigenvalues of K ab are less than unity and the equation of motion (5.1) implies that as the eigenvalues of the Kähler metric decrease, the second derivative of x is also smaller, and thus it needs to vary over a larger range to interpolate the same flux. This is indeed seen to be the case in our simulations (in addition, we have verified, as a check of their consistency, that while the 1D DWs become wider as is increased, the BPS DW tensions (2.15) remain the same, as they only depend on the boundary conditions and not on the Kähler metric). Second, as is seen from the fitted separation of the double strings given in (C.4), this widening of the BPS DWs, and subsequently of the double strings, is the main qualitative effect of the inclusion of W -boson loops. 41 At the same time, we find that including nonzero does not 40 The scaling of the string tension T in (6.10) in terms of parameters is usually written in a more familiar way that somewhat obscures the dependence: T = .  (9) f(9, k) = T (9, k) T (9,1) = 0 (L < L 1 ) L 1 ( = 0.09) L 2 3.6 L 1 ( = 0.12) Figure 15. The k-string tension ratio for SU (9), as a function of L at fixed Λ, for three different values of L, increasing from the top curve towards the bottom one. Clearly, as L increases, the k-string ratios are seen to decrease. Qualitatively, however, the curve near its maximum is not as flat as for = 0. This effect is entirely due to W -boson (and superpartner) loops. The values of the semiclassical expansion parameters corresponding to the chosen are discussed after eq. (6.11). lead to collapse of the double string for k > 1, or to a separation of the collapsed string for k = 1. In other words, the main conclusion from our results is that the qualitative behavior of double-string confinement is not changed by the inclusion of virtual W -boson effects.
Let us now discuss the results in some more detail. The string tensions T ( , Λ, N ) increase as (or L) is increased at fixed Λ: as (6.10) shows, increasing increases the prefactor, and numerics shows that the dimensionless partT ( ) also increases. The results forT ( ), shown in eq. (C.3) for up to N = 9, show thatT ( ) slightly increases as changes from .09 to .12. Thus, the increase of the dimensionful string tension at fixed Λ is dominated by the increase of the prefactor in (6.10). The effect of the quantum correction is small, as expected in the semiclassical regime. As our results show, this smallness holds even when the expansion parameter is not so small (e.g. for = .12). Thus, as expected from (6.10), the k-string tensions at fixed N and Λ increase with L, and are essentially proportional to L (i.e., they increase roughly 3.6 times), with small variations due to W -boson loop effects.
However, the behavior of the k-string tension ratio, f (N, k), as a function of L cannot be deduced without knowledge ofT ( ) of (C.3), as the prefactor in (6.10) cancels when taking the ratio. Even before discussing these data, we can make the following qualitative especially as N is increased. Perhaps, the results of (C.4) will be useful for a future model of the double string properties. (We stress, however, that the separations in (C.4) are based on only a few data points for the larger values of N , because the strings tended to collapse at smaller quark separations. For N = 9, = 0.12, for instance, the string is collapsed all the way to R = 40.)

JHEP01(2021)044
statement. If the double-string picture remains intact (as it is seen to), we expect that the string tension is approximately twice the BPS 1-wall tension. Then the Kähler potentialindependence of the latter implies that the approximately flat behavior of f (N, k) as a function of k will persist. This is indeed seen from the results of (C.3). The result which is difficult to predict without simulating is that f (N, k) varies even less with k as L increases between the two values shown in (6.11). In figure 15, we show the k-string tension ratio for the largest value N = 9, for the two values of L of (6.11) as well as for = 0. 42 One qualitative observation is that while the range of variation of f (9, k) is smaller for larger L, the f (9, k) appears to have more curvature near its maximum.
To conclude this section, we find that the qualitative double-string picture is not changed by including the leading virtual W -boson correction, upon increase of L within the semiclassical regime. The main qualitative effect of the one-loop correction is the widening of the double string. At the same time, the range of variation of the k-string tension ratio f (N, k) decreases upon increase of L.

Domain wall interactions in one dimension
In order to elucidate the behavior of the double-string configurations observed in two dimensions, we now perform a study of domain wall interactions in one spatial dimension. The morphology and behavior of individual BPS solitons in 1D was studied extensively in [28]. Here, we extend these results with a numerical study of the interaction of multiple DWs. In this section, we set the Kähler metric to K ab = δ ab .

Interaction energy of two BPS domain walls
As observed in our 2D numerical study, quarks of N -ality k = 1 are confined by a single string for N > 4, with a string tension lower than that of a BPS 1-wall. We may verify that this behavior is consistent with the 1D picture by considering the interaction energy of two parallel 1-walls as a function of their spatial separation. To do so, we first numerically solve for two BPS 1-walls with electric fluxes Φ 1,I k = 2πρ/N − 2πw k and Φ 1,II = 2πρ/N according to (5.2), and using the notation for fluxes from (4.6). Note that these match the fluxes of our double-string configuration with weight w k . We then sum the two walls, separating them by a certain distance d as measured from the centers of the kinks. Here, the "center" of a 1-wall is taken as the peak of the φ fields. An example of this configuration for N = 5 and k = 1 is shown in figure 16, for two different values of d. The total energy of this configuration as a function of the separation d can then be computed via (4.5). This gives an idea of the long-distance (compared to their size) interactions of these DWs, as in the "sum ansatz" used to study instanton-anti-instanton interactions, see e.g. [59] for a review.
For several different values of N and k, we computed the energy of these interacting DWs for kink separations d between 0 and 10. These energies are depicted for N = 4, 5 in figure 17. Here, we observe that DWs of N -ality k = 2 have a tendency to repel at small kink separations, as expected from our 2D simulations. In contrast, for N -ality k = 1, DWs exhibit either weak repulsion 43 in the case of N = 4, or slight attraction in the case of N = 5. For all N and k, the energy converges asymptotically to twice the BPS 1-wall tension as the kink separation is taken to infinity. These results hold for all other gauge groups considered in this study: N -ality k = 1 remains attractive at small kink separations for N > 5, and repulsive for N < 4; similarly, for all N -alities k > 1 and all N , the DWs are repulsive.
These results allow us to put the behavior of the 2D confining strings into greater context. As noted above, the repulsion of the DWs for k > 1 is consistent with the emergence of the double-string configurations observed in the 2D simulations. Moreover, the attraction of the DWs for k = 1 is consistent with the observed collapse into a single-string configuration for N ≥ 5 and k = 1. Notably, the behavior for N = 4, k = 1 is more complicated. Although these DWs seem to be slightly repulsive at small kink separations, in our 2D numerical simulations we observed these walls collapse into a single string at smaller quark separations. It is possible that since these walls are merely weakly repulsive relative to the JHEP01(2021)044 Twice the energy of a BPS 1-wall, as given by the dimensionless form of (2.15), is also marked by a dashed line. At large kink separations, the energy of the summed configuration approaches twice the BPS 1-wall energy asymptotically. At small kink separations, we can see from the slope of the energy curves that DWs of N -ality k = 2 tend to repel, whereas DWs of N -ality k = 1 are either only weakly repulsive in the case of N = 4, or slightly attractive in the case of N = 5. Note that these results hold for larger and smaller studied gauge groups: N -ality k = 1 continues to be attractive for N > 5 and repulsive for N < 4, whereas all other N -alities are repulsive for all N studied. k = 2 configuration, there exists some other property of the confining string for which a single string is energetically preferable at small quark separations. For example, morphological properties such as the additional string length and/or curvature of the double-string configuration may make a single string the preferred solution at these small quark separations.

Attraction and repulsion of domain walls
We further investigated the attraction and repulsion of DWs by numerically relaxing configurations consisting of the sum of two DWs similar to those depicted in figure 16 (this way of studying DW interactions is motivated by the "streamline" or "valley" method of [60], see also [59]).
An example of this for N = 5 and k = 1, 2 is shown in figure 18. Here, a configuration of two DWs of N -ality k = 1 is seen to attract to form a single string after relaxation, whereas a similar configuration of N -ality k = 2 is seen to repel after a similar procedure. In general, we observe that the attraction or repulsion of DWs configured as described above is consistent with the interaction energy of the walls studied in the previous section by the "sum ansatz." For example, configurations of N -ality k = 1 are seen to repel for N < 4 and attract for N ≥ 5.
For those DW configurations which are observed to repel, the solution resulting from numerical relaxation is seen to closely correspond with the simple sum of the φ and σ fields of the individual 1-walls taken a distance apart. For those configurations which attract,  figure 16. These initial conditions were relaxed for up to 100,000 iterations of the relaxation algorithm, with the usual error tolerance of 10 −9 . The φ fields for these solutions are depicted in the right column. In the top row, an N -ality k = 1 configuration of two DWs separated by 10 units is seen to attract to form a single string. The resultant φ field is comparable withbut not identical to -the sum of the φ fields of the individual 1-walls, having a smaller magnitude and an additional central divot as compared to the summed field (compare with the lower-left panel of figure 16). Although not depicted here, the σ field of the relaxed solution appears as the sum of the σ fields of the individual 1-walls, as depicted in the lower-right panel of figure 16. In the bottom row, an N -ality k = 2 configuration with an initial separation of 0 is seen to repel to form a double string. In this case, the resultant φ and σ fields both appear as the simple sum of the individual 1-walls, reminiscent of the top row of figure 16. although the σ fields of the 1-walls appear to simply sum together, the φ fields behave differently. The φ fields of these "collapsed" solutions are reminiscent of sum of the two BPS DWs, but have a smaller magnitude and posses an additional "divot" at their centers (for instance, compare the lower-left panel of figure 16 with the upper-right panel of figure 18). We make no attempt to explain this phenomenon in this work, other than to point out the obvious fact that these DW interactions are highly nonlinear. We are fortunate, therefore, that the σ fields of these DWs appear to interact linearly or almost linearly. 44 44 Scalar theories with a complex potential equal to our superpotential W , the affine Toda theories studied JHEP01(2021)044

Comparison with confining strings
We also compare the profile of these one-dimensional DWs with the y = 0 (midpoint between the quarks) cross sections of our confining string configurations. Although these cross sections are discontinuous at z = 0 due to the imposition of the monodromy, we can remove this discontinuity by translating the σ field in the region z > 0. In particular, for a quark-antiquark pair of weight ±w k , translating the σ field by 2πw k makes the field continuous. This allows us to directly compare the cross sections of the confining strings with the 1D solutions from section 7.2 above.
An example of such a comparison for N = 5 and k = 2 is shown in figure 19. Here, the relaxed sum of two 1D DWs is compared with a similar cross section to that shown in figure 6, although translated as described above to impose continuity. It is clear from this figure that the cross section closely replicates the 1D solution aside from a small difference in the string separations, for which the double string obeys the log-growth property described in section 6.3.
In general, we observed that the y = 0 cross sections for all of our confining strings were very similar to the corresponding one-dimensional configurations of two summed and relaxed BPS solitons. This held true for both the double-string configurations and the collapsed single-string configurations. Indeed, features of the collapsed 1D configurations -such as the nonlinear interaction of the φ fields like that depicted in the upper-right panel of figure 18 -are also replicated in the corresponding confining-string cross sections. As such, we believe that many properties of the confining strings, such as string tension and the tendency to collapse, can be predicted from the behavior of 1D DWs. This allows us to consider properties such as the behavior of higher gauge groups, or of weights in Z N orbits different from that of w k , without the significant computational demands of the two-dimensional configurations.

String tensions of differing Z N orbits
In this section, we study the tensions of double-string configurations with N -ality k weights belonging to different Z N orbits. Recall that our work up to this point has focused primarily on the string tensions of the orbits belonging to the highest weight, w k , of the k-index antisymmetric representation. As discussed in section 4, there exist double-string configurations of N -ality k belonging to differing Z N orbits, which can be confined by up-to BPS k-walls. However, if the string tensions of these configurations are higher than that of the highest-weight orbit, we would expect that such weights would be screened by the pair-production of W -bosons down to w k , i.e., such that the configuration is confined by BPS 1-walls. As discussed above, we can study the string tensions of these Z N orbits in a one-dimensional environment with the expectation that the corresponding 2D string tensions will be similar. In this section, therefore, we confirm that the orbits of N -ality k = 2 are expected to decay to the highest-weight orbit.
in [61], are integrable. We do not know if there is any connection between this fact and the "almost linear" property of the σ fields alluded to above.  Figure 19. One-dimensional DW versus the y = 0 cross section of a double string for N = 5 and N -ality k = 2. The top row shows the field of the 1D DW, which was produced by relaxing a configuration consisting of the sum of two 1-walls, as described in section 7.2. The bottom row shows a y = 0 cross section similar to that depicted in figure 6, although here the σ field has been translated by 2πw 2 in the region z > 0 such that it is continuous.
Recall from section 4.2.3 that the Z N orbits of N -ality k = 2 can be labeled by representative weights ν 1,1+p = ν 1 + ν 1+p , where p = 0, . . . , N − 1. We may estimate the string tension of these weights by numerically relaxing a 1D configuration with total flux Φ = 2πν 1,1+p . In particular, we set up our field with boundary values σ(−∞) = 0 and σ(∞) = 2πν 1,1+p , and with initial conditions linearly interpolating between these boundary values. When relaxing such a configuration, we found that the algorithm converged to the expected DW configuration of the given flux, and in particular replicated the y = 0 cross section of the corresponding 2D confining string whenever the 2D solution was available. Therefore, the tensions of such configurations can be used to approximate the string tensions of the 2D confining strings. We stress that these 1D tensions are only approximations of the 2D tensions, and are subject to some variation depending on the initial conditions used. In particular, these solutions may relax to some local minimum of energy instead of the global minimum. The string tensions reported herein represent the minimum energy configurations observed in our tests. The weight of each point is ν 1,1+p = ν 1 + ν 1+p . Also marked are twice the BPS 1-wall and 2-wall tensions for each N . The lowest tension orbit belongs to the highest weight, ν 1,2 = w 2 , and aligns closely with the tension of two BPS 1-walls. Notably, orbits with p > 1 and N > 5 appear to form a bound state with string tension lower than twice the BPS 2-wall tension, despite the expected double-string configuration for these weights being composed of 2-walls. This may be related to the N -ality k = 1 collapsing effect which led to the single-string configurations. Figure 20 plots the approximate string tensions of each N -ality k=2 orbit for N = 4, . . . , 9, computed from the one-dimensional setup described above. Also plotted are twice the BPS 1-wall and 2-wall tensions for comparison. As expected, the p = 1 orbit belonging to the highest weight consistently exhibits the lowest string tension, and closely follows twice the BPS 1-wall string tension for all N . As such, the assumption that strings will decay to the highest weight appears to be valid for N -ality k = 2. Note that the string tensions of the p = 1 orbit differ slightly from twice the BPS 1-wall tensions due to the effect of the binding energy of the two walls, although this effect is too small to be visible on the plot for N < 9. It is also notable that the Z N orbits with p > 1 do not strictly match the tension of two BPS 2-walls, despite the expectation that the double-string configurations for these orbits would consist of 2-walls. Instead, for N > 5 some of these orbits appear to form a lower-energy bound state. This may be related to the collapse of the confining string for N -ality k = 1 studied earlier.
We also confirmed that the p = 2 orbit is confined by 2-walls for N = 4, 5 by considering the two-dimensional double-string configurations with weight ν 1,3 . Figure 21 gives an example of this for N = 4. Here, the configuration is confined by a "magnetless" 2-wall, i.e., a 2-wall where the φ field vanishes identically. Indeed, it was shown in [28] that such magnetless k-walls can exist only for even N and for k = N/2. Thus, the vanishing of the φ field in figure 21 is consistent with confinement by a BPS 2-wall.  Figure 21. N = 4 double-string configuration with weight ν 1,3 . As a member of the p = 2 orbit, this configuration is confined by two BPS 2-walls, as expected from eq. (4.9). Note that these walls are "magnetless," i.e., the φ field vanishes everywhere (see [28]). It is expected that, for a large enough separation, pair production of W -bosons will cause this double-string configuration to decay to the energetically more favorable one made of two BPS 1-walls.

"Gluon confinement" and triple strings
As a final note, we also considered the behavior of a confining string with a weight belonging to the root lattice, corresponding to replacing the quark static sources by gluons (Wbosons). Such a configuration has N -ality 0, and accordingly does not belong to any of the Z N orbits described earlier. The potential energy density of such confining strings for N = 5 and with weight α 1 is shown in figure 22. Unlike all previous confining strings, this configuration possesses three distinct flux tubes: a single tube with a strong flux forming a straight line between the gluons, as well as two weaker flux tubes above and below this line. In fact, the uppermost and lowermost flux tubes are BPS 1-walls with fluxes Φ 1,I 1 = 2πρ 5 − 2πw 1 , while the thick wall is a BPS 3-wall with flux Φ 3,I 1,3,4 = 3 2πρ 5 − 2π(w 1 + w 3 + w 4 ), using the notation from (4.6). The two vacua inside the triple strings are the k = 1 and k = 4 vacua (2.14) and the non-Cartan gluon sources are embedded in the k = 0 vacuum. We took the signs of the fluxes such that the total flux is the sum of the fluxes of the three DWs: 2Φ 1,I 1 + Φ 3,I 1,3,4 = −2πα 1 . This effect is consistent with the one-dimensional configuration of equivalent flux when compared with the y = 0 cross section of the confining string. As explained earlier, at large enough separation between the root-lattice sources our abelian EFT breaks down and pair production of W -bosons is expected to break this triple flux tube.  (3) and N-ality k = 1 is tested separately in a 30 by 30 grid and a 60 by 60 grid. The fact that they agree even when the quarks are right at the edge of the smaller grid (R = 28) shows that the boundary effect in this confining theory is minimal, and that taking the finite distance boundary to approximate spatial infinity is justified.
The answer is that unlike standard electromagnetism, our theory is confining, that is, the field lines do not naturally spread to infinity and so the edge effect is weak. In fact, it is easy to show that if we turn off the superpotential W in our problem, the desired logarithmic growth of the electrostatic energy as a function of quark separation R quickly gets distorted by the edge effect.
To dispel all doubts that the edge effect is minimal here, we tested the energy dependence on R separately in a grid of size 30 by 30 and another of size 60 by 60. We push the quarks right up to the edge in the size 30 grid (with R = 28) and show that the energy still reproduces the "true" value, that of the larger grid; see figure 23.

B Methods for speed up
There are various standard methods to speed up the Gauss-Seidel convergence process. We will not describe them in any detail, except for mentioning that while using the successive overrelaxation method [53], we found that the overrelaxation parameter can be set to as high as 1.96 in most cases tested, while a value any larger might prevent convergence.
We now briefly describe two speed-up methods specific to the present problem (a third one is the BPS initial condition, as described in section 5.2. Half-grid method. We can exploit the fact that (5.1) and its boundary condition are symmetric under y → −y to reduce the computational time by a factor of 2. We first cut the standard grid in half along the y-direction. Along the three original sides of the JHEP01(2021)044 resultant half-rectangle, the boundary condition remains the same, that is, x is still equal to one of the vacua there. Those three sides continue to be untouched by the Gauss-Seidel iterations, i.e. they respect a Dirichlet boundary condition. The new boundary that comes from slicing the grid now respects a Neumann boundary condition: it must have a vanishing derivative in the y-direction. We implement this by setting the values along this boundary to be equal to its neighboring value in the y-direction, at the end of every iteration. At the end, we reflect the solution to restore the full grid. Unfortunately, this trick cannot be repeated in the z-direction because ∂ z δ(z) is odd.
Simplifying the EOM. We are often interested in the regime where the one-loop correction can be ignored, in which case the Kähler metric is just the Kronecker delta: K ab = δ ab . Then, (5.1) becomes Explicitly writing out the first term and using the identity α i · α j = 2δ ij − δ i,j+1 − δ i,j−1 , where the index is taken modulo N , we have This reduces a double nested for-loop to a single for-loop, which is a significant improvement in speed as this term must be evaluated at every point for every iteration. In addition, the 1D EOM (5.2) admits the same simplification.

C.1 String tensions T (N, k) for N ≤ 10
In eq. (C.1), we summarize the numerical results for string tensions T (N, k) and their ratios f (N, k) = T (N, k)/T (N, 1), as well as the uncertainties arising from the least squares fit.

JHEP01(2021)044
The string tension ratios f (N, k) for the largest number of colors we studied, N = 10, are plotted in figure 1 alongside other scaling laws; see also the discussion in section 6.2.
N kT f