Precision computation of hybrid static potentials in SU(3) lattice gauge theory

We perform a precision computation of hybrid static potentials with quantum numbers $\Lambda_\eta^\epsilon = \Sigma_g^-,\Sigma_u^+,\Sigma_u^-,\Pi_g,\Pi_u,\Delta_g,\Delta_u$ using SU(3) lattice gauge theory. The resulting potentials are used to estimate masses of heavy $\bar{c} c$ and $\bar{b} b$ hybrid mesons in the Born-Oppenheimer approximation. Part of the lattice gauge theory computation, which we discuss in detail, is an extensive optimization of hybrid static potential creation operators. The resulting optimized operators are expected to be essential for future projects concerning the computation of 3-point functions as e.g. needed to study spin corrections, decays or the gluon distribution of heavy hybrid mesons.


Introduction
The success of the quark model, following the realization of the importance of SU(3) flavor symmetry in the context of the eightfold way, led to understanding the properties of a large number of mesons and baryons. However, the quark model does not contain gluons. In the framework of QCD it is, thus, of utmost importance to investigate and to understand what kind of additional hadronic states or resonances can appear, when gluons are allowed to be in excited states. In this work we are particularly interested in heavy hybrid mesons, i.e. in mesons composed of heavy c or b quarks, where gluons contribute to the quantum numbers J P C in a non-trivial way.
With the discovery of the first of the so-called XY Z mesons around fifteen years ago, the X(3872), an entirely new chapter of hadronic physics was opened and flourished. At present there are about thirty such exotic states, which have been observed (for a theoretical summary cf. e.g. [1][2][3][4], for an experimental review cf. e.g. [5]). Many of these exotic states are believed to be tetraquark resonances, but some of them are also being considered as candidates for hybrid mesons. It is very challenging to understand the internal structure of exotic hadrons and, even though there is little doubt that hybrid mesons and baryons exist, not much else is known about them.
It is a notable feature of hybrid mesons that, due to their excited gluonic degrees of freedom, part of them have J P C quantum numbers, which are forbidden in the quark model 1 . In this sense, observing a meson with J P C = 0 +− , 0 −− , 1 −+ , 2 +− , . . . indicates an exotic structure, possibly the presence of excited gluons.
Obtaining solid results for exotic hadrons is highly non trivial, both on the theoretical and on the experimental side. There is currently a lot of experimental activity in the field of exotic hadrons, but even the exact attribution of the J P C quantum numbers is difficult for many of the experimentally observed XY Z states. Moreover, all up to now experimentally observed heavy candidates for hybrid mesons exhibit non-exotic quantum numbers, which makes theoretical investigations of their properties even more important. There are a few exotic states, which could be heavy hybrid mesons, the most prominent candidate being the Y (4260) with quantum numbers J P C = 1 −− , but there are also arguments disfavoring a hybrid identification (cf. e.g. the discussions in [6,7]). Several existing and future experiments will be taking data in the next couple of years (for example the GlueX and the PANDA experiment) and, thus, many more candidates for hybrid mesons are likely to be discovered. On the theoretical side models devised to explain the properties of exotic hadrons typically possess a limited applicability and there is no overall coherent theoretical picture of these hadrons. Lattice field theory, however, which is a non-perturbative first principles approach, is an ideal method to predict masses or to explore properties of hybrid mesons. Such results might also be useful as input for effective theories like pNRQCD or to calibrate or devise improved theoretical models.
In this work we carry out a precise computation of several hybrid static potentials using SU (3) lattice gauge. Gluonic excitations are included by considering trial states containing a static quark-antiquark pair and gluons, which are characterized by non-trivial quantum numbers, i.e. orbital angular momentum, parity or charge conjugation. Our aim is to improve on existing sim-ilar lattice field theory computations ) by providing results with smaller statistical errors and at finer spatial resolution (section 5) and by discussing all technical details of the optimization of creation operators and the computation of the potentials (section 2 to section 4). The latter is a necessary and important preparatory step for the computation of 3-point functions, which we recently started, and which we briefly discuss in our conclusions in section 7.
We also use some of the resulting hybrid static potentials to estimate masses of heavy-quark hybrid mesons, where the quarks are eithercc orbb (section 6). This is done in the Born-Oppenheimer approximation [34], where effects from the quark spins are neglected, by numerically solving an appropriate Schrödinger equation. This is expected to be a good approximation, because the time scales of the gluons and of the heavy charm or bottom quarks are significantly different and, thus, their dynamics decouples almost completely. In this context also effective theory approaches like potential Non Relativistic QCD (pNRQCD) are extremely useful, for example when parameterizing discrete lattice field theory results for hybrid static potentials by continuous functions (for recent pNRQCD articles on hybrid mesons cf. e.g. [7,35,36]).

