String breaking by light and strange quarks in QCD

The energy spectrum of a system containing a static quark anti-quark pair is computed for a wide range of source separations using lattice QCD with $N_\mathrm{f}=2+1$ dynamical flavours. By employing a variational method with a basis including operators resembling both the gluon string and systems of two separated static mesons, the first three energy levels are determined up to and beyond the distance where it is energetically favourable for the vacuum to screen the static sources through light- or strange-quark pair creation, enabling both these screening phenomena to be observed. The separation dependence of the energy spectrum is reliably parameterised over this saturation region with a simple model which can be used as input for subsequent investigations of quarkonia above threshold and heavy-light and heavy-strange coupled-channel meson scattering.


Introduction
QCD is believed to be responsible for confinement, the experimental observation that quarks are never seen as asymptotic states [1]. A full understanding of why QCD confines remains elusive. The simplest theoretical probe of the phenomenon is provided by the potential energy V (r) of a system made of a static quark and anti-quark pair immersed in the QCD vacuum in a colourless combination. This energy depends on the distance between the sources, and in the Yang-Mills theory of gluons alone it grows linearly at asymptotically large separations. The rate of increase is the well-known string tension, σ. If the static sources interact with the full QCD vacuum including light-quark dynamics, pair-creation of light quarks means the system can also resemble two separately-colourless static-light mesons allowing the potential energy to saturate at large separations. This phenomenon, induced by light-quark pair creation, is termed 'string breaking'.
Arising from purely non-perturbative effects, the static potential requires a robust method for direct determination from the QCD Lagrangian. Lattice QCD provides such a framework, but studying string breaking on the lattice is a technical and numerical challenge. The simplest approach would be to evaluate the expectation value of large Wilson loops in the QCD vacuum and observe deviation from an area law. However, Monte Carlo determinations of large loops suffer from poor signal-to-noise properties. Also, viewed in terms of eigenstates of the Lattice QCD Hamiltonian, the creation operator comprising the Wilson line forming one edge of the rectangular Wilson loop has very small overlap onto the two-meson system dominating the ground state at large separations. As a consequence, the saturation effect proves impossible to resolve from Monte Carlo studies of Wilson loops alone. Instead, a mixing analysis including two-meson 'broken' string states is needed [2][3][4][5][6].
In this work, by building a suitably diverse basis of creation operators, mixing between the state made by a gluonic flux tube and the broken string state resembling two static-light or staticstrange mesons is investigated fully in the N f = 2 + 1 theory on the lattice for the first time. This enables us to compute reliably the energies of the lowest three Hamiltonian eigenstates up to separations where string breaking saturation occurs. A simple parameterisation of the resulting spectrum in the breaking region is given, which should provide invaluable first-principles input into models of coupled-channel scattering of heavy-light mesons and the decays of quarkonia near threshold.

Methodology
We compute the potential energy of a system containing a heavy quark Q at spatial position x and a heavy anti-quarkQ at y in the static approximation. In this limit, the quarks remain separated by r = y − x, where x and y are conserved quantum numbers. To determine the energies of the ground state as well as the first-and second-excited state arising from mixing in QCD, a variational technique is employed. The input for this mixing calculation is a matrix of temporal correlation functions between the interpolator for a Wilson line O W , the two-static-light O BB and the twostatic-strange meson states O BsBs . It is essential that the interpolators transform irreducibly under the appropriate symmetry group.

Interpolators
A suitable interpolator O W creating a gluon string connecting sites x and y at time t is given by where γ is the three-vector of spatial Dirac matrices and the Wilson line W(y, x, t) is a product of spatial lattice links with time argument t connecting x and y. Wilson lines are constructed using a variant of the Bresenham algorithm [7] to approximate the shortest connection between x and y. The heavy-quark spins are coupled symmetrically via γ · r/r, which has a zero component of angular momentum projected along r. In the static limit, both the symmetric and antisymmetric combinations give an energy level in the Σ + g irreducible representation of the rotation group around r after the heavy spins are decoupled. Decoupling the heavy spins in the two-static-meson system similarly yields a composite operator in the Σ + g irrep. Details of this construction can be found in Ref. [6]. For mass-degenerate quark flavours q i , i = 1, ..., n f , the interpolators for a two-staticlight-or two-static-strange-meson state are given by For the two-static-light-meson state, the sum projects onto the isospin-zero channel.

