The Static Quark Potential from the Gauge Independent Abelian Decomposition

We investigate the relationship between colour confinement and the gauge independent Cho-Duan-Ge Abelian decomposition. The decomposition is defined in terms of a colour field $n$; the principle novelty of our study is that we have defined this field in terms of the eigenvectors of the Wilson Loop. This establishes an equivalence between the path ordered integral of the non-Abelian gauge fields with an integral over an Abelian restricted gauge field which is tractable both theoretically and numerically in lattice QCD. We circumvent path ordering without needing an additional path integral. By using Stokes' theorem, we can compute the Wilson Loop in terms of a surface integral over a restricted field strength, and show that the restricted field strength may be dominated by certain structures, which occur when one of the quantities parametrising the colour field $n$ winds itself around a non-analyticity in the colour field. If they exist, these structures will lead to a area law scaling for the Wilson Loop and provide a mechanism for quark confinement. We search for these structures in quenched lattice QCD. We perform the Abelian decomposition, and compare the electric and magnetic fields with the patterns expected theoretically. We find that the restricted field strength is dominated by objects which may be peaks a single lattice spacing in size or extended string-like lines of electromagnetic flux. The objects are not isolated monopoles, as they generate electric fields in addition to magnetic fields, and the fields are not spherically symmetric, but may be either caused by a monopole/anti-monopole condensate, some other types of topological objects or a combination of these. Removing these peaks removes the area law scaling of the string tension, suggesting that they are responsible for confinement.


Introduction
An enduring problem in QCD is to find the mechanism which causes quark confinement, which is known to be non-perturbative in its origin. Although several models have been proposed -for example, center vortices [1], and a dual Meissner effect due to magnetic monopoles [2,3,4,5]there has not yet been a convincing demonstration that any of them are correct. In this initial study, we investigate the Cho-Duan-Ge (CDG) Abelian decomposition (sometimes referred to as the Cho-Faddeev-Niemi decomposition) [6,7,8,9,10], and concentrate on one of the signals of confinement, a linear static quark potential. Unlike Dirac and 't Hooft (Maximum Abelian Gauge) monopoles, the CDG decomposition allows for monopole solutions while respecting the gauge symmetry and does not require a singular gauge field or an additional Higgs field. The decomposition is constructed from a colour field, n, which may be built from a SU(N C ) matrix θ, where SU(N C ) is the gauge group of QCD. Each n defines a different decomposition, and an important question is the best way of choosing this field. Recent work [11,12,13,14] has demonstrated that the magnetic part of the field strength dominates the confining string, using a decomposition constructed from one particular choice of θ ∈ SU (N C )/U (N C − 1); however in this case only one of the possible N C − 1 types of monopole is visible.
Here we consider a different choice of θ ∈ SU (N C )/(U (1)) NC −1 . Our initial goal was to investigate whether monopoles apparent in this construction may also lead to confinement 1 . Concentrating on the Wilson Loop, an observable used to measure the string tension, we show that the path ordering may be removed by diagonalising the gauge links along the Wilson Loop by a SU (N C )/(U (1)) NC −1 field θ. This, in principle, allows us to use Stokes' theorem to express the Wilson Loop in terms of a surface integral over an Abelian restricted gauge field strength tensor. We use the CDG decomposition, constructed from a colour field n j = θλ j θ † for a diagonal Gell-Mann matrix λ j , to find a consistent construction of this restricted field across space-time, and not just for those gauge links along the Wilson loop where it was originally defined. Our relationship for the string tension in terms of this restricted field is exact: we do not require any approximations or additional path integrals. We search for topological structures in this field strength, and find that discontinuities in the colour field related to the winding of one of the parameters describing it around a particular location may lead to large 'electric' and 'magnetic' fields (deriving the terms from an obvious analogy to electromagnetism) in both the restricted field strength and the part of the field strength derived from the colour field; in particular the component of the field strength in the plane of the Wilson Loop may exhibit structures, decaying with a power of the distance in all directions, which would lead to an area law scaling of the Wilson Loop. We label these structures merinthons. We also examine the other components of the electromagnetic field strength and find that if certain simplifying assumptions are reliable they will be dominated either by similar peaks or lines of electromagnetic flux. We do not show that this winding occurs in practice, but if it does it provides an explanation for the confining potential. We also discuss how string breaking might arise in this picture at large distances.
The second part of this article checks whether these structures are present using SU(3) quenched lattice QCD. We have generated 16 3 32 and 20 3 40 gauge field configurations, and used them to test this theory. We have concentrated on gauge invariant observables to avoid certain ambiguities with the gauge fixing: the restricted field strength is gauge invariant, but our choice of θ field is not, making it harder to directly observe merinthons in the underlying field. We perform the Abelian decomposition and calculate θ, n, the field strength tensor and the portion of the field strength tensor due to the structures in the θ field. Our goal in this initial study is to show that the Electromagnetic field contains the expected structures, and that the merinthons can account for the string tension. Our simulation is limited because although in principle we should recalculate θ and use a different Abelian decomposition for each Wilson Loop, we do not do so due to limited computer resources; nor have we considered the deconfinement transition, different representations of the gauge group, flux tubes of the confining string, or the relationship between this picture of confinement and the center vortex or dual Meissner models. These would have to be considered in later works.
Preliminary results were presented in [16], [17] and [18]. In section 2 we discuss the Abelian decomposition and its relation to the Wilson Loop and thus quark potential. We discuss the differentiability of the colour field in section 3, and then how confinement may arise in this picture in section 4. After a brief note on string breaking in section 5, we present numerical evidence in section 6 and conclude in section 7. There are four appendices outlining the details of some of the calculations and our numerical methods.