Quantum numbers and trial states
A hybrid static potential is a potential of a static quark and a static antiquark, where the gluons form non-trivial structures and contribute to the quantum numbers. We compute such hybrid static potentials from Wilson loop-like correlation functions using SU(3) lattice gauge theory. The gluonic excitations are realized by replacing the straight spatial Wilson lines of the Wilson loops by parallel transporters, which have a less trivial structure.
Hybrid static potentials are characterized by the following quantum numbers: • Λ = 0, 1, 2, . . ., the absolute value of the total angular momentum with respect to the axis of separation of the static quark-antiquark pair, i.e. with respect to the z axis.
• η = +, −, the eigenvalue corresponding to the operator P •C, i.e. the combination of parity and charge conjugation.
• = +, −, the eigenvalue corresponding to the operator P x , which denotes the spatial reflection along the x axis, which is perpendicular to the axis of separation of the static quark-antiquark pair.

Angular momentum Λ
We start in the continuum and consider hybrid static potential creation operators and trial states where |Ω is the vacuum and R(ϕ) denotes a rotation by an angle ϕ around the z axis. Moreover, where Q(+r/2) andQ(−r/2) are operators creating a spinless quark-antiquark pair and U (−r/2, r 1 )S(r 1 , r 2 )U (r 2 , +r/2) is a parallel transporter connecting the quark and the antiquark in a gauge invariant way. U (−r/2, r 1 ) and U (r 2 , +r/2) denote straight parallel transporters along the z axis (in the simplest case r 1 = −r/2 and r 2 = +r/2, i.e. U (−r/2, r 1 ) = U (r 2 , +r/2) = 1), while the operator S(r 1 , r 2 ) is different from a straight line and, thus, generates a gluonic excitation. It is easy to show that the trial state (1) has definite angular momentum Λ (see appendix A).
The corresponding lattice expression is where the angle of rotation is restricted to multiples of π/2 and U (−r/2, r 1 ), S(−r/2, +r/2) and U (r 2 , +r/2) are products of gauge links. For example for Λ = 1 i.e. one has to compute Wilson loops, where each of the straight spatial Wilson lines is replaced by a sum over the four rotations of the operators O with weight factors +1, +i, −1 and −i.

P • C and P x quantum numbers η and
It is straightforward to show where S P•C (−r 2 , −r 1 ) is the charge conjugated spatial reflection of S(r 1 , r 2 ) with respect to the center of the separation axis. Consequently, one has to include both S and S P•C in the final operator, to obtain a trial state with definite η. Similarly, where S Px (r 1 , r 2 ) is the spatial reflection of S(r 1 , r 2 ) along the x axis.
To construct a trial state, which has definite quantum numbers Λ η , we take the state (3), which has angular momentum Λ, and project that state onto the subspace of eigenstates of the operators P • C and P x characterized by η and , respectively: with projectors and Notice that not every operator S(r 1 , r 2 ) is suited to construct trial states for any given set of quantum numbers Λ η = Λ η , i.e. for some Λ η the trial state defined in (7) is zero, |Ψ hybrid S;Λ η = 0. Also note, that even though sectors with Λ ≥ 1 (i.e. in this work the Π and the ∆ sectors) are degenerate with respect to the quantum number , creation operators constructed via eq. (7) with either = + or = − might be different, i.e. yield non-identical correlation functions. In such cases it is important to consider the = + and = − operators as separate creation operators, when identifying an optimal set of creation operators for each correlation matrix (cf. section 4.3; an example is the Π g sector, when using S IV,1 ). In practice, we use eq. (7) to automatically generate creation operators with definite quantum numbers Λ η from all considered operators S(r 1 , r 2 ) (cf. Figure 1 for a graphical illustration of an example).

Correlation functions
We determine hybrid static potentials with quantum numbers Λ η , which we denote by V Λ η (r), from the asymptotic exponential behavior of temporal correlation functions Expressing eq. (10) in terms of a path integral and performing the integration over the static quarks leads to where U (r; t 1 , t 2 ) denotes a straight line of temporal gauge links at r from time t 1 to t 2 and . . . U is the average on an ensemble of gauge link configurations distributed according to e −S . The right hand side of this equation can be computed using standard techniques from lattice field theory as briefly summarized in section 3.