Correlation matrix
For the N f = 2 + 1 theory, a 3 × 3 matrix can be constructed from the pair-wise correlations of the two isospin-zero meson-pair creation and annihilation operators in combination with the string interpolation operator When mixing occurs, the physical energy eigenstates are not natural basis vectors and off-diagonal elements of the correlation matrix are non-vanishing. Following the method presented in [8], all gauge-links are smeared using HYP2 parameters [9, 10] α 1 = 1.0, α 2 = 1.0, α 3 = 0.5, where the smearing of the temporal links amounts to a modification of the Eichten-Hill static action and propagator [11] which reduces the divergent mass renormalisation and improves the signal-to-noise ratio at large Euclidean times [10]. As a second step, we construct a variational basis for the string state using 15 and 20 levels of HYP-smeared spatial links with parameters: Some correlation functions in C(r, t) include multiple light-quark field insertions. The lightquark fields must be integrated analytically prior to Monte Carlo evaluation of C(r, t) and the resulting Wick contractions involve numerically challenging disconnected contributions. In order to calculate these contributions, and to reduce statistical variance by exploiting translational invariance of the lattice, propagators between all space-time points are needed. For this, we employ the stochastic LapH method [12]. Based on distillation [13], the method facilitates all-to-all quark propagation in a low-dimensional subspace spanned by N v low-lying eigenmodes of the threedimensional gauge-covariant Laplace operator, which is constructed using stout-smeared gauge links [14] with parameters (ρ, n ρ ). This projection onto the so-called LapH subspace amounts to a form of quark smearing, where N v increases in proportion to the spatial volume for a fixed physical smearing radius. Introducing a stochastic estimator in the LapH subspace in combination with dilution [15,16] helps to reduce the rise in computational costs as the volume increases. It was shown [12] that for certain dilution schemes, the quality of the stochastic estimator remains approximately constant for increasing volume, while maintaining a fixed number of dilution projectors. The triplet b = (T, S, L) specifies a dilution scheme with time T , spin S and LapH eigenvector projector index L, where F indicates full dilution and In the interlacing of dilution projectors in index n. The dilution scheme and other stochastic-LapH parameters used are given in table 1.
id N ev n ρ × ρ line type dilution scheme N r light,strange source time t s /a N200 192 36 × 0.1 fixed (TF,SF,LI8) 5 , 2 32, 52 relative (TI8,SF,LI8) 2 , 1 - Table 1: Number of eigenvectors, stout-smearing parameters, dilution schemes, number of noise sources, and source times employed in this work. Quark lines starting and terminating at the same time slice t s = t f are referred to as 'relative', while quark lines with t s = t f are called 'fixed'. For definition of the LapH subspace and specification of dilution schemes, see [12].

