The quark-mass dependence of the potential energy between static colour sources in the QCD vacuum with light and strange quarks

The low-lying energy spectrum of the static-colour-source-anti-source system in a vacuum containing light and strange quarks is computed using lattice QCD for a range of different light quark masses. The resulting levels are described using a simple model Hamiltonian and the parameters in this model are extrapolated to the physical light-quark masses. In this framework, the QCD string tension is found to be $\sqrt{\sigma}=445(3)_{\rm stat}(6)_{\rm sys}$ MeV.


Introduction
The energy levels of a static quark anti-quark pair V n (r) as a function of the interquark distance r probe Quantum Chromodynamics (QCD) at all scales.For very small distances r ≲ 0.1 fm the potential can be computed by perturbation theory.At large distances the behaviour of the potential depends on the matter content.In puregauge theory at large separations, the potential can be modelled by an effective bosonic string.With dynamical quarks, the ground-state potential V 0 (r) flattens at large r due to the formation of a pair of static-light mesons ("light" refers here to the dynamical quarks).It was noted already in [1] that "the ground state potential therefore be called a static quark potential or a static meson potential".A first estimate of the distance where saturation of the ground state potential sets in was provided in [2].Later, a model for string breaking as a mixing phenomenon between a string and a two-meson state was formulated in Ref. [3], which predicted for lighter quark mass the energy gap becomes larger and the width of the mixing region also expands.First observations of string breaking were made in the three-dimensional [4] and fourdimensional [5,6] SU (2) gauge-Higgs model with scalar matter fields.Evidence for string breaking in QCD with N f = 2 quark flavours was shown in [7] for a sea quark mass slightly below the strange quark mass corresponding to m π ≈ 650 MeV.String breaking occurred for staticsource separations of r c ≈ 1.25 fm.Our previous work [8] presented a computation of the three lowest potential energy levels in QCD with N f = 2 + 1 quark flavours at quark masses corresponding to m π ≈ 280 MeV and m K ≈ 460 MeV.The computation was performed on the N200 gauge ensemble generated by the CLS consortium [9].We found two string breaking distances r c ≈ 1.22 fm and r cs ≈ 1.29 fm corresponding to the saturation of the potential levels by two static-light and two static-strange meson levels, respectively.
This work extends the analysis of Ref. [8] to study the quark-mass dependence by including two additional gauge ensembles with the same lattice spacing a ≈ 0.063 fm and with one ensemble at a lighter quark mass and one heavier ensemble, covering a range of pion masses m π ∈ [200,340] MeV.Ref. [8] showed the energy levels are modelled reliably by a simple three-state Hamiltonian with six parameters.One model parameter can be interpreted as a string tension.In this work we determine the quark-mass dependence of these model parameters.We introduce one more parameter in the model to include a curvature term proportional to 1/r as seen in the Cornell potential [10].We extrapolate values of the seven parameters to the physical point following a trajectory where the sum of the bare quark masses is kept constant and close to its physical value.
The string tension enters as a parameter into models describing the fragmentation of partons into hadrons [11,12,13,14].The so called "Lund model" [15] is implemented to describe hadronisation in the "Pythia" Monte-Carlo event generator [16].In such models the string tension corresponds to a constant energy per unit length of the flux tube containing the colour field of a quark and anti-quark pair.The result of the present work provides input from first principles which could be used to refine these models by adding information on the excited states, for example.
The paper is organised as follows.The methodology of the lattice calculation, closely following the techniques used in our earlier work, is reviewed in Sec. 2. Section 3 introduces the simple model describing the dependence of the lowest energy levels on the static-colour-source separation, including the addition of the 1/r term.Having extracted the model parameters from the data at three values of the light quark mass, Section 4 establishes their dependence on the light quark dynamics and extrapolates to the physical quark-mass value.This is used to present the energy levels in the physical theory as the final result of the paper in Fig. 9. Our conclusions are presented in Section 5.