Lattice setup
All computations presented in this work have been performed using SU(3) lattice gauge theory. The gauge link configurations have been generated with the standard Wilson gauge action (cf. standard textbooks on lattice field theory, e.g. [40]) and the Chroma QCD library [41]. Since we are considering purely gluonic observables, we expect that there is little difference between our SU(3) Yang-Mills results and corresponding results in full QCD (cf. also the discussion of systematic errors at the end of section 6.2 and [21]).
In this work we use a single ensemble with lattice extent 24 3 × 48 and gauge coupling β = 6.0 corresponding to lattice spacing a ≈ 0.093 fm and spacetime volume ≈ (2.22 fm) 3 × 4.44 fm, when identifying r 0 with 0.5 fm (for a determination of r 0 cf. section 6.1). The gauge link configurations are separated by 20 lattice updates, where each update comprises a heatbath and four over-relaxation steps. We have performed standard binning analyses with bins containing either 1, 2 or 4 gauge link configurations. We have found that the statistical errors of the Σ + g potential are essentially independent of the bin size, which indicates that performing 20 lattice updates largely eliminates correlations in Monte Carlo time. For the final results for hybrid static potentials presented in section 5 and Table 8 we have generated more than 5 500 gauge link configurations. During the time-consuming optimization of hybrid static potential creation operators and trial states discussed in section 4, we use a subset of 100 gauge link configurations, to reduce the computational effort to an acceptable level.
To improve the signal quality, standard smearing techniques are applied to the gauge links of the Wilson loop-like correlation functions (10). The temporal gauge links in U (r; t 1 , t 2 ) are HYP2 smeared gauge links [42][43][44], which lead to a reduced self energy of the static quarks and, consequently, to smaller statistical errors. The spatial gauge links in a S;Λ η (r 1 , r 2 ; t) are APE smeared gauge links (for detailed equations cf. e.g. [45]), where the parameters are tuned to optimize the ground state overlaps (cf. section 4.2) and, thus, allow to extract the potentials at smaller temporal separations.
All statistical errors shown and quoted throughout this paper, e.g. for the hybrid static potentials in section 5 or the potential parameterizations and hybrid meson masses in section 6, are determined via an evolved jackknife analysis starting at the level of the correlation matrices.
To exclude statistical correlations between gauge link configurations, which are close in Monte Carlo simulation time, we perform a suitable binning of these configurations.

Optimization of hybrid static potential creation operators and trial states
Since the signal-to-noise ratio of correlation functions (10) decreases exponentially with respect to the temporal separation, it is essential to identify hybrid static potential creation operators, which generate trial states with large ground state overlap. This allows to extract hybrid static potentials at rather small temporal separations, where the signal-to-noise ratio is favorable.
The starting point is a large set of quite distinct operators S, some of them simple, others of more complicated shape, which are shown in Figure 2. All these operators extend over regions, which are of the same order as the quark-antiquark separation. This is quite different from our previous exploratory study in SU(2) Yang-Mills theory [29], where we used local chromoelectric and chromomagnetic field strength insertions. While the latter are theoretically easier to handle and thus are quite common in analytical studies, e.g. based on pNRQCD [7], the extended operators we are using here are much better suited for numerical lattice field theory studies, because they lead to trial states with larger ground state overlaps and, thus, to results with significantly smaller statistical errors.
The operators S can be categorized into planar operators, Some of these operators, e.g. S I,1 , S I,3 , S III,1 , S III,2 , S III,4 , S III,5 , S IV,1 , S IV,2 , S IV,3 , S IV,5 and S V,2 were already used in previous lattice field theory studies of hybrid mesons, e.g. in [13], while other operators are explored for the first time in this work. From these operators we construct a large number of different trial states |Ψ hybrid S;Λ η using eqs. (7) to (9).
To check, whether a trial state |Ψ hybrid S;Λ η has large ground state overlap, we compute the effective mass at small temporal separations, in particular at t = a, where contributions of excited states are most prominent. Small effective masses indicate trial states with large ground state overlaps, while operators leading to large effective masses can be discarded. In the following we discuss in detail, how we identify and optimize a small set of relevant operators S for each hybrid static potential sector Λ η .
Figure 2: Operators S used to generate trial states |Ψ hybrid S;Λ η according to eqs. (7) to (9). Each arrow represents a straight path of gauge links. Arrows with the same color (red, green or blue) have the same length, i.e. represent the same number of gauge links. Dotted arrows can have length zero, while solid arrows represent at least one gauge link. If the starting point and the end point of S are marked by black dots, U (−r/2, r 1 ) and U (r 2 , +r/2) in eq. (2) have the same length, i.e. r 1 − (−r/2) = +r/2 − r 2 , else their length can be different. n ∈ {1, 2, 4} is the number of differently oriented operators (i.e. operators, which cannot be transformed into each other by rotations around the z axis) obtained by applying P • C, P x and (P • C)P x .

Optimization of the extents of the operators S
In a first step we consider each of the operators S shown in Figure 2 separately and optimize their extents for each hybrid potential sector Λ η . In other words, for each arrow in Figure 2 we determine the number of gauge links it represents, such that the ground state overlap of the corresponding trial state is maximal.
As an example we briefly discuss the optimization of the operator S I,1 for the Π u hybrid static potential. Variations of S I,1 are denoted by S Ex,Ez , where E x and E z are the operator extents in units of the lattice spacing in the x direction and the z direction, respectively. To keep the computational cost of the optimization on a feasible level, we first determine the optimal value for E x and after that the optimal value for E z . Figure 3 shows that the optimal value for E x weakly depends on the separation of the quark-antiquark pair r. For r/a ≤ 3 the operator extent E x = 2 minimizes the effective mass V eff;S Ex,Ez I,1 ;Πu (r, t = a) and, thus, the corresponding trial state has better overlap to the ground state in the Π u sector than trial states generated with E x ∈ {1, 2, 4}. Similarly, for r/a ≥ 4 the operator with extent E x = 3 minimizes V eff;S Ex,Ez I,1 ;Πu (r, t = a) and, thus, maximizes the ground state overlap. Since operator extents ;Πu (r, t = a) for any of the considered quark-antiquark separations r, they are discarded. In Figure 4 we show an analogous comparison of effective masses V eff;S Ex,Ez I,1 ;Πu (r, t = a) for different E z and the previously optimized E x ∈ {2, 3}. We find that the optimum is E z = r/a independent of r and E x , i.e. the z extent of operator S I,1 should be identical to the quark-antiquark separation, when used to compute the ground state hybrid static potential in the Π u sector. ;Πu (r, t = a) on E x . Red spheres, red arrows, and black dots represent quarks, gauge links and lattice sites, respectively. ;Πu (r, t = a) on E z (E x ∈ {2, 3}, the optimum according to Figure 3). Red spheres, red arrows, and black dots represent quarks, single gauge links and lattice sites, respectively.
All 19 operators S shown in Figure 2 are optimized for each sector Λ η and each separation r in a similar way. In the majority of cases more than two extents have to be optimized.