The Abelian Decomposition and Stokes' theorem
The Wilson Loop in an SU(N C ) gauge theory is defined as for a closed curve C s of length L which starts and finishes at a position s, where P represents path ordering and the gauge field, A µ , can be written in terms of the Gell-Mann matrices, λ a , as 1 2 A a µ λ a . We will use the summation convention that the superscripts a, b, . . . on a Gell-Mann matrix implies that it should be summed over all values of a, λ a A a ≡ N 2 C −1 a=1 λ a A a , while the indices j, k, . . . are restricted only to the diagonal Gell-Mann matrices, so, in the standard representation, A j λ j ≡ j=3,8,...,N 2 C −1 λ j A j . We shall often leave the gauge field dependence of W and W L implicit.
The Wilson Loop when C s is an R × T rectangle, with spatial extent R and temporal extent T , can be used to measure the confining static quark potential, V (R), [19] where . . . denotes the vacuum expectation value. It is expected, and observed in lattice simulations, that for intermediate distances the confining potential is linear in R, so V (R) ∼ ρR + c, where ρ is the string tension, and c is a constant. At very small distances, the potential is expected to be Coulomb, while in the presence of fermion loops at very large distances the string is broken and the potential becomes independent of R [20]. The main focus of this work is on the intermediate regime, so we expect that the expectation value of the Wilson Loop will scale with the spatial and temporal extents of the Wilson Loop as This is known as the area law scaling of the Wilson Loop. As discussed in Appendix B, this is satisfied, if on each individual configuration, the Wilson Loop scales as where F is randomly distributed from configuration to configuration (according to one of a certain set of distributions which includes those discussed in this work) with a mean value proportional to the area contained within the curve C s . The eventual goal (and this work is intended as a step towards that goal) is to demonstrate from first principles that equation (4) is satisfied in pure gauge QCD (and later full QCD), and thus that the quarks are linearly confined. A difficulty with evaluating equation (1) is the path ordering. If the fields A µ (x) at different x and µ commuted with each other, then we could ignore the path ordering, and use Stokes' theorem to convert the line integral to a surface integral. The problem would then reduce to showing that there was some flux flowing through the surface so that the surface integral was proportional to the area (perhaps by counting lines of flux). However, A µ are non-Abelian fields: we cannot immediately do this. Our approach is then first to construct an Abelian fieldÂ µ (x) so that W L [C s ,Â] = W L [C s , A], and the calculated string tension does not depend on which of the two fields we use. We may then remove the path ordering, apply Stokes' theorem to replace the line integral with a surface integral over some field strength, and then show that the surface integral is proportional to the area enclosed within the loop.
First we shall define what is meant by path ordering. We split C s into infinitesimal segments of length δσ, and define the gauge link as U σ ∈ SU (N C ) = P[e −ig σ+δσ σ Aσ dσ ] ∼ e −igδσAσ . 0 ≤ σ ≤ L represents the position along the curve and we write A σ ≡ A µ(σ) (x(σ)). We have assumed and will require throughout this work that the gauge field, A, is differentiable. This limits us to only a certain subset of gauges, and once we have found a suitable gauge we are restricted to continuous gauge transformations, i.e. only those gauge transformations which can be built up from repeatedly applying infinitesimal gauge transformations, A µ → A µ + g −1 ∂ µ α + i[α, A µ ], where α ≡ α a λ a and ∂ µ α are both infinitesimal. We also neglect the effects of the corners of the Wilson Loop; the discontinuity in A σ at the corner can, for example, be avoided by using a rounded corner. The path ordered integral over gauge fields is defined as the ordered product of the gauge links around the curve in the limit δσ → 0. W We wish to now replace U σ by an equivalent Abelian field. Previous work [21,22,11] has achieved this by inserting the identity operator, as parametrised below, between each pair of gauge links, where θ ∈ SU (N ), e NC is a unit vector in colour space, and dθ is a suitable Haar measure for θ. Proof that this is an identity operator is given in [11]. This allows one to replace U σ witĥ . The exponent is Abelian, so there is no need for path ordering. In effect, this replaces the path ordering with a path integral over θ. It can be shown that the exponent is related to the CDG Abelian decomposition, which allows for certain monopole solutions which reside in the field θ, and these monopoles may lead to confinement via a dual Meissner effect.
We believe that this approach contains a number of difficulties which have to be resolved before it can be used practically: 1. It introduces an extra gauge degree of freedom to the system -contained within the variable θ; 2. It is necessary to evaluate an additional path integral, making calculations harder; 3. In gauge groups larger than SU (2), only a small subset of the possible CDG monopoles are visible to the restricted field (i.e. only the monopoles contained within the n N 2 C −1 colour direction. n will be defined in equation (9)); 4. There is no obvious relation between the physical gauge field and the unphysical θ field which contains the monopoles and leads to confinement. We expect topological features within the gauge field to contribute to confinement; and this is only possible if the monopoles in the θ field are related to the monopoles in the gauge field. However, given that we integrate over every possible θ, it is difficult to see how the gauge field itself can affect the confining potential if both the CDG monopole picture and this model are correct; 5. It is possible to introduce a gauge transformation which leaves θ constant across space time, destroying any monopole structure; 6. The θ field at each location along the curve is independent (as each needs to be integrated separately), and therefore e † NC (θA σ θ † + iθ∂ σ θ † )e NC will not be a smooth function of the position for a general configuration of θ. This makes the conversion lim δσ→0 δσ → dσ problematic -this is required in the conversion from a product of gauge links to the exponential of a line integral of A σ around the Wilson Loop -and there are also difficulties defining ∂ σ θ. Some regularisation procedure would be required to ensure that θ is smooth. Normally, in a path integral in Euclidean space time, discontinuities if the fields are suppressed by the differential terms in the action, so we can rightly assume that the fields are differentiable. However, in this case θ does not contribute to the action, so there is no similar suppression of discontinuities (see [23] for a more detailed criticism along these lines).
In [11], many of these challenges are resolved by fixing the additional degrees of freedom by a procedure analogous to gauge fixing in QCD. This, in effect, selects one particular θ field and uses that for the subsequent analysis (an alternative way of expressing the approach of [11] is to say that their procedure fixes various components of the gauge field A µ ; but this breaks gauge invariance). However, this breaks the identity in equation (6) and thus the relationship between the Wilson Loop of the restricted gauge fieldÛ and the Wilson Loop of the original gauge field U . In practice, to pursue a program along these lines would require sampling numerous different choices of θ-fixing.
In view of these challenges, we propose that different approach should be taken when trying to express the Wilson Loop in terms of an Abelian field.
We introduce a field θ σ ≡ θ(x(σ)), which, for the moment, we shall take to be an element of U(N C ), at each location along C s and insert the identity operator θ σ θ † σ between each of the gauge links. θ is chosen so that θ † σ U σ θ σ+δσ is diagonal. There are L/δσ gauge links along the path, and we introduce L/δσ θ fields, so there is no obvious reason why the system cannot be solved. In fact, it is easy to construct a solution: it is easy to show (see Appendix A.1) that θ s contains the , for some real α j . As the phases of the eigenvectors are arbitrary, this definition only determines θ up to a (U (1)) NC transformation θ → θχ. χ makes no difference to any physical observable, but for practical purposes it is useful to select the phases and ordering of the eigenvectors by some arbitrary fixing condition to give a unique choice of θ ∈ SU (N C )/(U (1)) NC −1 . We define SU (N C )/(U (1)) NC −1 by considering the following parametrisation of a U (N ) matrix: a i , c i , d 0 and d j are real parameters, and there are N C (N C − 1)/2 Givens matrices (i.e. one for each of the possible ways of embedding a 2 × 2 matrix into a N C × N C matrix) parametrised by one particular a and c. An SU (N C )/(U (1)) NC −1 matrix is parametrised in the same way, but without the final (U (1)) NC term (i.e. by setting d j and d 0 to some arbitrary fixed value, most conveniently d j = d 0 = 0). Under a gauge transformation, U σ → Λ σ U σ Λ † σ+δσ for Λ = e iα a λ a ∈ SU (N C ), θ → Λθχ, where the (U (1)) NC −1 factor χ depends on the fixing condition (in this case Λθ ∈ SU (N ) so there is no contribution to χ from d 0 ). This follows from the definition of θ as containing the eigenvectors of the Wilson Loop. The Wilson Loop transforms under a gauge transformation as W [C s , U ] → Λ s W [C s , U ]Λ † s , so the operator which diagonalises it transforms according to θ s → Λ s θ s ; although we also need to reselect the U (1) NC −1 factor so that the fixing condition remains satisfied. With θ † σ U σ θ σ+δσ = e i λ j diagonal δσu j λ j for real u, removing the non-Abelian structure and the path ordering without introducing an additional path integral.
Our goal is to apply Stokes' theorem to convert this line integral into a surface integral, and this requires extending the definition of θ and u j across the surface bounded by C s . In practice, we construct these fields across all of space time. To generalise θ, we construct nested curves, C i , in the same plane as C s and then stack these rectangles on top of each other in the other dimensions, so that every location in Euclidean space-time is contained within one and only one curve. We then define θ so it diagonalises the gauge Links (and only the gauge links) which contribute to We cannot naively extend u j across all of space-time, because its definition requires that all the gauge links U are diagonalised by θ, not just those that contribute to the Wilson loop, and in general this cannot be satisfied. Instead, we replace the gauge links U with a fieldÛ , defined in a consistent way so that it is both diagonalised by θ across all of space time, and equal to U along the path C s . The first of these conditions means that for each diagonal λ j , which can be re-written in the form This condition is satisfied across all of space-time and for all directions µ. Note that n j is independent of the choice of χ. As we shall see later, the θ-dependence of the restricted field strength F µν [Â] only appears within n j , and objects contributing to the restricted field strength drive confinement. This is the justification of our earlier statement that the choice of χ does not affect the physical observable, which is the restricted field strength. To give thisÛ field a physical meaning we need to relate it to the gauge field U , and we do so via a second fieldX defined according tô For later convenience (equation (50)), we restrictX µ by imposing the condition Under a gauge transformation n transforms as n x → Λ x n x Λ † x (which follows from the transformation rule for θ) and the requirement that equations (9) and (11) are satisfied in every gauge leads to the transformation rulesÛ (9) and (11) are the lattice versions of the defining equations of the gauge independent CDG decomposition [6,7,8,9,10], which in the continuum is described by 2 tr(n j X) =0 (14) withÛ Equation (12) is the naive continuum limit of equation (10); equations (13) and (15) are the naive continuum limit of equation (9), and equation (14) is the continuum limit of equation (11). Proof that equation (16) solves the continuum decomposition has previously been given in (for example) [6,7,11], and is repeated in Appendix A. 2 The original discoverers of this decomposition prefer to call it a gauge independent decomposition rather than gauge invariant decomposition. In these earlier models, the choice of θ was left arbitrary. The procedure was to fix to some arbitrary gauge, and then select some θ. Since the decomposition only depends on the choice of θ, and proceeds regardless of which gauge was originally chosen, the decomposition was referred to as gauge independent to distinguish it from other approaches which required fixing to particular gauges, such as the decomposition based on the Maximum Abelian Gauge. Our work differs in philosophy from the original approach because we do not need to perform any gauge fixing to extract our observables, since our key quantities (such as the field strength and Wilson Loop) are gauge covariant and thus the corresponding observables are gauge invariant. However, our particular choice of θ is gauge dependent, and thus the quantities which we use to parametrise it can only be examined after gauge fixing to some arbitrary gauge. The authors of [6,7,8,9,10] gauge fixed to an arbitrary gauge and then performed an Abelian decomposition; we decompose and then (if required, which isn't the case for any physical observable) gauge fix.
The condition (11) may then be interpreted as requiring that tr(n jÂ ) = tr(n j A), the component of the gauge field A parallel to n is fully contained withinÂ. In the continuum, the solution forÂ and X is unique, but on the lattice we found that there were sometimes several distinct solutions to equations (9) and (11). In this case, we choose the solution which had the largest value of tr(X), a condition which is both gauge invariant and satisfied along C s whereÛ = U and thuŝ X = 1.
The continuum defining equations ensure that the field strengthF µν associated withÂ, defined by for any field α in the adjoint representation of the gauge group, satisfiesF µν [Â] = β j µν n j for some real scalars β j . The proof of this follows by substituting each of the n j fields in turn in place of α, using equation (13) to show that [F , n j ] = 0 for all the n j , and noting that the only objects which commute with each of the n j are the other n j fields.
We express the restricted field asÛ µ,x ≡ θ x e iλ j δσû j µ,x θ † x+μδσ for realû, and sinceÛ = U along the curve C s , we see that whereF j (likeû) is gauge invariant, Σ the (planar) surface bound by the curve C s , and dΣ an element of area on that surface. Note thatû does depend on the fixing condition, althoughF is independent of it. Equation (20) can be easily proved by considering for smooth θ andû which implies that and which follows from the definition n j = θλ j θ † .
We can now consider the Wilson Line around an infinitesimal plaquette p, which for a smootĥ u field gives which leads to and building the integral over the surface bounded by C s from the product of integrals over these small plaquettes, using that the exponent is Abelian, gives equations (19) and (20). Proof that the definitions ofF given in equations (18) and (23) are equivalent is given in Appendix A, along with a proof that the definition ofÂ σ in equation (22) is equivalent to that given in equation (16).
Finally, as alluded to earlier, we note thatF µν [Â] can be written in the formŝ as is proved in Appendix A, equations (A.25) and (A.27). These functions solely depend on n, and θ only indirectly through n, and since n is independent of χ andF is the physical observable we want to study, we conclude that the choice of χ will not affect any of the physical observables we need. Equation (19) is only valid ifû is differentiable. Equation (19) is also similar to what we see in QED, which is, of course, not confining. In analogy to QED, we may expect the contribution of those portions of space time whereû is continuous to have little contribution to the string tension. However, we must also add to this equation the effects of discontinuities inû. We do so by only extending the area integral over those areas whereû is continuous, and add additional line integrals around the areas where it is discontinuous. The linear string tension will, at least in part, arise from these discontinuities. Sinceû j is built from the gauge field, and we are working on a gauge where this is assumed to be differentiable, and θ, we must therefore consider whether θ is differentiable.

SU(2)
We wish to investigate how θ varies over space, which requires evaluating the derivative of the eigenvectors of the Wilson Loop operator. Suppose that an operator M has two eigenvectors (which is the case for SU(2)), ψ 1 and ψ 2 , with eigenvalues λ 1 and λ 2 and after a small change in M → M + δM , the eigenvectors change to ψ ′ 1 and ψ ′ 2 , with eigenvalues λ 1 + δλ 1 and λ 2 + δλ 2 . For normal M (so that the eigenvectors are orthogonal, left and right eigenvectors the same, and δM is (anti-)Hermitian or an anti-Hermitian matrix multiplied by M if M is unitary), we may parametrise the new eigenvectors as and the eigenvalue equations read which, after applying ψ † 1 and ψ † 2 to these equations gives At |λ 1 − λ 2 | ≫ |ψ † 1 δM ψ 1 − ψ † 2 δM ψ 2 |, this agrees with the result from first order perturbation theory [24] In this case, the change in the eigenvector is proportional to δM . However, if the eigenvalues are near-degenerate (i.e. the difference between them is of O( δM )), then β may be large even as δM → 0, and the eigenvectors would evolve discontinuously. As we will be interested in the spatial derivatives of the eigenvectors of the closed Wilson Loop, I will set M to be the Wilson Loop operator W [C r s ], and consider the derivatives in two directions, one parallel to the Wilson Line, and one perpendicular to it in the plane of the Wilson Loop. Any small shift in location in an arbitrary direction can be constructed in terms of a shift parallel to the Wilson Loop followed by a shift perpendicular to the Wilson Loop. r is an index denoting the size of the Wilson Loop; an increase in r means that we move to the next of the nested Wilson Loops used to define θ, so both the spatial and temporal extents of the curve C r s would increase by a small amount 2δr. For an R × T planar Wilson Loop with T > R we have r = R/2. We investigate the change of the Wilson loop both when we have a small shift along the curve, s → s + δσ, and a small shift perpendicular to the curve, r → r + δr.
Firstly, there will clearly be discontinuities in θ when the gauge field is discontinuous, but as we require a continuous gauge field for the decomposition to be valid we neglect this case. We wish to consider if there can be discontinuities in the eigenvectors of the Wilson line even for a continuous gauge field.
In the direction parallel to the line, we have and with U s = e −igδσA r s , we have for a smooth gauge field This gives (if the eigenvalues of the Wilson Loop are not degenerate) β is proportional to δσ, so, baring changes in the phase, the eigenvectors evolve smoothly around the Wilson Line.
Considering changes perpendicular to the Wilson Line, we have which gives, since λ 1 λ 2 = 1 for W ∈ SU (2), We see that, expect when the eigenvalues are near degenerate, β will be proportional to δrthough not δσ as s,s ′′ δσ 2 . . . ∼ O(1). This means that the θ field is discontinuous when the Wilson Loop has degenerate eigenvalues -and we shall investigate the consequences of this later. For now, we just consider the other situations with infinitesimal β, and ask if it is possible that there can be a large change in θ without a large change in β. We obtain where the phase δ is chosen so that whatever fixing condition we have chosen is satisfied. We denote B r s as the appropriate transformation when performing a small translation in the direction of increasing s and B s,r when translating in the direction of increasing r. In section 5, we shall use these relations in our discussion of string breaking.
There is one further discontinuity as we evolve θ over s and r. We may parametrise θ as θ = cos a i sin ae ic i sin ae −ic cos a with 0 ≤ a ≤ π/2, c ∈ R and d ∈ R. Note that the parameters a, c and d are not gauge invariant, so this parametrisation is only defined after we have fixed to a suitable gauge. The phase d makes no physical contribution, and should be fixed by some suitable condition, but for completeness we still give it explicitly. We perform a small translation in either the r or s direction, and obtain θ ′ which is parametrised by a ′ , c ′ and d ′ : which allows us to read off the new parameters a ′ , c ′ and d ′ , cos a ′ = cos 2 a cos 2 β + sin 2 a sin 2 β − 2 cos a cos β sin a sin β sin(c + γ − 2d) cos a sin β cos(2d − γ) + sin a cos β sin(c) − cos a sin β sin(2d − γ) + sin a cos β cos(c) .
It is clear that when cos a and sin a ≫ β, we have a ′ = a + O(β), c ′ = c + O(β) and d ′ = d + O(β), which means that θ is differentiable with respect to the spatial coordinates at least once. However, if a ∼ π/2 we have d ′ = c + γ − d + δ + (ν + 1 2 )π and c ′ = 2c − 2d + γ + (ν + 1 2 )π, where ν is an integer, and, at a ∼ 0, c ′ ∼ 2d − γ + (ν + 1 2 )π. At both a ∼ 0 (where θ is diagonal) and a ∼ π/2 (where θ is off-diagonal), the parameters which describe θ are discontinuous even when β is small. a itself, however, remains continuous and the eigenvalues of θ are dependent only on a once we fix the phase d. δ can be chosen according to our rule for fixing the U(1) phase d. For example, if we enforce d = d ′ = 0, then δ = − arctan sin a sin β cos(c+γ) cos a cos β−sin a sin β sin(c+γ) . This discontinuity of θ may be in principle be at a single point, or there may lines in the surface spanned by r and s where these discontinuities occur (if they occur), and these lines will either be open, form closed loops or will be infinite in extent. If they are lines, then we require that a ′ = a = π/2 (for example), which means that θ will be unchanged as we make a small displacement in space in some arbitrary direction. If θ is unchanged, then the Wilson Loop must also be unchanged as we move across space time. However, baring a remarkable gauge field configuration, it is unlikely that this will occur for any more than an infinitesimal distance. We may as well then neglect the possibility that these discontinuities are in lines, and focus on the case that they are at isolated points.

