Hybrid static potential flux tubes from SU(2) and SU(3) lattice gauge theory

We compute chromoelectric and chromomagnetic flux densities for hybrid static potentials in SU(2) and SU(3) lattice gauge theory. In addition to the ordinary static potential with quantum numbers $\Lambda_\eta^\epsilon = \Sigma_g^+$, we present numerical results for seven hybrid static potentials corresponding to $\Lambda_\eta^{(\epsilon)} = \Sigma_u^+, \Sigma_g^-, \Sigma_u^-, \Pi_g, \Pi_u, \Delta_g, \Delta_u$, where the flux densities of five of them are studied for the first time in this work. We observe hybrid static potential flux tubes, which are significantly different from that of the ordinary static potential. They are reminiscent of vibrating strings, with localized peaks in the flux densities that can be interpreted as valence gluons.


Introduction
The majority of mesons, i.e. hadrons with integer total angular momentum, are quark-antiquark pairs. It is, however, expected that some mesons, so-called exotic mesons, have a more complicated composition in terms of quarks and gluons. An important example are hybrid mesons, where gluons contribute to the quantum numbers J P C (J: total angular momentum; P : parity; C: charge conjugation) in a non-trivial way. In the quark model, where mesons are quark-antiquark pairs, quantum numbers are restricted to P = (−1) L+1 and C = (−1) L+S with spin S = 0, 1 and orbital angular momentum L = 0, 1, 2, .... Thus, mesons with J P C = 0 +− , 0 −− , 1 −+ , 2 +− , . . ., which are not allowed in the quark model, are obvious candidates for exotic mesons like hybrids. Moreover, a higher density of states than obtained by the quark model might also indicate hybrid mesons.
Experimentally observed examples, which could be hybrid mesons, are the J P C = 1 −+ states π 1 (1400) and π 1 (1600). They could, however, also be tetraquarks, i.e. two quarks and two antiquarks without excited glue. For heavy-heavy mesons the situation seems to be even less clear. There are several exotic candidates, which could be hybrid mesons, but for none of them such an interpretation seems to be likely (see. e.g. the experimental review of exotic hadrons [1] and the discussion in section VII.A of ref. [2]). Thus, the search for gluonic excitations is an important part of the research program of current and future experiments, e.g. the GlueX experiment at the JLab accelerator or the PANDA experiment at the FAIR accelerator.
Also on the theoretical side there are many open questions concerning hybrid mesons (see e.g. the theoretical reviews [3][4][5][6]). They are difficult to study, because in QCD total angular momentum J and parity P are not separately conserved for gluons on the one hand and for the quarkantiquark pair on the other hand. Only the overall J P are quantum numbers. For heavy hybrid mesons, e.g. composed of a b and ab quark and gluons, a simplification and good approximation is to study the static limit. In that limit the quark positions are frozen, which allows to separate the treatment of gluons and quarks.
In this work we use SU (2) and SU(3) lattice gauge theory to study heavy hybrid mesons in the static limit. Since quite some time hybrid static potentials haven been computed by various groups, mainly with the intention to compute masses of heavy hybrid mesons using the Born-Oppenheimer approximation (see refs.  and the recent review article [32]). We focus on a different problem, the computation of the gluonic flux densities for hybrid potential states, i.e. the structure of the flux tube, for several hybrid channels. While such flux tubes have been studied for the ordinary static potential using lattice gauge theory for quite some time (see refs. [33][34][35][36][37][38][39][40][41][42][43][44][45][46][47][48]), this is a rather new direction for hybrid static potentials, where first results appeared only recently [49][50][51][52]. In this paper we substantially extend existing work by performing computations for seven hybrid static potential sectors characterized by quantum numbers Λ ( ) Five of these sectors are studied for the first time, where preliminary results have been presented at a recent conference [52].
The paper is structured as follows. In section 2 we discuss theoretical basics, including quantum numbers for hybrid static potentials, the construction of corresponding trial states and the computation of chromoelectric and chromomagnetic flux densities. Section 3 contains a brief summary of our lattice setup. In section 4 we present our numerical results. We start with a discussion of systematic errors and symmetries, before showing and interpreting our main results, the chromoelectric and chromomagnetic flux densities for the seven hybrid static potential sectors Λ ( ) η = Σ + u , Σ − g , Σ − u , Π g , Π u , ∆ g , ∆ u . In section 5 we conclude with a short summary and an outlook.
2 Hybrid static potentials and flux tubes 2.1 Hybrid static potential quantum numbers and trial states A hybrid static potential is the potential of a static quark Q and a static antiquarkQ, where the gluons form non-trivial structures and, thus, contribute to the quantum numbers. Such potentials can be computed from temporal correlation functions of hybrid static potential trial states. After replacing the static quark operators by corresponding propagators, these correlation functions are similar to Wilson loops. Instead of straight spatial Wilson lines there are, however, parallel transporters with more complicated spatial structures. For a detailed discussion of such correlation functions see e.g. our recent work [31], where we have carried out a precision computation of hybrid static potentials using SU(3) lattice gauge theory.
Hybrid static potentials can be characterized by the following quantum numbers: • Λ = 0, 1, 2, . . ., the absolute value of the total angular momentum with respect to the QQ separation axis, 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 (an axis perpendicular to the QQ separation axis).
In [31] we discussed hybrid static potential creation operators and trial states both in the continuum and in lattice gauge theory in detail and performed a comprehensive optimization of these operators in SU(3) lattice gauge theory. In this paper we use the information obtained during this optimization to define suitable hybrid static potential creation operators both for SU(2) and SU(3) lattice gauge theory. These operators are important building blocks of the 2-point and 3-point functions, which need to be computed for the investigation of hybrid static potential flux tubes (see section 2.2).
U (−r/2, r 1 )S(r 1 , r 2 )U (r 2 , +r/2) has been optimized in ref. [31], such that the overlap of |Ψ Λ η (r) to the ground state in the Λ η sector is rather large. In contrast to ref. [31] we use in this work only a single operator S for each Λ η sector, not a linear combination obtained by a variational analysis. This reduces computation time to a feasible level, while the suppression of excited states in the 2-point and 3-point functions is still sufficiently strong. U (−r/2, r 1 )S(r 1 , r 2 )U (r 2 , +r/2) is different for each of the eight Λ η sectors as well as for the two QQ separations r = 6 a, 10 a considered. For the Σ + g sector, i.e. the ordinary static potential, it is just a straight line, while for Table 1. For each Λ η sector we take that operator from the set of three to four operators we optimized in ref. [31], which minimizes the effective potential at t = a. Thus Table 1 contains that part of the information shown in Table 1 to Table 7 of ref. [31], which is relevant in the context of this work.