Optimization of APE smearing parameters
To further improve the ground state overlap of the trial states, we use APE smeared spatial gauge links in a S;Λ η (r 1 , r 2 ; t) in (11). For detailed equations cf. e.g. [45]. We set α APE = 0.5, which is a common choice in the literature and we investigate the dependence of V eff;S;Λ η (r, t = a) on the number of APE smearing steps N APE . An example plot for operator S 1,r/a I,1 and the Π u hybrid static potential is shown in Figure 5. While there is a significant increase of the ground state overlap, when increasing N APE from 0 to around 20, there is no further gain, when using N APE > 20. This behavior is observed for various quark antiquark separations r = 2a, . . . , 8a. Similar findings are obtained also for the other operators S shown in Figure 2 and for all sectors Λ η . Therefore, we use N APE = 20 for all computations presented throughout this paper.

Selecting optimal sets of trial states
To further improve the ground state overlaps of the trial states, we resort to variational techniques for our final analyses in section 5. For each sector Λ η we use an "optimal set" of operators S, compute the corresponding correlation matrix (11) and solve generalized eigenvalue problems (see e.g. [46,47]). In this way the static potentials are determined using an optimized linear combination of creation operators. To keep the computational effort on an acceptable level, we have restricted these variational analyses to the three or four most promising operators S for each sector Λ η and separation r.
To select these operators, we have first performed an optimization of the extents of each operator as discussed in section 4.1. We have then taken those three or four operators, which yield the smallest effective masses (12) at t = a. Results are collected in the Table 1 to Table 7. Each table corresponds to another hybrid static potential Λ η . The operators are sketched in the left column of each table. In two cases (Λ η = Σ + u and Λ η = Π g , i.e. Table 2 and Table 4) not only the extents, but also the operators S change with the separation r, indicated by "-". The operators are also defined mathematically in the tables. For example the left hand side Table 1) represents U (−r/2, r 1 )S(r 1 , r 2 )U (r 2 , +r/2) from eq. (9), while the right hand side denotes 2 links in positive x direction, 2 links in positive y direction, E z links in positive z direction (where E z as a function of r is listed directly below), 2 links in negative y direction, 2 links in negative x direction.
r/a 2 3 4 5 6 7 8 9 10 11 12 E z 2 3 4 5 6 7 8 9 10 11 12   Table 6: Optimized creation operators for V ∆g (r). Note that, even though the ∆ g hybrid potential is degenerate with respect to , the construction of creation operators via eq. (7) is not independent of ; the optimized set of creation operators corresponds to = + as indicated by ∆ + g in the first line of the table.

