# A Novel Gauge Invariant Multi-State Smearing Technique

## I Introduction

Smearing is a technique that has been used for some time in Lattice QCD simulations  in order to improve the quality of the results obtained for a given computational expense by reducing the extent of excited state contamination. The choice of operator used to represent a state in any given simulation is essentially free within the constraint that it must have the appropriate quantum numbers of the desired state, which (usually) guarantees non-zero overlap with the groundstate of those quantum numbers. In Euclidean space, all such operators asymptotically behave as the lowest lying state with quantum numbers matching those of the operator <sup>*</sup><sup>*</sup>*Ignoring questions such as scalar glueball mixing in quenched simulations, and as we shall see special cases where the overlap is fine tuned to be near zero.. We consider the two point correlation function of a set of operators $`𝒪_k`$ with some definite $`J^{PC}`$ in a Euclidean lattice of temporal extent $`T`$. Inserting a complete set of physical intermediate states of the same quantum numbers we find the same asymptotic behaviour for all operators

$`𝒪_k(x)𝒪_k^{}(0)`$ $`=`$ $`{\displaystyle \underset{n}{}}{\displaystyle \frac{0|𝒪_k|nn|𝒪_k^{}|0}{2E_n}}e^{E_n\frac{T}{2}}\mathrm{cosh}E_n(t{\displaystyle \frac{T}{2}})`$ (1)   
$`=`$ $`{\displaystyle \underset{n}{}}C_{kn}C_{kn}^{}{\displaystyle \frac{e^{E_n\frac{T}{2}}}{2E_n}}\mathrm{cosh}E_n(t{\displaystyle \frac{T}{2}})`$ (2)   
$`\stackrel{t\mathrm{}}{}`$ $`|C_{k0}|^2{\displaystyle \frac{e^{E_n\frac{T}{2}}}{2E_0}}\mathrm{cosh}E_0(t{\displaystyle \frac{T}{2}})`$ (3)   
This somewhat egalitarian situation is the principal reason why lattice calculations have traditionally had difficulty extending their scope to radially excited states. However, some choices of operators turn out to be more equal than others, depending on the relative values of the $`C_{nk}`$. For ground state phenomenology we would like $`C_{0k}C_{nk}n1`$. To this end it is common to construct extended operators which project more predominantly onto the groundstate of the desired system; it is this extension procedure that has been dubbed smearing. Typically we require that the smearing procedure be parity and charge positive, and possess at least cubic symmetry

There are two different schools of thought on how to construct smeared operators - gauge invariant smearing and gauge fixed smearing.

### A Gauge Invariant Smearing

Explicit gauge invariance is one of the most attractive features of non-perturbative simulation techniques, and remains to this day a necessary (but sadly insufficient) condition for correct operation of simulation code. Further the gauge fixing process may be multiply defined due to Gribov ambiguities of some gauge fixing conditions , and one has to assume that the ambiguity does not cause any bias in the quantities being measured. For these reason many people are unwilling to use gauge fixing, and proceed with one of a number of gauge invariant techniques.

The most commonly used gauge invariant techniques are Wuppertal , Jacobi , and fuzzed smearing . Wuppertal smearing covariantly smears the quark fields with a spatial distribution which is the solution of the scalar Klein Gordon equation on a given gauge configuration, while Jacobi smearing is a numerically efficient way to approximate this. In order to preserve gauge invariance both the first two methods iteratively apply a gauge invariant nearest neighbour operator to smear out a quark field, with limiting behaviour to weight sites in a “bell shaped” distribution. Given the extended nature of the physical states, and the wavefunction picture in the non-relativistic limit one would expect smeared operators to more accurately project onto the ground state, and indeed this is borne out by simulations. However these gauge invariant techniques give little control over the precise functional form of the smearing beyond a basic radius parameter.

The cost of iteratively applying, say, the Jacobi operator $`N`$ times on a propagator is somewhat prohibitive. Lacock and Michael  suggested the use of a covariantly transported “cross” using low noise *fuzzed* links. The fuzzing prescription iterates