Expectation values of squared field strength components
The energy density of the gluon field is where E a j (x) and B a j (x) denote the components of the chromoelectric and chromomagnetic field strengths with spatial indices j and color indices a (a = 1, . . . , 3 for gauge group SU(2) and a = 1, . . . , 8 for gauge group SU(3)). The main goal of this work is to compute the expectation values of the six gauge invariant terms F 2 or F a j = B a j ) contributing to eq. (3) for states with a static quark-antiquark pair and quantum numbers Λ η . These chromoelectric and chromomagnetic flux densities provide information about the shapes of hybrid static potential flux tubes and the gluonic energy distributions inside heavyheavy hybrid mesons.
To compute the flux densities, we need the following quantities: The notation in the caption of each of the tables follows ref. [31]. Note that, even though the Π η and the ∆ η hybrid potentials are degenerate with respect to , the construction of creation operators via eq. (2) is not independent of . One can obtain an optimized Π − η operator from an optimized Π + η operator by applying a π/2 rotation with respect to the z axis. For ∆ η operators there is no analogous simple prescription. Therefore, we provide four different optimized ∆ η operators.
• Vacuum expectation values: These quantities can be combined to expressions for the expectation values of E 2 j (x) and B 2 j (x) for static potential states with quantum numbers Λ η with the vacuum expectation value subtracted: (see also ref. [33], where this quantity was first defined and used to study flux densities for the ordinary static potential with quantum numbers Λ η = Σ + g ). The right hand side of eq. (7) can be evaluated using Euclidean lattice gauge theory path integrals, where . . . U denotes the path integral expectation value and W (r, t 2 , t 0 ) = = Tr a S;Λ η (−r/2, +r/2; t 0 )U (+r/2; t 0 , t 2 ) a S;Λ η (−r/2, +r/2; (for Λ η = Σ + g , i.e. the ordinary static potential,W (r, t 2 , t 0 ) is the standard Wilson loop). P µν is a symmetrized plaquette in the µ-ν plane, also denoted as clover leaf. Eqs. (8) and (9) as well as the clover leaf are illustrated in Figure 1. Figure 1: Graphical illustration of lattice field theory quantities needed to determine ∆E 2 j,Λ η (r; x) and ∆B 2 j,Λ η (r; x). Red spheres, black dots, black solid lines and black dashed lines represent static (anti)quarks, lattice sites, gauge links and operators a S;Λ η , respectively. (a)W (r, t 2 , t 0 ) · P 0x (x, t 1 ) for r = 3 a, t 2 − t 0 = 4 a, t 1 = (t 2 − t 0 )/2 (grey dashed lines parallel to the t axis and the z axis are drawn to guide the eye). (b) The corresponding Wilson loop W (r, t 2 , t 0 ). (c) The symmetrized plaquette P µν .
One can show that under rotations around the z axis with angle α the states |0 Λ η (r) transform according to (see appendix A), while the field strength components transform as R(α) denotes the corresponding standard 3 × 3 rotation matrix, i.e.
and we have defined x −α = R(−α)x. Now we consider the rotated flux densities 0 Λ ± η (r)|R † z (α)F 2 j (x)R z (α)|0 Λ ± η (r) − Ω|F 2 j |Ω and rewrite them in two different ways, using first eq. (11), then eq. (12), with the shorthand notation c α = cos(α) and s α = sin(α), and where the index j on the right hand side indicates the j-th component of [. . .]. Equating eqs. (14) and (15) relates the flux densities ∆F 2 j,Λ η (r; x) and ∆F 2 j,Λ η (r; x −α ), i.e. yields their transformation law with respect to rotations around the z axis 1 . Clearly, one cannot expect that the flux densities ∆F 2 j,Λ + η (r; x) and ∆F 2 j,Λ − η (r; x) are invariant under rotations, nor that they appear to be identical, in particular not for Λ ≥ 1, even though the corresponding potentials are degenerate. Numerical computations confirm that these flux densities are not invariant under rotations and that they are different from each other (see the discussion in section 4.2 and the example plots in Figure 5 and Figure 6).
Instead of quantum numbers Λ η one can also use quantum numbers λ η to label hybrid static potential states with Λ ≥ 1, where λ = . . . , −2, −1, +1, +2, . . . is the total angular momentum with respect to the z axis, i.e. Λ = |λ|. Of course, there are again the same pairs of degenerate potentials, i.e. V +λη (r) = V −λη (r) = V Λ + η (r) = V Λ − η (r). In this case the behavior of the corresponding states and flux densities under rotations is different, and eq. (14) simplifies, while eq. (15) remains essentially unchanged (one just has to replace Λ η by λ η ). Consequently, the transformation law with respect to rotations around the z axis and the angular dependence of ∆F 2 j,Λ η (r; x) and ∆F 2 j,λη (r; x) is different, even though the corresponding hybrid static potentials are identical.
To eliminate this somewhat arbitrary angular dependence, which is a consequence of (when using quantum numbers Λ η ) or the sign of λ (when using quantum numbers λ η ), but not related to Λ = |λ| or η (Λ and η fully characterize hybrid static potentials V Λη (r) for Λ ≥ 1), we define for Λ ≥ 1 This quantity represents the average over an ensemble of states with fixed Λ and η, but arbitrary . After the last equality the projector to the corresponding 2-dimensional space of states has been introduced. This projector shows explicitly that ∆F 2 j,Λη (r; x) is independent of the basis used for that 2-dimensional space, i.e. independent of whether we use use or the sign of λ.
The transformation law with respect to rotations around the z axis for ∆F 2 j,Λη (r; x) is where the left hand side can be obtained by combining eqs. (17), (18) and (19) and the right hand side is essentially eq. (15). To simplify this even further it is convenient to define as e.g. also done in a similar way in [51]. This quantity as well as ∆F 2 z,Λη (r; x) are invariant under rotations around the z axis, i.e.

Lattice setup
The computations presented in this work have been performed using SU(2) and SU(3) lattice gauge theory. The corresponding gauge link configurations have been generated with the standard Wilson gauge action (see textbooks on lattice field theory, e.g. ref. [53]). Since we are considering purely gluonic observables, we expect that there is little difference between our pure gauge theory results and corresponding results in full QCD. This expectation is supported by ref. [20], where hybrid static potentials were computed both in pure gauge theory and QCD and no statistically significant differences were observed.
For the SU(2) simulations we have used a standard heatbath algorithm. To eliminate correlations in Monte Carlo time, the gauge link configurations are separated by 100 heatbath sweeps. For the SU(3) simulations we have used the Chroma QCD library [54]. There, the gauge link configurations are separated by 20 update sweeps, where each update sweep comprises a heatbath and four over-relaxation steps. Details of our simulated ensembles are collected in Table 2, including the gauge coupling β, the lattice extent (L/a) 3 × T /a and the number of gauge link configurations used for the flux tube computations. We also list the lattice spacing a in fm, which is obtained by identifying r 0 with 0.5 fm (see refs. [31,55]). For the majority of computations for gauge group SU (2)   To improve the signal quality, standard smearing techniques are applied, when samplingW appearing in eqs. (8) and (9) and defined in eq. (10): • Spatial gauge links, i.e. links in a S;Λ η (−r/2, +r/2; t 0 ) and a S;Λ η (−r/2, +r/2; t 2 ) (defined in eq. (2)), are APE smeared gauge links (for detailed equations see e.g. ref. [56]), where the parameters α APE = 0.5 and N APE = 20 have been optimized in ref. [31] to generate large overlaps with the ground states |0 Λ η (r) 2 . This allows to identify plateaus in ∆F 2 eff;j,Λ η (r, t 2 , t 0 ; x, t 1 ) at smaller temporal separations t 2 − t 1 and t 1 − t 0 (see eq. (7)).
All statistical errors shown and quoted throughout this paper are determined via jackknife. We perform a suitable binning of gauge link configurations to exclude statistical correlations in Monte Carlo time. We have determined ∆F 2 j,Λ η (r; x) by fitting a constant to the lattice result for ∆F 2 eff;j,Λ η (r, t 2 , t 0 ; x, t 1 ) at sufficiently large t 2 −t 1 and t 1 −t 0 , where the data points are consistent with a plateau (see eqs. (7), (8) and (9)). For even (t 2 − t 0 )/a we use t 1 = (t 0 + t 2 )/2, while for odd (t 2 − t 0 )/a we use t 1 = (t 0 + t 2 + a)/2, i.e. equal or similar values for t 2 − t 1 and t 1 − t 0 . Example plots of ∆F 2 eff;j,Λ η (r, t 2 , t 0 ; x, t 1 ) as a function of t 2 − t 0 for gauge group SU(2), all investigated Λ η sectors, quark-antiquark separation r = 10 a and x = 0 are shown in Figure 2. We have performed an uncorrelated χ 2 -minimizing fit of a constant corresponding to ∆F 2 j,Λ η (r; x) in the region t min ≤ t 2 − t 0 ≤ t max . Since the statistical errors of ∆F 2 eff;j,Λ η (r, t 2 , t 0 ; x, t 1 ) rapidly increase with increasing t 2 − t 0 , the results are almost independent of t max . We have taken the largest t max , where the signal is not lost in noise. t min has been chosen such that χ 2 /dof < ∼ 1 for the majority of fits. This results in t min ≈ 3 a . . . 4 a and t max ≈ 5 a . . . 8 a for hybrid static . . 6 a and t max = 10 a for the ordinary static potential with Λ η = Σ + g . As an additional check that t min is chosen sufficiently large, i.e. that excited states are strongly suppressed, we have repeated the computation of ∆F 2 j,Λ η (r; x) for gauge group SU(2), Λ η = Π u , r = 6 a and x = (x, 0, 0) using a creation operator S (see eq. (2)), which has a structure significantly different from that shown in Table 1, namely S I,1 as defined in ref. [31], Figure 2. Within statistical errors we find identical flux densities ∆F 2 j,Λ η (r; x), which we interpret as confirmation, that we indeed measure the flux densities of the ground states in the Λ η sectors and not flux densities, which depend on the creation operators and trial states we are using.

Discretization errors and smearing
Until now we have performed computations only at a single value of the lattice spacing a. Therefore, we are not yet able to study the continuum limit. Strong discretization errors are expected, when either r = |r Q − rQ|, |x − r Q | or |x − rQ| is small, where r Q = (0, 0, +r/2) and rQ = (0, 0, −r/2) are the positions of the static charges and x is the spatial argument of the flux density ∆F 2 j,Λ η (r; x). These discretization errors are expected to be even more pronounced, when using HYP2 smeared temporal links inW (see eq. (10)), which can be interpreted as increasing the radii of the static charges. We, therefore, compare results for ∆F 2 j,Λ η (r; x) obtained with and without HYP2-smeared temporal links.
In Figure 3 we show results for Λ η = Σ + g , Σ − u and QQ separation r = 10 a on the separation axis, x = (0, 0, z). For |x−r Q | ≤ a or |x−rQ| ≤ a drastic discrepancies between unsmeared and HYP2smeared results can be observed, while for |x−r Q |, |x−rQ| ≥ 4 a and ∆F 2 j,Λ η (r; x) = ∆E 2 z,Λ η (r; x) as well as for |x − r Q |, |x − rQ| ≥ 3 a and all other field strength components there is agreement within statistical errors. When using HYP2-smeared temporal links, the pronounced peaks at the positions of the charges, which are present in the unsmeared results, are essentially gone. This is expected and can be observed in a qualitatively similar way also in much simpler theories, for example in classical electrodynamics, when smearing the charge density of a point charge. Analogous plots for other Λ η sectors look very similar and are not shown. Therefore, for the computations of ∆F 2 j,Λ η (r; x) in a plane containing the separation axis (see section 4.3) we do not use HYP2-smearing. Note, however, that even unsmeared results within a radius of ≈ 2 a around either of the two static charges will exhibit sizable discretization errors and should be considered as crude estimates only. In other words, instead of the poles related to the infinite self energy of the static charges, ∆F 2 j,Λ η (r; x) will exhibit pronounced but finite peaks.
We also study the effect of HYP2-smearing on ∆F 2 j,Λ η (r; x) in the mediator plane defined by z = 0 for various QQ separations r. For r ≥ 6 a we find agreement within statistical errors for all Λ η sectors and all field strength components with exception of ∆E 2 z,Λ η (r; x = (x, y, 0)). For ∆E 2 z,Λ η (r; x = (x, y, 0)) there is agreement for r ≥ 10 a for Λ η = Σ + g and for r ≥ 8 a for all other Λ η sectors. Example plots for Λ η = Σ + g , Σ − u and x = (0, 0, 0) are shown in Figure 4. Therefore, for computations of ∆F 2 j,Λ η in the mediator plane, which we show for r = 10 a in section 4.3, we use HYP2-smearing, which reduces statistical errors significantly.

Finite volume corrections
Finite volume corrections are rather mild for static potentials, when r < L/2, where L is the spatial lattice extent. In particular for pure gauge theory, where the lightest particle (the J PC = 0 ++ glueball) is very heavy, finite volume corrections should be almost negligible. A comparison of flux densities ∆F 2 j,Λ η (r; x) for gauge group SU(2) on the two gauge link ensembles with (L/a) 3 × T /a = 18 4 and (L/a) 3 × T /a = 24 4 (see Table 2) supports this expectation.

Angular dependence and symmetrization of hybrid static potential flux densities
In section 2.2 we have discussed, how hybrid static potential flux densities transform under rotations around the z axis. On a hypercubic lattice the relevant eqs. (14) and (15) are exact only for cubic rotations, i.e. for rotations with angle α, which is a multiple of π/2. For α = ±π/2 they become for even Λ, i.e. Λ = Σ and Λ = ∆, and for odd Λ, i.e. Λ = Π, We have verified our numerical computation of flux densities using these equations, i.e. we have checked that all our results are consistent with these equations within statistical errors. In a second step we have used these equations to reduce the statistical errors of our results by averaging related flux densities.
In section 2.2 we have also discussed that hybrid static potential flux densities ∆F 2 and ∆F 2 j,Λ − η (r; x) with Λ ≥ 1 are not expected to be identical, even though the corresponding potentials are degenerate (see eqs. (14) and (15)). This expectation is confirmed by the plots in the upper row of Figure 5 and Figure 6, where we show two examples, the flux densities ∆F 2 j,Λ η (r; x) for Λ η = Π + u , Π − u and for Λ η = ∆ + g , ∆ − g in the mediator plane z = 0.
In the plots at the bottom of Figure 5 and Figure 6 we show the flux densities ∆F 2 j,Λη (r; x) defined in eq. (18), again for Λ η = Π u and for Λ η = ∆ g . As discussed in section 2.2 these are ensemble averages over states with fixed Λ and η, but indefinite . Note that ∆F 2 z,Λη (r; x) is invariant under cubic rotations, while ∆F 2 x,Λη (r; x) and ∆F 2 y,Λη (r; x), even though quite similar, are related by rotations with angle α = ±π/2 (see eq. (20)). Averaging ∆F 2 x,Λη (r; x) and ∆F 2 y,Λη (r; x) according to eq. (21) would lead to another quantity invariant under cubic rotations. From now on we always show the flux densities ∆F 2 j,Λη (r; x) for Λ ≥ 1, i.e. not anymore ∆F 2 j,Λ η (r; x).  x) and ∆F 2 j,∆g (r; x) in the mediator plane (z = 0) for gauge group SU(2) and QQ separation r = 10 a.