SU(3)
In SU(3), we introduce the following operators which form a complete basis of the Hermitian traceless generators of SU(3) θ 31 =i cos a 3 cos a 1 sin a 2 e −ic2 − sin a 3 sin a 1 e −i(c1+c3) θ 32 =i sin a 3 cos a 1 e −ic3 − cos a 3 sin a 1 sin a 2 e i(c1−c2) Note that d 3 and d 8 can easily be constructed from the phase of the 11 and 33 components of θ, c 2 and a 2 from the 13 component once the phase has been removed, and a 3 and c 3 and a 1 and c 1 from the 23 and 12 components of θ. We can again construct the differential of θ with respect to a coordinate by multiplying θ by a matrix B constructed in the same parametrisation as θ: with β ij taken from equation (45) but with the substitution Again, except in the neighbourhood of degenerate eigenvalues of the Wilson Loop, β will be infinitesimal and we can neglect terms of O(β 2 ) and higher. With the shifted θ ′ = θB, we have where, for the d = 0 fixing condition, δ is tuned to ensure that θ ′ 11 and θ ′ 33 are real. Once again, we see that a i , c i and d i will be differentiable at least once where β is continuous except where either a 1 = π/2, a 2 = π/2, a 3 = π/2, a 1 = 0, a 2 = 0, or a 3 = 0, where the parameters c i are discontinuous (possibly because the δ required to fix d is not infinitesimal).

SU(2)
Now suppose thatû j contains a non-analyticity. We integrate the field around a loopC parametrised byσ surrounding the discontinuity, bounding the surface integral by an additional line integral C dσû j σ , far enough away from the discontinuity thatû j σ is analytic alongC. We define {C n } as the set of curves surrounding all these discontinuities, andΣ the area bound within these curves. We can write and sinceû is continuous onC, after fixing the gauge we can expand U = 1 − i 1 2 gδσA a λ a and θ † x θ x+δσ = 1+δσθ † ∂σθ. We define X 0 ≡ 1 2 θ † (X +X † )θ. For example, in SU(2) we may parametrisê X asX thenX +X † = 2 cos x cos wI = 2X 0 , where I is the identity operator (this expression will obviously not be exact in higher gauge groups). If everything is analytic alongC s , and we consider infinitesimal displacements, we expectX andÛ to be close to the identity operator, and x and w will both be O(δσ). X 0 is thus I + O(δσ 2 ). This gives Using (11) the first term in equation (50) gives zero, while the second term will not contribute to an integral aroundC if U and X 0 are continuous and the area of the loop is small enough. We therefore concentrate on the contribution from the final term. Equation (20) is then replaced by where dΣ is an element of area.
In SU(2), after fixing to a suitable gauge (where the gauge field U is smooth), we parametrise θ as in equation (39), with c ∈ R, 0 ≤ a ≤ π/2 and d 3 determined by the fixing condition. As this contains the eigenvectors of W [C s ], it is differentiable except where W [C s ] has degenerate eigenvalues and those locations where a = 0 or a ≈ π/2, where c is ill-defined. The parameter c may wind itself around those points where a = π/2 or a = 0. We may parametrise the plane of the Wilson Loop in polar coordinates (r, φ), with the origin at the location where a = π/2 or a = 0. The continuity of θ away from r = 0 implies that at infinitesimal but non-zero r, c(r, φ = 0) = c(r, φ = 2π) + 2πν n for integer winding number ν n , and if c is ill defined at r = 0 then we may find that ν n = 0. It is straight-forward to calculate We integrate around a curve at fixed a, and we may choose a fixing condition which keeps d constant. This leads to where the last equality is only valid if X 0 is close to the identity operator, as we expect. In SU(2) (but not higher gauge groups, where the expression depends on more parameters) we discover that were we to use an infinitesimal loop around the discontinuity, then the area integral gives θ † s W [C s ]θ s ∼ e 2πiνλ 3 for integer ν, which cannot contribute. However, if we integrate at a larger radius, where a has fallen to some average background value a 0 , we do see a contribution to the Wilson Loop from the discontinuity, proportional to 2πν n sin 2 a 0 , and we may still find a large contribution to the string tension. The structure in the field is extended over a small region of space rather than being a δ-function spike. This conclusion may also be reached by considering gauge invariance: whileF is gauge invariant, θ is not and therefore the precise location of the discontinuity in θ will depend on the gauge (although it should not depend too strongly on the gauge). The peak inF cannot therefore be precisely at the discontinuity in θ in every gauge, so it cannot be a δ-function peak around that centre. We shall discuss gauge invariance in more detail in section 4.4. We will refer these objects in the field strength as CDG merinthons. They are topological objects in the homotopy group π 1 (distinct from monopoles which have the homotopy group π 2 , and instantons which are of the group π 3 ). They can be identified by a large value ofF , a = π/2 or a = 0, and a non-zero winding of the parameter c around the merinthon.
It is also important to note that the merinthons at a = 0 and a = π/2 with the same winding contribute with opposite signs. If we compare the results for the integrals of the restricted potential around fixed curves at a = δ and a = π/2 − δ, for some small δ, we find that the term proportional to the winding number gives When this is multiplied by i and exponentiated, e i2πν = 1 can't contribute to the Wilson Loop, so we can neglect this part of the expression and write, Thus with the same winding number, the two topological objects contribute with opposite signs to the Wilson Loop. This is important as we consider how the merinthons with positive and negative winding numbers can be distributed across the lattice (see figure 1). Assuming that c is continuous apart from at the merinthons, then we must have a positive winding number object next to a negative winding number object. Thus two a = π/2 objects next to each other (with opposite winding numbers) will give contributions to the Wilson Loop which approximately cancel (ultimately leading to a total contribution proportional to the square root of the number of objects, i.e. a perimeter law scaling for the Wilson Loop), and the same can be said for two a = 0 objects neighbouring each other. However, if we have a a = π/2 merinthon with positive winding number situated as a spatial neighbour to an a = 0 merinthon with negative winding number, then the contributions of the two objects would be additive (ultimately leading to a total contribution proportional to the number of objects, i.e. an area law scaling). 3 The total Wilson Loop will be the product of a boundary term and contributions from all the merinthons contained within the planar Wilson Loop. As we can expect the number of merinthons to be proportional to the area of the loop, this leads to an area law scaling for the Wilson Loop and a linear string tension.

SU3
The SU(3) system is considerably more cumbersome than the SU(2) theory. We start by constructing the operators n 3 = θλ 3 θ † , n 8 = θλ 8 θ † and θ † ∂ µ θ, where we parametrise θ = Figure 1: An Illustration of how topological objects can be distributed across a plane. The cross indicates the centre of the topological object. The circles and arrows indicate the winding of the parameter c around each topological object. In this simplified picture, it is clear that if c is continuous outside the topological object, the arrows of neighbouring objects must point in the same direction at the point of contact between them. In this simple picture, each object of positive winding number (clockwise arrows) must have objects with negative winding number neighbouring to it. While this can be avoided by having the direction of flow of c doubling back on itself, such a scenario would require that the flow should be three times as large on another part of the circle, meaning that the problem is not eliminated but just transferred elsewhere. e ia3φ3 e ia2φ2 e ia1φ1 e id3λ 3 +id8λ 8 . We find that where explicit expressions for n 3 , n 8 and S 1µ , . . . , S 4µ are given in Appendix C.
Close to the merinthon we can write the Abelian decomposition as which gives The two components of the field strength will be given by the integral ofû 3 orû 8 around a loop surrounding the merinthon at an a i that have fallen to some background level. If the fieldsX and U are smooth, and with θ differentiable along this path, we havê which implies that (using equations (C.4) -(C.8)) The components of the field strength tensor are then given by (excluding terms which will not contribute either because they are multiples of 2π or zero) where dΣ is an element of area which covers the merinthon andC the curve which surrounds the merinthon. Only those c i s which wind around the merinthon will vary on this infinitesimal loop; all the other parameters can be treated as constant. These field strengths are non-zero if one of the c i has a non-zero winding number around the merinthon, and they will not always be an integer multiple of 2π. In particular, for the terms proportional to ∂ µ c 2 and ∂ µ c 3 we may pick infinitesimal loops where the a parameters do not vary significantly. The merinthons thus will lead to a non trivial structure in the restricted field strength, and if the number of merinthons within an area is proportional to that area, an area law for the expectation value of the Wilson loop and confinement.