$$U_j^{}(x)=cU_j(x)+\underset{ij}{}\left[U_i^{\mathrm{Staple}}(x,x+\widehat{j})+U_i^{\mathrm{Staple}}(x,x+\widehat{j})\right],$$
(4)   
$$U_j(x)=𝒫U_j^{}(x),$$
(5)   
where $`𝒫`$ is a projection onto SU(3) (not closed under addition) by Cabibbo Marinari maximisation of $`\mathrm{Tr}[U_j^{}U_j^{}]`$ with six hits on each SU(2) subgroup. Typically $`c`$ is set to 2.0, and there are 5 iterations performed. These fuzzed links are used to transport the quark field by some number of sites $`N`$ along each of the spatial axes.

$$\psi ^{}(x)=\underset{i=x,y,z}{}\left[\underset{n=1}{\overset{r}{}}U_i(x+(n1)\widehat{i})\right]\times \psi (x+r\widehat{i})+\left[\underset{n=1}{\overset{r}{}}U_i^{}(xn\widehat{i})\right]\times \psi (xr\widehat{i})$$
(6)   
The fuzzing radius, $`r`$, is chosen to minimise the contamination of the ground state operator.

### B Gauge Fixed Smearing

Gauge fixed smearing (either Coulomb or Landau gauge) allows complete freedom for the functional form of the smearing; in a fixed gauge we are freed from the need to iterate local covariant transport to space fill, and convolve quark fields with arbitrary smearing functions using fast Fourier transforms. In particular, the use of smearing functions with nodes has allowed reliable extraction of excited state energies from NRQCD through the use of multi-exponential fits . Neither the cost of Fourier transforming the data set multiple times to implement this smearing on a massively parallel machine nor the loss of gauge invariance are desirable features, and there is clearly a need for a gauge invariant smearing technique with free radial functional form.

In order to proceed, we consider the non-relativistic interpretation of $`q\overline{q}`$ systems in terms of spatial wavefunctions. This is certainly appropriate to heavy-heavy mesons, however, we can apply this technique in heavy-light and light-light systems without needing to worry about its relevance.

## II Smearing In A Non-relativistic Potential

We consider smearing in the context of a non-relativistic potential model with spherical symmetry. The solutions may be written as $`\psi _{nlm}(r,\theta ,\varphi )=R_{nl}(r)P_{lm}(\theta ,\varphi )`$ where $`P_{lm}(\theta ,\varphi )`$ are the usual spherical harmonics. The orthogonality relation of the wavefunctions is then

$$𝑑V\psi _{mko}^{}\psi _{nlp}=\delta _{kl}\delta _{op}\underset{r=0}{\overset{\mathrm{}}{}}R_{ml}^{}(r)R_{nl}(r)r^2𝑑r$$
(7)   
$$\underset{r=0}{\overset{\mathrm{}}{}}R_{ml}^{}(r)R_{nl}(r)r^2𝑑r=\delta _{mn}$$
(8)   
We shall ignore the $`l`$ and $`m`$ indices from here, and consider spherically symmetric functions. Suppose we construct non-local meson operators

$$𝒪_n(x)=d^3y\overline{\psi }(y)S_n(y,x)\psi (x)$$
(9)   
where the smearing function $`S_n(y,x)=S_n(yx)`$ is thought of as a spatial wavefunction defined on lattice sites

$$S_n(y,0)S_n(y)\underset{i,j,k}{}\delta ^3(yi\widehat{x}j\widehat{y}k\widehat{z})$$
(10)   
Note we have dropped the colour indices since in our wavefunction approximation the gauge degrees of freedom are represented by a colour singlet potential.

The operator $`𝒪`$ will have an overlap with each physical state $`|\psi _m`$ of the relevant $`J^{PC}`$