Lattice field theory results for hybrid static potentials
We compute the ground state hybrid static potential for each of the sectors Λ η = Σ − g , Σ + u , Σ − u , Π g , Π u , ∆ g , ∆ u as well as the ground state and first excited static potential for the sector Λ η = Σ + g . The latter is in the same energy region as the ground state hybrid static potentials and of particular interest, since it is expected to become degenerate with the Π g hybrid static potential in the limit of small quark-antiquark separations r. For these computations we use correlation matrices • Σ − g , Σ + u , Σ − u , Π u : 3 × 3 correlation matrices with operators as specified in Table 1, Table 2, Table 3 and  Table 5.
We solve generalized eigenvalue problems with t 0 = a and n = 0, 1, . . . (we sort the resulting eigenvalues according to λ (0) (r, t, t 0 ) > λ (1) (r, t, t 0 ) > . . .; for details concerning the generalized eigenvalue problem cf. e.g. [47] and references therein). The resulting "effective potentials" are constant with respect to t for sufficiently large t within statistical errors. The plateau values of V (0) eff;Λ η (r, t, t 0 ) correspond to the ground state potentials V Λ η (r) and we extract them by fitting a constant to V (0) eff;Λ η (r, t, t 0 ) for each r in the range t min ≤ t ≤ t max . Similarly, we determine the first excitation V Σ + g (r) from V (1) eff;Λ η (r, t, t 0 ). In principle we could determine the first excitations of the hybrid static potentials in the same way, but since statistical errors for these excitations are quite large, we decided not to include the corresponding results in this work.
We choose t min sufficiently large to guarantee a strong suppression of excited states. On the other hand, statistical errors of effective potentials are smaller at smaller t. Consequently, the statistical error on V Λ η (r) will be smaller, when using smaller t min . We have implemented the following algorithm with the intention to automatically determine t min and t max for each of the static potentials and each r in a fair way.
• t max = 9 a, the maximum t, where correlation functions have been computed.
• Fit constants V Λ η (r) to V eff;Λ η (r, t, t 0 ) for all ranges t min . . . t max with t min ≤ t min , t max ≤ t max and t max − t min ≥ 2 a. Results with χ 2 red > 1.0 are discarded, where χ 2 red denotes the uncorrelated reduced χ 2 of the corresponding fit. If all fits yield χ 2 red > 1.0, keep that one with the smallest χ 2 red and discard all others. • As final result for V Λ η (r) take the fit result corresponding to the longest plateau, i.e. with maximum t max − t min . If there are several fit results with the same maximum t max − t min , take the fit result with the smallest t min .
We have checked each of the resulting fitting ranges and corresponding plateau fits and have found that in almost all cases the algorithm decided for a reasonable range t min . . . t max . Only in very few cases (less than 5%) one gets the impression that there is a slight mismatch between the range t min . . . t max determined by the algorithm and the plateau region. In these cases we have changed t min manually by either +a or −a.
To illustrate the quality of our numerical data and the automatic determination of fitting ranges, we show effective potentials for all eight sectors Λ η = Σ + g , Σ − g , Σ + u , Σ − u , Π g , Π u , ∆ g , ∆ u for separations r/a = 2, 5, 8 in Figure 6. The fitting ranges t min . . . t max are indicated by red lines.
The resulting hybrid static potentials are shown in units of r 0 2 together with the ordinary static potential and its first excitation in Figure 7. Separations r < 2a are not shown, because such data points are known to exhibit strong lattice discretization effects. The corresponding numbers V Λ η (r), Λ η = Σ + g , Σ − g , Σ + u , Σ − u , Π g , Π u , ∆ g , ∆ u and V Σ + g (r) are collected in Table 8 for future reference. This data might be of interest for similar recent or future lattice studies as a benchmark (cf. e.g. [32]) or as input for effective field theories like pNRQCD and mass determinations of heavy hybrid mesons (cf. e.g. [7,48]).
In the following we compare our results from Figure 6, Figure 7 and Table 8 to a previous computation of hybrid static potentials [15,16,18,19,22,23,25,28]. Even though this computation was performed more than 15 years ago, the resulting potentials are frequently used in current projects (cf. e.g. [7,48]) and seem to be the most accurate lattice field theory results for hybrid static potentials, which are currently available. Unfortunately, the above references correspond to rather short publications, where the resulting potentials are shown, but details are missing, e.g. tables containing numerical values, in particular error bars, or effective mass plots. More detailed information is, however, available on the webpage [49] of one of the authors of the listed publications. Thus, the following comparative discussion is to a large extent based on the data collected at [ R / a = 2 R / a = 5 R / a = 8 Figure 6: Effective potentials V    Figure 7: The ordinary static potential V Σ + g (r)r 0 and the corresponding first excitation V Σ + g (r)r 0 as well as the hybrid static potentials V Λ η (r)r 0 , Λ η = Σ − g , Σ + u , Σ − u , Π g , Π u , ∆ g , ∆ u as functions of the separation r/r 0 , where r 0 = 0.5 fm. To allow a straightforward comparison with results from the literature, e.g. with [13,23], the vertical scale has been shifted by an additive constant such that V Σ + g (2r 0 ) = 0.
• The maximal separations provided at [49] are r = 12 a ≈ 1.44 fm (Run A) and r = 10 a ≈ 1.90 fm (Run B ). In our computation we considered separations up to r = 12 a ≈ 1.12 fm.
• We use the ordinary Wilson gauge action, i.e. lattices with the same lattice spacing in all four spacetime dimensions. At [49] anisotropic lattices are employed, where the temporal lattice spacing is smaller by a factor of 3 (Run A) and 5 (Run B ) than the spatial lattice spacing.
• We use a 24 3 ×48 lattice, i.e. the same number of lattice sites in all three spatial directions. At [49] the number of spatial lattice sites is larger in the direction of the quark antiquark separation, than in the two other directions, i.e. (18 2 ×24)×54 (Run A) and (16 2 ×20)×80 (Run B ).
• We use HYP2 smeared temporal links, when computing correlation functions (11), to reduce the self energy of the static quarks. At [49] unsmeared temporal links are used.
Comparing our effective potentials ( Figure 6) to those at [49] in a meaningful way is difficult. Our effective potentials are obtained from correlation matrices by solving generalized eigenvalue problems, while at [49] optimized correlation functions are used. Also the extraction of the potentials is done in a different way. We fit constants at larger temporal separations, where effective potentials are consistent with a plateau, while at [49] sums of two exponentials are fitted to the optimized correlation functions including also data points at small temporal separations. Nevertheless, when comparing to [49]/Run A it seems that our effective potentials start to be consistent with plateaus at smaller temporal separations (in physical units) for Λ η = Σ − u , ∆ g , ∆ u , at similar temporal separations for Λ η = Σ + u , Σ − g , Π g and at larger temporal separations for Λ η = Π u . This is most likely a consequence of different operator sets used in this work and at [49]. While we have documented our operator optimization in section 4 in detail, equivalent information for the computation from [49] does not seem to be available. Another observation is that our plateau-like regions tend to be somewhat longer (in physical units), which we attribute to the smaller self energy of the static quarks due to the use of HYP2 smeared temporal links.
The temporal separations, where effective potentials start to be consistent with plateaus, are also reflected in the statistical errors of the hybrid static potentials (cf.
and larger by a factor ≈ 2.0 for Λ η = Π u .
[49]/Run B has slightly smaller statistical errors than [49]/Run A, but the overall picture is the same.
There is no obvious discrepancy for the majority of potentials concerning their shape (cf. Figure 7 of this work and e.g. FIG. 2 in [23]). Clearly visible differences can be observed for V Πg (r) and V ∆u (r), in particular at small separations r < ∼ 0.25 fm. Our results for these potentials are somewhat lower than those from [23] and exhibit the expected approximate degeneracy with V Σ + g (r) and V Σ + u (r), respectively (for a detailed discussion of these degeneracies and their relation to gluelump masses cf. e.g. [50]). Interestingly, we have found that the resulting potentials V Πg (r) and V ∆u (r) are quite sensitive to the creation operators used in the correlation matrices. In both cases the operator S IV,2 significantly increases the ground state overlap and, thus, is essential to observe the previously mentioned and expected degeneracies at short r. We interpret this as indication that our selected sets of operators are better able to isolate the groundstate potentials for short r in the Π g and ∆ u sectors than the operators used in [23]. It should, however, be noted that hybrid static potentials at sufficiently small separations can decay to the ordinary static potential and a glueball. Thus, extracting the potential from the exponential decay of a correlation function, as done in our work as well as in [23], might give contaminated results for small r. The lightest glueball has quantum numbers J P C = 0 ++ and mass m 0 ++ ≈ 4.21/r 0 [51]. Using this mass one can read off from Figure 7 that V Πg (r) can decay for r < ∼ 0.25 fm and V ∆u (r) for r < ∼ 0.5 fm. The lowest hybrid static potentials V Πu (r) and V Σ − u (r), which are used in section 6 to estimate masses of heavy hybrid mesons, can only decay for r < ∼ 0.12 fm.
There are further lattice field theory computations of hybrid static potentials, which are interesting to discuss or to compare with: • In [27] the Σ − u and Π u hybrid static potentials were computed in pure SU(3) gauge theory using several lattice spacings a ≥ 0.16 × r 0 ≈ 0.08 fm as well as off-axis separations. The focus of the paper is on the phenomenology of static sources and gluonic excitations at short separation. The results presented for the Σ − u and Π u hybrid static potentials seem to agree with our findings.
• Very recently color field densities of static potential flux tubes in the sectors Σ + g , Σ + u and Π u have been computed in pure SU(3) gauge theory [32]. As a byproduct the potentials in these three sectors have been obtained, which seem to agree with our results.

Masses of heavy hybrid mesons in the Born-Oppenheimer approximation
In this section we parameterize the ordinary static potential V Σ + g (r) and the two lowest hybrid static potentials V Πu (r) and V Σ − u (r) computed in section 5, to estimate masses of heavy hybrid mesons with quarksQQ =cc andQQ =bb for various J P C quantum numbers. This is done in the Born-Oppenheimer approximation [34], which is a two-step procedure commonly used e.g. in molecular physics. In the first step, which is the computation of hybrid static potentials using lattice field theory (cf. sections 5 and 6.1), the gluons are the only dynamical degrees of freedom, whereas the positions of the heavy quarksQ and Q are frozen. In the second step, which is discussed in section 6.2, this constraint is relaxed by solving the Schrödinger equation for the relative coordinate of theQQ pair using the hybrid static potentials computed in the first step.
6.1 Parameterization of lattice field theory results for the ordinary Σ + g and the hybrid Π u and Σ − u static potentials For the ordinary static potential (the ground state in the Σ + g sector) a common choice, which is able to parameterize lattice data for separations r > ∼ 0.15 fm, is σ is the string tension, α is a positive constant and V 0 is a physically irrelevant shift, which contains the self energy of the static quarks and, thus, depends on the lattice spacing. For a detailed recent discussion of this parameterization cf. e.g. [52].
In this section we focus on the two lowest hybrid static potentials with quantum numbers Π u and Σ − u . Our parameterizations are based on the pNRQCD prediction which is valid for small separations r 1/Λ QCD ≈ 0.5 fm [7,50].
o (ν f ) is the Renormalon Subtracted (RS) octet potential, ν f the subtraction scale and Λ H (ν f ) a constant. Λ H (ν f ) is the same for those hybrid static potentials, which are degenerate for r → 0, i.e. for V Πu (r) and V Σ − u (r). Eq. (17) suggests to use as fit function for small separations r, where both A 1 und A 2 are the same for Λ η = Π u and Λ η = Σ − u , while A 3 is different. We are interested in parameterizations of our lattice field theory results over the whole available range of separations r, i.e. up to r = 12a ≈ 1.1 fm. While eq. (18) is suited to parameterize V Πu (r) up to r = 12a, which is beyond the region of validity of the pNRQCD prediction, this is not the case for V Σ − u (r). Therefore, we use for V Σ − u (r) a fit function with additional degrees of freedom, which reduces to eq. (18) in the limit of small r. In [6] it was suggested to use for the difference of the two lowest hybrid static potentials. While this is a reasonable crude description of this difference, it is not sufficient to parameterize our precise lattice field theory results from section 5 in a consistent way, i.e. with reduced χ 2 < ∼ 1. Therefore, we extend eq. (19) by introducing another parameter B 2 , Altogether we parameterize V Πu (r) and V Σ − u (r) by where A 1 , A 2 , A 3 , B 1 , B 2 and B 3 are fit parameters. Note that (22) still reduces to the NRQCD prediction (18) in the limit of small r, i.e. the Π u and Σ − u hybrid static potentials are dominated by the same repulsive "octet-like" 1/r term and become degenerate (A 3 in eq. (18) is then equivalent to A 3 + B 1 in eq. (22)).
To determine the unknown parameters in eq. (16), eq. (21) and eq. (22), we perform uncorrelated χ 2 minimizing fits to our lattice data points in the region r min ≤ r ≤ r max .
For the Σ + g static potential we use r min = 3a (for r < 3a lattice field theory results obtained with HYP smearing typically exhibit non-negligible discretization errors) and r max = 12a. From a 3-parameter fit we obtain V 0 a = 0.1515(13) α = 0.2626(23) σa 2 = 0.04749 (17) , where χ 2 red = 0.80 indicates a consistent fit. The parameterization (16) with the parameters (23) is shown in Figure 8 (left) together with our corresponding lattice field theory results.
From these results we can determine r 0 , which is the typical length scale in lattice gauge theory and quite often used to set the scale. r 0 is defined via V (r 0 )r 2 0 = 1.65 (24) and we obtain r 0 a = 1.65 − α σa 2 1/2 = 5.405 (6).  (23)). (right) Π u and Σ − u hybrid static potentials (eq. (21) and (22) with parameters (26)). Vertical dotted lines indicate the data points, which have been considered in the χ 2 minimizing fits.
When identifying r 0 with 0.5 fm, which is a common choice in lattice gauge theory 3 , we find a = 0.0925(1) fm (cf. also section 3). This value is consistent with an independent simulation and scale setting analysis quoted in [54].
For the Π u (r) and Σ − u hybrid static potentials we use r min = 2 a and r max = 12 a. From a single 6-parameter fit we obtain where χ 2 red = 1.48 indicates a reasonable fit. The parameterizations (21) and (22) with the parameters (26) are shown in Figure 8 (right) together with our corresponding lattice field theory results.

Prediction of masses of heavy hybrid mesons
The Born-Oppenheimer approximation for heavy hybrid mesons was pioneered in [12,18,55] and is explained in detail in [6]. One has to solve the radial Schrödinger equation 2µr 2 + V Λ η (r) u Λ η ;L,n (r) = E Λ η ;L,n u Λ η ;L,n (r), (27) where r is the separation of the heavyQQ pair, V Λ η (r) is one of the static potential parameterizations from section 6.1, eq. (16), (21) or (22), and µ = mQm Q /(mQ + m Q ) is the reduced mass of the QQ pair. We use m Q = m c = 1628 MeV and m Q = m b = 4977 MeV from quark models [56]. The wave function of the relative coordinate of theQQ pair is ψ Λ η ;L,n,m L (r, ϑ, ϕ) = (u Λ η ;L,n (r)/r)Y L,m L (ϑ, ϕ). L ∈ {Λ, Λ + 1, . . .} is the quantum number corresponding to the operator L, the sum of all angular momenta excluding the heavy quark spins S, i.e. J = L + S, where J is the total angular momentum of the meson. In the limit r → 0 the gluon field configuration of a hybrid static potential is identical to that of a gluelump, where J Λ η is the gluon spin of this gluelump. 50]. The derivation of the Schrödinger equation (27) is based on the following approximations (cf. also [6]): • In the adiabatic approximation the gluon field is assumed to be in a stationary state in the presence of the heavyQQ pair, i.e. the gluon field configuration is one of the hybrid static potentials computed in section 5 labeled by quantum numbers Λ η . Errors are proportional to Λ QCD /m Q , i.e. the adiabatic approximation is suited for heavy quarks . It is consistent with using static potentials, where also 1/m Q corrections are neglected.
• The Schrödinger equation to determine masses of a heavy hybrid mesons with quantum numbers J P is a multi-channel equation including all hybrid static potentials V Λ η (r) consistent with J P (cf. [7] for a detailed derivation of a coupled channel Schrödinger equation). In the single channel approximation only a single component of this multi-channel Schrödinger equation is considered and couplings to other channels are ignored. The single channel approximation is good, if the resulting wave function is small for separations r, where the used hybrid static potential has avoided crossings with the other hybrid static potentials.
• Finally the gluon spin is approximated by the gluon spin of a gluelump, which is a good approximation for small separations r, where the system resembles a gluelump. Consequently, the approximation is good for resulting wave functions ψ(r, ϑ, ϕ), which are localized near r = 0.
The Schrödinger equation (27) can be solved numerically with standard techniques. We employ a 4th order Runge-Kutta shooting method combined with Newton's method for root finding. Note that the resulting energies E Λ η ;L,n contain the self-energies of the static quarks, which depend on the lattice spacing. To predict heavy hybrid meson masses, one has to eliminate these self-energies, which we do by subtracting E Λ η =Σ + g ;n=1,L=0 , the lowest energy from the ordinary static potential computed within the same setup. This energy E Λ η =Σ + g ;n=1,L=0 corresponds for QQ =cc to the η c (1S) and J/Ψ(1S) meson, which are degenerate in the static limit, and similarly forQQ =bb to the η b (1S) and Υ(1S) meson. Heavy hybrid meson masses are then given by as discussed in [6]. Our predicted heavy hybrid meson masses are collected in Table 9 and summarized in a graphical way in Figure 9. The errors are statistical uncertainties, which have been obtained via an elaborate jackknife analysis (cf. section 3). In Figure 10 we also show the probability density for the separation r, which is |u Λ η ;L,n (r)| 2 , for the Π u and the Σ − u hybrid static potentials and forQQ =cc andQQ =bb.  There are also systematic errors, which are difficult to quantify. The derivation of the Schrödinger equation is based on several approximations, as discussed above, most notably the neglect of 1/m Q corrections. In principle such corrections can be computed for hybrid static potentials, but this is expected to be very challenging, since it turned out to be difficult already for the ordinary static potential with quantum numbers Λ η = Σ + g [54,[58][59][60][61]. At the moment we crudely estimate the magnitude of this error by half the experimental mass differences of the non-exotic    We expect the corresponding corrections to be rather small based on our experience from previous projects concerned with the static potential, in particular [62]. Moreover, dynamical light quarks have been neglected. Note, however, in [21] the Π u hybrid static potential was computed in full QCD, i.e. with dynamical quarks corresponding to pion masses m π ≈ 500 . . . 1200 MeV.
No statistically significant differences were found, when comparing to an equivalent computation in pure SU(3) gauge theory. We have also investigated the dependence of the predicted heavy hybrid meson masses on the c and b quark masses used in the Schrödinger equation. This dependence has been found to be very weak. For example, when using quark masses m c = 1480 MeV and m b = 4890 MeV as in [6], which differ from our choice m c = 1628 MeV and m b = 4977 MeV by 2% and 9%, respectively, the mass differences E Λ η ;L,n − E Λ η =Σ + g ;n=1,L=0 (cf. eq. (28)) change for the majority of states only on the per mille level. Finally, these mass differences also depend on the scale setting procedure, i.e. on the value of r 0 used for scale setting. The corresponding relative systematic error on E Λ η ;L,n − E Λ η =Σ + g ;n=1,L=0 is roughly the same as the relative uncertainty on r 0 , which is around 4%.

