Simplicity of confinement in SU(3) Yang-Mills theory

Simplicit y of confinemen t in SU(3) Y ang-Mills theory Xa vier Crean , 1 , ∗ Jeffrey Giansiracusa , 2 , † and Biagio Lucini 3 , ‡ 1 Dep artment of Mathematics, F aculty of Science and Engine ering, Swanse a University, F abian W ay, Swanse a, SA1 8EN, UK 2 Dep artment of Mathematic al Sciences, Durham University, Upp er Mountjoy Campus, Durham, DH1 3LE, UK 3 Scho ol of Mathematic al Scienc es, Que en Mary University of L ondon, Mile End R o ad, L ondon, E1 4NS, UK (Dated: F ebruary 11, 2026) W e in tro duce a no vel observ able asso ciated to Ab elian monopole currents defined in the Maximal Ab elian Pro jection of SU(3) Y ang-Mills theory that captures the top ology of the current lo op. This observ able, referred to as the simplicity , is defined as the ratio of the zeroth ov er the first Betti n umber of the curren t graph for a giv en field configuration. A n umerical study of the exp ectation v alue of the simplicity p erformed in the framew ork of Lattice Gauge Theories enables us to determine the deconfinement temp erature to a higher degree of accuracy than that reached by con ven tional metho ds at a comparable computational effort. Our results suggest that Ab elian current lo ops are strongly correlated with the degrees of freedoms of the theory that determine confinemen t. Our in vestigation op ens new p erspectives for the definition of an order parameter for deconfinement in Quan tum Chromo dynamics able to exp ose the p oten tially rich phase structure of the theory . Intr o duction — The problem of colour confinement in the strong in teractions, i.e., the absence of fractionally c harged particles in asymptotic states, is broadly ac- kno wledged as one of the most urgent gaps to fill in or- der to gain a full understanding of Quan tum Chromo- dynamics (QCD) and, more in general, of the dynam- ics of non-Abelian gauge theories. It is kno wn that at high temperature confinemen t is lost. The widely ac- cepted scenario (for a recent review, see, e.g., Ref. [1]) is that ab o ve a critical temp erature the system b eha v es as a plasma of its elementary constituents, quarks and glu- ons. This picture has recently b een challenged, and the suggestion that an intermediate regime separates the con- finemen t regime from the deconfined quark-gluon-plasma one (see, e.g., Refs. [2–7]) is gaining increasing consen- sus. In the light of this developmen t, a full characteri- sation of confinement app ears to be ev en more crucial. T o study QCD, a first-principles approach that has b een used, often in conjunction with analytic guidance, is Lat- tice Gauge Theory , of whic h the demonstration of con- finemen t through the area law of Wilson lo ops has b een the very first application [8]. It is an old idea that confinemen t can b e under- sto od in terms of the dynamics of top ological exci- tations carrying a non-trivial magnetic charge [9, 10]. This has led to the proposal of defining Ab elian mag- netic monop oles through a partial gauge fixing pro cedure kno wn as Ab elian pr oje ction [11]. A num b er of studies w ere p erformed whic h explored these ideas from v ari- ous angles [12–14], from dominance of degrees of free- dom that are Ab elian in certain pro jections [15] to the construction of disorder parameters [16–18] 1 and effec- 1 Pathologies iden tified in the original disorder parameter con- struction [19] led to the improv ed proposal [20]. tiv e monop ole potentials [21]. Studies of the gauge- indep endence of monop ole condensation [22] and in vesti- gations of thermal monop ole prop erties [23, 24] provided imp ortan t insights into the connection b et w een monop ole dynamics and the confinement/deconfinemen t transition, with an approach based on an effective mo del for the con- densation of monop oles at the critical temp erature [25] enabling the computation of the latter in SU(2) Y ang- Mills [26]. More recent inv estigations [27, 28] reaffirmed the centralit y of monop oles in colour confinement. Despite this noticeable progress, we still lack a robust (dis)order parameter for confinement in Y ang-Mills the- ory based on features of these top ological excitations that is free from lattice artefacts, enables determining with sufficien t accuracy quantities that c haracterise the dy- namics of the transition, and can b e readily adopted in the presence of dynamical fermions. A key difficulty is the top ological nature of monop ole excitations, whic h is harder to exp ose in the lattice discretisation. Recently , T opological Data Analysis (TDA) has emerged in compu- tational top ology as a robust metho dology to characterise top ological prop erties of discrete sets of p oin ts. F or this reason, applications of TDA to field theory and statistical mec hanics are raising at a fast pace (see, e.g., [29–38]). Building on the intuition developed from previous studies and b orro wing rigorous to ols of TDA along the path laid in Refs. [39, 40], in this letter we are going to prop ose and test in SU(3) Y ang-Mills an observ able derived from mag- netic monop ole currents that quan titatively characterises the deconfinement phase transition. W e will demonstrate n umerically that this observ able, whose definition in the presence of fermions remains identical, enables us to ac- curately determine the critical v alue of the coupling and the order of the phase transition in the pure gauge sys- tem, henceforth providing evidence of the relev ance of Ab elian monop oles for deconfinement in Y ang-Mills. 2 SU(3) Y ang-Mil ls on the L attic e — W e consider the SU(3) lattice gauge theory describ ed b y the Wilson ac- tion, S = β X i X µ<ν  1 − 1 3 ℜ eT r U µν ( i )  , (1) with β ≡ 6 /g 2 and g the coupling of the theory . This action b ecomes the Y ang-Mills action in the contin uum limit. The quantit y U µν ( i ) ≡ U µ ( i ) U ν ( i + ˆ µ ) U † µ ( i + ˆ ν ) U † ν ( i ) (2) is the plaquette v ariable, with U µ ( i ) ∈ SU(3) the link v ariable, defined on the link ( i ; ˆ µ ) stemming from the p oin t i and ending at the point i + ˆ µ (with ˆ µ the unit v ector in the p ositiv e direction µ ) of a four-dimensional Euclidean lattice of dimension N t × N 3 s . Periodic b ound- ary conditions are imp osed in all directions. The tem- p erature T of the system is given by T = 1 a ( β ) N t , (3) with a ( β ) the lattice spacing, which is a monotonically decreasing function of β . The path integral is given by Z = Z   Y i,µ dU µ ( i )   exp {− S } , (4) with dU µ ( i ) representing the Haar measure. A t thermal equilibrium, the exp ectation v alue of a giv en observ able O is computed using ⟨ O ⟩ = 1 Z Z   Y i,µ dU µ ( i )   O exp {− S } . (5) By fixing N t and taking the infinite v olume limit N s → ∞ , we can study the system’s b eha viour in the thermo- dynamic limit. The thermo dynamic observ ables thus ex- tracted can then b e extrap olated to the contin uum limit b y taking progressively finer discretisations in N t . Magnetic monop oles in the Maximal Ab elian Gauge — The identification of Abelian magnetic monopole cur- ren ts in SU( N ) gauge theories is p erformed following a gauge-fixing procedure known as an A b elian pr oje ction . T o define an Ab elian pro jection, an operator F trans- forming in the adjoint representation is c hosen, and a partial gauge fixing is p erformed. The latter consists in diagonalising F at eac h point and ordering its eigen v al- ues λ 1 , . . . λ N in non-decreasing order. This partial gauge fixing leav es a residual U(1) N − 1 symmetry . The pro cess exp oses N − 1 sp ecies of magnetic monop oles, each corre- sp onding to one of the U(1) residual gauge factors. Mag- netic monop oles of the type j occur where λ j = λ j +1 . In SU(3) Y ang-Mills, we then hav e tw o types of monop oles. While the gauge-dep endence of the iden tification of monop ole currents might obscure the gauge-inv ariant pic- ture, physical prop erties of magnetic monop oles in some gauge can b e a useful starting p oint to b etter understand prop erties of more complex underlying top ological struc- tures that the gauge-dep enden t ob jects exemplify . In this resp ect, a particularly useful gauge is the Maximal Ab elian Gauge ( MA G ), since gauge-indep enden t excita- tions maximally ov erlap with Ab elian monop oles defined in this gauge [41]. F ollo wing the prescription of Ref. [42], for the SU(3) lattice Y ang-Mills theory , we define the adjoin t op erator ˜ X ( i ) as ˜ X ( i ) = X µ h U µ ( i ) ˜ λU † µ ( i ) + U † µ ( i − ˆ µ ) ˜ λU µ ( i − ˆ µ ) i , ˜ λ = diag(1 , 0 , − 1) . (6) The MAG is defined as the gauge in which ˜ X ( i ) is diagonal. This gauge choice is equiv alent to requiring that the op erator ˜ F MAG ( U, g ) = X µ,i tr  g ( i ) U µ ( i ) g † ( i + ˆ µ ) ˜ λg ( i + ˆ µ ) U † µ ( i ) g † ( i ) ˜ λ  (7) is maximised ov er g , i.e., the gauge fixing transformation { g } can b e deriv ed from the condition { ˜ g } = argmax { g } ˜ F MAG ( U, g ) . (8) In the MA G, the diagonal elements of the link matrices ˜ U ii read ˜ U ii = r i e iφ i , X i φ i = 2 π n + δ φ . (9) In general, δ φ  = 0 , as a consequence of the fact that even after gauge fixing the links are not fully diagonal. Angle v ariables ϕ i are then defined through the redistribution 3 of the excess phase as ϕ i = φ i − δ φ    ˜ U ii    − 1 P j    ˜ U j j    − 1 . (10) The tw o lattice Ab elian fields in the residual gauge are θ 1 = ϕ 1 and θ 2 = − ϕ 3 . The corresponding species of monop oles are defined follo wing the prescription of Ref. [43]. The truncation of the theory to the degrees of freedom represen ted b y these tw o angles pro vides a definition of Ab elian pr oje ction using the MAG. T op olo gic al simplicity — T o inv estigate the top ology of monop ole curren ts, it is conv enien t to define the dual lattice Λ ∗ , obtained by shifting each p oin t by half a lattice spacing in all p ositiv e directions. Eac h geometric element of the original lattice Λ of dimension d is pierced at its cen tre by an element of Λ ∗ of dimension 4 − d . These t wo elements are dual to eac h other. Being defined on cub es of Λ , on Λ ∗ magnetic charges are link v ariables and are more appropriately referred to as currents. F or eac h lattice configuration, each dual link ( j ; ˆ µ ) carries a curren t m µ ( j ) = 0 , ± 1 , ± 2 . A non-trivial current m on the dual link ( j ; ˆ µ ) can b e seen as a charge m moving from site j to site j + ˆ µ . Akin to Kirchoff ’s law, the total c harge entering a dual site j is equal to the total charge exiting from that site. F or the same reason, curren ts cannot end at sites but m ust form closed lo ops. The set of links G ≡ { ( j ; ˆ µ ) , m µ ( j )  = 0 } is the curren t graph asso ciated to the configuration on which the m µ ( j ) hav e b een computed. This ob ject is interpreted as a graph, with v ertices coinciding with the end p oin ts of links in G and edges coinciding with the links in G . The n umber of connected comp onen ts of G defines the Betti n umber b 0 , while the total n umber of lo ops defines the Betti num ber b 1 . In [40], we hav e shown that b oth quantities ρ 0 = ⟨ b 0 ⟩ N 3 s and ρ 1 = ⟨ b 1 ⟩ N 3 s (11) are sensitive to the phase transition, by pro viding nu- merical evidence that their susceptibilities, χ 0 and χ 1 , displa y a p eak at a v alue β c ( N s , N t ) of β that for fixed N t scales as β c ( N s , N t ) = β c ( N t ) + a/ N 3 s , (12) where β c ( N t ) is the critical v alue of β at fixed N t . While these results show a remark able connection b et ween topo- logical prop erties of current graphs and the deconfine- men t phase transition, suggesting that exp ectation v alues of the Betti num b ers b 0 and b 1 are key to understanding confinemen t, they fail short of providing an order param- eter for the deconfinement phase transition. In this w ork, w e use the information ab o v e and exp ec- tations developed from previous calculations to define a N t N s β min β max N β N meas 4 16 , 20 , 24 , 28 , 32 5 . 6600 5 . 7200 15 600 6 24 , 30 , 36 , 42 , 48 5 . 8100 5 . 9300 12 600 8 32 , 40 , 48 , 56 , 64 6 . 0350 6 . 0800 12 400 T ABLE I. This table sp ecifies the v arious lattice sizes, N t and N s , that we use in this study . F or N t = 4 , we draw a sample of N = 600 configurations, for a range of 15 β v alues. F or N t = 6 , w e draw a sample of N = 600 configurations, for a range of 12 β v alues. F or N t = 8 , w e draw a sample of N = 400 configurations, for a range of 12 β v alues. quan tity that b eha v es as an order parameter. This ob- serv able is the top olo gic al simplicity , or just simplicity , whic h is defined as λ =  b 0 b 1  . (13) The ratio that defines this observ able, b 0 b 1 , is the recip- ro cal of the num b er of lo ops p er connected comp onen t, whic h is som etimes called the c omplexity and is at least 1. While the simplicity and complexity are undefined for an empty netw ork, a typical configuration at non-infinite β will hav e a non-empty monop ole netw ork. In the con- fined phase at low temp erature (corresp onding at fixed N t to low β ), one typically sees a p ercolating current net- w ork with one component and many lo ops, and so the simplicit y approaches 0. In the deconfined phase at high temp erature (corresp onding at fixed N t to high β ), large fluctuations are suppressed b y the Boltzmann weigh t and dominan t configurations ha ve current netw orks consist- ing of a sparse gas of lo ops, so the simplicity is approx- imately 1. In the critical region, due to the breaking of the p ercolating lo op, the simplicity raises from 0 to 1. Ho wev er, this qualitative b eha viour is not sufficient to evidence the coupling of λ with the degrees of freedom that are relev ant for the phase transition. In order to test whether λ is an observ able that is sensitive to the deconfinemen t phase transition, we study the behaviour of its susceptibility χ λ = V *  b 0 b 1  2 + −  b 0 b 1  2 ! (14) in a range of β v alues in the critical region, for a range of spatial extensions N s and temp oral extensions N t = 4 , 6 , 8 , as specified in T ab. I. Larger v alues of N t cor- resp ond to finer discretisations of the direction associ- ated with the temperature (see Eq. (3)). V arying N s will instead enable us to prov e sensitivit y to the ex- p ected first-order deconfinemen t phase transition, which will be shown b y the detection of the scaling b eha viour in Eq. (12). W e rep ort in Fig. 1 the b eha viour of λ and of χ λ together with a determination for β c for N t = 8 with 4 6.050 6.055 6.060 6.065 6.070 6.075 6.080 0.56 0.57 0.58 0.59 0.60 0.61 0.62 0.63 0.64 L i t e r . c ( , N t = 8 ) e s t i m a t e L i t e r . c ( , N t = 8 ) s t d . e r r o r N t , N s = 8 , 3 2 N t , N s = 8 , 4 0 N t , N s = 8 , 4 8 N t , N s = 8 , 5 6 N t , N s = 8 , 6 4 6.050 6.055 6.060 6.065 6.070 6.075 6.080 10 20 30 40 50 60 70 80 90 100 L i t e r . c ( , N t = 8 ) e s t i m a t e L i t e r . c ( , N t = 8 ) s t d . e r r o r N t , N s = 8 , 3 2 N t , N s = 8 , 4 0 N t , N s = 8 , 4 8 N t , N s = 8 , 5 6 N t , N s = 8 , 6 4 FIG. 1. F or N t = 8 , λ and χ λ as functions of β zo omed in to the critical region with translucent lines to guide the eye. The vertical line and band sho w resp ectiv ely the central v alue and the statistical error for the extrap olated β c determined in Ref. [44]. con ven tional metho ds taken from Ref. [44]. λ shows a monotonically increasing b eha viour, with slop e that gets steep er near β c and becomes more pronounced as the v olume increases. As a consequence, the susceptibility χ λ has a p eak in the critical region with height grow- ing with the v olume. Figs. 2, 3 and 4 sho w the scaling of the p osition of the p eak of χ λ with the spatial vol- ume, resp ectiv ely at N t = 4 , 6 , 8 . The data hav e b een fitted with Eq. (12) and with a higher-order correction in 1 / N 3 s for v arious fitting ranges, and the results hav e b een combined with the Ak aike information principle in the implementation of Ref. [45] to determine a weigh ted extrap olation and the corresp onding errors (b oth also re- p orted in the figures). 0.0 0.5 1.0 1.5 2.0 2.5 N 3 s 1e 4 5.6920 5.6925 5.6930 5.6935 5.6940 5.6945 5.6950 L i t e r . c ( , N t = 4 ) e s t i m a t e L i t e r . c ( , N t = 4 ) s t d . e r r o r c ( N s , N t = 4 ) c ( , N t = 4 ) e s t i m a t e F its P olyn. degr ee: 1 N s v a l u e s : 2 4 , 2 8 , 3 2 P olyn. degr ee: 1 N s v a l u e s : 2 0 , 2 4 , 2 8 , 3 2 P olyn. degr ee: 1 N s v a l u e s : 1 6 , 2 0 , 2 4 , 2 8 , 3 2 FIG. 2. P osition of the p eak of χ λ for lattices with N t = 4 and N s = 16 , 20 , 24 , 28 , 32 (blue triangles). Dashed lines represen t the p olynomial regression fits used for the infinite v olume extrap olation of these p eak v alues, and the red circle is the resulting estimate of β c in the thermodynamic limit (see the supplemental material for additional details). The horizon tal line and band show resp ectiv ely the central v alue and the statistical error for the extrap olated β c determined in Ref. [44]. 0 1 2 3 4 5 6 7 8 N 3 s 1e 5 5.893 5.894 5.895 5.896 5.897 5.898 5.899 L i t e r . c ( , N t = 6 ) e s t i m a t e L i t e r . c ( , N t = 6 ) s t d . e r r o r c ( N s , N t = 6 ) c ( , N t = 6 ) e s t i m a t e F its P olyn. degr ee: 1 N s v a l u e s : 3 6 , 4 2 , 4 8 P olyn. degr ee: 1 N s v a l u e s : 3 0 , 3 6 , 4 2 , 4 8 P olyn. degr ee: 1 N s v a l u e s : 2 4 , 3 0 , 3 6 , 4 2 , 4 8 FIG. 3. P osition of the p eak of χ λ for lattices with N t = 6 and N s = 24 , 30 , 36 , 42 , 48 (blue triangles). Dashed lines represen t the p olynomial regression fits used for the infinite v olume extrap olation of these p eak v alues, and the red circle is the resulting estimate of β c in the thermodynamic limit (see the supplemental material for additional details). The horizon tal line and band show resp ectiv ely the central v alue and the statistical error for the extrap olated β c determined in Ref. [44]. Discussion and outlo ok — The smaller error bar of the extrap olated v alue with resp ect to the reference v alue, whic h has b een obtained by av eraging ov er around an or- der of magnitude more configurations, provides evidence of the relev ance of the top ological prop erties captured b y the simplicity for the dynamics of colour confinement. 5 0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 N 3 s 1e 5 6.060 6.062 6.064 6.066 6.068 6.070 L i t e r . c ( , N t = 8 ) e s t i m a t e L i t e r . c ( , N t = 8 ) s t d . e r r o r c ( N s , N t = 8 ) c ( , N t = 8 ) e s t i m a t e F its P olyn. degr ee: 1 N s v a l u e s : 4 8 , 5 6 , 6 4 P olyn. degr ee: 1 N s v a l u e s : 4 0 , 4 8 , 5 6 , 6 4 P olyn. degr ee: 1 N s v a l u e s : 3 2 , 4 0 , 4 8 , 5 6 , 6 4 FIG. 4. P osition of the p eak of χ λ for lattices with N t = 8 and N s = 32 , 40 , 48 , 56 , 64 (blue triangles). Dashed lines represen t the p olynomial regression fits used for the infinite v olume extrap olation of these p eak v alues, and the red circle is the resulting estimate of β c in the thermodynamic limit (see the supplemental material for additional details). The horizon tal line and band show resp ectiv ely the central v alue and the statistical error for the extrap olated β c determined in Ref. [44]. W e remark that the simplicity can b e adopted also in the presence of dynamical fermions. If a nov el intermediate- temp erature b eha viour exists in QCD, based on our stud- ies of spin systems, where observ ables constructed with similar metho dology as the simplicity hav e prov ed to b e sensitiv e to the order-disorder phase transition, we ex- p ect the susceptibility of the simplicit y to sho w a non- monotonic structure at the higher-temp erature change of regime. A cknow le dgements — W e thank M. D’Elia, A. Gonzalez-Arro yo and T. Sulejmanpasic for discussions. X C was supported by the A dditional F unding Pro- gramme for Mathematical Sciences, delivered by EPSRC (EP/V521917/1) and the Heilbronn Institute for Math- ematical Researc h. JG w as supp orted by EPSRC grant EP/R018472/1 through the Centre for TDA and the Er- langen Hub for AI through EPSRC grant EP/Y028872/1. The work of BL was partly supp orted by the EPSR C ExCALIBUR ExaTEPP pro ject EP/X017168/1 and b y the STFC Consolidated Gran ts No. ST/T000813/1 and ST/X000648/1. Numerical sim ulations ha ve been p erformed on the Sw ansea SUNBIRD cluster (part of the Sup ercomputing W ales pro ject) and AccelerateAI A100 GPU system. The Sw ansea SUNBIRD system and A ccelerateAI are part funded b y the European Regional Dev elopmen t F und (ERDF) via W elsh Gov ernment. This work used the DiRAC Data Intensiv e service (CSD3) at the Universit y of Cam bridge, managed by the Universit y of Cambridge Universit y Information Ser- vices on b ehalf of the STF C DiRAC HPC F acility (www.dirac.ac.uk). The DiRAC comp onen t of CSD3 at Cam bridge was funded by BEIS, UKRI and STF C capi- tal funding and STFC op erations gran ts. DiRA C is part of the UKRI Digital Research Infrastructure. Op en A ccess Statemen t — F or the purp ose of op en access, the authors hav e applied a Creativ e Commons A ttribution (CC BY) licence to any Author A ccepted Man uscript version arising. Researc h Data Access Statemen t — The data and analysis co de for this manuscript can b e do wnloaded from Ref. [46]. The Monte Carlo co de can b e found from Ref. [47]. ∗ 2237451@sw ansea.ac.uk † jeffrey .giansiracusa@durham.ac.uk ‡ b.lucini@qm ul.ac.uk [1] G. Aarts et al. , Prog. Part. Nucl. Phys. 133 , 104070 (2023), arXiv:2301.04382 [hep-lat]. [2] L. Y. Glozman, Acta Phys. Polon. Supp. 10 , 583 (2017), arXiv:1610.00275 [hep-lat]. [3] A. Alexandru and I. Horváth, Phys. Rev. D 100 , 094507 (2019), arXiv:1906.08047 [hep-lat]. [4] M. Cardinali, M. D’Elia, and A. P asqui, arXiv:2107.02745 [hep-lat] (2021), [5] M. Hanada, H. Ohata, H. Shimada, and H. W atanabe, PTEP 2024 , 041B02 (2024), arXiv:2310.01940 [hep-th]. [6] Y. F ujimoto, K. F ukushima, Y. Hidak a, and L. McLerran, Ph ys. Rev. D 112 , 074006 (2025), arXiv:2506.00237 [hep- ph]. [7] J. A. Mickley , C. Allton, R. Bignell, and D. B. Lein weber, Ph ys. Rev. D 112 , 054505 (2025), arXiv:2507.17200 [hep- lat]. [8] K. G. Wilson, Phys. Rev. D 10 , 2445 (1974). [9] S. Mandelstam, Phys. Rept. 23 , 245 (1976). [10] G. ’t Ho oft, in New Phenomena in Subnucle ar Physics , edited by A. Zichic hi (Editrice Comp ositori, Bologna, 1976) pp. 1225–1249, pro ceedings of the talk given at the 13th In ternational School of Subnuclear Ph ysics: New Phenomena in Subnuclear Physics, 11 July-1 August 1975. Erice, Italy. [11] G. ’t Ho oft, Nucl. Phys. B 190 , 455 (1981). [12] A. S. Kronfeld, G. Sc hierholz, and U. J. Wiese, Nucl. Ph ys. B 293 , 461 (1987). [13] A. S. Kronfeld, M. L. Laursen, G. Schierholz, and U. J. Wiese, Phys. Lett. B 198 , 516 (1987). [14] T. L. Iv anenko, A. V. Pochinsky , and M. I. P olik arpov, Ph ys. Lett. B 302 , 458 (1993). [15] Z. F. Ezaw a and A. Iwazaki, Phys. Rev. D 25 , 2681 (1982). [16] L. Del Debbio, A. Di Giacomo, G. Paffuti, and P . Pieri, Ph ys. Lett. B 355 , 255 (1995), arXiv:hep-lat/9505014. [17] A. Di Giacomo, B. Lucini, L. Montesi, and G. Paffuti, Ph ys. Rev. D 61 , 034503 (2000), arXiv:hep-lat/9906024. [18] A. Di Giacomo, B. Lucini, L. Montesi, and G. Paffuti, Ph ys. Rev. D 61 , 034504 (2000), arXiv:hep-lat/9906025. [19] J. Greensite and B. Lucini, Phys. Rev. D 78 , 085004 (2008), arXiv:0806.2117 [hep-lat]. [20] C. Bonati, G. Cossu, M. D’Elia, and A. Di Giacomo, 6 Ph ys. Rev. D 85 , 065001 (2012), arXiv:1111.1541 [hep- lat]. [21] M. N. Cherno dub, M. I. Polik arp o v, and A. I. V eselov, Ph ys. Lett. B 399 , 267 (1997), arXiv:hep-lat/9610007. [22] J. M. Carmona, M. D’Elia, A. Di Giacomo, B. Lucini, and G. Paffuti, Phys. Rev. D 64 , 114507 (2001), arXiv:hep- lat/0103005. [23] M. N. Cherno dub and V. I. Zakharov, Phys. Rev. Lett. 98 , 082002 (2007), arXiv:hep-ph/0611228. [24] A. D’Alessandro and M. D’Elia, Nucl. Phys. B 799 , 241 (2008), arXiv:0711.1266 [hep-lat]. [25] M. Cristoforetti and E. Shury ak, Phys. Rev. D 80 , 054013 (2009), arXiv:0906.2019 [hep-ph]. [26] A. D’Alessandro, M. D’Elia, and E. V. Shury ak, Phys. Rev. D 81 , 094501 (2010), arXiv:1002.4161 [hep-lat]. [27] M. Nguyen, T. Sulejmanpasic, and M. Ünsal, Phys. Rev. Lett. 134 , 141902 (2025), arXiv:2401.04800 [hep-th]. [28] J. Giansiracusa, D. Lanners, and T. Sulejmanpasic, Ph ys. Rev. Lett. 135 , 221901 (2025). [29] I. Donato, M. Gori, M. Pettini, G. Petri, S. De Nigris, R. F ranzosi, and F. V accarino, Phys. Rev. E 93 , 052138 (2016), [30] Q. H. T ran, M. Chen, and Y. Hasega wa, Phys. Rev. E 103 , 052127 (2021), arXiv:2004.03169 [cond-mat.stat- mec h]. [31] K. Kashiwa, T. Hirakida, and H. Kouno, Symmetry 14 , 1783 (2022), arXiv:2103.12554 [hep-lat]. [32] B. Olstho orn, J. Hellsvik, and A. V. Balatsky , Ph ys. Rev. Res. 2 , 043308 (2020), arXiv:2009.05141 [cond-mat.stat- mec h]. [33] A. Cole, G. J. Loges, and G. Shiu, Phys. Rev. B 104 , 104426 (2021), arXiv:2009.14231 [cond-mat.stat-mech]. [34] N. Sale, J. Giansiracusa, and B. Lucini, Phys. Rev. E 105 , 024121 (2022), arXiv:2109.10960 [cond-mat.stat- mec h]. [35] D. Sehay ek and R. G. Melko, Phys. Rev. B 106 , 085111 (2022), arXiv:2201.09856 [cond-mat.stat-mech]. [36] N. Sale, B. Lucini, and J. Giansiracusa, Phys. Rev. D 107 , 034501 (2023), arXiv:2207.13392 [hep-lat]. [37] D. Spitz, J. M. Urban, and J. M. Pa wlowski, Phys. Rev. D 107 , 034506 (2023), arXiv:2208.03955 [hep-lat]. [38] D. Spitz, J. M. Urban, and J. M. Pa wlowski, Phys. Rev. D 111 , 114519 (2025), arXiv:2412.09112 [hep-lat]. [39] X. Crean, J. Giansiracusa, and B. Lucini, SciPost Phys. 17 , 100 (2024), arXiv:2403.07739 [hep-lat]. [40] X. Crean, J. Giansiracusa, and B. Lucini, PoS LA T- TICE2024 , 395 (2025), arXiv:2501.19320 [hep-lat]. [41] C. Bonati, A. Di Giacomo, L. Lep ori, and F. Pucci, Phys. Rev. D 81 , 085022 (2010), arXiv:1002.3874 [hep-lat]. [42] C. Bonati and M. D’Elia, Nucl. Phys. B 877 , 233 (2013), arXiv:1308.0302 [hep-lat]. [43] T. A. DeGrand and D. T oussaint, Phys. Rev. D 22 , 2478 (1980). [44] B. Lucini, M. T ep er, and U. W enger, JHEP 2004 (01), 061, arXiv:hep-lat/0307017. [45] W. I. Ja y and E. T. Neil, Phys. Rev. D 103 , 114502 (2021), arXiv:2008.01069 [stat.ME]. [46] X. Crean, J. Giansiracusa, and B. Lucini, Simplicit y of confinemen t in S U (3) Y ang-Mills theory — Data and Analysis Co de Release (2026). [47] X. Crean, J. Giansiracusa, and B. Lucini, Simplicit y of confinemen t in S U (3) Y ang-Mills theory — Mon te Carlo Co de Release (2026). [48] N. Cabibb o and E. Marinari, Ph ys. Lett. B 119 , 387 (1982). [49] P . Dlotk o, in GUDHI User and R efer enc e Manual (GUDHI Editorial Board, 2015). [50] A. M. F erren b erg and R. H. Swendsen, Phys. Rev. Lett. 63 , 1195 (1989). [51] M. Lusc her, Comput. Ph ys. Commun. 79 , 100 (1994), arXiv:hep-lat/9309020. 1 Supplemen tal Material: Simplicity of confinemen t in SU(3) Y ang-Mills theory Xa vier Crean, 1 Jeffrey Giansiracusa, 2 and Biagio Lucini 3 1 Dep artment of Mathematics, F aculty of Scienc e and Engine ering, Swanse a University, F abian W ay, Swanse a, SA1 8EN, UK 2 Dep artment of Mathematic al Scienc es, Durham University, Upp er Mountjoy Campus, Durham, DH1 3LE, UK 3 Scho ol of Mathematic al Scienc es, Que en Mary University of L ondon, Mile End R o ad, L ondon, E1 4NS, UK T o p ological Data Analysis for Lattice Gauge Theories T opology is ab out quantifying the asp ects of shap es that are inv ariant under stretching and deforming. T opological data analysis (TDA) emerged from this domain as a set of computational to ols that pro duce numerical inv arian ts of the shap e of ob jects arising from data. These numerical inv ariants can provide insight in to the shap e of the distribution underlying a set of samples. W e take a different approach here, using TDA to ols to pro duce non-lo cal observ ables for lattice gauge theories that are sensitive to large scale structures app earing in configurations. In algebraic top ology , one frequently represen ts geometric ob jects with simplicial complexes. How ever, since we are w orking with lattice gauge theory on cubical lattices, it is more conv enien t to represent geometric ob jects as cubical complexes, which are simply ob jects built by gluing cub es together along their b oundaries. One of the fundamental algebraic inv ariants of top ology is homolo gy , which is ab out quan tifying the holes and voids of a shap e. Holes come in different types. A point missing from the plane is a 1-dimensional hole in the sense that it can be captured by a circle, which is a 1-dimensional manifold. Ho wev er, a circle cannot capture a p oin t missing from 3-space b ecause it could alw a ys slip off ov er the top or b ottom; this hole m ust b e wrapp ed by a sphere, whic h is a 2-dimensional manifold, so this is a 2-dimensional hole. Homology formalises this idea. The input is a simplicial or cubical complex X . The output is a vector space H n ( X ) for each natural num ber n . The dimension of H n ( X ) represen ts the count of n -dimensional holes. Here are some useful prop erties of homology . 1. If X is a p oin t (or con tractible), then H 0 ( X ) is 1-dimensional, and H n ( X ) = 0 (the trivial vector space) for all n > 0 . 2. In general, H 0 ( X ) has dimension equal to the num b er of connected comp onen ts of X. 3. Disjoint u nion X ∪ Y corresp onds to direct sum: H n ( X ∪ Y ) ∼ = H n ( X ) ⊕ H n ( Y ) . 4. H n ( X ) = 0 for all n larger than the dimension of X . 5. If X is a d -dimensional closed and orientable manifold, then H d ( X ) is 1-dimensional. 6. If X and Y are homotopy equiv alent, then H i ( X ) ∼ = H i ( Y ) for all i . In some situations, having a vector space is more information than just having its dimension. Ho wev er, in this pap er w e will only make use of the dimensions. The dimension of H i ( X ) is known as the i th Betti num ber of X , and it is denoted b i . Cubical complexes An n -dimensional cub e is simply a Cartesian pro duct of n copies of the unit interv al [0 , 1] . More systematically , giv en a finite set A , the asso ciated cub e C ( A ) is [0 , 1] A . Let ( x a ) a ∈ A b e co ordinates on the cub e. A fac e of this cub e is a subspace where some collection of co ordinates x a are either 0 or 1. Each face is canonically linearly homeomorphic to a cub e of lo wer dimension. An orientation of a cube C ( A ) corresponds to a c hoice of ordering of the set A up to even p erm utations. Since a cub e is a manifold with boundary and corners, an orien tation of C ( A ) induces an orientation of eac h co dimension 1 face. In terms of orderings, the induced orientation of the face at x a = 0 is given by representing the orientation as 2 an ordering of A with a last and then restricting this to an ordering of A ∖ { a } , and the induced orien tation of the opp osite face at x a = 1 is given by the opp osite of this. A cubic al c omplex is a top ological space X together with a collection of maps { f i : C i → X } of cubes C i in to X , satisfying the following conditions: 1. Each f i is a homeomorphism onto its image. 2. The union of the images is all of X . 3. If the images of f i and f j ha ve non-empty in tersection K , then the comp osition of f i (restricted to the preimage of K ) follow ed by f − 1 j is a linear homeomorphism of a face of C i on to a face of C j . A sub c omplex is a subspace of X that is the union of the images of a subset of the cub es. The example we are concerned with here is the cubical complex corresponding to a dual lattice Λ ∗ for a lattice discretisation of spacetime Λ and sub complexes corresp onding to a collection of links. Chain complexes and their homology A chain c omplex is a sequence of vector spaces and linear maps . . . ∂ d +1 − → V d ∂ d − → V d − 1 ∂ d − 1 − → V d − 2 − → . . . suc h that the comp ositions ∂ d ◦ ∂ d +1 are all equal to 0. Giv en a cubical complex ( X, { f i : C i → X } ) , whic h we abbreviate simply as X , one obtains a chain complex as follo ws. The degree d space V d is spanned by the pairs ( C i an d -cub e of X , o an orientation on C d ) mo dulo the relation that − ( C , o ) = ( C , − o ) . The b oundary map ∂ d sends ( C , o ) to the sum of the co dimension 1 faces of C with their induced orientations. Note that if all vector spaces are defined o ver the field Z / 2 , then one do esn’t need to keep track of orientations and signs, as ( C , o ) = ( C , − o ) . Homology Elemen ts in the kernel of ∂ d are called cycles in degree d , and the subspace of cycles is written Z d . Elements in the image of ∂ i + d are called b oundaries of degree d , and the space of b oundaries is B d . Every b oundary is a cycle, so B d ⊂ Z d , but the con verse is not necessarily true. The homology H d is defined to b e the quotient vector space Z d /B d . It measures the extent to whic h there are cycles that are not b oundaries. The abov e algebraic definition ma y not appear all that intuitiv e, but it has the adv antage that it can easily be turned in to an algorithm to pro duce a basis for the homology in any giv en degree and hence determine the Betti n umbers. Graphs A graph (without edges that start and end at the same vertex) is an example of a cubical complex. In fact, these are the only kind of cubical complexes that we need to consider in this pap er. If X is a graph, then the homology is nontrivial only in degrees 0 and 1. The degree 0 part tells us the num ber of connected comp onents, which is the zer oth Betti numb er , b 0 . If G is a connected graph and T ⊂ G is a spanning tree (a subgraph that is a tree with the prop ert y that adding any additional edge of G to T results in it no longer b eing a tree), then con tracting all the edges of T results in a single vertex v and a collection of lo op edges that start and end at v , one for each edge in the complement of T . See Fig. S1. Since the spanning tree T is a tree, and hence con tractible, the original graph G is homotopy equiv alen t to the result of contracting T . One can show that the first Betti numb er , b 1 , (the dimension of H 1 ( G ) ) is equal to the num b er of edges in the complemen t of a spanning tree, which is the num b er of lo ops in G/T . 3 a b c d a b c d FIG. S1. Left: A graph G and a spanning tree T shown in red. Right: the graph G/T resulting from contracting T to a single v ertex. Mon te Carlo simulations In this study , we draw a sample of SU(3) lattice gauge field configurations using MCMC imp ortance sampling metho ds, namely the heat-bath and ov errelax algorithms [48]. W e define one comp osite up date as 1 heat-bath up date follo wed b y 4 o v errelax updates. Note that imp ortance sampling reliably approximates the system’s Boltzmann distribution once the Mark ov chain has b een successfully equilibrated; th us, for every sim ulation, w e discard the first 10 , 000 comp osite up dates. F urther, in order to minimise the auto correlation betw een successive recorded configurations, w e separate each recording by 2 , 000 comp osite up dates. F or each resp ectiv e lattice size, we select an appropriate range of β v alues to study that are strategically lo cated to co v er the critical region of the deconfinement transition; this has b een previously sp ecified in the literature, e.g., in Ref. [44]. F or eac h β v alue in our selected range, w e record a sample of N meas configurations. Using our giv en computational allocation, w e ha ve b een able to record N meas = 600 configurations for lattices with N t = 4 , 6 and N meas = 400 configurations for lattices with N t = 8 . Thermal exp ectation v alues of observ ables, as p er Eq. (5), are then estimated as the sample mean as a function of β . W e estimate the standard error in the sample mean by using the bo otstrap method with N bs = 2 , 000 . F urther details on our computational pip eline can b e seen in this publication’s accompanying co de release in Ref. [46]. Data analysis In this w ork, we analyse the union of the Ab elian magnetic monopole currents – treating b oth species on an equal fo oting. W e hav e verified that considering each individual species pro duces compatible results to the analysis of the union of sp ecies. F urther, we observ e no statistically significant difference b et w een the t wo sp ecies – this is consistent with the picture seen in the literature, e.g., Refs. [4, 40, 42]. Our aim is to analyse readily computable top ological inv ariants of the monop ole current net works called Betti n umbers. By fixing an arbitrary orientation (the results are indep enden t of this choice), we map the netw ork of monop ole currents to a directed graph G . W e then compute b 0 ( G ) ≡ dim H 0 ( G ) representing the comp onen ts of the graph and b 1 ( G ) ≡ dim H 1 ( G ) representing the num b er of loops in the graph. In our computation of b k ( G ) , we lev erage a highly optimised implementation of a homology computing algorithm, namely Ref. [49], designed to deal with a cubical complex. This is p ossible since a directed graph is a 1 -dimensional cubical complex; therefore, a cubical complex is a conv enien t data structure for represen ting graphs. See Ref. [39], for further details. Our results for ρ i ( β ) and χ i ( β ) (as defined in the main b ody) at N t = 4 , 6 , 8 are plotted in Figs. S2, S3, S4, S5, S6 and S7 resp ectiv ely . A suitable range of β -v alues cov ering the critical region ha ve b een chosen with error bars computed using the b ootstrap metho d with N bs = 2 , 000 . One can see that the susceptibilities χ i p eak at the critical v alue with p eaks b ecoming larger and more concentrated as the spatial volume increases as exp ected for an extensiv e observ able at the critical v alue of a first order phase transition. Similarly , our results for the simplicity λ ( β ) and χ λ ( β ) for N t = 4 , 6 , 8 are plotted in Figs. S8, S9 and 1 resp ectiv ely . One can see that the susceptibility χ λ p eaks at the critical v alue with p eaks b ecoming larger and more concentrated as the spatial volume increases. In Fig. S10, we plot a wider range of β -v alues for N t , N s = 8 , 32 so that the global trend is made clear for finite N s size. 4 5.685 5.690 5.695 5.700 0.020 0.022 0.024 0.026 0.028 0.030 0.032 0.034 0 L i t e r . c ( , N t = 4 ) e s t i m a t e L i t e r . c ( , N t = 4 ) s t d . e r r o r N t , N s = 4 , 1 6 N t , N s = 4 , 2 0 N t , N s = 4 , 2 4 N t , N s = 4 , 2 8 N t , N s = 4 , 3 2 5.685 5.690 5.695 5.700 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0 L i t e r . c ( , N t = 4 ) e s t i m a t e L i t e r . c ( , N t = 4 ) s t d . e r r o r N t , N s = 4 , 1 6 N t , N s = 4 , 2 0 N t , N s = 4 , 2 4 N t , N s = 4 , 2 8 N t , N s = 4 , 3 2 FIG. S2. F or N t = 4 , a scatter plot of ρ 0 and χ 0 as functions of β with translucent lines to guide the eye. The vertical line and band show respectively the central v alue and the statistical error for the extrap olated β c determined in Ref. [44]. 5.685 5.690 5.695 5.700 0.20 0.22 0.24 0.26 0.28 1 L i t e r . c ( , N t = 4 ) e s t i m a t e L i t e r . c ( , N t = 4 ) s t d . e r r o r N t , N s = 4 , 1 6 N t , N s = 4 , 2 0 N t , N s = 4 , 2 4 N t , N s = 4 , 2 8 N t , N s = 4 , 3 2 5.685 5.690 5.695 5.700 0 5 10 15 20 25 30 1 L i t e r . c ( , N t = 4 ) e s t i m a t e L i t e r . c ( , N t = 4 ) s t d . e r r o r N t , N s = 4 , 1 6 N t , N s = 4 , 2 0 N t , N s = 4 , 2 4 N t , N s = 4 , 2 8 N t , N s = 4 , 3 2 FIG. S3. F or N t = 4 , a scatter plot of ρ 1 and χ 1 as functions of β with translucent lines to guide the eye. The vertical line and band show respectively the central v alue and the statistical error for the extrap olated β c determined in Ref. [44]. Rew eighting via density of states estimation Our aim is to estimate the lo cation of the p eaks of χ i ( i = 0 , 1 , λ ) so that we may extrap olate to the infinite volume limit N s → ∞ . W e use m ultiple histogram reweigh ting to extract a more precise location of the peak b y estimating χ i ( β ) at a high resolution of β -v alues cov ering the p eak of our simulated v alues. This pro cedure inv olv es learning the densit y of states of the system ρ ( E ) where the partition function may b e expressed Z β = X E ρ ( E ) e − β E . (S1) 5 5.888 5.890 5.892 5.894 5.896 5.898 5.900 0.0400 0.0405 0.0410 0.0415 0.0420 0.0425 0.0430 0.0435 0 L i t e r . c ( , N t = 6 ) e s t i m a t e L i t e r . c ( , N t = 6 ) s t d . e r r o r N t , N s = 6 , 2 4 N t , N s = 6 , 3 0 N t , N s = 6 , 3 6 N t , N s = 6 , 4 2 N t , N s = 6 , 4 8 5.888 5.890 5.892 5.894 5.896 5.898 5.900 0.05 0.10 0.15 0.20 0.25 0.30 0 L i t e r . c ( , N t = 6 ) e s t i m a t e L i t e r . c ( , N t = 6 ) s t d . e r r o r N t , N s = 6 , 2 4 N t , N s = 6 , 3 0 N t , N s = 6 , 3 6 N t , N s = 6 , 4 2 N t , N s = 6 , 4 8 FIG. S4. F or N t = 6 , a scatter plot of ρ 0 and χ 0 as functions of β zo omed into the critical region with translucent lines to guide the eye. The vertical line and band sho w resp ectiv ely the central v alue and the statistical error for the extrap olated β c determined in Ref. [44]. 5.888 5.890 5.892 5.894 5.896 5.898 5.900 0.106 0.108 0.110 0.112 0.114 0.116 0.118 0.120 1 L i t e r . c ( , N t = 6 ) e s t i m a t e L i t e r . c ( , N t = 6 ) s t d . e r r o r N t , N s = 6 , 2 4 N t , N s = 6 , 3 0 N t , N s = 6 , 3 6 N t , N s = 6 , 4 2 N t , N s = 6 , 4 8 5.888 5.890 5.892 5.894 5.896 5.898 5.900 0.4 0.6 0.8 1.0 1.2 1.4 1.6 1.8 1 L i t e r . c ( , N t = 6 ) e s t i m a t e L i t e r . c ( , N t = 6 ) s t d . e r r o r N t , N s = 6 , 2 4 N t , N s = 6 , 3 0 N t , N s = 6 , 3 6 N t , N s = 6 , 4 2 N t , N s = 6 , 4 8 FIG. S5. F or N t = 6 , a scatter plot of ρ 1 and χ 1 as functions of β zo omed into the critical region with translucent lines to guide the eye. The vertical line and band sho w resp ectiv ely the central v alue and the statistical error for the extrap olated β c determined in Ref. [44]. In practice, following Ref. [50], this is achiev ed implicitly by solving via iteration the equation Z β = R X i =1 N i X a =1 g − 1 i e − β E a i P R j =1 N j g − 1 j e − β j E a i − log Z β j (S2) where • R is the num b er of sim ulations each resp ectiv ely conducted at β i with sample size N i , • E a i is the energy measurement (in our case action) of the a -th recorded configuration in a simulation 6 6.050 6.055 6.060 6.065 6.070 6.075 0.0406 0.0408 0.0410 0.0412 0.0414 0.0416 0.0418 0.0420 0.0422 0 L i t e r . c ( , N t = 8 ) e s t i m a t e L i t e r . c ( , N t = 8 ) s t d . e r r o r N t , N s = 8 , 3 2 N t , N s = 8 , 4 0 N t , N s = 8 , 4 8 N t , N s = 8 , 5 6 N t , N s = 8 , 6 4 6.050 6.055 6.060 6.065 6.070 6.075 0.04 0.06 0.08 0.10 0.12 0.14 0 L i t e r . c ( , N t = 8 ) e s t i m a t e L i t e r . c ( , N t = 8 ) s t d . e r r o r N t , N s = 8 , 3 2 N t , N s = 8 , 4 0 N t , N s = 8 , 4 8 N t , N s = 8 , 5 6 N t , N s = 8 , 6 4 FIG. S6. F or N t = 8 , a scatter plot of ρ 0 and χ 0 as functions of β zo omed into the critical region with translucent lines to guide the eye. The vertical line and band sho w resp ectiv ely the central v alue and the statistical error for the extrap olated β c determined in Ref. [44]. 6.050 6.055 6.060 6.065 6.070 6.075 0.066 0.068 0.070 0.072 0.074 1 L i t e r . c ( , N t = 8 ) e s t i m a t e L i t e r . c ( , N t = 8 ) s t d . e r r o r N t , N s = 8 , 3 2 N t , N s = 8 , 4 0 N t , N s = 8 , 4 8 N t , N s = 8 , 5 6 N t , N s = 8 , 6 4 6.050 6.055 6.060 6.065 6.070 6.075 0.15 0.20 0.25 0.30 0.35 0.40 0.45 0.50 0.55 1 L i t e r . c ( , N t = 8 ) e s t i m a t e L i t e r . c ( , N t = 8 ) s t d . e r r o r N t , N s = 8 , 3 2 N t , N s = 8 , 4 0 N t , N s = 8 , 4 8 N t , N s = 8 , 5 6 N t , N s = 8 , 6 4 FIG. S7. F or N t = 8 , a scatter plot of ρ 1 and χ 1 as functions of β zo omed into the critical region with translucent lines to guide the eye. The vertical line and band sho w resp ectiv ely the central v alue and the statistical error for the extrap olated β c determined in Ref. [44]. • and g i = 1 + 2 τ i is a co efficient measuring the auto correlation time b et w een recorded configurations. W e ma y then estimate the exp ectation v alue of an observ able O using ⟨ O ⟩ β ≈ R X i =1 N i X a =1 O a i g − 1 i e − β E a i +log Z β P R j =1 N j g − 1 j e − β j E a i − log Z β j . (S3) Note that it can b e a c hallenge to keep exp onen tials of the form exp {− (∆ β ) S } in a stable numerical range (a voiding underflo w to 0) esp ecially as the action S scales with lattice spacetime volume N t N 3 s . Our algorithmic implemen tation in Ref. [46] uses long double precision and recentring to av oid numerical underflow. 7 5.685 5.690 5.695 5.700 0.08 0.10 0.12 0.14 0.16 0.18 L i t e r . c ( , N t = 4 ) e s t i m a t e L i t e r . c ( , N t = 4 ) s t d . e r r o r N t , N s = 4 , 1 6 N t , N s = 4 , 2 0 N t , N s = 4 , 2 4 N t , N s = 4 , 2 8 N t , N s = 4 , 3 2 5.685 5.690 5.695 5.700 0 5 10 15 20 25 30 35 L i t e r . c ( , N t = 4 ) e s t i m a t e L i t e r . c ( , N t = 4 ) s t d . e r r o r N t , N s = 4 , 1 6 N t , N s = 4 , 2 0 N t , N s = 4 , 2 4 N t , N s = 4 , 2 8 N t , N s = 4 , 3 2 FIG. S8. F or N t = 4 , a scatter plot of λ and χ λ as functions of β with translucent lines to guide the eye. The vertical line and band show respectively the central v alue and the statistical error for the extrap olated β c determined in Ref. [44]. 5.888 5.890 5.892 5.894 5.896 5.898 5.900 0.33 0.34 0.35 0.36 0.37 0.38 0.39 0.40 0.41 L i t e r . c ( , N t = 6 ) e s t i m a t e L i t e r . c ( , N t = 6 ) s t d . e r r o r N t , N s = 6 , 2 4 N t , N s = 6 , 3 0 N t , N s = 6 , 3 6 N t , N s = 6 , 4 2 N t , N s = 6 , 4 8 5.888 5.890 5.892 5.894 5.896 5.898 5.900 10 20 30 40 50 60 70 80 L i t e r . c ( , N t = 6 ) e s t i m a t e L i t e r . c ( , N t = 6 ) s t d . e r r o r N t , N s = 6 , 2 4 N t , N s = 6 , 3 0 N t , N s = 6 , 3 6 N t , N s = 6 , 4 2 N t , N s = 6 , 4 8 FIG. S9. F or N t = 6 , a scatter plot of λ and χ λ as functions of β with translucent lines to guide the eye. The vertical line and band show respectively the central v alue and the statistical error for the extrap olated β c determined in Ref. [44]. In order to estimate the systematic error in the reweigh ting pro cedure, we use the b ootstrap method with N bs = 2 , 000 . More sp ecifically , follo wing Ref. [50], giv en a sim ulation at β i of sample size N i , we estimate the auto correlation τ i b et ween recorded configurations in the sample. Appro ximating g i ≈ τ i , w e then randomly subsample to pro duce an effective sample of decorrelated configurations of size N eff i = ⌊ N i /τ i ⌋ , (S4) where ⌊·⌋ is the flo or function. T o minimise auto correlation, we use the RANLUX random n umber generator defined in Ref. [51]. Since configurations in this b ootstrap sample are effectively decorrelated, this allows us to iteratively solve Eq. (S2) with g i = 1 . W e take a conv ergence tolerance of 10 − 14 (smaller tolerances sho w compatible con v ergence). W e th us hav e a resultan t b o otstrap distribution consisting of 2 , 000 reweigh ting curves. 8 5.2 5.4 5.6 5.8 6.0 6.2 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 L i t e r . c ( , N t = 8 ) e s t i m a t e L i t e r . c ( , N t = 8 ) s t d . e r r o r L i t e r . c ( , N t = 8 ) e s t i m a t e L i t e r . c ( , N t = 8 ) s t d . e r r o r 6.00 6.05 6.10 0.50 0.55 0.60 0.65 0.70 5.2 5.4 5.6 5.8 6.0 6.2 0 5 10 15 20 25 30 35 L i t e r . c ( , N t = 8 ) e s t i m a t e L i t e r . c ( , N t = 8 ) s t d . e r r o r L i t e r . c ( , N t = 8 ) e s t i m a t e L i t e r . c ( , N t = 8 ) s t d . e r r o r 6.00 6.05 6.10 10 15 20 25 FIG. S10. F or N t , N s = 8 , 32 , a scatter plot of λ and χ λ as functions of β ov er a wider range of β – with inset plot zo omed into the critical region. T ranslucen t lines are included to guide the eye. One can see that λ ∈ [0 , 1] as exp ected. The red vertical line and band show respectively the central v alue and the statistical error for the extrap olated β c determined in Ref. [44]. W e extract β c = argmax { χ i ( β ) } , the lo cation of the p eak of each curve, for ev ery sample in the b ootstrap distri- bution. W e then take the mean as the estimate and standard deviation as the standard error. Fitting of numerical results In order to extrap olate the p eaks of our observ ables β c ( N s , N t ) , referred to as pseudo-critical v alues, to the thermo- dynamic limit N s → ∞ , we fit the data using a p olynomial ansatz in inv erse p o w ers of the spatial volume β c ( N s , N t ) = β c ( N t ) + k max X k =1 α k ( N 3 s ) − k (S5) where k max < ∞ defines the degree of the p olynomial. W e p erform a v ariety of fits v arying the range of the data and p olynomial degree. Note that w e discard mo dels with an extremal v alue within the range of datap oin ts since only subleading non-linear corrections are v alid in the thermo dynamic limit. The estimate for β c ( N t ) is the intercept of eac h fit and we estimate the standard error again using the b ootstrap metho d using N bs = 2 , 000 . Rather than select an individual fit, w e pro duce an estimate based on an ensemble of fits. F ollowing R ef. [45], w e use a statistic constructed from the Ak aik e Information Criterion (AIC) that considers go odness of fi t, mo del complexit y and degrees of freedom in the mo del. Roughly , the idea is to give strong weigh ting to mo dels that hav e a high go o dness of fit score, lo w mo del complexity and use all av ailable data p oin ts. The normalised w eights, interpreted as a probability , are calculated using the expression w i = 1 N exp {− 1 2 ( χ 2 + 2 n par − n data ) } (S6) where χ 2 is the standard chi-squared statistic, n par is the n umber of parameters in the mo del, n data is the n umber of data p oin ts used in the fit and N is the normalisation with resp ect to all fits. W e then model the distribution by a w eighted sum of Gaussian distributions with PDF p ( x ) = X i w i N ( x ; µ i , σ i ) (S7) with mean µ i and standard deviation σ i tak en to b e the statistical error in the estimate of the intercept (as estimated ab o ve via the b o otstrap metho d). T aking the CDF of this distribution, we are able to solve numerically for the 16% , 9 Observ able N t = 4 N t = 6 N t = 8 ρ 0 5 . 69247 +0 . 00005 − 0 . 00006 5 . 8942 +0 . 0004 − 0 . 0004 6 . 0645 +0 . 0014 − 0 . 0017 ρ 1 5 . 69257 +0 . 00010 − 0 . 00008 5 . 8942 +0 . 0006 − 0 . 0004 6 . 0625 +0 . 0006 − 0 . 0009 λ 5 . 69240 +0 . 00005 − 0 . 00005 5 . 8940 +0 . 0004 − 0 . 0004 6 . 0625 +0 . 0008 − 0 . 0010 Literature v alue 5 . 69236 +0 . 00015 − 0 . 00015 5 . 8941 +0 . 0012 − 0 . 0012 6 . 0625 +0 . 0018 − 0 . 0018 T ABLE S1. Performing a finite-size scaling analysis using the Betti num ber observ ables, simplicit y and their resp ectiv e susceptibilities, this table sp ecifies our estimates for β c in the thermodynamic limit N s → ∞ for lattices with temp oral size N t = 4 , 6 , 8 . Error bars are computed using the 68% confidence interv al of the weigh ted sum of Gaussians enco ding the v arious fits. As a comparison with the literature, the final row are results from Ref. [44]. This study used comparable MC metho ds to generate its configurations, whic h allo ws us to mak e a faithful comparison b et w een the sensitivit y (to the deconfinemen t phase transition) of our Ab elian monop ole observ ables with standard observ ables. In that study , the exp ectation v alue of the mo dulus of the av erage Poly ako v lo op ⟨| ¯ P |⟩ was used such that the p eaks of the resp ectiv e susceptibilities χ P ( β ) were used to giv e an estimate for β c as N s → ∞ . One can see that our results are compatible with those given in the literature. 0.0 0.5 1.0 1.5 2.0 2.5 N 3 s 1e 4 5.6920 5.6922 5.6924 5.6926 5.6928 5.6930 5.6932 5.6934 L i t e r . c ( , N t = 4 ) e s t i m a t e L i t e r . c ( , N t = 4 ) s t d . e r r o r c ( N s , N t = 4 ) c ( , N t = 4 ) e s t i m a t e F its P olyn. degr ee: 1 N s v a l u e s : 2 4 , 2 8 , 3 2 P olyn. degr ee: 1 N s v a l u e s : 2 0 , 2 4 , 2 8 , 3 2 P olyn. degr ee: 1 N s v a l u e s : 1 6 , 2 0 , 2 4 , 2 8 , 3 2 0.0 0.5 1.0 1.5 2.0 2.5 N 3 s 1e 4 5.6900 5.6905 5.6910 5.6915 5.6920 5.6925 L i t e r . c ( , N t = 4 ) e s t i m a t e L i t e r . c ( , N t = 4 ) s t d . e r r o r c ( N s , N t = 4 ) c ( , N t = 4 ) e s t i m a t e F its P olyn. degr ee: 1 N s v a l u e s : 2 4 , 2 8 , 3 2 P olyn. degr ee: 1 N s v a l u e s : 2 0 , 2 4 , 2 8 , 3 2 P olyn. degr ee: 1 N s v a l u e s : 1 6 , 2 0 , 2 4 , 2 8 , 3 2 P olyn. degr ee: 2 N s v a l u e s : 1 6 , 2 0 , 2 4 , 2 8 , 3 2 P olyn. degr ee: 3 N s v a l u e s : 1 6 , 2 0 , 2 4 , 2 8 , 3 2 FIG. S11. Finite-size scaling analysis for the observ ables ρ 0 (left) and ρ 1 (righ t) both for lattices with N t = 4 and N s = 16 , 20 , 24 , 28 , 32 . Blue triangles represent β c ( N s , N t = 4) , i.e., lo cation of the resp ectiv e p eaks of the reweigh ted susceptibility curv es, with error bars computed via bo otstrapping with N bs = 2 , 000 . Dashed lines represen t the polynomial regression fits used for the infinite volume extrap olation of these peak v alues. The red point is our calculated estimate of the β c in the thermo dynamic limit N s → ∞ using the w eighting pro cedure outlined in the b ody . The horizon tal line and band show resp ectiv ely the central v alue and the statistical error for the extrap olated β c determined in Ref. [44]. 50% , 84% centiles as our low er b ound, estimate and upp er b ound resp ectiv ely . This final estimate represen ts weigh ted con tributions from all fits and so we claim this is more robust than quoting estimates based on any single individual fit. Our estimates of β c are presented in T ab. S1. Plots of scaling analysis for ρ 0 and ρ 1 are presented in Figs. S11, S12 and S13. Plots of scaling analysis for λ are presen ted in Figs. 2, 3 and 4. 10 0 1 2 3 4 5 6 7 8 N 3 s 1e 5 5.893 5.894 5.895 5.896 5.897 L i t e r . c ( , N t = 6 ) e s t i m a t e L i t e r . c ( , N t = 6 ) s t d . e r r o r c ( N s , N t = 6 ) c ( , N t = 6 ) e s t i m a t e F its P olyn. degr ee: 1 N s v a l u e s : 3 6 , 4 2 , 4 8 P olyn. degr ee: 1 N s v a l u e s : 3 0 , 3 6 , 4 2 , 4 8 P olyn. degr ee: 1 N s v a l u e s : 2 4 , 3 0 , 3 6 , 4 2 , 4 8 P olyn. degr ee: 2 N s v a l u e s : 2 4 , 3 0 , 3 6 , 4 2 , 4 8 0 1 2 3 4 5 6 7 8 N 3 s 1e 5 5.890 5.891 5.892 5.893 5.894 5.895 L i t e r . c ( , N t = 6 ) e s t i m a t e L i t e r . c ( , N t = 6 ) s t d . e r r o r c ( N s , N t = 6 ) c ( , N t = 6 ) e s t i m a t e F its P olyn. degr ee: 1 N s v a l u e s : 3 6 , 4 2 , 4 8 P olyn. degr ee: 1 N s v a l u e s : 3 0 , 3 6 , 4 2 , 4 8 P olyn. degr ee: 1 N s v a l u e s : 2 4 , 3 0 , 3 6 , 4 2 , 4 8 FIG. S12. Finite-size scaling analysis for the observ ables ρ 0 (left) and ρ 1 (righ t) both for lattices with N t = 6 and N s = 24 , 30 , 36 , 42 , 48 . Blue triangles represent β c ( N s , N t = 6) , i.e., lo cation of the resp ectiv e p eaks of the reweigh ted susceptibility curv es, with error bars computed via bo otstrapping with N bs = 2 , 000 . Dashed lines represen t the polynomial regression fits used for the infinite volume extrap olation of these peak v alues. The red point is our calculated estimate of the β c in the thermo dynamic limit N s → ∞ using the w eighting pro cedure outlined in the b ody . The horizon tal line and band show resp ectiv ely the central v alue and the statistical error for the extrap olated β c determined in Ref. [44]. 0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 N 3 s 1e 5 6.060 6.062 6.064 6.066 6.068 6.070 6.072 6.074 L i t e r . c ( , N t = 8 ) e s t i m a t e L i t e r . c ( , N t = 8 ) s t d . e r r o r c ( N s , N t = 8 ) c ( , N t = 8 ) e s t i m a t e F its P olyn. degr ee: 1 N s v a l u e s : 4 8 , 5 6 , 6 4 P olyn. degr ee: 1 N s v a l u e s : 4 0 , 4 8 , 5 6 , 6 4 P olyn. degr ee: 1 N s v a l u e s : 3 2 , 4 0 , 4 8 , 5 6 , 6 4 0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 N 3 s 1e 5 6.060 6.062 6.064 6.066 6.068 L i t e r . c ( , N t = 8 ) e s t i m a t e L i t e r . c ( , N t = 8 ) s t d . e r r o r c ( N s , N t = 8 ) c ( , N t = 8 ) e s t i m a t e F its P olyn. degr ee: 1 N s v a l u e s : 4 8 , 5 6 , 6 4 P olyn. degr ee: 1 N s v a l u e s : 4 0 , 4 8 , 5 6 , 6 4 P olyn. degr ee: 1 N s v a l u e s : 3 2 , 4 0 , 4 8 , 5 6 , 6 4 P olyn. degr ee: 2 N s v a l u e s : 3 2 , 4 0 , 4 8 , 5 6 , 6 4 FIG. S13. Finite-size scaling analysis for the observ ables ρ 0 (left) and ρ 1 (righ t) both for lattices with N t = 8 and N s = 32 , 40 , 48 , 56 , 64 . Blue triangles represent β c ( N s , N t = 8) , i.e., lo cation of the resp ectiv e p eaks of the reweigh ted susceptibility curv es, with error bars computed via bo otstrapping with N bs = 2 , 000 . Dashed lines represen t the polynomial regression fits used for the infinite volume extrap olation of these peak v alues. The red point is our calculated estimate of the β c in the thermo dynamic limit N s → ∞ using the w eighting pro cedure outlined in the b ody . The horizon tal line and band show resp ectiv ely the central v alue and the statistical error for the extrap olated β c determined in Ref. [44].