$$C_{nm}=\psi _m^{}(x)S_n(x)𝑑V=\underset{i,j,k}{}\psi _m^{}(i\widehat{x}+j\widehat{y}+k\widehat{z})S_n(i\widehat{x}+j\widehat{y}+k\widehat{z})$$
(11)   
Now if $`S_n(x)`$ is chosen to well approximate the state $`\psi _n(x)`$ and the lattice is sufficiently fine that the sum (11) is close to the integral

$$C_{nm}=\psi _m^{}S_n(x)𝑑V\psi _m^{}(x)\psi _n(x)𝑑V=\delta _{nm}$$
(12)   
Thus with a well chosen set smearing functions $`S_n`$ we can write

$$C_{nm}=\delta _{nm}+\eta _{nm}$$
(13)   
where $`\eta _{nm}`$ is a small contamination to the signal, giving a set of correlation functions corresponding to each of the physical radial excitations. In fact $`\eta _{nm}0m<n`$ guarantees the lowest state significantly contributing to the correlation function of the $`n`$th radially excited operator is in fact the $`n`$th state, and the effective mass will not decay to the ground state of the $`J^{PC}`$ on a finite lattice.

If we choose a non-space-filling smearing function of the form

$$S_n(\stackrel{}{x},0)\underset{r=0}{\overset{N}{}}r^2\varphi _n(r)\underset{\widehat{\mu }=\widehat{x},\widehat{y},\widehat{z}}{}\left(\delta ^3(\stackrel{}{x}r\widehat{\mu })+\delta ^3(\stackrel{}{x}+r\widehat{\mu })\right)$$
(14)   
where $`\varphi _n`$ is some arbitrary modulating function that is chosen to approximate the true solution $`\psi _n`$, and the factor $`r^2`$ compensates for the non-space-filling nature of the smearing, then the overlap is

$`C_{nm}`$ $`=`$ $`{\displaystyle \psi _m^{}(\stackrel{}{x})S_n(\stackrel{}{x},0)𝑑V}`$ (15)   
$`=`$ $`{\displaystyle \underset{r=0}{\overset{N}{}}}r^2\varphi _n(r)\psi _m^{}(r)`$ (16)   
$`\stackrel{a0}{}`$ $`{\displaystyle \underset{r=0}{\overset{\mathrm{}}{}}}r^2\psi _m^{}\varphi _n𝑑r`$ (17)   
$``$ $`\delta _{nm}.`$ (18)   
Thus if we choose the $`\varphi _m`$ correctly (i.e. to be $`\psi _n`$) then we establish the orthogonality relation required without the need for space filling. The new smearing technique proposed is a generalisation of the above form to incorporate colour:

$$\psi (x)\underset{r=0}{\overset{N}{}}\left(r+{\scriptscriptstyle \frac{1}{2}}\right)^2\varphi _n(r)\underset{i=x,y,z}{}\left\{\left[\underset{n=1}{\overset{r}{}}U_i(x+\left(n1\right)\widehat{i})\right]\psi (x+r\widehat{i})+\left[\underset{n=1}{\overset{r}{}}U_i^{}(xn\widehat{i})\right]\psi (xr\widehat{i})\right\}$$
(19)   
where the links are fuzzed links, c.f. Section I A. We use the factor $`(r+\frac{1}{2})^2`$ as opposed to $`r^2`$ so that there is a non-zero contribution from the local current, since this was observed to improve the statistical noise. The procedure is very similar in cost to fuzzing with a radius $`N`$, since all the terms can be formed as part of the process of fuzzing to the largest radius. The technique therefore eliminates the two main problems with implementing gauge-invariant multi-state smearing, namely that gauge invariance requires expensive iterative space filling techniques, and that it restricts the functional form. Here we have a low-cost gauge invariant technique with an arbitrary functional form, allowing the insertion of both ground and radially excited wavefunctions.

We shall consider the case when the physical solutions are approximately hydrogenic wavefunctions.

## III Optimising the Smearing