The Electromagnetic field
We saw in equation (20) that the field strength can be written as, In [6], the potential is split into two parts A E = tr n j A µ and A M = 1 2g itr n j θ∂ µ θ † . The field strength, is often split into two components, 4 The magnetic current k µ can be defined through Maxwell's equation and while 1 2 ǫ µνρσ ∂ νÊρσ = 0, as in electromagnetism, in general 1 2 ǫ µνρσ ∂ νĤρσ = 0. This means that this field strength, coming from the colour field rather than the gauge field, is associated with a magnetic source. Note, however, that it is not in general a pure monopole (except for certain specific choices of θ), since the electric current, defined through ∂ µĤ µν = 0. NeitherĤ µν norÊ µν are by themselves gauge invariant, but onlyF =Ê +Ĥ. We expect that close to the merinthon solutionsF will be dominated byĤ since this is the part ofF which is sensitive to the winding of c around the centre of the merinthon. In analogy to QED, we shall refer to the components of F (and alsoĤ) as the electromagnetic fields.

SU(2)
In SU (2), we see that (using equations (53) and (A.25) and setting d = 0), The Yang-Wu monopole [25,26,27] is one classic magnetic monopole solution for winding number ν = 1. Here we parametrise the coordinates as (with the origin at the location of the monopole) The reason we choose this parametrisation is so that we can easily respect the symmetries of the theory. Since we are considering point-like objects, it is most natural to use a spherical polar coordinate system. QCD, after averaging over the gauge links, possesses a spherical symmetry, which we have broken by inserting the Wilson Loop in the xt plane. However, we should still be free to rotate the coordinate axes in the yz plane. This particular parametrisation allows us to absorb such rotations into a shift of ψ 1 . If, however, we parametrised the coordinates differently (for example by exchanging t and z in equation (73)), a rotation in the yz plane would require some non-trivial change of all the angular variables. The only changes we can easily make to the coordinate parametrisation are to interchange the t and x axes or rotate in the yz plane (we shall only consider rotations of π/2 since we are primarily interested in comparing with a lattice theory). For the Yang-Wu monopole we select the solution with c = ψ 1 and a = ψ 2 /2. 5 This leads to the solution for the electric and magnetic fields E = 0, B = 1 2g rs r 3 s , where r s = (x, y, z); in other words the true magnetic analogue of a point like electric charge. Note that there is no time dependence in the monopole field strength: these objects will reveal themselves as stings of electromagnetic flux extending in the time direction.
The literature concentrates on the possibility of finding these monopoles in the colour field, which are then assumed to condense and provide confinement through a dual Meissner effect. We will later, therefore, in addition to the merinthon solution suggested above, consider what structures may appear in the fields if these Yang-Wu monopoles dominate.
In our case, we have something a little different. We will start by considering the coordinate system given in equation (73). We consider the objects close to a = π/2 here; the objects at a = 0 can be treated in a similar way. We can then write the variables, at all locations away from the singularity, as a = π/2 − G(r, ψ 1 , ψ 2 , ψ 3 ) where G is zero at r = 0 and c = R(r, ψ 1 , ψ 2 , ψ 3 ); and for a winding number one merinthon we require R(ψ 3 + 2π) = R(ψ 3 ) + 2π. The functions R and G will depend on the gauge, but in a way which allowsF 3 to remain gauge-invariant. This gives and .
The simplest possibility is to set G = G(kr) for some k and R = R(ψ 3 ), neglecting the ψ 1 and ψ 2 dependence. We shall consider a more general solution later. Using equation (72), we obtain, where R ′ and G ′ indicate the first derivative of these functions and E 3 is the electric field. Similarly, For the magnetic field B 3 we find, Here, the merinthon generates a pure electric field. This is, of course, not the only way we can parametrise the coordinates. For example, we could have chosen the coordinate system which in effect swaps t and x. This gives a different solution for the field strength tensor with the other components ofĤ 3 zero. We could also use any other rotation of the coordinates consistent with the space time symmetries. For a general function G(r, ψ 1 , ψ 2 , ψ 3 ) and R(ψ 1 , ψ 2 , ψ 3 ) (neglecting any radial dependence of c, which may in practice only be valid for small enough distances and certain choices of gauge), with ∂ i ≡ ∂/∂ψ i and r 2 yz = y 2 + z 2 and r 2 xyz = x 2 + y 2 + z 2 , we find for the parametrisation (73) We must also consider every permutation and rotation of the coordinates to obtain other possible solutions.
At small r, each of these terms gives a finite and non-zero result if the derivatives of G and R are non-zero and both G and ∂ ψi G are proportional to r. Our ultimate goal is to measure the electromagnetic fields on the lattice to try to identify these structures, and this means evaluating how the fields behave at large and intermediate distances. This is complicated by the unknown functions G and R (which are also not gauge invariant and therefore uniquely defined). We require that at least one of the ∂ ψi R = 0 (and, in practice, on average take an integer value). We can, however, study the explicit coordinate dependence within equation (80) for G ∼ r (obviously this is in practice only an approximation, as we only consider part of the expression for the field strength, but nonetheless we may expect these features to be visible to a certain extent). We can see three types of behaviour at large distance. Firstly, there are terms such as ( yx ) which fall away at least as 1/z i for large spatial distances, while increasing at large t if ∂ i G is linear in r -we find a possible solution where the field grows with distance, although we may hope that these are excluded in practice. Secondly, we have terms which form a line of large electromagnetic field in one of the directions, for example inĤ 3 xz we find the term and the maxima of this function fall off at large distances from the origin as 1/z 3 , 1/y 2 and 1/t but remain roughly independent of x. This will correspond to a line of large magnetic field (which may be broken or change sign) in the x-direction, and we will label it as a µ-string, where µ indicates the direction (in practice, the line of field strength will only have a finite extent because of the G and R dependence). We do not mean anything by this terminology beyond that the large field may extend in the µ-direction and falls of rapidly in the other directions. We require that the leading order r dependence in G should be linear if we wish to avoid zero or infinite field strengths at r = 0. Alternatively, by permuting the coordinates, we also find these strings may exist in the other electromagnetic fields. Finally, we have what we will label as points where the electromagnetic flux falls as at least 1/r in each direction at large distances, r 1/k, for example the term sin 2G∂ r G∂ 1 R 1 r inĤ 3 yz . We will concentrate on solutions with ∂ 1 R = 0 and ∂ 2 R = 0. Our reason for doing so is that we expect the discontinuity to occur in a point; thus there is no obvious reason why the objects should wind around this discontinuity in the yz plane and not the xt plane; if the origin is at the discontinuity and c proportional to the angle we might expect c(−δx µ ) = c(δx µ ) + π, and only setting c = ψ 3 accomplishes this. If there is no winding number related to variations in ψ 1 and ψ 2 , then we might expect any fluctuations in the fields in these directions to cancel out as we integrate around a loop surrounding the discontinuity. Also, only winding in the xt plane can contribute to confinement. In our numerical simulations, we only use Wilson Loops in the xt-plane, so we restrict ourselves to these loops here. This means that we should expect a symmetry between the structures in the x and t components of the electro-magnetic field, and also between the y and z components of the field, but not necessarily between (for example) the x and y components. This all suggests that it is safe to neglect this angular dependence in R. With this simplification, we obtain,Ĥ In this parametrisation of the coordinates, the electromagnetic fields will therefore appear either as points (not necessarily spherically symmetric, but falling in every direction) or strings in either the x or t direction. We also consider the three other parametrisations of the coordinates which arise from the exchanges x ↔ t and y ↔ z. We summarise the results in table 1. With these choices of coordinate parametrisations, we see thatÊ x will only have points,Ê y and E z strings along the x axis,B x may have a string along the x or t axis, whileB y andB z strings along the time axis. Of course, in practice this picture will be obscured by fluctuations in the fields and effects from the unknown functions G and R; we also cannot say how long the strings will be. It also depends on that G has some angular dependence, and in particular has some variation in the yz plane (as parametrised by the angle ψ 1 ). Equally, the strings depend on different angular differentials of G; if G has no angular dependence they will not be present, and if it depends on ψ 2 but not on ψ 1 only some of these strings will appear. We will search for these patterns of electromagnetic fields in in lattice QCD in section 6.4. We can also consider the patterns we can expect from the Wang-Yu monopole and other solutions for the fields. We have already seen that the monopole solution with the coordinate parametrisation given in equation (73) gives t-strings in the magnetic field; we may also study the other options when we permute the coordinates. Additionally, we list a few similar simple alternative forms for the functions R and G. These are summarised in table 2. It should be observed that if a is linked to an angular variable then eitherÊ x = 0 orB x = 0; while if c is linked to the angular variable ψ 2 rather than ψ 3 with G ∼ r, then the Electric fields have t-strings and the magnetic fields x-strings: the opposite of our expected model in table 1. We note that a combination of a = ψ 3 /2; c = ψ 1 monopoles and a = ψ 3 /2; c = ψ 2 structures may have the same pattern of points and strings as the configuration of table 1; the difference is that there is in this case no correlation between the spatial locations of the strings and the points in theÊ 3 x field. That the structures in the field strength are not point like but contain vortices is important because it means that we cannot evade the large field strength by distorting the surface Σ. For example, we may create a bump in Σ to avoid the merinthon and the structure in theF xt field: however, this will require integrating over a surface in (for example) the yt plane whereF yt is large, and the yx plane where F yx is large. This means that it is not naively ruled out that the surface integral bound by the Wilson loop is not affected by the choice of surface to integrate over.
gives a similar picture as SU (2); here we only outline the argument about which structures will emerge in the field strength. The starting point is equation (64), and by differentiating this we can construct the two field strengthsF 3 andF 8 . For example, we find that The first two terms are similar to the quantities already discussed in SU (2), and the same analysis applies. The third term, mixing the differential of c 3 with that of a 2 , is new and we need to consider this as well.F 3 µν is similar toF 8 µν , with similar expressions arising, plus additionally terms proportional to the product of the differential of to of the c i s, for example ∂ µ c 1 ∂ ν c 2 . We therefore expect to see the same structures as in SU (2), plus possibly some additional structures from these extra terms. For example, inF 8 it with c 3 ∼ ψ 3 , we find the contribution xit ∂xi , while inF 8 ij we find the terms xit r 2 rxzy ∂a2 ∂xj − xj t r 2 rxzy ∂a2 ∂xi . In general a 2 , unlike a 3 , will not be at its minimum or maximum value at the merinthon . There is therefore no reason to require that ∂a 2 /∂ψ 1,2 will be proportional to r at small r (we can discount the possibility that c 2 and c 3 both have winding around the same point). We can also neglect terms such as ∂ µ c 1 ∂ ν c 2 because we do not expect both c 1 and c 2 to wind around the same point. This means that these additional objects in the electromagnetic field strength will just be points. The same argument applies to all the additional terms in theF 3 µν . Therefore, the electromagnetic structure in SU(3) will be the same as in SU(2) but there might be additional point-like structures appearing in all the components of the electromagnetic field; this applies equally for both of the electromagnetic field tensors,F 3 µν andF 8 µν .