Hybrid static potential flux densities for all Λ η sectors
In this section we show and discuss the main numerical results of this work, the flux densities ∆F 2 j,Λ ( ) η (r; x), j = x, y, z, ⊥ for the eight sectors Λ ( ) We decided to perform computations for two QQ separations, r = 6 a and r = 10 a. This allows to compare results for two significantly different r, i.e. to see, how the shapes of the hybrid static potential flux tubes change, when the quark and the antiquark are pulled apart. We did not study separations r < 6 a, because for such small separations flux densities exhibit sizable discretization errors in the region between the two charges (see the discussion in section 4.1.2). Since the signal for a Wilson loop decays exponentially with its area, we also refrained from performing computations for r > 10 a, which are very costly in terms of CPU time.
Since the resulting flux densities in the mediator plane for r = 6 a and r = 10 a are very similar, we only present them for r = 10 a. In Figure 7 these flux densities ∆F 2 j,Λ ( ) η (r; x = (x, y, 0)), j = x, y, z are shown as 2D color maps. In the upper panel of Figure 8 we present similar results, the rotationally invariant ∆F 2 j,Λ ( ) η (r; x = (x, 0, 0)), j =⊥, z along the x axis, i.e. in the mediator plane as a function of the radial coordinate. In contrast to the 2D color maps, these 1D curves allow to also show statistical errors and, thus, provide information about the precision of our numerical results. In the lower panel of Figure 8 we present in the same style differences of hybrid static potential flux densities to those of the ordinary static potential, i.e. ∆F The flux densities of the ordinary static potential form a cigar-shaped flux tube with strong positive contributions to the energy density from the chromoelectric and smaller negative contributions from the chromomagnetic field strength components. The maxima are on the QQ separation axis, i.e. at x = y = 0. While this is known from previous lattice gauge theory investigations of the ordinary static potential (see e.g. ref. [39]), the corresponding flux densities for hybrid static potentials show a variety of different and interesting structures. For example chromomagnetic flux densities of hybrid static potentials are typically larger close to the center of the flux tube than those of the ordinary static potential, as can be seen in Figure 8, lower panel. Hybrid static potential flux tubes are also wider, i.e. have a larger extension in x and y direction (see e.g. Figure 8, Figure 9 and Figure 10). Another interesting difference is that some hybrid static potentials show a clear reduction of the chromoelectric flux densities close to the center, while the chromomagnetic flux densities exhibit peaks (most prominently for Λ η = Π u , ∆ g , but to some extent also for Λ η = Σ − u ). For other sectors, Λ ( ) η = Σ + u , Π g , ∆ u , the opposite is the case, i.e. there is a positive localized peak at the center for the chromoelectric flux densities and a corresponding negative contribution of the chromomagnetic flux densities. These peaks in either the chromoelectric or chromomagnetic flux densities can be interpreted as "valence gluons" gen- and QQ separation r = 10 a.     Figure 9: Flux densities ∆F 2 j,Λ η (r; x = (x, 0, z)), j = x, y, z in the separation plane for gauge group SU(2) and sectors Λ η = Σ + g , Σ + u , Σ − g , Σ − u . (left) QQ separation r = 6 a. (right) QQ separation r = 10 a.  Figure 10: Flux densities ∆F 2 j,Λη (r; x = (x, 0, z)), j = x, y, z in the separation plane for gauge group SU(2) and sectors Λ η = Π g , Π u , ∆ g , ∆ u . (left) QQ separation r = 6 a. (right) QQ separation r = 10 a. erating the hybrid quantum numbers, as discussed in models and phenomenological descriptions of hybrid mesons. The positive or negative peaks are surrounded by spherical shells, where flux densities are smaller or larger, respectively (see in particular the 2D color maps in Figure 9 and Figure 10, where these shells are visible as rings). These structures remind of and might indicate vibrating strings, which have either nodes or maxima at z = 0. Moreover, the transverse extent of the structures formed by the chromoelectric or chromomagnetic flux densities as almost the same for QQ separation r = 6 a and r = 10 a, which is consistent with a string interpretation.
It is also interesting to compare the resulting flux densities to the gluonic excitation operators for hybrid static potentials at leading order in the multipole expansion of pNRQCD (see e.g. refs. [2,60]). Similar operators were also used in lattice gauge theory computations of hybrid static potentials as local insertions in Wilson loops (see e.g. ref. [28]), but numerically it turned out that they generate less ground state overlap than optimized non-local operators (like those discussed in ref. [31] and in section 2.1 of this work) and are, thus, less suited for computations in lattice gauge theory. The leading order gluonic excitation operators of pNRQCD are listed in Table 3, where the QQ separation axis is again the z axis. For certain Λ ( ) η sectors the flux densities we have obtained by our lattice computation closely resemble the pNRQCD operators. For example in the lower panel of Figure 8 one can clearly see that the chromomagnetic flux densities for Π u and ∆ g are significantly larger than for the ordinary static potential Σ + g , in particular the x and y components. The corresponding pNRQCD operators include B x and B y as well as D x B x − D y B y and D x B y + D y B x . Similarly, for Σ − u the z component of the chromomagnetic flux density is rather large, where one of the corresponding pNRQCD operators is B z . Further interesting structures are the double peaks in the the chromoelectric flux densities for Σ − u and Π u as shown in the upper panel of Figure 8. The pNRQCD operators for these sectors contain derivatives in x direction of the corresponding chromoelectric field operators, D x E y − D y E x and D x E z − D z E x , respectively. Again this is consistent, because from lattice gauge theory it is known that such derivative operators generate nodes in the corresponding wave functions. Table 3: Gluonic excitation operators at leading order in the multipole expansion of pNRQCD (the QQ separation axis is the z axis; D j denotes the covariant derivative; see e.g. refs. [2,60]).
Finally we compare and discuss our results in the context of a recent and similar lattice computation of hybrid static potential flux densities [51]. There the flux densities for two hybrid sectors Λ ( ) η = Σ + u , Π u were computed for gauge group SU(3), for Λ η = Π u not only for the ground state, but also for the first excitation. We have computed the flux densities for the ground states of the seven hybrid sectors Λ ( ) for gauge groups SU(2) as well as SU (3).
Lattice spacings, spacetime volumes and QQ separations are similar in both works. Compared to ref. [51] our presentation of results is different in the following aspects: • We show flux densities ∆F 2 j,Λ ( ) η (r; x) for j = x, y separately, while in [51] only the average of the x and the y component is shown, i.e. ∆F 2 ⊥,Λ ( ) η (r; x) (cf. also eq. (21)).
• In contrast to ref. [51] we do not show the flux densities on the QQ separation axis x = y = 0 as curves, because several of the hybrid static potentials have small flux densities on the separation axis, but large flux densities on spherical shells rather far away.
Since the latter information is lost in such 1D curve plots, we prefer to show 2D color maps including the separation axis (see Figure 9 and Figure 10 for SU (2) and Figure 13 and Figure 14 for SU (3)).
Comparing the ground state flux densities for Λ ( ) η = Σ + u , Π u and gauge group SU (3), which were computed both in ref. [51] (see Figure 7, Figure 10, Figure 11 and Figure 12 in ref. [51]) and this work (see Figure 11 to Figure 14), we find fair agreement. A detailed comparison is, however, difficult, because of the different QQ separations considered. Concerning statistical errors, our results are more precise by a factor of up to five. An obvious reasons for this is the larger number of gauge link configurations we have been using (4 500 compared to 1 199). Moreover, we have improved our statistical precision by averaging data points, which are related by symmetries, e.g. rotational symmetry as discussed in detail in section 2.3 and section 4.2. Such a symmetrization was not done in [51] as indicated by various plots presented in ref. [51].

