Constructing static quark-anti-quark creation operators from Laplacian eigenmodes

We investigate static quark anti-quark operators based on trial states formed from eigenvectors of the covariant three-dimensional lattice Laplace operator. We test the method by computing the static quark-anti-quark potential and comparing results to standard Wilson loop measurements. The new method is efficient not only for on-axis, but also for many off-axis quark-anti-quark separations when a fine spatial resolution is required. We further improve the ground-state overlap by using multiple eigenvector pairs, weighted with Gaussian profile functions of the eigenvalues, providing a variational basis. The method presented here can be applied to potential functions for all possible excitations of a gluonic string with fixed ends, hybrid or tetra-quark potentials, as well as static-light systems and allows visualization of the spatial distribution of the Laplace trial states.


I. INTRODUCTION
The potential of a static quark-anti-quark pair V 0 (r) has always played an important role in Quantum Chromodynamics (QCD). It can be computed via Wilson loops [1] and established an understanding of confinement and its interplay with asymptotic freedom, a central problem of particle physics, via the formation of a flux tube between quark-anti-quark static charges [2][3][4][5][6][7][8][9]. Confinement manifests itself in the linear rise of V 0 (r) at large r; the corresponding slope is known as the string tension. The static potential can be used in the Born-Oppenheimer approximation [10] to compute the spectrum of quarkonium [11][12][13][14][15][16]. It is also an important observable in setting the scale in lattice QCD. In quenched calculations, the scale has been set using the string tension, but in full QCD the string breaks at the pair-production threshold, making a precise definition difficult. The static energy allows determination of the strong coupling, α s , or, equivalently, Λ MS ; see Refs. [17,18] for recent reviews. Instead of the static energy, one can also use the force F (r) ≡ dV 0 (r)/dr, which is free of the self-energy linear divergence. The dimensionless product r 2 F (r) can be used to set the scale [19] at distances where statistical and systematic uncertainties are under good control, e.g., r 0 or r 1 , defined by r 2 i F (r i ) = c i , with c 0 = 1.65 [19], c 1 = 1 [20].
In this paper, we investigate a method for computing the static quark-anti-quark potential in lattice QCD not based on Wilson loops, but where trial states are formed from components of eigenvectors of the covariant lattice Laplace operator [21]. In this construction, the spatial Wilson lines in the Wilson loop are replaced by outer products of Laplacian eigenvectors. This idea was proposed in the context of adjoint string breaking [22] and of Polyakov loops and the static potential at finite temperature [23,24]. The main advantage is we can not only form straight lines (on-axis), but also offaxis paths very easily. These correspond to very complicated stair-like constructions of spatial link variables. It is important to compute the static potential for many off-axis separations whenever a fine resolution is required, e.g., for a detailed investigation of string breaking [25,26] or to determine the scale Λ MS via matching the perturbative and the lattice QCD static potential [27][28][29][30]. It is even mandatory to compute all possible on-and off-axis separations to determine the static potential in momentum space representation [31].
The implementation of [21] which uses only the eigenvector corresponding to the lowest eigenvalue can be significantly improved by summing over several eigenvectors, weighted by Gaussian profile functions of their corresponding eigenvalues. A similar method was successfully applied to hadronic correlation functions in [32] where an optimal smearing profile was introduced in the distillation framework [33], which can be equivalently expressed as an optimal creation operator for a meson. In the case of the static potential we get an improvement for the static energies, which reach their plateau values at earlier temporal distances, to be quantified below. The improved implementation can also be adapted to measure multi-quark potentials, hybrid static potentials of exotic mesons, where the gluonic string excitations can be realized by applying covariant derivatives to the Laplacian eigenvectors, as well as static-light potentials with insertions of light quark propagators. Further, we present a simple way to illustrate the flux tube between a static quark and antiquark pair using a Laplacian eigenvector pair as a 'test charge' scanning the chromo-electromagnetic field.
The article is organized as follows: First, we reintroduce the notation of Laplace trial states in section II. Next we reformulate the standard Wilson loop in terms of Laplace trial state correlators and discuss their improvement via Gaussian profile functions in section III, allowing us to formulate a generalized eigenvalue problem (GEVP) for the Laplace trial state correlation basis matrix of the static potential, resulting in optimal profile functions for ground and excited states. We test the new improved method on a dynamical fermion ensemble in section IV, presenting results for effective energies, static potentials as well as excited states. In section V we look at the spatial distribution of the optimal Laplace trial states which probe the ground and excited static potentials of a quark-anti-quark pair. We draw our conclusions and give a short outlook in section VI.