Gauge Invariancê
F j is gauge invariant, which is what we expect for the mechanism behind confinement. However θ, X and U are dependent on the gauge. This presents a conceptual difficulty, since we have stated that the gauge invariant field strength is influenced by objects which arise due to features in the gauge dependent parametrisation of θ. The key contribution, however, comes from the winding of c around the merinthon, and we here show that this is gauge-invariant (at least in SU (2)).
Since this formalism requires that the gauge field is continuous, we can only use continuous gauge transformations. However, to undo the winding in θ we require a discontinuous gauge transformation. Thus the discontinuity in θ at a = π/2 will survive any smooth gauge transformation, and the merinthons will remain. For example, in SU(2), we can parametrise an infinitesimal gauge transformation as with l 1 and l 3 infinitesimal and 0 < l 2 < 2π. Recall that under a gauge transformation θ → Λθ.
We fix d 3 = 0, so we remove the matrix on the right by a U (1) rotation. From the remaining term, we can read off e ic ′ as the phase of the top right component of the matrix on the left. In particular, we see that If a ∼ π/2 and a ∼ 0, we can read off, to lowest order in l 1 and l 3 and given that c = c ′ at l 1 = l 3 = 0 (i.e. there is no gauge transformation), The winding around the merinthon becomes λ 3 1 2 tr[X 0 ] ∂σc ′ dσ = λ 3 1 2 tr[X 0 ] ∂σ(c ′ − c)dσ + λ 3 πνtrX 0 , and since l 1 and l 3 are continuous functions and c only contributes to (c ′ − c) within a trigonometric function, the winding of θ around the merinthon is unaffected by a continuous gauge transformation. The location where a = π/2 or a = 0 may, however, be shifted by a small amount. The merinthons , of course, are not δ-function peaks, but extended over a small area.

Wilson Loop Dependence of the Merinthon Solution.
Each merinthon is constructed from a particular configuration of the θ field, and each θ field is given by the eigenvectors of a particular set of nested Wilson Loops. This means that in addition to the gauge field, the particular topological objects we find will depend on the choice of which set of Wilson Loops is being studied. So if we choose one particular set of Wilson Loops, we will uncover one distribution of merinthons, while if we choose a different set of Wilson Loops, we may well uncover a different distribution of the objects. An obvious question is whether this eliminates the usefulness of the merinthons as an explanation for confinement. After all, an instanton, Dirac magnetic monopole, or vortex are dependent only on the gauge field, and so can be measured without making any arbitrary choices concerning the θ field. The global topological charge of the field is an observable which is independent of any choices that we make. If this can be reduced to a set of local topological structures, then the natural thing is to say that these structures should also be independent of any arbitrary decisions we should make when parametrising the gauge field (such as our choice of Wilson Loop).
Obviously this issue would not be relevant if the distribution of merinthons was (to a good approximation) independent of our choice of Wilson Loop. However, our numerical tests indicate that this is not the case: we see no significant correlation between the topological field strength when we use one set of Wilson Loops to construct θ compared to when we use a different set of Wilson Loops.
The Wilson loop represents the propagation of a static quark/anti-quark pair through time. The choice of Wilson Loop is equivalent to choosing where the quark/anti-quark pair under study is located. So the question posed in the first paragraph of this subsection is asking whether the topological objects which explain the confinement of quarks (if quark confinement is explained by topological objects) should manifest themselves independently of the presence of those quarks. There is no obvious reason why that should be the case. Discussion of confinement makes no sense unless there are quarks present to confine; the study of confinement in general can be reduced to the study of why each particular pair of quarks is confined. As the explanation of the confinement of a particular pair of quarks, where the position of the quarks determines the Wilson Loop which then fixes the θ we require for the decomposition and uniquely determines the topological structures contained within θ, the topological objects under study here work as well as those which are present in the gauge field independently of any choice of Wilson Loop. There is no obvious reason why the key quantity that determines the strength of the confining potential (which is an observable independent of which quarks are being studied), the density (after averaging over the gauge fields) of merinthon -anti-merinthon pairs, should be dependent on which Wilson Loop is chosen. Thus we feel that the arbitrariness of the merinthon fields does not weaken this model's use as an explanation of quark confinement.

String breaking
So far, we have not discussed the effects of the discontinuity in θ when the Wilson Loop has degenerate eigenvalues. Although the main focus of our work is elsewhere, some discussion of this topic is required for completeness; however we will do no more than present a simple result leaving the questions arising for the physical consequences of this result unanswered. We will find that if the Wilson Loop has degenerate eigenvalues, then the confining potential will weaken at larger distances. This is similar to the effect known in dynamical simulations (including quark loops) known as string breaking. It is usually expressed in terms of the shielding of the potential by the quark loops. Once the quarks become separated far enough, then the binding energy becomes sufficiently large that it is advantageous to pull a quark anti-quark pair out of the vacuum, and have two shorter 'strings' plus an extra pair of quarks rather than one really long 'string'. Obviously, our numerical results are in a pure gauge theory, and we do not expect to observe string breaking. If our model is correct, then there must be something suppressing the occurrence of degenerate eigenvalues of the Wilson Loop in pure gauge theory, with this suppression weakened in the dynamical theory. However, this statement remains just a hypothesis, and a dedicated study, well beyond the scope of this work, is required to either confirm or deny it.
We expect that the same model of confinement will be valid in both quenched QCD (Pure Yang-Mills theory) and full QCD. This means that if confinement is explained by some particular mechanism just involving topological objects in the gauge field, then there should be some explanation in full QCD in terms of those objects or the theory surrounding them why the string is broken at large distances. But given that the explanation itself does not depend on the presence of dynamical quarks -only the weighting given to various gauge field configurations -then there should be some allowance for string breaking within the model, and this will apply for both quenched and full QCD. String breaking would then be activated by certain gauge fields which are suppressed in quenched QCD but present in full QCD. That our model has an allowance for string breaking even though we only consider the pure Yang Mills theory is thus not surprising.
We wish to evaluate the change of the individual gauge links of the Wilson line as we gradually increase the spatial extent of the Wilson Loop. Again, we consider a T × R rectangular Wilson Loop, of length L = 2(R + T ) and T > R. We parametrise the position along the Loop using 0 ≤ s < L and the spatial extent of the loop with r = R/2. We again neglect the effects of the corners of the Loop.
The building block of the Wilson Loop is writing θ = e iaφ (neglecting the U(1) phase), we can use equation (A.30) to writê where we use the notation a and φ were chosen so that ∆ r 1s is Hermitian, traceless and diagonal. We now consider the evolution ofû with r, using the machinery developed in section 3 and equation (37). We write θ r+δr s = θ r s e iβ r,s Γ r,s , with We want to calculate e iδσû r+δr,j s λ j = e −iβ r,s Γ r,s e −iaφ e iδσ(A µ(s) (s+ δs 2 ,r)+δr∂rA µ(s)) (s+ δs 2 ,r) e iaφ e iβ r,s+δσ Γ r,s+δσ Using equation (A.30), we can again write this in the form We can evaluate the last term in equation (96) using equation (A.33), which gives Given that ∆ 1 ∝ λ k (i.e. one of the diagonal Gell-Mann operators; λ 3 in SU (2)), we can see that trΓ∆ 1 = 0 and tr(λ j [Γ, ∆ 1 ]) = 0. This means that from these two terms only the 1 2 tr λ j cos(2β)∆ 1 part contributes to u r+δr,j s . In most situations, β is infinitesimal, and this just leads toû r+δr,j s = u r,j s + . . . where the . . . comes from the remaining terms in equation (96) and is of O(δr).û r+δr s isû r s plus a small amount, and therefore its expected absolute value increases (as we have seen in the previous sections) with r once it is clear of the Wilson Loop eigenvalue degeneracy at r = 0. However, where there is a degeneracy in the eigenvalues of the Wilson Loop, β will be large, and at some point close to the degeneracy cos 2β will cross zero. While the first terms in equation (96) will no longer all be infinitesimal,û r+δr,j s will lose its direct dependence onû r,j s . Effectively, the value of the Wilson Loop is reset and loses all knowledge of what occurred at R < 2r.
Let us suppose that for a particular configuration, these eigenvalue degeneracies in the Wilson Loop occur at distances R = 0, R 0 , R 1 , . . .. In this case, once we are clear of R = 0 we may expect u (on average) to rise linearly with R until we reach R 0 , where it will be reset to some valueũ. Subsequently, it will rise linearly again with R − R 0 until R = R 1 and it is reset again before rising with R − R 1 . The location of R 0 , R 1 etc. will vary from configuration to configuration, which means that when we average over configurations, we can consider some mean separation between those Wilson Loops with degenerate eigenvalues R 0 . The pattern we may expect is, apart from some non-linear behaviour at small R as the loop escapes the degeneracy at R = 0 the potential will increase linearly while R ≪ R 0 . For R ≫ R 0 we may expect the potential to be flat, since on each configuration it will be somewhere between ũ and ũ + ρ R 0 , and on the average over configurations this will become a constant as long as the distribution of R 0 , R 1 , . . . is not so sharply peaked that they occur in the same place on each configuration. There will obviously be some intermediate region between these two regimes where we transit from the linear to constant behaviour.
This, of course, precisely as we expect in dynamical QCD: at large distances, the quark loops screen the static quark potential leading to a breaking of the string. This has been observed numerous times in lattice QCD, for example in [20], although in our lattice study which excludes quark loops we do not see it.

Numerical results
We generated 16 3 32 and 20 3 40 quenched (i.e. pure Yang Mills) lattice QCD (SU(3)) configurations with a Tadpole Improved Luscher-Weisz gauge action [28,29,30,31,32] using a Hybrid Monte Carlo routine [33] (see table 3). The lattice spacing was measured using the string tension ρ ∼ (420MeV) 2 . We fixed to the Landau gauge and applied ten steps of improved stout smearing [34,35] with parameters ρ = 0.015 and ǫ = 0. θ andÛ were extracted from the gauge field numerically by solving equations (7), (9) and (11), as outlined in Appendix D. We constructed θ from planar Wilson Loops in the xt plane.F xt ≡Ê x is therefore the component of the field strength tensor responsible for confinement. All our error analysis was performed using the bootstrap method (except for the parameters obtained from fits which estimated the errors from the χ 2 distribution).
The main difficulty we have faced while investigating this proposed mechanism of confinement is the non-gauge invariance of θ. We need to work in a gauge where the gauge field is continuous, while on the lattice there is no such gauge. However, to completely establish the proposed mechanism of confinement, we would need to calculate the variables a and c used to parametrise θ. Extracting the components a and c from θ is straightforward, and we presented some results for the winding of c around the peaks in [16], with some encouraging results. They are, however, gauge dependent, and there is therefore no unambiguous definition of these quantities. While the winding number is protected from continuous gauge transformations in the continuum, on the lattice it is difficult Lattice size L (fm) β a (fm) # 16 3 Table 3: Parameters for our simulations: the lattice size, measured in units of the lattice spacing, the spatial extent of the lattice, L, in physical units, the inverse gauge coupling β, the lattice spacing a, and the number of configurations in each ensemble #.
even to ensure that you are in the right gauge. Therefore for this paper, we have concentrated on gauge invariant quantities such as the static potential and the restricted field strength.
Our aim is to show that the string tension is reproduced by that of the restricted field, and the restricted field dominated by the contribution from θ, or the merinthon term in the field strength (constructed from the potential trλ j θ † ∂ µ θ). We should show that the restricted field strengthF xt is dominated by points (which would appear on the lattice as peaks a single lattice spacing across), while the other components ofF contain structures extended in the directions specified in table 1.
Measuring the string tension associated with the restricted potential is straight-forward, since it is gauge invariant, so it can be extracted from the Wilson Loop using standard methods. Equally, the restricted field strength can be measured using the plaquette definition and as this is gauge invariant it can be measured without any ambiguities or difficulties. Extractinĝ F j µν from e iF j λ j , however, can only be done modulo 2π for F 3 and 4π/ √ 3 for F 8 . For smooth gauge fields, whereF is small (as is the case in most lattice applications), this is not a concern. We, however, as interested in those cases whereF is not small, and we shall find that in some applications |F j | may be larger than π. For example, ifF 3 ∼ 2π it is impossible to distinguish this large value from a fluctuation around zero. This is not of significance for the physics, since these shifts inF will not affect the value of the Wilson Loop, but it may be of importance in some of the following as we attempt to identify peaks in the field strength.
However, constructing the merinthon contribution to the restricted field strength is more challenging because there is no clean gauge-invariant way to observe this on the lattice. In the continuum, the observable is lim a→0 2 aμ ∂ µ θ x+ 1 2 aμ ) for lattice spacing a. Since this is not gauge invariant, a gauge transformation which would be discontinuous in the continuum would lead to additional discontinuities appearing in the observable, or the removal of the discontinuities we want to observe. Fixing the gauge would not help us, because we could well be fixing to a discontinuous gauge. We solved this problem by introducing a gauge invariant quantity which measures the merinthon contribution to the potential.
We introduce a fieldŨ , which is the gauge field U subjected to a large number of coarse stout smearing sweeps. A small amount of smearing is useful to remove discontinuities on the order of the lattice spacing; a large amount is usually seen as dangerous as it removes the physical features of the gauge field; however we can use this to our advantage. Each Stout smearing sweep replaces U x,µ → U ′ x,µ = e iQx U x,µ where Q is a Hermitian operator constructed from closed loops of gauge links starting and finishing at x containing the gauge link U x,µ (the original smearing algorithm [34] used plaquettes; we also used 2 × 1 rectangles following [35]). After a gauge transformation U x,µ → Λ x U x,µ Λ x+μ , the smeared link transforms in the same way, U ′ x,µ → Λ x U ′ x,µ Λ x+μ . The second effect of stout smearing is to smooth the gauge field, and reduce the value of the plaquette to zero and thus the potential to a constant, and in practice the potential becomes zero. Eventually, after repeating the smearing a very large number of times, the over-smeared linkŨ becomes, in effect, a gauge transformation of something close to the identity operator. 6 Thus θ † xŨ x,µ θ x+μ is both gauge invariant and in a gauge where U is continuous would measure θ † ∂ µ θ, the observable we need. We therefore use this observable to measure the merinthon portion of the field strength. In effect, we are replicating the calculation for the restricted field, but with the gauge field replaced by something close to zero, just leaving the merinthon contribution to the restricted field. To extract the Abelian component of θ † xŨ x,µ θ x+μ we perform another Abelian decomposition and extract the restricted fieldÛ x,µ which satisfies [Û x,µ , λ j ] = 0 and tr(λ j (X −X † )) = 0 with We can then compare the field strength from this restricted field, which represents the merinthon portion of the restricted field strength, with that of the restricted field U . Our expectation is that the observables calculated from the merinthon field and the restricted field should be similar: the string tensions should be in agreement, and the field strength should contain a similar pattern of peaks. We denote H j µν as the field strength tensor constructed usinĝ U .