Variational analysis
A correlation matrix is evaluated for a set of static-source separations r and time separations t. The energy spectrum of the system for each spatial separation can be extracted by solving a generalised eigenvalue problem (GEVP) for each r [17,18], where λ n and v n are the eigenvalues and eigenvectors respectively. After solving the GEVP, energies are extracted using a correlated-χ 2 minimisation to a two-parameter, single-exponential fit ansatz. We solve the GEVP for a fixed pair of reference time separation t 0 and another separation t d at which the correlation matrix is diagonalized [19]. Interpolators with maximal overlap onto a single eigenstate are defined asĈ where the parentheses denote an inner product over GEVP indices. A potential source of systematic error is introduced, as off-diagonal elements ofĈ ij are not exactly zero. We control this by assessing the stability of the GEVP against varying the operator basis and using different pairs (t 0 , t d ). Comparing the results for some sample distances to the GEVP as given in Eq. (2.4) shows agreement between both methods. We are interested in the difference between the energies V n (r) of the ground (n = 0), first (n = 1) and second excited state (n = 2) and twice the energy of the static-light meson 2E B , which can be directly extracted from fits to the ratio  [20]. Open temporal boundary conditions are imposed on the fields. The imprint on observables is expected to fall exponentially with distance from the boundary [21] so measurements are made on the central half of the lattice only.
Previous studies found the signal quality to be limited by the Wilson-loop correlation functions [6], which a preliminary analysis on a subset of our data corroborated. We mitigate this issue by measuring Wilson loops on 1664 configurations, while diagrams containing light or strange propagators are evaluated on a subset of 104 samples. The Wilson loops are then averaged into 104 bins containing 16 configurations each, with the centre of the bin aligned with one entry in the 104-configuration subset. For each separation r = (r 1 , r 2 , r 3 ), we exploit cubic symmetry and average over spatial rotations to increase our statistics. As expected, no dependence on the direction is observed in our data. The analysis uses a Jupyter notebook 1 adapted from Ref. [22]. Statistical uncertainties are estimated using 800 bootstrap resamplings [23,24] and the uncertainty quoted is given by 1σ bootstrap errors. The covariance matrix entering the single exponential fits to Eq. (2.6) is estimated once and kept constant on every bootstrap sample.
In Fig. 1, the potential energy relative to 2E B is shown; this subtraction removes the divergent mass renormalisation of a static quark. The grey line corresponds to twice the static-strange mass, with an energy difference 2E Bs − 2E B = 0.014(2)a = 42 (7)MeV. As expected due to the choice of light and strange quark masses, the difference is smaller than the physical energy difference between B and B s mesons.
The avoided level crossing between the ground-and first-excited states is clearly visible and the expected second avoided crossing due to the formation of two static-strange mesons is also evident for the first time. As the distance over which this phenomenon occurs is small, this effect can not be resolved using only on-axis separations [25]. For distances beyond the string breaking scale, the ground state tends rapidly towards the mass of two non-interacting static-light mesons.
The inset in Fig. 1 shows the string breaking region in more detail. Both avoided crossings are visible and the energy gap between the ground state and first level is larger than the gap between first and second level. Qualitatively, the first mixing region appears to be broader, but it is not possible to determine the difference between the first string breaking distance r c and the second string breaking distance r cs by eye. The quantification of string breaking involving three levels is more complex in comparison to the two-level situation. For the N f = 2 vacuum, the string breaking distance r c can be defined by the minimum of the energy gap ∆E [6]. When the strange quark is included, an alternative definition of the two string breaking distances r c and r cs is needed as there is not necessarily a minimum energy gap.