We consider the correlation function with a smeared source denoted $`R`$ and a local sink denoted $`L`$, and restrict the argument to one contaminating excited state, and ignore the periodicity of the lattice

$$C_t=Ae^{E_1t}+Be^{E_2t}$$
(20)   
where $`AC_{L1}C_{R1}^{}`$ and $`BC_{L2}C_{R2}^{}`$. The effective mass of this correlator is loosely the negative of the derivative of the log:

$$M_t\frac{d}{dt}\mathrm{log}C_t=\frac{E_1+\frac{B}{A}E_2e^{(E_2E_1)t}}{1+\frac{B}{A}e^{(E_2E_1)t}}.$$
(21)   
By inspection it can be seen that there are two distinct cases for the approach to the plateau:

* $`\frac{B}{A}>0`$ approach from above
* $`\frac{B}{A}<0`$ approach from below, vertical asymptote at $`t=\frac{\mathrm{log}|\frac{B}{A}|}{E_2E_1}`$

Suppose the physical wave functions are well described by the hydrogenic wavefunctions, $`\psi _n`$, with some Bohr radius, $`r_1`$. We wish to approximate these solutions with a trial wavefunction, which we choose to be the hydrogenic solutions $`\varphi _m`$ with Bohr radius $`r_2`$. Consider the overlap of the trial excited state with the physical groundstate,

$`C_{10}`$ $`=`$ $`{\displaystyle r^2\psi _0\varphi _1𝑑r}`$ (22)   
$``$ $`{\displaystyle r^2(1\frac{r}{2r_2})e^{\frac{r}{2r_2}}e^{\frac{r}{r_1}}𝑑r}`$ (23)   
$`=`$ $`c^3\left[\mathrm{\Gamma }(3){\displaystyle \frac{c}{2r_2}}\mathrm{\Gamma }(4)\right]`$ (24)   
$``$ $`2c^3{\displaystyle \frac{\delta }{r_2}}`$ (25)   
where $`c=\frac{2r_1r_2}{r_1+2r_2}`$, $`\delta =r_2r_1`$, and the last step is made for small $`\delta `$. Thus the overlap of the excited state smearing with the physical state is proportional to the error in the Bohr radius. The overlap is positive for smearing radii that are too large, and negative for radii that are too small, and so the effective mass of the smeared - local correlation function will have a vertical asymptote for radii that are too small, and a smooth form for radii that are too large. This dramatic change of behaviour may be used to tune the Bohr radius on very small statistical samples.

In fact it should always be possible to orthogonalise the first radially excited smearing with respect to the physical groundstate by tuning the Bohr radius parameter even when the physical states are not hydrogenic. One choice for the radius will not, however, simultaneously minimise the overlaps $`C_{10}`$ and $`C_{01}`$. The better the ansatz for the smearing functions, the closer together we expect the optimal radii for the groundstate and excited state correlators to lie. Where the ansatz is better we also expect there to be less contamination from the higher excited states, since the physical wavefunctions are all mutually orthogonal.

In order to demonstrate the applicability of this smearing technique, a feasibility study was carried out on a small sample, with a large number of smearing radii, in the light-light, heavy-light and heavy-heavy sectors.

## IV Feasibility Study

We use the notation $`R_{nl}`$ to denote smearing with the corresponding hydrogenic wavefunction with some Bohr radius. Fourteen quenched configurations of a $`24^3\times 48`$ lattice at $`\beta =6.2`$ were used with non-perturbatively $`O(a)`$ improved quarks (i.e. $`C_{sw}=1.61377`$). Kappa values 0.13460 and 0.12300 were used, corresponding to quarks near the strange and charm masses respectively. The lighter quark propagators had local sources, and both local and $`R_{10}`$ sources were generated for the heavy propagator with Bohr radii of 2.4, 2.5, 2.6, 2.75, 2.8, 2.9, 3.0, 3.1, 3.25, 4.0, and 5.5 in lattice units. This allowed the smearing combinations in Table I to be generated for each Bohr radius. The correlation functions were analysed using uncorrelated fits, because the small statistical sample caused the correlation matrix to be too noisy for reliable inversion in correlated fits. As such the estimates of $`\chi ^2`$ are probably unreliable.