Methodology
To investigate the dependence of the static energies on quark mass, we extend our previous analysis to include three ensembles of gauge configurations generated by the CLS consortium [9,17].These ensembles include the dynamics of N f = 2 + 1 flavours of non-perturbatively O(a) improved Wilson fermions [20] with a tree-level O(a 2 ) improved Lüscher-Weisz gauge action [21] and are described in Table 1.In addition to ensemble N200 analysed in Ref. [8], we consider one ensemble at a smaller value of the pion mass, D200 and one at a larger pion mass, N203.These ensembles have the same lattice spacing a = 0.0633(4)( 6) fm (corresponding to an inverse bare gauge coupling 6/g 2 0 = β = 3.55) [22].The quark masses m up = m down = m l and m strange = m s vary along a chiral trajectory with constant sum of the bare quark masses.The bare quark masses were chosen such that the chiral trajectory approximately passes through the physical point [19].Since the bare parameters corresponding to the physical point are only determined after completing the analysis of the trajectory, some discrepancy can remain.For the ensembles considered in this work, these quark mass mistunings have been computed in [19,22].We are confident the effects these mistunings have on our analysis are small and subsequently neglect them.We introduce the quark-mass parameter, where m π , m K are the pion and kaon masses.We take these meson masses in lattice units from Ref. [19] to compute µ l .To leading-order in chiral perturbation theory µ l ≈ 3m l 2m l +ms and so along a trajectory which holds the sum of the quark masses 2m l + m s fixed, µ l is proportional to the light quark mass m l .At the N f = 3 flavoursymmetric point µ l = 1 holds.For physical quark masses and correcting for isospin-breaking effects, m π = 134.8MeV, m K = 494.2MeV [19] and consequently µ phys l = 0.1076.To mitigate against topological freezing, open boundary conditions in time [23] are used.Reweighting factors [24] are included in the analysis of observables to modify the action to the appropriate one.Statistical uncertainties are determined by the Γ-method [25,26].There are no measurable autocorrelations in our data and we do not add a tail to the autocorrelation function.
We compute the allowed energies of a system comprising a static quark at spatial position ⃗ x and a static antiquark at position ⃗ y as a function of separation ⃗ r = ⃗ y − ⃗ x.
The calculation involves the computation of a matrix of time-correlation functions C(⃗ r, t) between states created by a flux tube operator or states consisting of a static meson and anti-meson pair with isopin zero.In a theory with degenerate dynamical up and down light quarks and a strange-quark, two sets of mesons are formed by combining a static source with either a light or strange quark, giving the static-light or static-strange mesons.C(⃗ r, t), the relevant matrix correlation function needed to study this spectrum is shown schematically in Fig. 1. expressions for the interpolation operators of string and two-meson states are given in [8].All gauge links in the interpolating operators are smeared using HYP2 parameters [27,28,29,30,31] α 1 = 1.0, α 2 = 1.0, α 3 = 0.5.For the string operators, 15 and 20 levels of spatial HYP smearing with parameters α 2 = 0.6, α 3 = 0.3 [31] are used.To increase spatial resolution of the energies in the string breaking regions we use on-and off-axis distances ⃗ r.For off-axis displacements, gauge-link paths which follow the straight line connecting the static source and anti-source as closely as possible are constructed [32], using the Bresenham algorithm [33].Following this closest path significantly enhances the overlap between creation operator and physical string states.
Two distinct computational techniques to compute quark propagators (represented by the wiggly lines in Fig. 1) are needed.We refer to propagators starting and terminating on the same time slice t s = t f as "relative", while quark propagators with t s ̸ = t f are called "fixed".Evaluation of quark line diagrams uses the stochastic LapH method [34], based on the distillation quark-smearing technique [35].Distillation projects the quark fields on a time-slice into the space spanned by the lowest N v eigenmodes of the three-dimensional gauge-covariant Laplace operator, constructed from stout-smeared gauge-links [36] with parameters ρ = 0.1, n ρ = 36.Propagators between distillation spaces are estimated using stochastic LapH, with the variance reduced by dilution [37,38,39].Fixed quark propagators are evaluated on two source times t s /a = {32, 52} using full time and spin dilution with interlace-8 LapH eigenvector dilution, denoted collectively by (TF,SF,LI8) [34].Relative quark propagators are evaluated on all source times t s /a ∈ {32, 33, . . ., 95} using interlace-8 time and full spin dilution with interlace-8 LapH eigenvector dilution, labelled (TI8,SF,LI8).A total of S  [9,17] used in this work.N conf is the number of configurations on which fermionic observables were measured.N W conf is the number of Wilson loop measurements, binned into N conf bins.The third column gives the flow scale t 0 [18] in lattice units.The lattice sizes are T = Nta in time and L = Nsa in space.The last three columns list the pion mass, the kaon mass and the value of mπL from [19].
N r (fixed+relative) solutions of the Dirac equation (inversions) are required per gauge configuration for each quark flavour.Here N r denotes the number of random vectors in the space defined by the dilution projectors.For our choice of dilution schemes the number of inversions per quark flavour is 4 These parameters are summarised in Table 2: The number of Laplacian eigenmodes Nv spanning the distillation space and the inversion costs for all ensembles used in this calculation.For each type of quark line and flavour, the number of stochastic vectors in the space defined by the dilution projectors and the corresponding total number of inversions is given.The source times and the dilution schemes employed are specified in the text.
We start by extracting energies of the static-light meson E Ql and static-strange mesons E Qs .The thresholds for the avoided level crossings of string breaking are given by twice these energies.When an energy level reaches one of these thresholds, the state formed by a meson anti-meson pair becomes energetically more favourable.The staticlight and static-strange energies are obtained from single exponential correlated fits to the time-correlation function C Ql (t) and C Qs (t) of a static-light meson and a staticstrange meson respectively.For all ensembles fits are carried out in the time interval [15,21] in lattice units.Fig. 2 shows the effective masses aM eff (t + a/2) = ln[C(t)/C(t + a)] (points) as a function of t/a for ensemble D200.The top and bottom panels show the static-strange and staticlight mesons respectively.Horizontal bands indicate the values of aE Ql and aE Qs extracted from the correlator fits extending over the fitted time range.
For each fixed inter-quark separation r = |⃗ r| we solve a generalized eigenvalue problem (GEVP) using the corre- to extract the potential energies for n = 0, 1, 2 and t > t 0 .We fix t 0 = 5a for all ensembles.To ensure numerical stability of the analysis we first prune [40,41,42] the correlation matrix from size 4 × 4 down to 3 × 3, replacing C(t) with Ū contains the three most significant eigenvectors at t 0 +a.The bar in Ū and C(t 0 ) indicates these matrices are determined on the average over gauge configurations.Then we solve the GEVP in Eq. ( 2) at time t = t d = 10a, which is fixed for all the ensembles.Using the three eigenvectors v i (t 0 , t d ), i = 0, 1, 2 we project to a single correlation function ("fixed GEVP") by computing where the parentheses denote the inner product over the components of the eigenvectors.The statistical precision of our analysis is improved by exploiting the beneficial covariance between Ĉ and the static-light meson correlator C Ql (t).The correlation function ratios for n = 0, 1, 2 are then computed and fits to a single exponential with time dependence exp(−t(V n − 2E Ql )) are performed.In these ratios, energy levels become renormalised by the subtraction of twice the divergent staticlight energy since the additive self-energy contribution of the static-quarks exactly cancels.The fit ranges [t min /a, t max /a] are carefully chosen in each case by inspection.Any remaining systematic uncertainty mainly comes from t min /a dependence and we choose values of t min /a such that these effects are within the limits of statistical precision.Finally we set t max = t min + 6a.ples of effective masses aE eff (t+a/2) = ln[R n (t)/R n (t+a)] for n = 0, 1, 2 for two on-axis inter-quark separations r = 15a and r = 19a on the ensemble D200.Note the effective masses derived from ratios in Eq. ( 5) can approach their plateau values from below.The bands indicate the energy values from the exponential fits of R n and extend over the fitted time ranges.The distance 15a (bottom plot in Fig. 3) is the smallest distance down to which we can reliably determine the excited states.The separation 19a (top plot in Fig. 3) is in the middle of the first string-breaking region where the ground state and first excited state potential levels would cross.Notice all three energy levels are very close at r = 19a, and yet the GEVP is able to disentangle the three levels very efficiently.Fig. 4 shows the three lowest potential energy levels a(V n − 2E Ql ) in lattice units as a function of the distance r/a in the region where string breaking occurs for each of the three ensembles described in Table 1.Along our trajectory of decreasing light quark mass and increasing strange quark mass the energy gap between the first and second string breaking increases.In order to have a better resolution of the avoided level crossings we include off-axis distances.Some fluctuations in the size of the statistical errors can be seen especially for the off-axis distances.For on-axis distances we profit from self-averaging due to the cubic symmetry.