A model of the string breaking spectrum
We describe the potential-energy spectrum in the breaking region using a simple Hamiltonian that extends the model for N f = 2 given in [26]. Consider a three-state system with Hamiltonian: (4.1) The diagonal elements are a functionV (r) describing the unbroken string andÊ 1 ,Ê 2 , the energies of a noninteracting pair of static-light and static-strange mesons, respectively. As withV , these energies are measured relative to 2E B . g 1 and g 2 are two coupling constants describing the strength of the mixing between the gluon flux and the separated colour-screened static sources. The offdiagonal term that would mix the two-static-light-meson state with the two-static-strange-meson state is set to zero. There is no constraint on this mixing in the energy spectrum as we assume no dependence on the source separation for this effect. A basis rotation shows any non-zero value of this parameter yields an equivalent spectrum to the Hamiltonian of Eq. (4.1). A suitable choice for the function representing the string state is the Cornell potential [27]. Since we are modeling the string breaking region and not the potential at small distances, only the with string tension σ, is constrained by our data. The eigenstates of the Hamiltonian are mixtures of the unbroken string and the two staticlight-and two static-strange-meson states while the eigenvalues of H correspond to the three extracted energy levels. After diagonalizing H, we perform an uncorrelated six-parameter fit to the spectrum, whose result is shown in Fig. 2. We find for our fit parameters: The model in Eq. (4.1) assumes a three state system. While it is possible that the physical eigenstates receive contributions from higher lying states Fig. 2 shows that our data is welldescribed by the fit parameters.
We now turn to a quantitative definition of the string breaking distance in the N f = 2+1 case. To investigate the dependence of the string breaking distance on the sea quark masses, a robust definition is needed. By using the asymptotic states of our model, which for large distances r are given by the diagonal entries of our Hamiltonian, we extract two distinct string breaking distances corresponding to the light and strange mixing phenomenon. The string breaking distance, r c associated with the formation of two static-light quarks is given by the crossing distance wherê V (r c ) =Ê 1 and a corresponding definition can be employed to define r cs usingV (r cs ) =Ê 2 . The quoted errors for the physical units take into account the uncertainty of a = 0.06426(76)fm [20]. Fig. 3 shows that the energy gap ∆E 1 (r) = V 1 (r) − V 0 (r) between the ground and firstexcited state does not exhibit a minimum, thus making it impossible to use the minimal gap distance to define r c . It can also be observed that in the string breaking region, ∆E 1 (r) is twice as large as ∆E 2 (r) = V 2 (r) − V 1 (r) , the energy gap between the first-excited and second-excited state.
Only one previous study [6] of string breaking for N f = 2 QCD on the lattice exists with m π = 640MeV. However, the differing definitions and quark content of the vacuua mean it is not possible to make a statement on quark mass dependence. Using the position of the minimal energy gap, they findr c ≈ 1.248(13)fm. Even though in this study the pion mass was relatively heavy,r c is of the same order of magnitude as our result and falls between our values for r c and r cs .

Phenomenology
To make a first, simple comparison between our data and the experimentally observed quarkonium spectrum, we use the ground-state energy computed from our model Eq. (4.1) as input potential in the Schrödinger equation and extract the bound state energies of quarkonia in the Born-Oppenheimer approximation. We use input bottom-and charm-quark masses m b = 4977MeV, m c = 1628MeV from quark models [28,29].  The number of bottomonium bound states is listed in the first row of Tab. 2. We observe three S-wave states in agreement with the physical spectrum of Υ mesons. We do not find bound states beyond l > 4. For the states closest to the threshold for each angular momentum, the size of the wave function as determined from the root mean-square radius is between 0.63 r c and 0.79 r c . We also solved the Schrödinger equation when setting V (r) = min{σr +V 0 , 0} and find changes in energies of less than 10 MeV, within the errors from the fit parameters, so the bound-state spectrum is not sensitive to the width of the mixing region. The second row of Tab. 2 shows the number of bound states we find for charmonium for comparison.
This calculation is only a check that counting the bound states leads to sensible results. Our model Eq. (4.1) provides input for more refined calculations of the coupled-channel scattering of heavy-light mesons and the decays of quarkonia around threshold.

Conclusions
This work presents a first calculation of string breaking with N f = 2 + 1 dynamical quarks from lattice QCD on a single ensemble with light/strange quark masses that are heavier/lighter than in nature, corresponding to m π = 280 MeV and m K = 460 MeV. We compute the three lowest energy levels of a static quark and anti-quark pair and observe the avoided level crossings due to pair-creation of light and strange sea quarks. Our main result can be summarized in the model of Eq. (4.1) which provides a very good description of our spectrum, shown in Fig. 2. In particular we extract the values of the parameters g 1 = 47.2(1.4) MeV and g 2 = 24.5(1.6) MeV which describe the mixing between the gluonic flux tube and the broken string.
Further anaylsis to study the dependence of string breaking on the quark masses is in progress.