Of primary interest is the simultaneous two exponential fit to the $`R_{10}`$ and $`R_{20}`$ source smeared - local sink correlation functions. We show the fitted effective masses for various Bohr radius for the heavy-heavy, heavy-light and light-light sectors in Figures 1 through 6.

The effective mass for the radially excited smearing demonstrates the expected transition in behaviour, and appears to be forming a plateau for $`r_0=2.6(3.0,3.1)`$ in the heavy-heavy (heavy-light, light-light) sector. While the formation of an excited state plateau is clearly strongly dependent on the selection of $`r_0`$ to within about $`0.1`$ (this is possible to do using only one configuration by looking for vertical asymptote in the effective mass plot), we do not wish the fitted values to be dependent on the smearing. To investigate this we study the stability of the fitted masses with the input parameters to the fitting procedure, namely the Bohr radius, and the timeslice range of the fit.

### A Stability Analysis of Fitted Masses

The fitted values for the ground state and excited state masses are given in Table II to Table V, and are found to be remarkably stable in the fit range typically agreeing within statistical error for all values of $`t_{min}3`$.

The level of statistical error on the excited state mass is also particularly pleasing. The error on the excited light-light state is typically between only 5 and 10 times larger than that on the ground state.

In order to study the stability of the fits with respect to the chosen Bohr radius, I tabulate the fitted mass values corresponding to $`t_{min}=4`$ for all the simulated Bohr radii in Table VI, Table VII, and Table VIII. They are stable within the statistical error over almost the entire range of radii simulated.

### B Selecting the Bohr Radius

From inspection of the effective mass plots presented, we can see that the optimal values of the Bohr radius parameter (at least for the $`R_{20}`$ correlator) were about 2.6 for the heavy-heavy state, 3.0 for the heavy-light and 3.1 for the light-light. Using the discussion in Section III we plot the fitted amplitudes of the local-smeared correlation functions versus Bohr radius, looking for the zero of the overlap of the $`R_{20}`$ correlator with the groundstate.

Figure 7, Figure 8 and Figure 9 present the dependence of the four fitted amplitudes on the Bohr radius used in the heavy-heavy, heavy-light and light-light systems respectively. The linear fits are a good approximation, producing the optimal values for the smearing radii given in Table IX.

While the effective mass plots presented were on fourteen configurations, it is clear that following the pole on the plot gave very accurate optimisation of the Bohr radius, as demonstrated by later fits to the amplitudes. It was anticipated in Section III that, since our smearing function basis is not the physical one, the radius for minimal contamination of the groundstate by the excited state would differ slightly from the radii which gave minimal contamination of the excited state. This is not apparent from the data presented in Figure 7, Figure 8 and Figure 9, where there is no minimum of the $`R_{10}`$ \- excited state amplitude. We must bear in mind however that the contamination of the $`R_{10}`$ correlator is measured at early times, whereas the contamination of the $`R_{20}`$ correlator is found from its asymptotic time dependence, leaving the determination of the optimal Bohr radius from the $`R_{10}`$ correlator much more subject to contamination from higher excited states. For this reason, it is dangerous to look only at the local-smeared amplitude, since cancelling contributions from the higher states can cause a seemingly flat effective mass.

Instead we study the smeared-smeared correlators, which have positive definite contaminating contributions, and minimise the amplitude for the contamination of the $`R_{10}`$ \- $`R_{10}`$ correlator. These correlation functions were only generated for the heavy-heavy combination, and Figure 10 presents a sample fit to the pseudoscalar data. Here the contamination below timeslice 5 is apparent, as well as a weaker drift to the plateau beyond this point. It is this ultimate drift we must study. Figure 11 plots the fitted contaminating amplitude for this correlator versus radius; we expect from Section III a quadratic dependence on the smeared-smeared correlator. We obtain a minimising radius of $`r_0=2.1(4)`$ with a quadratic fit.