Model
The model used in Ref. [8] is extended to fit the data at smaller separations by the inclusion of a Cornell 1/r term.The three-state model Hamiltonian becomes After introducing a parameter γ for the 1/r term, the new model has seven parameters, labelled { Ê1 , Ê2 , g l , g s , σ, V0 , γ}.The diagonal entries in Eq. ( 6), V (r), Ê1 , Ê2 are the asymptotic energy levels for r → ∞ up to O(r −1 ).When there are N l degenerate light flavours in the mixing matrix elements between string and two-meson states of the model Hamiltonian, a flavour factor √ N l multiplies the mixing coefficient g l , in analogy to the corresponding flavour factor in Fig. 1.For N l = 2 the factor √ 2 in the three-state model Hamiltonian Eq. ( 6) can be derived from a four-state Hamiltonian with non-degenerate up and down quarks by taking the limit of degenerate light quarks.The dependence of the three eigenvalues e n (r), n = 0, 1, 2 of the model Hamiltonian Eq. ( 6) on the distance r is fitted to the three lowest potential energy levels V n (r) determined from the data.These levels are normalised by subtracting twice the energy of the static-light meson 2E Ql .Hence the fitted eigenvalues inherit this normalisation, which removes the divergent self-energy contribution from the temporal static-quark lines.The model fit parameters in Eq. ( 6) are obtained by minimisation of the correlated χ 2 , where C−1 is the inverse of the covariance matrix for the potential levels V n .The matrix C is determined by the Γ-method.For each pair of potential levels we compute the squared statistical errors of their sum and their difference.The covariance between two levels is obtained from the difference of these two squared errors.In order to obtain a positive matrix we switch off the autocorrelations between the potential levels when computing C. The distance ranges [r 1 (n), r 2 (n)] used to fit potential levels V n on individual ensembles are given in Table 3.The minimal distance for the ground state V 0 is set to 4a in order to be sensitive to the Cornell 1/r-term.For the first and second excited states, V 1 and V 2 respectively, the minimal distance is chosen so the energy gap to the ground state is smaller than 2m π .By this choice we avoid fitting our excited state levels at too small distances where there could be unresolved intermediate energy levels in our analysis.For ensembles N200 and N203 we set the maximal distance for the fit equal to L/2.Potential values at distances larger than L/2 suffer from finite volume effects and are not considered.We impose a cut on the covariance matrix C of the potential levels using a singular value decomposition.Eigenvalues of C less than 10 −4 times the maximal eigenvalue are removed from the inverse covariance which we denote by C−1 in Eq. (7).We investigated the dependence of the fit parameters on this cut parameter and found this a reasonable choice; cuts of 10 −5 or 10 −3 times the maximal eigenvalues give consistent results.
In Fig. 4 the bands show the seven-parameter model fits to Eq. ( 6) for each of the three ensembles analyzed in this work and including statistical errors.We see that the model describes our data very well in the region where string breaking happens.