Conclusions
We have computed hybrid static potentials V Λ η (r) for Λ η = Σ − g , Σ + u , Σ − u , Π g , Π u , ∆ g , ∆ u in SU(3) lattice gauge theory. Compared to results from the literature we use a rather fine lattice spacing and statistical errors are quite small. In contrast to the majority of existing publications, technical aspects of the computation and the analysis are discussed in detail. This offers the possibility of direct and meaningful comparison with similar computations and of methodological improvement. Moreover, we provide the numerical values of the discrete lattice data points for all computed hybrid static potentials V Λ η (r), to allow straightforward usability of our results in future effective field theory or phenomenological work by other authors.
We also estimate masses of heavy hybrid mesons in the Born-Oppenheimer approximation, where we follow closely the approach discussed in [6]. The resulting spectra are, however, rather crude estimates, mainly because 1/m Q corrections, e.g. from the quark spins, are neglected at the moment.
An essential point of this work is the extensive optimization of hybrid static potential creation operators. We plan to use these optimized operators in follow-up projects concerned with the computation of 3-point functions. Such 3-point functions might allow to drastically reduce systematic errors in the above mentioned prediction of heavy hybrid meson masses, e.g. by computing quark spin corrections or by studying possible decays to ordinary quarkonium states and glueballs. Furthermore, 3-point functions will provide interesting insights concerning the gluon distribution inside heavy hybrid mesons. We have presented first corresponding results at recent conferences [63,64].

A Angular momentum of the trial states
In this appendix we show that the trial state (1), has definite total angular momentum Λ with respect to the z axis.
Any state can be rotated by an angle α around the z axis using the total angular momentum operator J z . For example the rotated trial state |Ψ hybrid S;Λ is given by R(α) |Ψ hybrid S;Λ = exp(−iJ z α) |Ψ hybrid S;Λ .
This proves that the trial state |Ψ hybrid S;Λ is an eigenstate of the total angular momentum operator J z with eigenvalue Λ.