String tension -Fixed Wilson Loop
Firstly, in figure 2 and table 4, we extract the string tension, ρ, for the original gauge field U , the restricted fieldÛ and the merinthon component of the restricted field M . We have calculated this in two ways: firstly by calculating the expectation value of the R × T Wilson Loop for one of the fields, and fitting its logarithm firstly to V T + a + b/T and then V to ρR + c + d/R for additional fit coefficients a, b, c and d; and secondly by fitting the whole data set according to ρRT + aR + bT + c + d/R + e/T + f /(RT ) + gR/T + hT /R. The results from these two approaches were in good agreement. We tabulate the results from the combined RT fit, while plot the results of the single fits. The errors on the expectation value were calculated using the bootstrap method; we calculated the errors on the fitted parameters by measuring the surfaces of constant χ 2 . In principle we should recalculate a new θ field for each Wilson Loop; however the numerical cost to do so was prohibitive for an initial study, and we firstly show results obtained by using the same θ field calculated from one set of Wilson Loops in the xt plane for our entire simulation [18]. We discuss this simplification a little more in the next subsection. We average over all Wilson Loops for the restricted field, and not just those with the correct θ field. Were we to recalculate a different θ for each Wilson Loop, the string tension extracted from the U andÛ fields would be identical. Because we have not done this, we expect the two quantities to differ by a small amount.
We calculateŨ after 100, 300, 500, 600, 800 and 1000 sweeps of stout smearing with parameters ǫ = 0, ρ = 0.015 following the notation of [35], and extract the merinthon portion of string tension by constructing Wilson Loops from the Abelian decomposition of θ †Ũ θ. We also show the string tension forŨ , and can confirm that it is much smaller than that of the original gauge field and decreases as we increase the level of smearing. The observable we are using to measure the merinthon portion of the string tension is unaffected by the smearing, suggesting that we have indeed measured the contribution from θ∂ µ θ † rather than any remnant of the gauge field remaining after the smearing. We see a good agreement between the merinthon and restricted field string tensions, within 10% on the smaller lattices (although the difference is larger on the larger lattice where we have less statistics), suggesting that it is indeed structures within the θ field that lead to confinement. There is no obvious change in this pattern as we change the lattice spacing.

String Tension -Varied Wilson Loop.
We generated the results presented in the previous section as a quick and dirty initial study, largely to demonstrate to ourselves and the wider community that this method was worth pursuing, and we first published these results in [18] (these are the 10 smearing sweeps results of table 5).  Since the initial drafts of this paper, we have commenced a more accurate study removing the simplification of the previous section, this time recomputing the Abelian decomposition for each Wilson Loop under consideration. This obviously maintains the identity between the confining potential of the full and restricted gauge fields. Our results so far are only tentative, and far from being finalised. We present early results here mainly to update our initial presentation of this study in [37]. In those proceedings, we found that although the restricted and full QCD potentials were in perfect agreement, the topological contribution to the string tension seemed to be relatively small. If this result is genuine, there are two possible explanations: firstly, that the Maxwell term in the restricted potential, contrary to expectations, also has a large contribution to confinement. The second possibility is that our means of extracting the topological part of the string tension is flawed, and what we were measuring was not the topological part of the string tension but something else. Since our initial publication, however, we have expanded our study, and have noted another issue. We apply a small number of smearing sweeps to the gauge field initially, before applying the Abelian decomposition, to remove various unphysical fluctuations of the order of the inverse lattice spacing. This is necessary to allow us to obtain a clear result. How many smearing sweeps should be used is difficult to assess, and for a full presentation of the results we should measure the string tension at different levels of smearing to gain some idea of the systematic error caused by this arbitrary choice. As we did so, we observed that the measured topological contribution to the smearing depended strongly on the number of initial smearing sweeps. Furthermore (unlike in the results of the previous section), we also found that at a larger number of initial smearing sweeps, there was a strong lattice spacing dependence on the ratio between the restricted and topological potentials: as the lattice spacing grew smaller, the restricted potential became increasingly dominated by the topological part. This means that our full analysis is considerably more difficult to do accurately than we had first envisaged, and, given our limited computing resources, will take longer to complete.
Some preliminary results (which are not finalised, and could change when we make our final publication) are shown table 5. As we vary the number of smearing sweeps, the string tension of the original gauge field first increases as the unphysical dislocations are removed, and then slowly decreases as the smearing gradually distorts the gauge field. This obviously begs the question of which value we should use, and our intention has been to use the number of smearing sweeps where the string tension is at its maximum (at the present time, we have kept the smearing parameters fixed at ǫ = 0, ρ = 0.015 for this study). So far, we have only conducted this analysis for the 16 3 × 32 configurations. Table 5 lists the β value of the configuration, the number of smearing sweeps where the maximum value of the string tension was reached (out of those which we have sampled so far), and the measured string tension for the QCD (and restricted) gauge fields, and the proportion of that string tension which we have measured as coming from the topological part of the restricted potential (after 2500 smearing sweeps).
It can be seen that firstly the maximum string tensions for the β8.3 and β8.52 configurations are higher than the results given above (they were calculated following a different amount of initial smearing). Secondly, the proportion of the string tension attributable to the topological part increases as we decrease the lattice spacing. This pattern is not merely a matter of the results being taken after a different number of smearing steps, but is apparent at each value of the initial smearing, although it is more pronounced as we increase the smearing. Our analysis is less developed on the larger lattice, but the results we have suggest that there is no significant volume dependence.
This strong dependence on the lattice spacing is awkward because it necessitates an extrapolation to the continuum limit. This is challenging with the data we have gathered so far, firstly because we do not know of a good formula to guide the extrapolation, and secondly because we are still on a fairly coarse lattice. We need to finalise our results on these small lattices, and then expand our data set so that we have access to finer lattice spacings before it is useful to attempt to perform such an extrapolation; we will have to leave the completion of this study to a future publication.
However, it is not unreasonable to suggest on the basis of this data that the string tension will   be dominated by the topological part of the restricted potential (as we have measured it) in the continuum limit. This still needs to be confirmed by a complete and more careful analysis, which we are currently undertaking.

Presence of peaks inF xt
The next question is whether theÊ x component of the gauge invariant restricted field strength is dominated by point-like objects as expected in the theoretical model, and we plot the distribution ofF 3 xt in figure 3 on two neighbouring slices of the lattice, using a contour plot to display the data. The corresponding plot forF 8 xt has a similar set of structures. It can be seen thatF 3 xt is indeed dominated by these objects a single lattice spacing across (or, on occasion, two lattice spacingssee also figure 9). There is no relationship between the peaks on neighbouring slices of the lattice, suggesting that these are indeed points rather than strings or surfaces. Are these peaks from the merinthon part of the restricted potential? The background shading of figure 3 gives a similar plot for the field strength extracted fromÛ , and the is a strong correspondence between the location of the peaks in the merinthon field and the restricted field, although on occasion one of the peaks may be shifted by a lattice spacing (which is the limit of the resolution of the lattice, so we can expect errors of up to a half lattice spacing on each of the data sets), and a few objects in the merinthon field are not visible in the restricted field. Nonetheless, the similarity between the two sets of data indicates that these peaks are caused by the magnetic portion of the ensemble as predicted by the theory. Although figure 3 shows data from just one slice of one configuration for one ensemble, the pattern of peaks is similar across all our ensembles.
In figure 4 and table 6, we investigate whether these peaks are responsible for the string tension, i.e. if excluding the peaks would reduce or eliminate the confining potential. Usually, we measure the expectation value of the Wilson Loop by averaging over every planar loop in the configuration. However, it is a straight-forward exercise to restrict the average to loops in the xt plane and either only include or exclude those Wilson Loops which contain one of the peaks from the average. We show in figure 4 the result of excluding the Wilson Loops which contain peaks. We introduce a cut-off C, and do not include any Wilson Loops where |F | > C for one of the lattice sites contained within the loop. If the peaks are responsible for confinement, we would expect to see the string tension disappear as we exclude more of the peaks. If they are merely incidental, and confinement is dominated by some other mechanism than that proposed here, we might expect to see little difference in the string tension as we change the cut-off. The picture is, however, consistent with the idea that these peaks are responsible for confinement.   Table 6: TheÛ string tension, ρ , excluding Wilson loops containing peaks of height |Fxt| > C from the average. The β8.3L data refers to the largest lattice; there was insufficient data to get a reliable estimate of the string tension at the lowest cut-off on this lattice.