Mass dependence and physical point
With three ensembles generated with different light and strange quark masses, the mass dependence of the model and its parameters can be investigated and an extrapolation to the physical values of the quark masses attempted.To begin, the difference between the static-light and static-strange meson masses was determined and a linear extrapolation to the physical light quark mass performed.The result of this extrapolation can be seen in Fig. 5.The subtraction removes the divergent static source energy, yielding a splitting with a well-defined continuum limit.We performed a fit to a simple model assuming linear dependence of this splitting on the light-quark mass parameter, µ l away from the three-flavour symmetric point corresponding to µ l = 1.The data clearly shows agreement with this simple behaviour.Extrapolation to µ phys which gives E Qs −E Ql = 74(3) stat (1) sys MeV using √ t 0 = 0.1443(7)(13) fm [22].This compares with the experimentally determined value of the mass difference, averaged over isospin for the B s and B mesons of 87.42 (14) MeV [43].The discrepancy is most likely due to finite heavy-quark-mass effects, neglected in this calculation.Other possible effects include finite-lattice-spacing artefacts and systematic uncertainty in the light-and strange-quark mass determinations.
The two mixing parameters, g l and g s of Eq. ( 6) determined at the three light-quark values are displayed in Fig. 6.Without a robust prediction from theory, their light-quark mass dependence was again modelled assuming a simple linear behaviour, cf.[44] with the two straightlines constrained to a common value at the flavour-symmetric point, corresponding to µ l = 1.For this fit, χ 2 /n df = 16.8 and some tension between the data and the model is seen close to the three-degenerate-flavour point.In the absence of a more detailed model for this mass dependence, the phenomenon was not investigated further.As seen clearly in this figure however, the two mixing parameters, g l and g s differ by a factor of almost two for physical values of the quark masses.Assuming linear behaviour yields extrapolated values of g l,s at the physical light-quark mass ratio of Fig. 7 shows the light-quark mass dependence of the coefficient γ weighting the 1/r term in V (r) modelling the gluonic flux-tube in Eq. (6).|γ| is plotted for clarity.Significant mass dependence in this term is observed and the value extrapolated to the physical point is found to be γ phys = −0.434(11).
The uncertainty is purely statistical as there is no scale dependence on this dimensionless parameter.Understanding the mass dependence in detail would require closer study of the transition from the short-distance behaviour of the static potential into the string-breaking region and is beyond the scope of this first study.The dependence of the string tension on µ l is seen in Fig. 8, along with two simple models to describe this behaviour and extrapolate to the physical point.The simplest model asserts no dependence of σ on µ l while the second model assumes linear dependence.The change in σ over the full range of µ l is seen to be mild and is only a few percent from the three-flavour symmetric theory to the physical point.The mass dependence is not described very well by either model and the statistical uncertainties at the physical point from the two fits do not capture the where the uncertainty quoted includes both statistical and extrapolation uncertainties from our determination of σt 0 combined in quadrature with the scale-setting uncertainty from √ t 0 = 0.1443(7)(13) fm [22].lations between the model parameters after extrapolating to the physical quark masses.Some parameters are seen to be strongly correlated, in particular the three coefficients in V (r) have high statistical correlations, including a pair with a correlation of 92%.The right table in Table 4 lists the model parameters with their statistical errors after linear extrapolations in µ l (best fit) to the physical quark masses.
Using the model with parameters extrapolated to the physical point, the dependence of the energy spectrum of the static source-anti-source system on source separation in the presence of physical up, down and strange quarks can be computed and is presented in Fig. 9.The zero of the potential is fixed by subtracting twice the energy of the static-light system, which makes the ground-state energy at large separations vanish asymptotically.The excited states are not displayed at short separations, where the basis of operators may not couple efficiently to possible multi-meson states and so the determination of the spectrum in this range for this calculation is uncertain.
The string-breaking distances, r ⋆ l and r ⋆ s are indicated by two vertical bands in Fig. 9.These distances are derived from the model by determining the static source-sink separation where V (r) is equal to the asymptotic energy of two static-light (for r ⋆ l ) or static-strange (for r ⋆ s ) mesons.This yields r ⋆ l = 8.39(3) √ t 0 = 1.211 (7) stat (11) sys fm, r ⋆ s = 9.26(2) √ t 0 = 1.336 (7) stat (12) sys fm.(11) Again, the values quoted in physical units include estimates of the uncertainties arising both from limits to the statistical precision of our calculation and systematic uncertainty from setting the physical scale.