### C Use of Doubly Smeared Correlation Functions

In this section I demonstrate that the signal obtained from doubly smeared correlators is useful. The double smearing is created by constructing the meson propagator from two propagators which are each source smeared, while, for the purposes of this demonstrate, we use local sinks. The double exponential fit to the double $`R_{10}`$ and single $`R_{20}`$ smeared correlation functions are presented in Figure 12 and Figure 13 for the heavy-heavy and light-light systems respectively.

Clearly this smearing combination, while not tailored for any particular state constitutes a valid signal that can be included in multi-channel multi-exponential fits, which is always of use, especially for the p-states.

## V Comparison With Fuzzing

In order to fairly compare the level of statistical error with traditional techniques, I present double exponential fits to the $`R_{10}`$ sink and the fuzzed sink, both with the optimal radius in Figure 14, Figure 15 and Figure 16 for the heavy-heavy, heavy-light and light-light systems respectively, on the same set of local source propagators.

It can be seen that the statistical error from the $`R_{10}`$ correlator is much reduced with respect to the fuzzed data. Admittedly the fuzzed data does appear to plateau earlier, with a caveat however; the noise on the fuzzed data is such that there is quite possibly a “double approach” to the plateau, due to contributions from more than one excited state, as seen in the large statistics fuzzed data at $`\beta =6.2`$ in Figure 17, so that with higher statistics the plateau may well not survive. It has also been found that the fuzzed-fuzzed combination tends to plateau after the fuzzed-local, indicating opposing contributions from excited states, and meaning that the plateau at early times is a misleading balance between contaminations of opposite sign (at least until after the fuzzed-fuzzed correlator has reached a plateau), c.f. Figure 17.

## VI Conclusions

This smearing method provides, for the first time, a gauge invariant technique for inserting arbitrary radial wavefunctions in a smeared operator. So far results have only been presented using hydrogenic wavefunctions, though in principle any form could be used. Excellent stability of fits to both the ground and radially excited states was demonstrated, and the level of statistical error was significantly smaller than traditional techiques on the same ensemble. Tuning the radius of the wavefunctions allowed a clear plateau for the radially excited state to be isolated. An understanding of the dependence of the fitted amplitudes on the radius of the smearing was obtained, which can be used to determine more accurately the optimal radius.

The feasibility study did not include a sufficient number of quark masses to extrapolate to the physical spectrum, however a rudimentary comparison of the $`2S1S`$ splittings with experiment is made in Table X, and can be seen to be plausible. The main source of error, particularly in the light-light system, is expected to be finite volume effects, since the Bohr radius required suggests the extent of the wavefunction is significantly bigger than the lattice used.

Since the technique allows a free choice for the wavefunction it would be interesting to try to find a better basis than the hydrogenic wavefunctions. Harmonic oscillator wavefunctions have been tried and found to be inferior to hydrogenic. Other possible options are numerically solving the Cornell potential and inserting the solution as the smearing function, as has been performed by the SESAM collaboration in the gauge fixed context , or alternatively measuring the wavefunction within a lattice calculation and re-inserting this. It would be interesting to compute the $`R_{20}`$ smeared correlators on a large volume for light quarks since we expect finite volume effects on the current volume.

## VII Acknowledgements

The work of this paper was carried out under the supervision of Richard Kenway, Ken Bowler and Brian Pendleton, whom I wish to thank for many useful conversations. I also wish to thank Christine Davies for proof reading this work. The calculations were performed on the Cray T3D at the EPCC using UKQCD computer time. I wish to acknowledge the support of EPSRC grant GR/K41663 and PPARC grant GR/K55745. The author was funded by the Carnegie Trust for the Universities of Scotland while this work was carried out at the University of Edinburgh, and is grateful for PPARC grant PP/CBA/62, and to the University of Glasgow for support while writing this paper.