Conclusions
We have computed chromoelectric and chromomagnetic flux densities for hybrid static potential states for seven different sectors Λ ( ) (2) and SU(3) lattice gauge theory. These flux densities can be interpreted as flux densities inside heavy hybrid mesons and, thus, provide insights into the structure of such mesons. Five of these sectors, Λ ( ) , are investigated for the first time in this work, while our computation of the remaining two sectors, Λ ( ) η = Σ + u , Π u , confirm results recently published [51], now provided with improved precision. We find flux tubes with interesting structure, significantly different from that of the ordinary static potential with Λ = Σ + g and reminiscent of different vibrational modes of a string. There are also localized peaks in the flux densities, which can be interpreted as valence gluons. Moreover, we compared the resulting flux densities to local operators typically used to study such states, e.g. in pNRQCD.
Concerning future work a straightforward direction would be to consider smaller lattice spacings and larger spatial volumes, i.e. to study the continuum and infinite volume limit. However, we do not expect significant changes in the numerical results presented here, since we already have partly investigated discretization errors (by comparing results obtained with unsmeared and with HYP2 smeared static propagators) and finite volume corrections (by comparing our SU(2) main results to an identical computation with a smaller volume and lower statistics). A more interesting direction would be to extend the computation by including also dynamical light quarks. In principle, one could then study not just heavy-heavy hybrid mesons, but, more generally, heavy-heavy exotic mesons and explore their gluon and light quark distribution at the same time. In practice, however, this might be very challenging, because a hybrid static potential state might decay into the ordinary static potential and one or more light mesons.
B Hybrid static potential flux densities for all Λ η sectors: plots for SU(3) gauge theory Hybrid static potential flux densities for SU(3) gauge theory are shown in Figure 11 to Figure 14. Qualitatively, these plots are very similar to the corresponding SU(2) plots in Figure 7 to Figure 10. For a discussion see section 4.3.