Conclusions
This paper extends the analysis of Ref. [8] to investigate the light-quark mass dependence of the potential energy of a system of a static colour-anti-colour source pair.The three lowest energy levels in the resulting spectrum are determined up to the scale where mixing between states resembling a gluon string and two static-light or static-strange mesons is largest.A robust determination is achieved by computing correlations within a suitable variational basis of interpolating fields and solving the resulting generalised eigenvalue problem.A simple model of the spectrum is presented, starting from the Cornell potential and allowing mixing with the asymptotic static-meson states.
Using the new data from calculations with different light-quark masses, an extrapolation of the model parameters to the physical quark mass values is performed, enabling us to predict the physical potential in QCD with dynamical light and strange quarks and to give a simple parameterisation of this potential and its excitations.Perhaps the most studied parameter is the string tension where the value determined in this work is seen in Eq. (10).In comparison, the string tension in pure-gauge theory is in the range σt 0 ∈ [0.143, 0.159].These values are computed using data on the deconfining temperature T c /σ compiled in [45] combined with T c r 0 from [45] and r 0 / √ t 0 from [46].

Figure 1 :
Figure 1: Schematic representation of the string breaking mixing matrix.

Figure 2 :
Figure 2: Effective masses of the static-light (bottom) and staticstrange (top) mesons on the ensemble D200.The horizontal bands indicate the time intervals used and energies obtained from a single exponential fit to the corresponding correlators.