Electromagnetic field tensor
We now look at the full strength tensor to examine how closely it resembles our expectations from section 4.3, although for conciseness we concentrate only onĤ 8 µν (the results forĤ 3 µν and F j µν are indistinguishable). First of all, we plot histograms of each component of the field strength in figure 5, showing what proportion of lattice sites had a particular value of the electric or magnetic field. It can be seen that the distribution forÊ y ,B y ,Ê z andB z are all similar, while the distribution forÊ x andB x differ. Most of the plots show similar features: a large peak around H 8 µν = 0, meaning that most of the points have no field strength beyond a small fluctuation around zero, and a plateau for larger field strength up to the maximum value ofĤ 8 µν , meaning that on a minority of lattice sites there is a large contribution to the field strength, where the number of sites with that particular value of the field strength seems to be independent of the field strength. This is as we expected. We also see that the proportion of the points with small values of the field strength increases as we decrease the lattice spacing, which suggests that these objects have a fixed density per physical volume but the size of the objects decreases as the lattice spacing tends to zero, which is consistent with the idea that these are thin points or string-like objects. B x andÊ x are different; there are twice as many lattice sites with a large value ofB x compared to the fields in the y and z direction, and very few lattice sites with a large value ofÊ x . There also seems to be little volume dependence.
The picture for all the fields is consistent with having some (perhaps Gaussian distributed) fluctuations around zero, and then some additional objects at large field strength. ForÊ x , the number of the objects diminishes as the field strength increases, while for the remaining fields the proportion of lattice sites with a particular field strength seems to be roughly independent of the field strength, at least for those field strengths we can measure. This means that there may be peaks in the field with a field strength ∼ 2π forF 3 or ∼ 4π √ 3 forF 8 which are impossible for us to distinguish from zero. It should therefore be borne in mind that not every object in these fields will be visible to us.
This picture can be studied in more detail by plotting slices of the electromagnetic fields: figure  6 shows this in the xt plane; figure 7 in the yt plane, 8 in the zt plane.
These objects are not just isolated magnetic monopoles: there are peaks in the electric field, and the magnetic fields are not spherically symmetric. They may, of course, be some magnetic monopoles and some other structures, or possibly a signature of a monopole/anti-monopole condensate. The results are consistent with our expectations given in table 1, and the patterns seen on this configuration and slice of the lattice are duplicated across all our ensembles. The electric field is clearest (although not completely unambiguous):Ê x is a point, whileÊ y andÊ z form strings parallel to the x-axis. The magnetic field is a little less clear, and expected string like behaviour ofB y andB z is less obvious, although we have seen clear examples of t-strings.B x is, as expected, extended along both the x and t axes, though not the z or y axis. In particular, other    orientations of the fields with strings in directions perpendicular to the xt plane are inconsistent with this data.