II. LAPLACE TRIAL STATES
LetQ a ( x) denote a static color source with a = 1, 2, 3 at spatial position x. Wilson loops arise from correlations in time of trial statesQ( x)U s ( x, y)Q( y) for a static color anticolor source pair located at spatial positions x and y respectively 1 . Note the same Wilson loops are obtained when the static color sources are replaced by static quarks since the heavy quark spins decouple in the static limit and the trace over spin yields a constant, see [34]. The spatial Wilson line U s ( x, y) = exp(i y x A µ dx µ ) = U µ is a path-ordered product of link variables from x to y. We want to replace the spatial part of trial states in each time-slice with an alternative operator which respects the gauge transformation behavior of the spatial Wilson line, given by to ensure gauge invariance of the trial state. The three-dimensional gauge-covariant lattice Laplace operator ∆, acting on a field ψ( x) on a single time-slice of the four-dimensional lattice gives and has the required transformation behavior ∆ ( [35]. It follows, that we can write down a combination of eigenvector components for a given eigenvalue λ, namely v( x)v † ( y), which has the same behavior under gauge transformations as the spatial Wilson line U s ( x, y): At this point, inspired by the distillation operator [33] ab ( z, we introduce the more general operator by including a quark profile ρ i , which modulates contribution from different eigenmodes. Note is a projection matrix, 2 = onto V , the vector space spanned by {v i }, while˜ is no longer idempotent, it still has an image given by the span of v i . Next, we define the auxiliary field on each time-slice 1 We omit the time coordinate in this section since trial states exist on single time-slices only. χ a ( z| x) can be interpreted as an effective smeared colorelectromagnetic field over the whole time-slice induced by the static source at x. At first this seems contradictory to a 'static' color source, but it follows the notation of distillation. We stress the role of the 'smearing parameter' N v , the number of eigenvectors to be summed over in Eq. (3), behaves opposite to intuition. N v = 1 corresponds to the maximal smearing and in the limit where all eigenvectors are included N v → 3N 3 s with N 3 s the spatial lattice volume of a time slice, the smearing operator becomes the identity. This we have to keep in mind when constructing gauge invariant trial states for a color anti-color source pair located at spatial positions x and y, respectively, via where we used the orthonormality of the Laplacian eigenvectors, which ensures that the standard distillation operator is idempotent, i.e., 2 = . We denote Eq. (6) as a Laplace trial state, the positions x and y label the sector of the Hilbert space in which the static energies will be determined. Notice the sum over eigenvectors in Eq. (6) must be truncated at finite N v or a non-trivial profile ρ i must be applied to avoid the collapse of the Laplace trial state, or the annihilation of the quark-anti-quark pair. For example ρ i = δ ik corresponds to the choice of a single eigenvector v k . A simple truncation of the sum at some finite N v = k could be formulated via ρ i = Θ(k − i) and we can of course introduce multiple profile functions to define an operator basis i /4σ 2 k corresponds to a sum over eigenvectors weighted with Gaussian profiles in eigenvalue space with different Gaussian widths σ k , which turned out to be very efficient for meson operators in [32]. In the following section we will reformulate the usual Wilson loops in terms of Laplace trial state correlators and follow the same strategy as in [32] by introducing a set of Gaussian profile functions into the the correlators and solving a generalized eigenvalue problem (GEVP) for the Laplace trial state correlation matrix to extract optimal trial state profilesρ (n) i for ground and excited states of the static potential V n (R), (n = 0, 1, 2 . . .). We also tried other profile functions, e.g., δor Θ-functions to construct an N v × N v transfer matrix with individual eigenmode pair contributions or summing up different numbers of eigenmodes N v to construct a GEVP basis matrix like the ordinary construction using Wilson loops with different spatial smearing levels. Different profiles yield the same results, yet the Gaussian basis seems the most natural (vs. δor step-functions) and numerically stable choice.
i ( y, t) (right, eigenvector pairs to be read in anti-clockwise direction), to form Laplace trial state correlators via the two static perambulatorsτ ij ( x, t 0 , t 1 ) and τ ij ( y, t 0 , t 1 ).
can be rewritten using Laplace trial state correlators by replacing the spatial Wilson lines representing static time-like propagation for a color source at space point y from time t 0 to t 1 is sandwiched between eigenvectors at the corresponding start-and end-times v † i ( y, t 0 ) and v j ( y, t 1 ). Distinct eigenvector indices appear at the source and sink times, so this can be interpreted as the static perambulator at y of time extent T = |t 1 − t 0 |. Its expectation value τ ij ( y, t 0 , t 1 ) vanishes of course. When combined with another static perambulator τ ji ( x, t 1 , t 0 ) at x, it gives the Laplace trial state correlator using the same t S = 3a for all R/a and corresponding ground state energies aV 0 (R) from a cosh-fit, for more details see [32]. These fractional overlaps are listed in table I and demonstrate that a large number N v of eigenvector pairs gives better overlaps for small distances R/a, but with decreasing importance for large distances, where already N v < 100 shows better overlaps.
FIG. 2: The effective energies for R/a = 2 − 4 from Wilson loops and Laplace trial state correlators with increasing numbers of eigenvectors N v . The ground state overlap drastically improves by using more eigenvectors, we see earlier plateaus for larger N v , also quantified in table I. The lines connecting the measured points just help to guide the eye.
Next, instead of trivial quark profiles ρ i,j , we use Gaussian quark profile functions ρ  (GEVP) [36] to identify the optimal trial state profilesρ (n) R (λ) for various energy levels V n (R) (n = 0, 1, 2, . . .). First, we apply the strategy presented in [37,38] and prune L kl using the three most significant singular vectors u i from a singular value decomposition 2 (SVD) at a specific t G = 4 viã L mn = u † m L kl u n , which keeps a smaller set of distinct profiles which improves the stability of the GEVP. We perform the latter at the same t G , separately for all spatial distances R: From the eigenvalues or so-called principal correlators lim t→∞ µ (n) (t, t G ) = e −En(t−t G ) we get the effective energies for a fixed t G , by performing a cosh-fit in practice, due to periodic boundary conditions. From the generalized eigenvectors ν (n) k we can construct the optimal trial state profiles ρ (n) R for the energy states provided by the GEVP, which also depend on the quark separation R, obviously. First, we use the singular vectors u l to get the pruned (or most significant) profilesρ (k) R (λ i ) = l u k,l e −λ 2 i /2σ 2 l . Then we form the linear combination of pruned profiles using the generalized eigenvectors ν k to give the optimal trial state profiles figure 3 for the ground and excited states at R = 4a. The optimal profiles suggest a number N v < 100 of significant/important eigenvectors in the correlator, because each trial state comes with a profile and the combination falls off about twice as fast compared to figure 3. The fractional overlaps with the ground state in table I also favor the Laplace trial 2 L kl = U DV † with U = V (because L is Hermitian in our case) being a unitary matrix, whose column vectors u i form an orthonormal basis, and D being diagonal with non-negative real numbers on the diagonal.
states from a GEVP with optimal profiles in the 6th column, which are even better than standard Wilson loop results from a GEVP with different HYP smearing levels (col. 7).

IV. RESULTS FROM OPTIMAL LAPLACE TRIAL STATES
We performed all our measurements on 48 × 24 3 lattices with periodic boundary conditions except for anti-periodic boundary conditions for the fermions in the temporal direction. They were produced with the openQCD package [39] using the plaquette gauge action and two dynamical nonperturbatively O(a) improved Wilson quarks [40] with a mass equal to half of the physical charm quark mass. The bare gauge coupling is g 2 0 = 6/5.3 and the hopping parameter is κ = 0.13270. The scale r 0 /a = 4.2866(24) [19] and the flow scale [41] is t 0 /a 2 = 1.8477 (3). The corresponding lattice spacing is a = 0.0658(10)fm [42,43]. All measurements were performed by our C+MPI based library that facilitates massively parallel QCD calculations. A total of N v = 200 eigenvectors of the 3D covariant Laplacian were calculated on each time-slice of the lattices as described in [32]. A total of 20 3D APE smearing [44] steps with α AP E = 0.5 were applied on each gauge field before the eigenvector calculation so as to smooth the link variables that enter the Laplacian operator. When forming the correlations of the Laplace trial states, we apply one HYP2 smearing step to the temporal links [34,[45][46][47][48]. Standard Wilson loops were measured using the wloop package [49], also applying one HYP2 step to all gauge links, and 4 levels (0 10 20 30 steps) of spatial HYP smearing to form a variational basis. Wilson loops were measured on 4646 gauge configurations, while Laplace trial states were measured on every fourth configuration only (1160 measurements). The error analysis in this work was done using the Γ method [50,51] with a recent python implementation (pyerror) [52] with automatic differentiation [53].
We compare the effective energies using the improved Laplacian eigenvector approach with Gaussian profiles after solving the GEVP together with smeared Wilson loop results in figure 4. Results from Laplacian modes show better ground state overlaps and higher accuracy than those from Wilson loops with only a quarter of the statistics. The static potentials V n for the ground (n = 0) and excited (n = 1, 2) states. We compare with radially excited string states V 0 + (n + 1)π/R, the lowest 0 ++ isoscalar meson (possible glueball) V 0 + m G from [54] and two times the static-charm meson mass 2m Bc , also evaluated from Laplace trial states.
In figure 5 we present the static potentials V n for the ground (n = 0) and excited (n = 1, 2) states using the Laplace trial states with optimal quark profiles after solving the GEVP. The excited states are just included to show the potential of the method, we want to stress here, that we only have the Laplace trial states in the operator basis, which just like Wilson loops may not have a good overlap with multi-particle states. Note, that we only analyze the Σ + g state according to the nomenclature in [55,56] and its radial excitations, not the first-excited (hybrid) potential Π u , lying between Σ + g (V 0 ) and Σ + g ' (V 1 ), which will be investigated in a future work, using covariant derivatives of eigenvectors in the trial states. For comparison we plot the radially excited string states V 0 + (n + 1)π/R, as well as the lowest 0 ++ isoscalar meson (possible glueball) V 0 + m G from [54] and two times the static-charm meson mass 2m Bc . The latter was also evaluated using the new method, by combining our static perambulators τ ij ( x, t 0 , t 1 ) with a projector P + = (1 + γ 0 )/2 and charm-quark perambulators τ αβ [32], where the quark propagator D −1 includes the dependence on the mass of the quark.
The computational effort of this new method is less than the standard Wilson loop calculation, especially for off-axis separations. In fact, for our test ensemble on a 24 3 × 48 lattice the computation of on-axis Wilson loops using 4 spatial smearing levels (0, 10, 20, 30 HYP steps) is equally expensive as the calculation of 100 Laplacian eigenvectors and Laplace trial states with 7 Gaussian profiles including off-axis distances. The computational advantage of new method can be explained by the fact that the static perambulators can be computed first at each position, resulting in complex numbers, which then can easily be multiplied for arbitrary on-and off-axis separations without the need to compute spatial Wilson lines. In figure 6 we present the optimal static potential V 0 (R) for all on-and off-axis separations R/a from N v = 100 Laplacian eigenvectors compared to on-axis Wilson loop results, which agree well within errors. We also include a measurement of un-smeared Laplace trial state correlators for R/a ≤ 3 (no HYP smearing), showing the Coulomb behavior of the potential at small R. The green points in the plot are shifted vertically such that the un-smeared potential matches the potential with HYP2 smeared temporal links at R/a = 2, which corresponds to removing the free energy difference. Further, we want to note that contrary to Wilson loops, Laplace trial states have an exact symmetry of the potential around half the lattice extension (in a specific direction r), where in fact the force between QQ must vanish due to the periodic boundary conditions, i.e., the static potential should be flat.

V. THE SPATIAL DISTRIBUTION OF OPTIMAL LAPLACE TRIAL STATES
If we do not evaluate the spatial sum in the third line of the Laplace trial state in Eq. (6), we are left with an eigenvector pair v † ( z)v( z) which acts as a 'test-charge' in the original Laplace trial state Nv ijρ (n) which allows the scanning of individual contributions of the quark-anti-quark operator in a 3D time-slice via the free coordinate z. We average over the whole lattice ( x, t), which already gives a very smooth signal on a single configuration. Note that we include the optimal trial state profiles which in this case still depend on the two eigenvalues λ i and λ j , since we did not perform the sum over z in Eq. (6) and therefore did not get a δ ij . The singular vectors u k and generalized eigenvectors ν (n) come from the SVD and GEVP in the static potential calculations for specific quark separation distances R and allow us to look at the flux tube profiles for various energy states of V n (R).
In figure 7 we present the spatial distributions of the optimal Laplace trial states to measure the ground resp. first excited state potentials of a static quark-anti-quark pair at spatial distance R = 10a. The first excitation shows additional nodes in the spatial distribution along and perpendicular to the quark separation axis. The physical interpretation of these distributions in terms of the chromo-electromagnetic flux tube is not clear yet, the optimal profiles certainly contain some information of the ground and excited states of the static potential, the 'test-charge' v( z)v † ( z) however does not measure a specific color field component.

VI. CONCLUSIONS & OUTLOOK
Alternative creation operators for static-quark-anti-quark states based on Laplacian eigenmodes are investigated. The use of a large number of eigenvectors weighted with Gaussian profiles is found to improve performance. An operator basis can be defined via different Gaussian profiles which can be analyzed with the GEVP formalism to extract optimal profiles and Laplace trial states. Temporal correlations of the new operators are used to compute static quark-anti-quark ground and excited state potentials. We observe earlier plateaus in the effective masses compared to standard Wilson loops. One significant advantage of the approach is its efficiency for computing the static potential not only for on-axis, but also for many off-axis quark-anti-quark separations. Indeed the new method requires far less computing time in particular for the latter case, since the eigenvector components of the covariant lattice Laplace operator have to be computed only once and can then be used for arbitrary on-axis and off-axis separations without the need to compute stair-like gauge-link connections. Finally, we visualize the spatial distribution of the optimal Laplace trial states for ground and excited state creation operators of the quark-anti-quark pair. We are currently working on an adaptation of the method to compute hybrid static potentials of exotic mesons, where gluonic string excitations requiring gluonic handles in the standard Wilson loop approach can be realized with covariant derivatives acting on the Laplacian eigenvectors, and to static-light mesons, cf. [26]. First results were presented at the ConfinementXV [57], Lattice 2022 [58] and ExcitedQCD [59] conferences.