Fig. 3 Figure 3 :
Figure 3: Examples of effective masses of the ratios Eq. (5) on the ensemble D200 at two distances as indicated in the plots.The horizontal bands correspond to the energies Vn(r) − 2E Ql , n = 0, 1, 2 extracted from exponential fits to the ratios and extend over the fitted time interval.

Figure 5 :
Figure 5: Dependence of the mass difference between the static-light and static-strange meson masses on the light-quark mass parameter µ l of Eq. (1).

Figure 6 :
Figure 6: Dependence of the mixing parameters, g l and gs on the light-quark mass parameter, µ l of Eq. (1).The vertical line indicates the physical light-quark mass ratio.

Figure 7 :
Figure 7: Light-quark mass-dependence of the coefficient of the dimensionless 1/r Cornell potential term from fits to the model of Eq. (6).The vertical line indicates the physical light-quark mass ratio.

Figure 8 :
Figure 8: Dependence of the string tension on the light and strange quark masses.µ l , the dimensionless light quark mass parameter is defined in Eq. (1).The two bands indicate extrapolating fits using constant or linear dependence.The physical value of the mass ratio, µ phys l

Table 2 .
Note correlation functions for off-axis static source separations are easy to compute using relative and fixed quark propagators, unlike the off-axis gauge-link paths.

Table 3
(6),22]: The energy levels extracted from the GEVP analysis for separations in the range r/a ∈[17,22]for the three ensembles described in Table1.The data include off-axis separations.The bands indicate the energies extracted from a fit of data to the seven-parameter model of Eq.(6) Table 4 shows the corre-