Cluster analysis
To investigate the restricted electromagnetic fields more thoroughly, we employ a simple strategy to determine the size and shape of these structures. We first of all find clusters, connected regions of sign coherent high electric or magnetic field strength. We perform the analysis for each component of the field strength tensor separately. While we have performed this analysis forF 3 µν , F 8 µν ,Ĥ 3 µν andĤ 8 µν individually, there is no significant difference in the results, so here we just show the results forĤ 8 µν . We define a cluster as all connected sites with a value of |F µν | > 1 and identical sign ofF µν (i.e. except for clusters containing just one lattice site, each site in the cluster is the nearest neighbour of at least one more site in the cluster and all its nearest neighbours which satisfy the bound |F µν | > 1 are within the same cluster and each cluster is internally connected so it is possible to reach any site within the cluster from any other by crossing to a neighbouring site while remaining within the same cluster, while different clusters are not internally connected with each other). We have also studied (but do not show here) results with numerous values of this cut-off from 0.8 to 2.0, but we do not see any important difference as a result of varying this cut-off.
For all the fields exceptÊ x , the vast majority of the lattice sites where the restricted field satisfied our condition were contained within a large cluster which spanned the entire lattice. This is in certain respects similar to the structures found in the full field strength tensor [38,39]. The theoretical picture we were testing described a number of strings in either the x or t direction. This large structure may be formed if the density of these strings is sufficiently large that every string has at least one neighbour on an adjacent lattice site in the y or z directions. It may, however, relate to a different picture of the vacuum. There were a few clusters not attached to this large structure; and most of our analysis is dominated by results from these smaller structures, as our averaging with respect to the number of structures weights against the smaller number of structures with a large number of lattice sites. We also present a few results which specifically target these larger structures. Here we see the expected string behaviour, with a large number of objects at just a single point. We emphasise thatÊ x (orF xt ) differed, in that the field strength did not contain this large structure, but only single lattice site structures. To attempt to analyse the shape of this large cluster, we have also performed a similar analysis, but where the clusters are constructed from connected sites with a two dimensional slice of the lattice.
When we show our results, we will label those plots where we have included all clusters as unfiltered, those which only average over clusters with at least four lattice sites as filtered, and those which average over clusers with at least two hundred lattice sites as Large Clusters.
The theoretical expectation was that ifÊ x is dominated only by points of large field, then the remaining fields should be dominated by objects which are either points or strings in the x or t direction. The only expected difference in the results for theÊ y ,Ê z ,B z andB y fields is the orientation of the string. A point will show itself as a cluster containing just a few lattice sites (ideally only one lattice site, but due to lattice artefacts it may be smeared across two or three sites). A string would display itself as a line of flux, i.e. each site within a cluster would have only two nearest neighbours within a cluster if it is in the centre of the string and one if it is at the end of the string, and it will only be extended in one direction. In practice, the situation will be messier than this. Neighbouring strings may be counted as the same cluster, and there may be ambiguities caused by the poor resolution of the lattice leading to a thickness of the string being more than one lattice spacing in places, two different strings may overlap (and the expectation is only accurate if the functions G and R ′ have no strong angular dependence). Nonetheless, we would hope to see some signal of these strings when considering an average size and shape of the cluster.
We first plot the distribution of the number of lattice sites in each cluster in figure 9. The picture is similar across all our ensembles, so we have only plotted the data for the β = 8.52 ensemble. It can be seen that all the components of the field strength contain a large number of  points with only one or two sites contributing to the cluster. However, theÊ x field only contains these objects, while the remaining fields also contain a substantial number of structures extended across several lattice sites. This means that if the picture outlined in the earlier sections is correct, then theB x ,Ê y ,Ê z ,B y andB z fields can only contain points and the strings described in table 1. We also note that the data for theÊ y ,Ê z ,B y andB z fields are indistinguishable, in line with our expectations.
We next study the number of nearest neighbours within a cluster in figures 10, 11 and 12. The first of these plots shows the data across all the clusters, while the second filters only those clusters containing four or more lattice sites, so eliminating the points. We have studied the effects of changing this filter from all clusters with at least three sites to all clusters with at least six lattice sites, and there are no substantial changes in the results. The third plot focusses on the very large clusters. A point should have no nearest neighbours (or possibly on occasion one), a string one or two neighbours, a two dimensional surface one, two three or four neighbours (with the proportion depending on the size and shape of the surface), and higher dimensional objects would have a larger number of nearest neighbours within the cluster. The data forÊ x is again unambiguous: these are points. For the remaining fields, it is less clear. Focussing on figure 11, we see that firstly there is no substantial difference between the results for the different fields: the same structures appear in the fields with the same frequency. The majority of the lattice sites within the larger clusters had two or three neighbours, with a small number having four. We can rule out that these are extended in more than two dimensions. This may indicate strings where two strings are on neighbouring lattice sites and thus counted in the same cluster (and giving objects with three or more nearest neighbours), or it may indicate that there are also some two dimensional surfaces with a large electromagnetic field. There is a small but statistically significant increase in the proportion of sites with two nearest neighbours as we decrease the lattice spacing, and no noticeable difference in the distribution as we increase the volume of the lattice. For the very large clusters, which are extended across the whole lattice, we again see that the sites have two or three nearest neighbours in the cluster. This is consistent with, but does not demonstrate, that these clusters are constructed from densely packed one dimensional objects.
We predicted that the strings in theB x field could be in either the x or t direction, while the electric fieldsÊ y andÊ z would extend in the x direction only and the magnetic fieldsB y andB z in the time direction only. In figure 13, we plot the distribution of the spatial extent in each direction of each cluster (for example, a cluster that stretched from (t, x, y, z) = (3, 4, 5, 7) to (6,9,8,13)             directions. This is not fully realised in our data: there is a noticeable number of structures which are extended for up to four or five lattice sites in the 'unwanted' directions. However, for each field, the clusters extend considerably more in the directions of the expected strings. We once again restrict the results we present to those clusters containing at least four lattice sites, although once again we have studied numerous values of this filter without seeing any significant difference in the results. The shape of the curve changes as we only include larger clusters (larger extents become more likely as we increase the cutoff for the cluster size), although the maximum in the X and T directions remains at four lattice spacings. This may be consistent with our expectations, with strings overlapping each other, or it may indicate that there are additionally other types of objects in the electromagnetic field.
One question that can be raised is whether that some clusters have a large thickness in the unwanted directions and some sites within the clusters challenge whether these large fields are arranged in strings one lattice spacing thick. Clearly, in the context of this picture, this can be explained by having strings on neighbouring lattice sites, and some degree of the breakdown of the naive behaviour should be expected. To better judge whether this is indeed the case, we have generated configurations where the field strength is distributed as expected. We have counted the number of lattice sites where F j µν > 1 and F j µν < −1, and generated two different sets of configurations of the restricted field strength. In each case, we ensured that the same number of lattice sites had the large F j as on the actual quenched QCD configuration. In the first random configuration, these sites were allocated completely random positions around the lattice. In the second random configuration, one quarter of the sites were positioned randomly, and the rest were arranged into strings in either the x or t direction as appropriate for the field. The principle ambiguity is what to use as the distribution for the length of these strings, since the distribution given in figure 13 will contain examples from multiple strings. We examined several Poisson and Gaussian distributions, and the data we present here is taken from a Poisson distribution with the mean value extracted from the data presented in figure 13. Only the analogue of the curve in 13 in the direction of the string depends strongly on the distribution of the length of the string. We cannot predict the shape of this curve for this simple model. We then simply repeat the cluster analysis on these configurations.
Some sample results are plotted in figures 14 and 15. It can be seen that the data for figure 14 is qualitatively similar to the actual QCD data for E x , while there are disagreements when comparing the model and actual data sets for the other fields. The second model, with strings, compares well with the numerical data forB y ,B z ,Ê y andÊ z ; in particular the plots for the number of nearest neighbours and the extent of cluster in the directions away from the string are close to the QCD data. The extent of the cluster in the direction of the string strongly depends on the distribution for the string length, so we do not expect a perfect agreement here. In particular, that the 'thickness' of the string agrees with the data, and seems to be largely model independent, suggests that this cannot be used to rule out the predicted string model of flux. However, the broad agreement between our model calculation and the data implies that the QCDÊ x field is dominated by points and the other fields by a mixture of points and strings.
Finally, we check to see whether the locations of the clusters in the different fields are correlated. If our model is correct, then the space-time location of many of the peaks in E x should be correlated with high field strengths in the other fields. To do this, we take each of the peaks inÊ x (i.e. each point where |Ĥ xt | > 1) and plot a histogram of the value of the other fields around that site. For comparison, we produce a similar plot measuring the field strength around random lattice sites rather than the peaks of theÊ x field. We would expect the plot for the actual data to have a higher concentration of large field strengths than the random data, and a smaller concentration of lower field strengths. There is one complication: we measure the field strength at different locations in space-time for each field. For example, we measureF xt (t + 1 2 , x + 1 2 , y, z) ≡Ê x (t, x, y, z) and F yz (t, x + 1 2 , y, z + 1 2 ) ≡B y (t, x, y, z). This means that there may be a shift of up to one lattice spacing in the location of the peak within each of the measured fields. We therefore take the maximum value of |B y (t, x, y, z)|, |B y (t + 1, x, y, z)|, |B y (t, x, y, z − 1)| and |B y (t + 1, x, y, z − 1)|     Figure 15: The distribution of the spatial extent of the clusters for where the field strength is selected according to the distribution given in figure 9 for the β = 8.52 ensemble for theBy (top left) andBx (top right) fields, but with the lattice sites randomly distributed with three quarters of the larger field strengths arranged into strings; the distribution of the number of nearest neighbours (bottom left), and the filtered distribution of the number of nearest neighbours (bottom right) for this random ensemble.
(for example). This creates a bias towards larger field strengths for both the actual locations of the peaks inÊ x and the random data. This is particularly pronounced forB x , where there are 16 lattice sites which we need to consider, and coupled with the higher density of lattice sites with a large value of theB x field this means that actual data fromB x does not differ greatly from the randomised data. We plot the results in figure 16 for the β = 8.52 ensemble. There is a statistically significant difference between the curves for theÊ x peaks and the random data for all the fields, with a larger proportion of sites with a large value of theB y ,B zÊy andÊ z fields and a smaller proportion of sites with a small value of these field surrounding the peaks inÊ x than for the random data. There is, however, still a large number of sites where there does not seem to be a great deal of correlation between the peaks in (for example)Ê x and E y ; partly this can be explained by those sites where in practiceÊ y ∼ ± √ 3π but we measureÊ y = 0, but it still seems as though not every peak inÊ x corresponds to a peak in every one of the other fields. This is not inconsistent with equation (82); for example if ∂ 1 G is small then only three of the other fields are large when there is a contribution inÊ x , and if both ∂ 1 G and ∂ 2 G are small then only two of them do. This data is then qualitatively but not conclusively in line with the expectations: there is a clear correlation between the fields, however it is not large enough to say that all the the topological structures in the different components ofF µν are correlated with each other.
This analysis has concentrated on the smaller, isolated structures. We still have to consider the larger structure which extends throughout the lattice, which we have seen from the nearest neighbour analysis is either one or two dimensional. This is consistent with our expectations if it is constructed from neighbouring strings along the X or T directions (so we may have an x-string at y = z = t = 0 and extending from x = 0 to x = 5, and another string at y = 1, z = t = 0 extending from x = 5 to x = 10, and so on; if these are closely enough packed then our analysis will place them into the same cluster). We investigate this structure by modifying the cluster algorithm. We firstly select two directions; and we construct the cluster by only considering nearest neighbours in those two directions. This restricts us to isolated objects. By varying the two directions, we can see if the strings predominantly extend in particular directions.
We display the extent of the strings in figures 17, 18 and 19. We have shown the xt, xy and yt planes as examples. There was again little difference between our ensembles. Again, we restrict ourselves to clusters with at least four lattice sites, so there must be an extent of on average two or three lattice sites in at least one direction. We again see that the magnetic fieldsB y andB z are extended far more parallel to the t axis than in the other directions,Ê x andÊ z along the x-axis, whileB x is extended along both the t and x axes. None of the fields show similar features parallel to the y and z axis.
The picture we see is consistently that the restrictedÊ x field is dominated by points, while theB x field by strings in both directions in the xt plane, and the other electric fields by strings in the x direction, while the other magnetic fields by strings in the t direction. These are not magnetic monopoles (or not only monopoles): the behaviour of the electric fields andB x rules this out (and note that we presented results for the supposedly 'magnetic' trn[∂ µ n, ∂ ν n] part of the field strength; although the picture for the full field strength is the same). The results for the restricted electric and magnetic fields are consistent with the expectation (though, of course, they do not prove that the model of confinement we presented in the first part of this paper is correct: they may also be consistent with other models of the vacuum).

Conclusions
We have investigated whether the confining quark potential in quenched QCD is caused by CDG gauge-invariant merinthons. Our main novelty is to construct the CDG colour field n j = θλ j θ † from the eigenvectors of the Wilson loop, which allows us to access the full symmetry group of the merinthons and permits a theoretical discussion. From this analysis, we have shown that the Wilson Loop for the actual QCD gauge field is identical to the Wilson Loop constructed from the restricted Abelian field, and this allows us to express the Wilson Loop in terms of an area integral over its surface. We then find that the construction of θ allows for certain topological defects to    appear in the restricted field strength. In the component of the restricted field strength in the plane of the Wilson Loop, these manifest themselves as point-like objects; in the other restricted electric and magnetic fields they manifest themselves as strings parallel to one of the axes of the Wilson Loop. The net result is that the restricted field strength is dominated by a large structure which is constructed from these one dimensional strings. The full QCD field strength contains the restricted field, so it may be expected that these global one dimensional structure may extend to the full QCD field strength (albeit with O(4) symmetry suggesting that F µν [X] will also contain strings in other directions). If these strings are sufficiently dense, this picture is consistent with lattice simulations [38,39]. In particular, these objects are not magnetic monopoles; at least for this particular choice of the colour field. The peaks in theÊ x field lead to an area law for the Wilson Loop and thus confinement. We have also given an argument within the context of this decomposition which might explain string breaking at large spatial distances.
Our numerical results focus on searching for these strings and peaks in the restricted field strength, and we have uncovered firm evidence that the electric field in the direction of the Wilson Loop, which is responsible for the confining potential, is dominated by objects which are just present on a single lattice spacing. The remaining components of the electric and magnetic fields contain one dimensional objects which extend in the directions suggested by our theoretical analysis. While this does not prove that our analysis is correct, they are certainly consistent with it; and whatever these peaks are, they are responsible for the confining potential. Our numerical analysis is consistent with the vacuum related to the restricted field being dominated by a combination of magnetic monopoles and the π 1 topological objects we have labelled merinthons. It is also possible that the binding of monopoles and anti-monopoles into pairs (similar to the Cooper pairs in superconductivity) might lead to a monopole condensate with a different pattern in the field strength than one would expect from a naive picture of isolated monopoles and anti-monopoles, and we cannot rule out this possibility. We also show that in an initial test, making an important approximation, the string tension can be accounted for by the second, 'magnetic', term within the restricted field strength. The results of the full calculation (without any approximation) are less clear: all we have shown so far is that we will need to make a careful study of the continuum limit. The objects responsible for confinement are embedded in the colour field, and only indirectly in the gauge field.
We believe that this evidence is enough to merit further, more detailed, investigation of whether confinement can be explained within the context of this particular formulation of the CDG Abelian decomposition.
We have not yet investigated the effects of changing the representation of the group, of incorporating quark loops into our lattice study, or studied the de-confinement transition, and hope to investigate at least some of these questions in future work. which means that which, unless W [C s ] has degenerate eigenvalues, is only possible if θ † s U s θ s+δσ is diagonal, as required by the lemma.
There always exists a θ field (unique baring the U (1) NC phase factors and ordering of the eigenvectors) which diagonalises a W [C s ] ∈ SU (N C ) as long as W does not have degenerate eigenvalues. From lemmas Appendix A.1 and Appendix A.3, we see that this field is also a unique (up to the same caveats) solution for a field that diagonalises U . Thus we have shown that there is a SU (N C )/U (1) NC −1 field θ which diagonalises the gauge links along the closed Wilson Line.

Appendix A.2. The Abelian Decomposition in the continuum
Here we construct an explicit form forÂ. The argument here is based on that presented in [11].
Lemma Appendix A.4. For any Field Y in the adjoint representation of a SU(N C ) lie algebra where λ j are the diagonal elements of the lie algebra normalised so that tr(λ j λ k ) = 2δ jk .
Proof Y can be decomposed as Y = λ a Y a , where λ a /2 are the generators of the gauge group. λ a are Hermitian, traceless matrices which satisfy trλ a λ b = 2δ ab . We adopt the convention that the diagonal Gell-Mann matrices are referred to with the index j, k, . . ., while λ a , λ b , . . .  diag(1, 1, 1, 1, . . . , −(N C − 1)). The other elements can be expressed in terms of the basis E ±α , where E +|α| consists of one element of 1 above the diagonal and its symmetric counterpart below it, while E −|α| has −i above the diagonal and i below it. For example, Now, as λ j is diagonal, [λ j , E α ] must have the same non-zero components as E α , and thus be a linear combination of E α and E −α . Given that trE α [λ j , E α ] = 0, we can write [λ j , E α ] = a j α E −α for some coefficient a which needs to be determined (a could be 0). Similarly, [E α , E −α ] is diagonal and traceless, so [E α , E −α ] = k b k α λ k , for some coefficients b. The following results then proceed immediately using the trace theorem tr(A[B, C]) = tr(C[A, B]): The first equality of the last equation follows from the specific form of E α and E −α : [E α , E −α ] has two non-zero terms on the diagonal, one 2i and the other −2i. We now have the necessary ingredients to prove equation (A.7). Expanding Y = c j λ j + d α E α , the right hand side reads Remark In the Lemma above, we used a natural representation of λ j and E α in terms of the Gell-Mann matrices. However, it is easy to see that the same commutation relations between the generators apply if we rotate the basis λ j → θλ j θ † , E α → θE α θ † , with θ ∈ SU (N C ), while the traces of powers of the operators are also untouched. As the proof of (A.7) just depends on the commutation relations and the traces, it is clear that a form of this equation is also valid in the rotated basis. In particular, this means that we can replace λ j in equation (A.7) with the rotated basis n j , leaving us with From equation (13), The second field X µ may be constructed from X µ = A µ −Â µ .
Theorem Appendix A.5. The definitions of the restricted gauge field given in equations (16) and (22), i.e.

Lemma Appendix
Lemma Appendix A.9. For non-commuting objects X and Y in some Lie algebra, where X is differentiable with respect to a variable s and we can neglect terms of O(δs 2 ) e X(s) e δsY e −X(s+δs) = e δsY −δs∂sX+δs which is the result we wanted to prove.

Appendix D. Numerical Methods
Appendix D.1. Algorithm for fixing θ It is simplest to find θ by solving for the eigenvectors of the Wilson Loop; however we employed a different technique (the method discussed below has the advantage that it can be easily adapted to other rules for choosing θ discussed in the literature than just the one we use in this paper and to other problems such as gauge fixing and solution of the defining equations. It is also reasonably quick and efficient.).
We wish to find the θ that solves and λ j = λ 3 or λ 8 and θ ∈ SU (3) (in fact, the solution we want is at F [θ] = 0). The trace is over both lattice sites and SU(3) indices and V is the lattice volume. We solve this minimisation problem by combining two different algorithms; one based on molecular dynamics, which is an extension of the steepest descent method, and the second a Newton-Raphson algorithm. The goal is to have the molecular dynamics solve the system to a sufficient accuracy that the Newton-Raphson algorithm will converge, and then Newton-Raphson rapidly finds the solution to a high precision. The combination of these two methods, while it may not be the most efficient algorithm, worked well on all the lattice volumes we used. We also employed similar algorithms to gauge fix and solve the defining equations of the Abelian decomposition.
Having found θ, one must subsequently order the columns of θ and fix the U(1) phases according to the chosen fixing condition. We ordered the columns according to decreasing real part of eigenvalue, and either fixed θ using d = 0 or by having the same value for eachû around the Wilson Loop, depending on the observable.
for each direction µ. EachX µ can be solved independently, although in practice we set up our