format: text

text:
# A lattice study of interaction mechanisms in a heavy-light meson-meson system

## I Introduction

The term “lattice hadron physics” has been coined for strong interaction physics based on first principles, i.e. quantum chromodynamics (QCD). In recent years hadron physics has emerged as a field in its own right Capstick et al. (2000). The desire to explain strong interaction phenomena in terms of the underlying dynamics of quarks and gluons sets the field apart from traditional nuclear physics, which emphasizes an effective field theory point of view. Among the important issues faced by hadron physics are baryon and meson spectroscopy, and structure, as well as the mechanism of their strong interaction. Those issues have in common the need to deal with excited states. In lattice QCD, which affords the most direct access to hadron physics at nuclear energy scales, excited states spectroscopy is just now moving into reach due to the use of anisotropic lattices, advanced analysis techniques, and powerful computing facilities.

This work is concerned with learning about hadronic interaction mechanisms. We believe that much of the physics of hadronic interaction can be understood by investigating heavy-light systems. In those the relative distance of hadrons is a well defined quantity. The lattice ‘data’ can be interpreted in terms of intuitive pictures, like potentials. Insight into mechanisms of the strong interaction flows from looking at excitations due to hadron-hadron operators, say $`\mathrm{\Phi }(t)`$, at various relative distances $`r`$. Different choices for the structure of those operators in terms of their composition from quark and gluon fields may potentially point at interesting physics of the system, such as the importance of quark versus gluon exchange degrees of freedom as a function of $`r`$. To extract information of this kind the computation of matrix elements $`n|\mathrm{\Phi }(t_0)|0`$ between the vacuum $`|0`$ and ground and excited states $`|n`$, $`n>0`$, will be useful.

From the vantage point of a numerical lattice simulation this can be a notoriously difficult problem. Hadron-hadron operators are prone to produce very noisy correlation functions. Extracting spectral information in the standard fashion, i.e. trying to identify a plateau in an effective mass function, may not be practical.

An alternative analysis method that is just being ‘discovered’ by lattice practitioners is application of the maximum entropy method (MEM) or an otherwise constrained form of Bayesian inference Nakahara et al. (2000); Lepage et al. (2002); Yamazaki et al. (2002). From the Bayesian perspective the parameters of a model for the Euclidean time correlation function are viewed as random variables drawn from a certain probability distribution function. The latter, known as the posterior probability $`𝒫[\rho C]`$, is the conditional probability for a certain parameter set $`\rho `$ given a data set $`C`$. In our case the data set $`C`$ is the lattice-measured time correlation function $`C(t,t_0)`$ and the parameter set $`\rho `$ is the spectral density function in the model

$$F(\rho |t,t_0)=𝑑\omega \rho (\omega )e^{\omega (tt_0)}.$$
(1)
Discretization understood, in this approach the number of parameters is allowed to exceed the number of data points without causing conceptual problems. In the MEM the posterior probability is constructed from the $`\chi ^2`$-distance between the lattice data $`C`$ and the model $`F`$, and the information content of $`\rho `$ measured by the entropy $`S=𝑑\omega \rho (\omega )\mathrm{ln}\rho (\omega )`$. Usually, the result inferred from the Bayesian approach is the most likely parameter set $`\rho `$.

Analytically, the spectral density function is a sum of discrete $`\delta `$-peaks

$$\rho (\omega )=\underset{n0}{}\delta (\omega \omega _n)|n|\mathrm{\Phi }(t_0)|0|^2.$$
(2)
From a numerical viewpoint computational constraints on $`C`$ render the peak widths finite. Physical quantities are contained in each peak as low $`\omega `$ moments. Among those are the peak volume $`|n|\widehat{\mathrm{\Phi }}(t_0)|0|^2`$ and the peak energy $`E_n`$, i.e. the mean value of $`\omega `$.

The Bayesian spectral analysis of the lattice data is an interesting problem in itself. It leads to the discussion of a number of computational strategies in a general context. In order to keep this presentation focused, we refer the reader to a separate paper Fiebig (2002a) where selected aspects of Bayesian spectral analysis are discussed, using the same lattice simulation data as a testing ground.

Starting a little more than a decade ago lattice work on hadronic interaction has followed mostly two tracks; investigation of heavy-light systems for varying hadron relative distance Richards et al. (1990), and the computation of scattering phase shifts from energy spectra of two hadrons in a finite box Lüscher (1991). We will not review these subjects here, but rather mention a few leads to facilitate following the literature.111For reviews see Fiebig and Markum (2002); Green (2002) Extraction of scattering phase shifts has been successful in terms of structureless particles Goe . On the other hand, scattering of composite hadrons within Lüscher’s framework is considerably more difficult. Scattering phase shifts have been obtained for the S-wave interaction of the $`\pi \pi `$ system in the $`I=2`$ channel Fiebig et al. (2000), and by the CP-PACS collaboration Aoki et al. (2002). Considerable computing power was brought to bear by the JLQCD collaboration to extracting scattering lengths for the $`\pi \pi `$, $`\pi `$N, and NN systems Fukugita et al. (1995). While lattice volume limitations hinder a realistic treatment of the NN system, the $`\pi `$ scattering lengths of Ref. Fukugita et al. (1995) stand out as the only quantitative results for hadronic interaction from the lattice to date. Systems with heavy, even static, quarks will not yield quantitative results, but are useful to gain a deeper understanding of the mechanisms that lead to the phenomenology of nuclear forces. Green and coworkers have been studying energies of static four-quark systems in various geometric configurations, most recently with two static and two light quarks Green et al. (2000). Their results shed light on the importance of many-body versus two-body forces in hadrons. Other studies of heavy-light systems Mihály et al. (1997) focus on adiabatic potentials, i.e. the ground state energies as functions of the relative hadron-hadron separation. Along these lines, most of the hadronic interaction work by members of the UKQCD collaboration can be traced from Michael and Pennanen (1999). There, heavy-light meson-meson systems are classified according to isospin and spin configurations of the light quarks ($``$ mesons). Depending on the channel, interaction potentials at around $`0.5\mathrm{fm}`$ come out attractive or repulsive, with magnitudes rarely exceeding $`50\mathrm{M}\mathrm{e}\mathrm{V}`$. Though this is typical for nuclear physics, from the point of view of a lattice simulation these are very small energy differences to be measured. This is a generic problem for hadronic interaction physics on the lattice.

In the present work we report on a study of a system of two heavy-light mesons based on two operators. The first one, $`\mathrm{\Phi }_1`$, is the standard product of two local pseudoscalar heavy-light meson fields at distance $`r`$. The other one, $`\mathrm{\Phi }_2`$, is similar, but nonlocal, probing color rearrangement. In some sense $`\mathrm{\Phi }_1`$ and $`\mathrm{\Phi }_2`$ test the importance of quark and gluon exchange degrees of freedom for the interaction. They enter into a $`2\times 2`$ time correlation matrix. Our goal here is to learn about the interaction mechanisms represented by those operators as a function of the relative meson-meson distance. The current work goes beyond previous studies of heavy-light meson-meson systems mainly in that (i) nonlocal operators are considered and (ii) operator mixing is built into time correlation functions, allowing the system to ‘dynamically pick’ the dominant excitation mechanism as the meson-meson separation $`r`$ varies. Moreover, due to the spectral analysis used in this work we are able to (iii) extract the actual excitation strengths of the ground and excited two-meson states.

Preliminary results reported in Fiebig (2002b, c) were based on an analysis of the diagonal correlator elements, with no mixing. There, averaging over annealing start configurations had not been done, the spectral density functions came from single annealing runs. As it turns out this is a source of systematic error that can not be tolerated in the light of the smallness of the extracted energy shifts.

## II Lattice action

The meson-meson operators employed in this simulation lead to somewhat massive states. The resulting steep drop of time correlation functions, particularly for the excited states, makes it very difficult to analyze the lattice signal. To deal with this situation we use an anisotropic lattice action. If the aspect ratio $`\xi =a_s/a_t`$ of the lattice constants in space and time directions, respectively, is made larger than one the number of usable time correlation function data is increased before the signal ‘vanishes’ into noise. Anisotropic lattices have been essential for computing the glue ball mass spectrum Morningstar and Peardon (1999). The current simulation of hadronic interaction has in common the need for extracting excited states.

The gauge field part of the lattice action has the form

$`S_G[U]`$ $`=`$ $`\beta {\displaystyle \underset{\mathrm{}}{}}c_{\mathrm{}}\mathrm{\Omega }_{\mathrm{}}\mathrm{with}`$ (3a)
$`\mathrm{\Omega }_{\mathrm{}}`$ $`=`$ $`{\displaystyle \underset{𝒞𝒮_{\mathrm{}}}{}}{\displaystyle \frac{1}{3}}\mathrm{ReTr}[𝟙𝕌(𝒞)].`$ (3b)
Here $`\beta =6/g^2`$ in terms of the gauge coupling $`g`$, $`\mathrm{}`$ labels sets $`𝒮_{\mathrm{}}`$ of closed lattice contours $`𝒞`$ and $`U(𝒞)`$ denotes the path ordered product of gauge field link variables along $`𝒞`$. We adopt a tree-level tadpole improvement scheme with four classes of loops: all oriented spatial elementary plaquettes $`\mathrm{}=ss`$, temporal elementary plaquettes $`\mathrm{}=st`$, spatial planar rectangles $`\mathrm{}=sss`$, and short temporal planar rectangles with two spatial and one temporal link $`\mathrm{}=sst`$. Specifically

$`S_G[U]`$ $`=`$ $`\beta \{{\displaystyle \frac{5}{3\xi u_s^4}}\mathrm{\Omega }_{ss}+{\displaystyle \frac{4\xi }{3u_s^2u_t^2}}\mathrm{\Omega }_{st}`$ (4)
$`{\displaystyle \frac{1}{12\xi u_s^6}}\mathrm{\Omega }_{sss}{\displaystyle \frac{\xi }{12u_s^4u_t^2}}\mathrm{\Omega }_{sst}\},`$
where $`u_s`$ and $`u_t`$ are spatial and temporal link renormalization factors Lepage and Mackenzie (1993), and

$`\mathrm{\Omega }_{ss}`$ $`=`$ $`{\displaystyle \underset{x}{}}{\displaystyle \underset{1\mu <\nu 3}{}}{\displaystyle \frac{1}{3}}\mathrm{ReTr}[𝟙𝕌_\mu (𝕩)\times `$ (5)
$`U_\nu (x+\widehat{\mu })U_\mu ^{}(x+\widehat{\nu })U_\nu ^{}(x)]`$
$`\mathrm{\Omega }_{st}`$ $`=`$ $`{\displaystyle \underset{x}{}}{\displaystyle \underset{1\mu 3}{}}{\displaystyle \frac{1}{3}}\mathrm{ReTr}[𝟙𝕌_\mu (𝕩)\times `$ (6)
$`U_4(x+\widehat{\mu })U_\mu ^{}(x+\widehat{4})U_4^{}(x)]`$
$`\mathrm{\Omega }_{sss}`$ $`=`$ $`{\displaystyle \underset{x}{}}{\displaystyle \underset{1\mu \nu 3}{}}{\displaystyle \frac{1}{3}}\mathrm{ReTr}[𝟙𝕌_\mu (𝕩)𝕌_\mu (𝕩+\widehat{\mu })\times `$ (7)
$`U_\nu (x+2\widehat{\mu })U_\mu ^{}(x+\widehat{\nu }+\widehat{\mu })U_\mu ^{}(x+\widehat{\nu })U_\nu ^{}(x)]`$
$`\mathrm{\Omega }_{sst}`$ $`=`$ $`{\displaystyle \underset{x}{}}{\displaystyle \underset{1\mu 3}{}}{\displaystyle \frac{1}{3}}\mathrm{ReTr}[𝟙𝕌_\mu (𝕩)𝕌_\mu (𝕩+\widehat{\mu })\times `$ (8)
$`U_4(x+2\widehat{\mu })U_\mu ^{}(x+\widehat{4}+\widehat{\mu })U_\mu ^{}(x+\widehat{4})U_4^{}(x)]..`$
This action is the same as in Morningstar and Peardon (1999), the chosen parameters are $`\beta =2.4`$ and $`\xi =3`$. The spatial link renormalization is computed self consistently from the average spatial plaquette,

$$u_s=\left(1\mathrm{\Omega }_{ss}/3L^3\right)^{1/4},$$
(9)
while the temporal renormalization factor is set to one, $`u_t=1`$. According to Morningstar and Peardon (1999) this results in a spatial lattice constant of about $`a_s0.25\mathrm{fm}`$, or $`a_s^1800\mathrm{M}\mathrm{e}\mathrm{V}`$.

The fermion action is

$$S_F[U,\psi ,\overline{\psi }]=\underset{x,y}{}\overline{\psi }_{fA\mu }(x)Q_{fA\mu ,gB\nu }(x,y)\psi _{gB\nu }(y),$$
(10)
with $`f,g`$=flavor, $`A,B`$=color, and $`\mu ,\nu `$=Dirac indices. The fermion matrix $`Q`$ is assumed flavor diagonal

$$Q_{fA\mu ,gB\nu }(x,y)=\delta _{fg}Q_{A\mu ,B\nu }^{(f)}(x,y),$$
(11)
and has a Wilson and a clover term Wilson (1974); Sheikholeslami and Wohlert (1985)

$$Q^{(f)}=𝟙\kappa ^{(𝕗)}(𝕄𝕔_{\mathrm{SW}}𝕂).$$
(12)
In detail

$`M(x,y)=`$
$`{\displaystyle \frac{1}{u_s}}{\displaystyle \underset{\mu =1}{\overset{3}{}}}\left[(r_s\gamma _\mu )U_\mu (x)\delta _{x+\widehat{\mu },y}+(r_s+\gamma _\mu )U_\mu ^{}(y)\delta _{x,y+\widehat{\mu }}\right]`$
$`+{\displaystyle \frac{\xi }{u_t}}[(r_t\gamma _4)U_4(x)\delta _{x+\widehat{4},y}+(r_t+\gamma _4)U_4^{}(y)\delta _{x,y+\widehat{4}}].`$
The critical hopping parameter in the anisotropic case is $`\kappa _c=(2\xi r_t+6r_s)^1`$. The Wilson parameters are chosen as $`r_s=r_t=1`$. Because the lattice is coarse in the space directions and, more importantly, because we wish to avoid problems with ghosts (unphysical branches in the lattice-quark dispersion relation Alford et al. (1997)) we use clover improvement Sheikholeslami and Wohlert (1985) only in the spatial planes,

$$K(x,y)=\delta _{x,y}\frac{1}{u_s^4}\underset{1\mu <\nu 3}{}\sigma _{\mu \nu }\frac{1}{2i}\left(P_{\mu \nu }(x)P_{\mu \nu }^{}(x)\right).$$
(14)
Here $`\sigma _{\mu \nu }=\frac{i}{2}[\gamma _\mu ,\gamma _\nu ]`$ and $`P_{\mu \nu }(x)=\frac{1}{4}_{k=1}^4P_{k,\mu \nu }(x)`$ is made from four transport operators along oriented 4-link paths in the $`\mu `$$`\nu `$ plane, starting and terminating at $`x`$, for example $`P_{1,\mu \nu }(x)=U_\mu (x)U_\nu (x+\widehat{\mu })U_\mu ^{}(x+\widehat{\nu })U_\nu ^{}(x)`$, see Sheikholeslami and Wohlert (1985).

The strength coefficient is $`c_{\mathrm{sw}}=1`$. At $`\kappa ^{(f)}=0.0679`$ we have $`m_\pi /m_\rho 0.75`$ on a $`L^3\times T=10^3\times 30`$ lattice. Using $`m_{\mathrm{eff}}a_t=0.4676`$ from Fiebig (2002a) a crude estimate for the quark mass puts it within 15% above the strange mass scale. We have used a hybrid molecular dynamics algorithm (HMC) Duane and Kogut (1986) to generate $`N_U=708`$ quenched gauge configurations.

## III Operators

At this time we wish to address the physical mechanisms responsible for the features of hadronic interaction rather than making quantitatively precise predictions. Toward this end we employ, for a two-meson system, a set of operators meant to excite different QCD degrees of freedom. For two heavy-light pseudoscalar mesons, for example, the operator

$$\mathrm{\Phi }_1(t)=\underset{\stackrel{}{x},\stackrel{}{y}}{}\delta _{\stackrel{}{r},\stackrel{}{x}\stackrel{}{y}}\overline{Q}_A(\stackrel{}{x}t)\gamma _5q_A(\stackrel{}{x}t)\overline{Q}_B(\stackrel{}{y}t)\gamma _5q_B(\stackrel{}{y}t)$$
(15)
is built from separable products of two color singlets, describing two mesons at relative distance $`\stackrel{}{r}`$. In (15) $`Q`$ and $`q`$ are the heavy and light quark fields, respectively, and $`A,B`$ are color indices. Color contractions are done between quark fields which spatially coincide at $`\stackrel{}{x}`$ and $`\stackrel{}{y}`$, respectively. So far, only local operators of this kind have been employed for hadronic interaction studies Fiebig et al. (2000); Michael and Pennanen (1999); Mihály et al. (1997); Richards et al. (1990). At small values of $`r`$ color coupling schemes involving quarks in different mesons may become dynamically possible. An operator testing excitations of that nature is

$`\mathrm{\Phi }_2(t)`$ $`=`$ $`{\displaystyle \underset{\stackrel{}{x},\stackrel{}{y}}{}}\delta _{\stackrel{}{r},\stackrel{}{x}\stackrel{}{y}}U_{P;AA^{}}(\stackrel{}{x}t,\stackrel{}{y}t)U_{P^{};B^{}B}^{}(\stackrel{}{x}t,\stackrel{}{y}t)`$ (16)
$`\overline{Q}_A(\stackrel{}{x}t)\gamma _5q_B(\stackrel{}{x}t)\overline{Q}_B^{}(\stackrel{}{y}t)\gamma _5q_A^{}(\stackrel{}{y}t).`$
It involves link products $`U_P(x,y)`$ along spatial paths $`P`$ within a fixed time slice. The operator $`\mathrm{\Phi }_2(t)`$ interpolates fields that are still products of two color singlets, but those are now formed from quarks at different locations $`\stackrel{}{x}`$ and $`\stackrel{}{y}`$.

In order to simplify the notation we will tacitly assume dependence of all subsequent quantities on the relative distance $`\stackrel{}{r}`$, but suppress $`\stackrel{}{r}`$ in most of the expressions below.

We are thus lead to computing the elements

$$C_{ij}(t,t_0)=\widehat{\mathrm{\Phi }}_i^{}(t)\widehat{\mathrm{\Phi }}_j(t_0),i,j=1,2,$$
(17)
of $`2\times 2`$ time correlation matrices, one for each $`\stackrel{}{r}`$, where $`\widehat{\mathrm{\Phi }}=\mathrm{\Phi }\mathrm{\Phi }`$, means vacuum-subtracted operators. The correlation matrix $`C(t,t_0)`$ built from the operators (15) and (16) can be worked out by way of Wick’s theorem, symbolically

$`\overline{Q}q\overline{Q}q\overline{q}Q\overline{q}Q`$ $`=`$ $`\stackrel{4}{\text{}\overline{Q}}\stackrel{3}{\text{}q}\stackrel{2}{\text{}\overline{Q}}\stackrel{1}{\text{}q}\stackrel{1}{\text{}\overline{q}}\stackrel{2}{\text{}Q}\stackrel{3}{\text{}\overline{q}}\stackrel{4}{\text{}Q}+\stackrel{4}{\text{}\overline{Q}}\stackrel{3}{\text{}q}\stackrel{2}{\text{}\overline{Q}}\stackrel{1}{\text{}q}\stackrel{1}{\text{}\overline{q}}\stackrel{4}{\text{}Q}\stackrel{3}{\text{}\overline{q}}\stackrel{2}{\text{}Q}+\stackrel{4}{\text{}\overline{Q}}\stackrel{3}{\text{}q}\stackrel{2}{\text{}\overline{Q}}\stackrel{1}{\text{}q}\stackrel{3}{\text{}\overline{q}}\stackrel{2}{\text{}Q}\stackrel{1}{\text{}\overline{q}}\stackrel{4}{\text{}Q}+\stackrel{4}{\text{}\overline{Q}}\stackrel{3}{\text{}q}\stackrel{2}{\text{}\overline{Q}}\stackrel{1}{\text{}q}\stackrel{3}{\text{}\overline{q}}\stackrel{4}{\text{}Q}\stackrel{1}{\text{}\overline{q}}\stackrel{2}{\text{}Q},`$ (18)
where pairs $`nn`$ of numbers $`n=1\mathrm{}4`$ denote contractions. This leads to the following expression

$`C_{ij}(t,t_0)=2\delta _{\stackrel{}{r},\stackrel{}{r}^{}}^{(+)}{\displaystyle \underset{\stackrel{}{x},\stackrel{}{y}}{}}\delta _{\stackrel{}{r},\stackrel{}{x}\stackrel{}{y}}^{(+)}𝒰_{i;AA^{},BB^{}}(t,\stackrel{}{x}\stackrel{}{y})𝒰_{j;DD^{},CC^{}}(t_0,\stackrel{}{x}\stackrel{}{y})`$ (19)
$`H_{B^{}\nu ^{},D^{}\lambda ^{}}^{}(\stackrel{}{y}t,\stackrel{}{y}t_0)H_{A\nu ,C\lambda }^{}(\stackrel{}{x}t,\stackrel{}{x}t_0)[G_{A^{}\nu ^{},C^{}\lambda ^{}}(\stackrel{}{y}t,\stackrel{}{y}t_0)G_{B\nu ,D\lambda }(\stackrel{}{x}t,\stackrel{}{x}t_0)G_{A^{}\nu ^{},D\lambda }(\stackrel{}{y}t,\stackrel{}{x}t_0)G_{B\nu ,C^{}\lambda ^{}}(\stackrel{}{x}t,\stackrel{}{y}t_0)].`$
In (19) the symmetrized Kronecker symbol

$$\delta _{\stackrel{}{r},\stackrel{}{r}^{}}^{(+)}=\delta _{+\stackrel{}{r},\stackrel{}{r}^{}}+\delta _{\stackrel{}{r},\stackrel{}{r}^{}}$$
(20)
is related to $`O(3,)`$ symmetry, projecting heavy-quark distances to absolute lengths $`r=|\stackrel{}{x}\stackrel{}{y}|`$. The heavy anti-quark propagator is employed in the static approximation, i.e. the leading term of the hopping parameter expansion222The factor $`(2\kappa )^{tt_0}`$ only produces an irrelevant mass shift. It has been dropped from all subsequent correlator functions.

$$H_{A\nu ,C\lambda }^{}(\stackrel{}{x}t,\stackrel{}{y}t_0)=\delta _{\stackrel{}{x},\stackrel{}{y}}(2\kappa )^{tt_0}\frac{1}{2}(𝟙+\gamma _\mathrm{𝟜})_{\lambda \nu }𝕌_𝔸(\stackrel{}{𝕩}𝕥_\mathrm{𝟘},\stackrel{}{𝕩}𝕥),$$
(21)
where $`U(\stackrel{}{x}t_0,\stackrel{}{x}t)`$ is the path-ordered link variable product along a straight line from $`\stackrel{}{x}t_0`$ to $`\stackrel{}{x}t`$. Finally, depending on whether $`\mathrm{\Phi }_1`$ or $`\mathrm{\Phi }_2`$ is involved in the correlator matrix element, the contour operators

$`𝒰_{1;AA^{},BB^{}}(t,\stackrel{}{x}\stackrel{}{y})`$ $`=`$ $`\delta _{AB}\delta _{A^{}B^{}}`$ (22)
$`𝒰_{2;AA^{},BB^{}}(t,\stackrel{}{x}\stackrel{}{y})`$ $`=`$ $`U_{P;AA^{}}(\stackrel{}{x}t,\stackrel{}{y}t)U_{P^{};BB^{}}^{}(\stackrel{}{x}t,\stackrel{}{y}t)`$ (23)
$`+U_{P^{};AA^{}}(\stackrel{}{x}t,\stackrel{}{y}t)U_{P;BB^{}}^{}(\stackrel{}{x}t,\stackrel{}{y}t)`$
are trivial, or involve path-ordered link variable products along purely spatial paths $`P,P^{}`$ from $`\stackrel{}{x}t`$ to $`\stackrel{}{y}t`$. On the present lattice we consider straight on-axis paths of lengths $`r=1,2,3,4`$.

Diagrammatic likenesses of the correlator matrix elements (19) are shown in Fig. 1. The simplest one, the graph of $`C_{11}`$, is the standard quark exchange diagram that is usually considered Fiebig et al. (2000); Michael and Pennanen (1999); Mihály et al. (1997); Richards et al. (1990), while the rest involve explicit gluon degrees of freedom.

In (19) light-quark propagator elements, e.g. $`G_{A^{}\nu ^{},D\lambda }(\stackrel{}{y}t,\stackrel{}{x}t_0)`$, between arbitrary lattice sites $`\stackrel{}{x}t,\stackrel{}{y}t_0`$ emerge. The source time slice $`t_0`$ may be kept fixed, however, all-to-all spatial propagator elements are needed. More precisely, after working out one of the space-site sums in (19) the number of propagator columns needed may be minimized by using translational invariance of the gauge field average $`\mathrm{}`$. Then, the number of propagator columns becomes equal to the number of relative distances $`\stackrel{}{r}`$. On the other hand, employing all-to-all propagator elements has the advantage of greater flexibility with regard to varied choices for $`\stackrel{}{r}`$, like off-axis distances (future studies), and improved statistics due to space-site averaging in (19). Because the gauge link contour operators (21) and (23) tend to make the correlation matrix (19) quite noisy, we follow the latter strategy.

Random-source estimation is a proven technique for generating all-to-all propagators Scalettar et al. (1986); Duane and Kogut (1986); Gottlieb et al. (1987); Fiebig et al. (2000); Michael and Peisa (1998). Consider the following linear equation

$`{\displaystyle \underset{\stackrel{}{y}y_4}{}}{\displaystyle \underset{B\nu }{}}Q_{A\mu ,B\nu }(\stackrel{}{x}x_4,\stackrel{}{y}y_4)X_{B\nu }^{(A^S\mu ^Srx_4^S)}(\stackrel{}{y}y_4)`$ $`=`$ (24)
$`\delta _{AA^S}\delta _{\mu \mu ^S}R^{(A^S\mu ^Srx_4^S)}(\stackrel{}{x})\delta _{x_4x_4^S}.`$
On the right-hand side $`R^{(A^S\mu ^Srx_4^S)}`$ denotes a vector of length $`L^3`$ of complex random deviates. Indices superscripted with S characterize the source with respect to color, Dirac, and time slice. We choose sources that are nonzero on only one fixed time slice $`x_4^S=t_0`$, whereas for each of the 12 color-Dirac combinations $`A^S=1\mathrm{}3`$, $`\mu ^S=1\mathrm{}4`$, in turn, $`r=1\mathrm{}N_R`$ different random vectors are generated. In the simulation complex $`_2`$ distributed random deviates were used Dong and Liu (1994). With $`N_R=8`$ this results in $`96`$ statistically independent sources per gauge configuration. Writing $`_{<r>}`$ for the ensemble average, of which $`\frac{1}{N_R}_{r=1}^{N_R}`$ is a truncation, we have

$$\underset{<r>}{}R^{(A^S\mu ^Srx_4^S)}(\stackrel{}{x})R^{(A^S\mu ^Srx_4^S)}(\stackrel{}{y})=\delta _{\stackrel{}{x}\stackrel{}{y}}$$
(25)
for all $`A^S,\mu ^S,x_4^S`$ with appropriately chosen normalization. Application of (25) to (24) yields

$$G_{B\nu ,A\mu }(\stackrel{}{y}t,\stackrel{}{x}t_0)=\underset{<r>}{}X_{B\nu }^{(A\mu rt_0)}(\stackrel{}{y}t)R^{(A\mu rt_0)}(\stackrel{}{x})$$
(26)
as a stochastic estimator for spatial all-to-all propagator matrix elements.

We have used the BiCGStab algorithm Frommer et al. (1994); van der Vorst (1992) for solving (24).

Operator smearing Alexandrou et al. (1994) is a technique for enhancing the amplitude of the ground state in a correlation function. We smear the light quark fields only, defining iteratively

$`q_A^{\{0\}}(\stackrel{}{x}t)`$ $`=`$ $`q_A(x)`$ (27a)
$`q_A^{\{m\}}(\stackrel{}{x}t)`$ $`=`$ $`{\displaystyle \underset{B}{}}{\displaystyle \underset{\stackrel{}{y}}{}}𝒦_{AB}(\stackrel{}{x},\stackrel{}{y})q_B^{\{m1\}}(\stackrel{}{y}t),`$ (27b)
with $`m`$, and the matrix

$`𝒦_{AB}(\stackrel{}{x},\stackrel{}{y})`$ $`=`$ $`\delta _{AB}\delta _{\stackrel{}{x},\stackrel{}{y}}+\alpha {\displaystyle \underset{\mu =1}{\overset{3}{}}}[\text{}U_{\mu ,AB}(\stackrel{}{x}t)\delta _{\stackrel{}{x},\stackrel{}{y}\widehat{\mu }}`$ (28)
$`+U_{\mu ,AB}^{}(\stackrel{}{y}t)\delta _{\stackrel{}{x},\stackrel{}{y}+\widehat{\mu }}].`$
The real number $`\alpha `$ and the maximum value $`M`$ for iterations $`m=0\mathrm{}M`$ are parameters. We have also used APE type Albanese et al. (1987) fuzzy link variables $`USU(3)`$ in (28) with the same parameter set. The above iterative prescription translates directly to the random source and solution vectors $`R`$ and $`X`$, i.e. replacing $`RR^{\{M\}}`$ and $`XX^{\{M\}}`$ in (26) yields the propagator $`G^{\{M\}}`$ for smeared fermion fields. For more technical details see Fiebig et al. (2000).

Smearing was employed as a matter of course. In fact the operator (16) is designed with excited states of the two-meson system in mind. Thus light-quark field smearing and link variable fuzzing may be of limited value for the task at hand. We have thus chosen somewhat conservative parameter values $`\alpha =1.4`$ and $`M=4`$ for the current simulation.

## IV Spectral density analysis

We now turn to the problem of extracting spectral information from the correlation matrix (17,19). It has the decomposition

$`C_{ij}(t,t_0)`$ $`=`$ $`{\displaystyle \underset{n0}{}}\varphi _{ni}\varphi _{nj}^{}e^{E_n(tt_0)}`$ (29)
$`\varphi _{ni}`$ $`=`$ $`0|\mathrm{\Phi }_i^{}(t_0)|n,`$ (30)
where $`|n`$ denotes a complete set of states with energies $`E_n`$, some of them may be negative in a lattice simulation. The sum in (29) is truncated in practice, $`n=1\mathrm{}N`$. The truncation $`N`$ is determined by the physics of the system (ultimately the lattice action), the simulation parameters (most importantly the lattice energy cutoff), and by the choice of operators. The matrix elements (30) can be interpreted as components of $`N`$ $`M`$-dimensional vectors $`\varphi _n`$, where $`M`$ is the size of the correlation matrix $`C(t,t_0)`$. Following Lüscher and Wolff Lüscher and Wolff (1990) we diagonalize the correlation matrix separately on each time slice, say

$`{\displaystyle \underset{j=1}{\overset{M}{}}}C_{ij}(t,t_0)v_{mj}(t,t_0)`$ $`=`$ $`C_m(t,t_0)v_{mi}(t,t_0)`$ (31a)
$`{\displaystyle \underset{i=1}{\overset{M}{}}}v_{mi}^{}(t,t_0)v_{m^{}i}(t,t_0)`$ $`=`$ $`\delta _{mm^{}},`$ (31b)
denoting the eigenvectors by $`v_m`$ and the eigenvalues by $`C_m`$, $`m=1\mathrm{}M`$. Now assume the following:

1. The energies $`E_n`$ in (29) are non degenerate, and ordered $`E_1<E_2<\mathrm{}<E_N`$.
2. The vectors $`\varphi _n`$ are linearly independent, this implies $`NM`$.
3. There is a $`t_C`$ such that for all $`tt_C`$ the eigenvalues are ordered $`C_1C_2\mathrm{}C_M`$.

Under those conditions a theorem proven in Lüscher and Wolff (1990) states that for all $`n=1\mathrm{}N`$

$$\underset{t\mathrm{}}{lim}C_n(t,t_0)=Z_ne^{E_n(tt_0)}[1+O(e^{\mathrm{\Delta }E_n(tt_0)})],$$
(32)
where $`Z_n>0`$ and $`\mathrm{\Delta }E_n=\mathrm{min}_{n^{}n}\{|E_n^{}E_n|\}`$ is the distance to the energy closest to $`E_n`$. We are interested in the structure of the spectral representation of $`C_n(t,t_0)`$. Toward this end, applying (31b) to (31a) and then inserting (29), one obtains

$$C_m(t,t_0)=\underset{n0}{}|\underset{i=1}{\overset{M}{}}v_{mi}^{}(t,t_0)\varphi _{ni}|^2e^{E_n(tt_0)}.$$
(33)
In addition to items 1.–3. above we will also assume:

1. In the large-$`t`$ limit the first $`N^{}`$ of the eigenvectors of (31a,31b) converge in the sense that in

$$\underset{t\mathrm{}}{lim}v_{mi}(t,t_0)=\eta _m(t,t_0)v_{mi},m=1\mathrm{}N^{},$$
(34)
the vectors $`v_m`$ are constant, being multiplied by a $`t`$-dependent phase factor, $`|\eta _m(t,t_0)|=1`$.

Note that $`N^{}NM`$. Thus, using (34) and (32) in (33) we arrive at

$$\underset{n0}{}|\underset{i=1}{\overset{M}{}}v_{mi}^{}\varphi _{ni}|^2e^{E_n(tt_0)}=Z_me^{E_m(tt_0)},$$
(35)
for $`m=1\mathrm{}N^{}`$. Since all $`E_n`$ are different the exponentials are linearly independent functions of $`t`$, hence

$$|\underset{i=1}{\overset{M}{}}v_{mi}^{}\varphi _{ni}|^2=Z_n\delta _{mn},m=1\mathrm{}N^{}.$$
(36)
The square root of this is

$$\underset{i=1}{\overset{M}{}}v_{mi}^{}\varphi _{ni}=\zeta _m\sqrt{Z_m}\delta _{mn},m=1\mathrm{}N^{},$$
(37)
with $`|\zeta _m|=1`$. The eigenvectors of (31a,31b) satisfy a completeness relation in $`M`$-dimensional space. We split it into two parts

$$\underset{m=1}{\overset{N^{}}{}}v_{mj}(t,t_0)v_{mi}^{}(t,t_0)+\underset{k>N^{}}{\overset{M}{}}v_{kj}(t,t_0)v_{ki}^{}(t,t_0)=\delta _{ji}.$$
(38)
According to the assumption (34) all terms in the first sum will individually converge in the large-$`t`$ limit. The individual terms in the second sum might not, however, it must of course become $`t`$-independent as a whole,

$$\mathrm{\Pi }_{ji}=\underset{t\mathrm{}}{lim}\underset{k>N^{}}{\overset{M}{}}v_{kj}(t,t_0)v_{ki}^{}(t,t_0).$$
(39)
Clearly the projector $`\mathrm{\Pi }=\mathrm{\Pi }^2=\mathrm{\Pi }^{}`$ is orthogonal on the space defined by the span of the $`v_m`$,

$$\underset{i=1}{\overset{M}{}}\mathrm{\Pi }_{ji}v_{mi}=0,\mathrm{for}m=1\mathrm{}N^{}.$$
(40)
Thus, as $`t\mathrm{}`$, equation (38) assumes the form

$$\underset{m=1}{\overset{N^{}}{}}v_{mj}v_{mi}^{}+\mathrm{\Pi }_{ji}=\delta _{ji}.$$
(41)
Finally, operating with $`_{m=1}^N^{}v_{mj}\mathrm{}`$ on both sides of (37) and then using (41) gives

$$\varphi _n=\zeta _n\sqrt{Z_n}v_n+\mathrm{\Pi }\varphi _n,\mathrm{for}n=1\mathrm{}N^{}.$$
(42)
This result relates the matrix elements $`\varphi _{ni}=0|\widehat{\mathrm{\Phi }}_i^{}(t_0)|n`$ to the solutions of the $`t`$-dependent eigenvalue problem (31a, 31b) in the large-$`t`$ limit. An immediate consequence, derived by taking the square of (42), is $`\varphi _n^{}\varphi _n=Z_n+\varphi _n^{}\mathrm{\Pi }\varphi _n`$, or

$$\underset{i=1}{\overset{M}{}}|n|\mathrm{\Phi }_i(t_0)|0|^2=Z_n+\mathrm{\Pi }\varphi _n^2,\mathrm{for}n=1\mathrm{}N^{}.$$
(43)
Thus $`Z_n`$ is a lower bound on the total probability for (incoherent) excitations by a set of operators $`\mathrm{\Phi }_i,i=1\mathrm{}M`$, into a certain state $`|n`$. In principle the value of $`N^{}`$ can be computed from (34). In practice this is hard to accomplish, because the components of the eigenvectors fluctuate strongly. If the set of $`M`$ operators couples to all $`N`$ available physical states it seems reasonable to expect that $`N^{}=N`$. In case that $`N<M`$ we have used more operators than quantum states are available in the system. Then the projector term compensates for over counting the physical degrees of freedom.

Aside from (43) there is an alternative way of interpreting the $`Z_n`$. Based on (34) define the meson-meson fields

$$\mathrm{\Psi }_m(t)=\underset{i=1}{\overset{M}{}}v_{mi}\mathrm{\Phi }_i(t)\mathrm{for}m=1\mathrm{}N^{},$$
(44)
and consider the $`N^{}\times N^{}`$ correlator matrix

$$D_{mm^{}}(t,t_0)=\widehat{\mathrm{\Psi }}_m^{}(t)\widehat{\mathrm{\Psi }}_m^{}(t_0).$$
(45)
Inserting (44) into (45), then using (29) and (37), it is straightforward to show that

$$D_{mm^{}}(t,t_0)=\delta _{mm^{}}Z_me^{E_m(tt_0)}.$$
(46)
On the other hand, starting from (45), the diagonal elements have the standard decomposition

$$D_{mm}(t,t_0)=\underset{n0}{}|n|\mathrm{\Psi }_m(t_0)|0|^2e^{E_n(tt_0)}.$$
(47)
Comparison of (46) and (47), using linear independence of the exponentials again, then gives

$$|n|\mathrm{\Psi }_m(t_0)|0|^2=\delta _{nm}Z_n\mathrm{for}n,m=1\mathrm{}N^{}.$$
(48)
Thus $`Z_n=|n|\mathrm{\Psi }_n(t_0)|0|^2`$ is the excitation probability of the state $`|n`$ due to an operator $`\mathrm{\Psi }_n(t)`$ that is optimal, with regard to $`|n`$, within the linear space of the original set $`\mathrm{\Phi }_i(t),i=1\mathrm{}M`$.

## V Discussion of results

The significance of the above is that it suggests an analysis strategy for the spectral features of the two-meson system: Diagonalize the time correlation matrix (19) on each time slice,

$`\left(\begin{array}{cc}C_{11}(t,t_0)& C_{12}(t,t_0)\\ C_{21}(t,t_0)& C_{22}(t,t_0)\end{array}\right)\left(\begin{array}{c}v_{m1}(t,t_0)\\ v_{m2}(t,t_0)\end{array}\right)`$ $`=`$ (56)
$`C_m(t,t_0)\left(\begin{array}{c}v_{m1}(t,t_0)\\ v_{m2}(t,t_0)\end{array}\right),`$
$`m=1,2`$. Then view each eigenvalue as a separate time correlation function $`C_m(t,t_0)`$ subject to spectral analysis. In particular, seek to extract a spectral representation of the form $`𝑑\omega \rho (\omega )e^{\omega (tt_0)}`$, see (1). This can be done by Bayesian inference Jarrell and Gubernatis (1996); Box and Tiao (1973); Fiebig (2002a). The expected structure of the spectral density $`\rho `$ is a linear combination of peaks, see (2). In an actual numerical simulation those will have finite widths. More importantly, because the main results of Sec. IV, including (43) and (48), require the limit $`t\mathrm{}`$ only the lowest-energy peak of each $`C_m(t,t_0)`$ is significant. Physical information that can be extracted from the peak includes the energy $`E_n`$ of the state $`|n`$, and the strength $`Z_n`$ of excitations by means of the set of operators employed.

The eigenvalue correlators $`C_m(t,t_0),m=1,2`$, of (56) are displayed in Fig. 2 for four meson-meson relative distances $`r=1\mathrm{}4`$. Due to the non-local structure of the two-meson operators, leading to loop-loop correlations, somewhat noisy data are unavoidable. Link variable fuzzing and operator smearing was used to be able to use ‘earlier’ time slices, recall Sect. III. Not surprisingly, the lower eigenvalue $`C_2(t,t_0)`$ is more susceptible to noise than $`C_1(t,t_0)`$. Backward going propagation is present, though suppressed by four to five orders of magnitude.333Note that (anti)periodic boundary conditions in the lattice action do not in general translate into similar behavior of the correlation functions. Here, these are some combination of ‘$`\mathrm{cosh}`$’ and ‘$`\mathrm{exp}`$’ behavior akin to local meson and loop like operators. By inspection of Fig. 2 it is apparent that the condition $`C_1(t,t_0)C_2(t,t_0)`$ is fulfilled for $`tt_0`$. This is directly obvious for most of the forward temporal range, say $`0<t<20`$, and for $`20t<30`$ it can be inferred from the global fits to $`C_2(t,t_0)`$. Indicative of a violation would be the possibility of being able to smoothly connect two sets of consecutive data points in such a way that the two resulting curves would intersect.444For example $`c_1(t,t_0)=a_1e^{b_1(tt_0)}`$ and $`c_2(t,t_0)=a_2e^{b_2(tt_0)}`$ intersect at $`t_Ct_0=\mathrm{ln}(a_2/a_1)/(b_2b_1)`$. The observation simply means that $`t_Ct_0`$. We thus observe that the above ordering of eigenvalues, as stated in Sec. IV, is satisfied.

For asymptotic times individual components of the eigenvectors in (56) show large statistical fluctuations. In Fig. 3 a typical example of the (complex) values of $`z=v_{n1}(t,t_0)v_{n2}^{}(t,t_0)`$ is displayed. This quantity is relevant to the assumption (34). In the region of forward propagation, say $`t20`$, the real part of $`z`$ exhibits increasing stochastic errors as $`t`$ becomes large, but stabilizes within those bounds. A similar statement can be made for the imaginary part of $`z`$, adding that it is consistent with zero. This justifies (34) within the limitations given by the quality of the data.555Note that (34) is stated for the forward propagation region of the correlation functions. It may have to be modified for backward propagation.

In Fig. 4 we show the magnitudes $`|v_{11}(t,t_0)|^2`$ of eigenvector components versus $`t`$ for four relative distances. Substantial noise notwithstanding we have strictly applied the prescription of Sect. IV and computed the eigenvalue correlators following (31a,31b) verbatim

$$C_m(t,t_0)=\underset{i,j=1}{\overset{2}{}}v_{mi}^{}(t,t_0)C_{ij}(t,t_0)v_{mj}(t,t_0),$$
(57)
for $`m=1,2`$. The corresponding spectral analysis is discussed below.

Stochastic fluctuations of the eigenvector components are enhanced by diagonalizing (56) separately on each time slice. An attractive alternative is to replace the eigenvector components $`v_{mi}(t,t_0)`$ in (57) with time averaged components $`\overline{v}_{mi}`$ taken at asymptotic times. Within the present data set the available time window is $`10t20`$, excluding backward propagation. We have made fits to $`|v_{11}(t,t_0)|^2`$ with a constant model $`s`$ in the time slice range $`11t17`$, and also with a control set in $`11t19`$. The results are listed in Tab. 1. Specific values for $`s`$ give rise to, time independent, vectors

$$\overline{v}_1=\left(\begin{array}{c}\sqrt{s}\\ \sqrt{1s}\end{array}\right)\overline{v}_2=\left(\begin{array}{c}\sqrt{1s}\\ \sqrt{s}\end{array}\right),$$
(58)
up to an arbitrary phase. Thus we have also performed a spectral analysis based on the correlators

$$\overline{C}_m(t,t_0)=\underset{i,j=1}{\overset{2}{}}\overline{v}_{mi}^{}C_{ij}(t,t_0)\overline{v}_{mj},$$
(59)
for $`m=1,2`$. Besides smoothing out statistical fluctuations of the eigenvector components the advantage of (59) is that the asymptotic form of the correlator (57) is now used on all time slices. This should improve the signal derived from the Bayesian spectral analysis which makes use of data on all times slices.

The values of $`s`$ inform us about operator mixing as the relative distance $`r`$ changes. They are a measure of how strongly the operator $`\mathrm{\Phi }_1`$ couples to a meson-meson system in the ground state. As Fig. 5 shows this measure distinctly decreases from about $`1.0`$ as $`r`$ becomes smaller. Since $`\mathrm{\Phi }_1`$ is designed to test quark exchange degrees of freedom we see that those gradually become less important at smaller relative distances. By the same token $`1s`$ measures the ground state coupling strength of $`\mathrm{\Phi }_2`$. Thus we learn that quark exchange is the dominant interaction mechanism at large $`r`$, while gluon exchange gradually takes over as $`r`$ decreases. A glance at Fig. 5 reveals that the mechanisms become balanced ($`s0.5`$) at distances $`r`$ somewhat less than $`0.5`$, or $`0.25\mathrm{fm}`$ in physical units.

The Bayesian analysis of time correlation functions has been extensively discussed in Fiebig (2002a), using the same lattice (raw) data. We here briefly state the main points for coherence of presentation, but otherwise refer the reader to Fiebig (2002a). We expect the lattice data to fit the model (1), or rather its discretized form

$$F(\rho |t,t_0)\underset{k=K_{}}{\overset{K_+}{}}\rho _ke^{\omega _k(tt_0)}$$
(60)
with $`\rho _k=\mathrm{\Delta }\omega \rho (\omega _k)`$ and $`\omega _k=\mathrm{\Delta }\omega k`$.666 In Ref. Fiebig (2002a) backward ($`\omega <0`$) and forward ($`\omega >0`$) going exponentials were normalized differently, giving rise to a modified spectral function $`\rho _T`$. However, $`\rho _T(\omega )=\rho (\omega )`$ for $`\omega 0`$. Since we here only present $`\omega 0`$ data there is no need to distinguish between $`\rho _T`$ and $`\rho `$. The objective is to compute the spectral density function $`\rho `$. Toward this end consider the functional

$$W[\rho ]=\frac{1}{2}\chi ^2[\rho ]\alpha 𝒮[\rho ].$$
(61)
The first term involves the usual $`\chi ^2`$-distance between the lattice data $`C_m(t,t_0)`$ and the model $`F(\rho |t,t_0)`$, computed with the full covariance matrix derived from gauge field configuration statistics, and

$$𝒮[\rho ]=\underset{k=K_{}}{\overset{K_+}{}}\left(\rho _km_k\rho _k\mathrm{ln}\frac{\rho _k}{m_k}\right),$$
(62)
is the (Shannon-Jaynes) entropy Shannon and Weaver (1949); Jaynes (1957). It plays the role of a Bayesian prior Jarrell and Gubernatis (1996). The configuration $`m=\{m_k:K_{}kK_+\}`$ is called the default model. Another parameter in (61) is the entropy strength $`\alpha `$. The optimization problem $`\chi ^2[\rho ]=\mathrm{min}`$ is ambiguous because in practice a reasonable resolution $`\mathrm{\Delta }\omega `$ will result in the number $`K_+K_{}+1`$ of fit parameters $`\rho _k`$ being much larger than the number $`T`$ of simulation data. However, $`W[\rho ]`$ has a unique absolute minimum Jarrell and Gubernatis (1996). From the viewpoint of Bayesian statistics $`\rho `$ is interpreted as a random variable subject to a certain probability distribution (posterior probability). The most likely $`\rho `$ is the one that minimizes $`W[\rho ]`$. Finding the minimum within this framework is known as the maximum entropy method (MEM). It is designed to minimize the information not supported by the data. Loosely speaking, we seek to minimize $`\chi ^2[\rho ]`$ while assuming as little information as possible about the spectral density $`\rho `$. To solve the optimization problem $`W[\rho ]=\mathrm{min}`$, probabilistic methods seem closest in spirit to the Bayesian stochastic interpretation of $`\rho `$. We have thus employed simulated annealing Kirkpatrick (1984), or cooling, based on the partition function

$$Z_W=[d\rho ]e^{\beta _WW[\rho ]}.$$
(63)
In Fiebig (2002a) we have studied the dependence of the resulting spectral density on (i) the entropy strength parameter $`\alpha `$, (ii) the default model $`m`$, and (iii) the annealing start configuration. It was found that $`\rho `$ was essentially independent of $`\alpha `$ and of $`m`$, and that the expectation values of low $`\omega `$ moments were independent of the annealing start within errors native to the lattice data set.

Results of the MEM analysis are presented in Fig. 6 for the eigenvalue correlators (57) and in Fig. 7 for the asymptotic stabilized correlators (59). As discussed in Fiebig (2002a) the parameter choices are not critical. Specifically, the entropy strength is $`\alpha =5.0\times 10^7`$, and the default model is constant with $`m_k=10^{12},K_{}kK_+`$. The annealing schedule is given in Fiebig (2002a). All spectral densities are averages over eight random annealing starts. The $`\omega `$ discretization is set by $`\mathrm{\Delta }\omega =0.02`$, and $`K_{}=100`$, $`K_+=+200`$. Those numbers reflect the lattice design, like the energy cutoff $`a_{t}^{}{}_{}{}^{1}`$, and other considerations, see Fiebig (2002a). The $`\omega `$ interval is larger and the resolution much finer compared to the discretization used in Fiebig (2002a). Note that the entire spectral mass range $`2.0\omega +4.0`$, including backward going propagation, is utilized in the spectral analysis. This is true of both the eigenvalue correlators and the asymptotic stabilized correlators. Most of the spectral structure, however, is invisible on the linear scales used in Figs. 6 and 7, particularly for $`\omega 0`$ (backward propagation).

As discussed in Sect. IV, in the limit $`t\mathrm{}`$ only the lowest $`\omega `$ peak from each correlator $`C_m(t,t_0),m=1\mathrm{}M`$, should be used to extract physical information. We will refer to those as primary peaks. Suppose primary peaks are seen in the spectral densities belonging to $`C_n(t,t_0)`$, for $`n=1\mathrm{}NM`$. We loosely characterize those by $`\delta _n=\{\omega :\omega \mathrm{peak}\mathrm{\#}\mathrm{n}\}`$. In Fiebig (2002a) it was argued that low $`\omega `$ moments of the spectral densities $`\rho (\omega )`$ can be reliably extracted. Specifically, these are the peak volume $`Z_n`$, the mean energy $`E_n`$, and the width $`\mathrm{\Delta }_n`$ of the peak

$`Z_n`$ $`=`$ $`{\displaystyle _{\delta _n}}𝑑\omega \rho (\omega )`$ (64)
$`E_n`$ $`=`$ $`Z_n^1{\displaystyle _{\delta _n}}𝑑\omega \rho (\omega )\omega `$ (65)
$`\mathrm{\Delta }_n^2`$ $`=`$ $`Z_n^1{\displaystyle _{\delta _n}}𝑑\omega \rho (\omega )\left(\omega E_n\right)^2.`$ (66)
Note that the theorem (32) establishes the peak volume (64) as identical with the factors $`Z_n`$ introduced in Sect. IV, the caveat being that numerically extracted peaks have finite widths. In particular, the peak volumes (64) have the interpretation given by (43), or by (48) as optimal excitation probabilities.

In Figs. 6 and 7 the (lowest mass) primary peaks clearly dominate both ground state spectral density functions, $`C_1`$ and $`\overline{C}_1`$, for each $`r=1\mathrm{}4`$. With reference to (43) the secondary peaks at larger mass may indicate that $`\mathrm{\Pi }\varphi _1^2>0`$, their volumes are smaller though. It is much harder to make out a distinct peak structure in the spectral density functions for the excited state eigenvalue correlators $`C_2`$. We attribute this fact to strong statistical fluctuations of the eigenvector components $`v_{mi}(t,t_0)`$, spoiling the signal. The asymptotic stabilized excited state correlators $`\overline{C}_2`$ alleviate that problem, see Fig. 7. The spectral peaks of $`\overline{C}_2`$ are broad, for $`r>1`$. This is an indication that the lattice data support them with only a few consecutive points in the time correlation function. In other words, there is not enough information in the data for distinct narrow peaks to develop against the entropy background. Although the peaks are wider, lower mass primary peaks are clearly distinguishable for $`r=1\mathrm{}4`$ in Fig. 7. In Tabs. 2 and 3 we list the volumes, energies, and widths of all primary peaks of Figs. 6 and 7. For the excited states the $`C_2`$ data do not always clearly define a primary peak. The corresponding numbers in Tab. 2 use $`\omega `$ cuts of 1.32, 1.58, 3.20 and 2.80, for $`r=`$ 1,2,3 and 4, respectively. On the other hand, the $`\overline{C}_2`$ data originating with the asymptotic stabilized correlators provide us with a much improved picture. There, by inspection of Fig. 7, the $`\omega `$ cuts are 1.28, 1.72, 2.20 and 2.46, for $`r=1\mathrm{}4`$, respectively. Most of the subsequent discussion therefore uses results of the $`\overline{C}`$ analysis.

With an appropriate annealing schedule cooling fluctuations at the final ‘temperature’ can be made negligible. The size of the statistical error from $`N_U=708`$ gauge configurations is comparable to uncertainties from the MEM analysis. This was looked at in Fiebig (2002a). Those uncertainties are introduced through using different start configurations $`\rho `$ in the annealing process. There is no magic way to eliminate those, nor would it be desirable, because they test a property of the data set. Although theoretically $`W[\rho ]`$ has a unique absolute minimum sole knowledge of its location, say $`\rho _{\mathrm{min}}`$, is deficient. It should be supplemented by having some notion about the shallowness, or as the case may be, the distribution of local minima close to $`\rho _{\mathrm{min}}`$ with values $`W[\rho ]`$ not much different from $`W[\rho _{\mathrm{min}}]`$ which the annealing algorithm may settle into. Those manifest themselves in micro fluctuations (on a scale of $`\mathrm{\Delta }\omega `$) of the spectral density functions. Local minima close to $`\rho _{\mathrm{min}}`$ appear to be numerous. The numbers in Tabs. 2 and 3 are from averages over eight annealing start configurations, and the uncertainties are the corresponding standard deviations. They are indicative of the spread of local minima of $`W[\rho ]`$ in the vicinity of $`\rho _{\mathrm{min}}`$ and, ultimately of the uncertainty of the results.

In view of the broader peaks in Figs. 6 and 7 it may be argued that the widths $`\mathrm{\Delta }_2`$ and $`\overline{\mathrm{\Delta }}_2`$ supersede the standard deviations derived from annealing start configurations as a useful indicator for the uncertainties of the average energies $`E_2`$ and $`\overline{E}_2`$, see Tabs. 2 and 3. In fact there are at least three types of indicators: (i) the gauge field configuration statistical error, (ii) the annealing start configuration standard deviation, and (iii) the peak width. They all point at different aspects of the uncertainty of the $`E_n`$. For example, the $`\mathrm{\Delta }_n`$ convey the aspect of information content of the data in the sense of Shannon and Weaver (1949). A deeper discussion is beyond the scope of this work. For the current results the gauge field configuration statistical error is typically smaller than 2% and thus much less than the other two measures of uncertainty. Then, the annealing start configuration standard deviation in most cases is much less than the peak width. We here adopt the point of view that, since Bayesian inference is built on the information content of the data, the peak width is the appropriate measure of uncertainty. However, for reference to (i) refer to Fiebig (2002a), and we include (ii) in the results below as appropriate.

The ground state energies $`\overline{E}_1`$ of the meson-meson system are plotted in Fig. 8 together with the mass $`2m`$ and the error band for non-interacting mesons. The mass $`m=0.4676`$ and the gauge configuration statistical error $`\mathrm{\Delta }m=0.0075`$ are those from $`m_{\mathrm{eff},0}`$ in Fiebig (2002a). Error bars on the data points are the peak widths $`\overline{\mathrm{\Delta }}_1`$. Although the errors extend well below $`V=0`$, a systematic tendency for repulsion, averaging about $`90\mathrm{M}\mathrm{e}\mathrm{V}`$ above $`2m`$, is apparent. These features are consistent with previous lattice studies of the pseudoscalar meson-meson system in the same isospin channel Fiebig et al. (2000).

We also display in Fig. 8 results from Michael and Pennanen (1999) which relate to the current discussion. In Michael and Pennanen (1999) operators are characterized by the isospin and spin symmetry of the light quarks. Because in the present work all the light quarks have the same flavor, only the symmetric isospin combination is relevant; in the notation of Michael and Pennanen (1999) this is $`I=1`$. Figure 8 shows their $`S=0`$ and $`S=1`$ results from two lattices. The uncertainties are gauge configuration statistical errors and are thus much smaller than our peak width error bars. The energy values populate a $`\pm 50\mathrm{M}\mathrm{e}\mathrm{V}`$ band around $`V=0`$. It should be noted that the authors of Michael and Pennanen (1999) go through extraordinary length to obtain the ‘best’ ground state energies possible by means of diagonalizing a matrix of operators (similar to $`\mathrm{\Phi }_1`$) with different fuzzing levels of the gauge fields. This goes beyond the widely used practice to employ a few iterative steps, say $`N`$, of gauge link fuzzing Albanese et al. (1987) and quark field smearing Alexandrou et al. (1994). A reasonable choice for the number of iterations is to make $`Na_s`$ about the radius of the hadron considered. (In our case, with $`a_s=0.25\mathrm{fm}`$, this would be $`N=4`$, for example.) The result is a spatially extended operator, spreading across a hadron volume, that makes correlation functions assume asymptotic behavior at ‘earlier’ time slices. This is the procedure adopted in the current work. This practice also has the side effect of lowering the ground state energy obtained from the simulation because, numerically, contaminations from excited states are reduced. This effect is of course enhanced if the ground state is extracted from diagonalizing a correlation matrix from operators of different fuzzing levels, as employed in Michael and Pennanen (1999). The size of the effect can be seen in Fig. 8. We note that the physics goal of the authors of Michael and Pennanen (1999) is to investigate if the $``$ system can possibly be bound. In that context even small effects on the ground state energy are vital, so the elaborate matrix fuzzing procedure is justified. However, the physics goals of this work are completely different. Studying the interaction mechanism rests on a comparison of ground and excited state energies, as $`r`$ changes. Because the excited states energy levels are substantially larger, tiny shifts in the ground state energies are irrelevant.

This situation is reflected in Fig. 9 where we show excited state and ground state energies together, versus the relative meson-meson distance. The energy $`\overline{E}_2`$ drops considerably as $`r`$ decreases. At large $`r`$, according to Tab. 1 and Fig. 5, the operator $`\mathrm{\Phi }_2(t)`$ is mainly responsible for the values of $`\overline{E}_2`$. Bearing in mind that $`\mathrm{\Phi }_2(t)`$, see (16) and Fig. 1, involves color charges separated by a distance $`r`$, which should be confined, we have fitted the data with the model $`y=a/r+br+c`$, where the parameter $`b`$ was fixed to the string tension, $`\sqrt{b}=0.44\mathrm{GeV}`$ or $`b=0.968\mathrm{GeV}/\mathrm{fm}`$. The remaining fit returns $`a=0.30(14)\mathrm{GeV}\mathrm{fm}`$, $`c=1.09(52)\mathrm{GeV}`$, where the uncertainties are variances. The resulting curve is shown in Fig. 9. Around $`r0.2\mathrm{fm}`$ a level crossing between the ground and excited state energies of the meson-meson system apparently occurs. From the viewpoint of an adiabatic potential this distance defines the transition point (avoided level crossing) between a weakly repulsive interaction mostly mediated by quark exchange, and strongly attractive interaction from gluon exchange degrees of freedom. The transition distance is consistent with Fig. 5, where at a point somewhat smaller than $`r0.5`$ marks equal $`s`$ values for the $`\mathrm{\Phi }_1`$ and $`\mathrm{\Phi }_2`$ components of the asymptotic eigenvectors $`v_1`$ and $`v_2`$. The picture emerging from the current results thus is that of heavy-light pseudoscalar mesons having a radius of about $`R0.1\mathrm{fm}`$, with a sharp boundary, as probed by their mutual interaction. Note that the dramatic change of the adiabatic potential, at the distance of the avoided level crossing nevertheless comes from a smooth transition between the interaction mechanisms, as manifest in Fig. 5,

By extrapolation it appears that below $`r0.2\mathrm{fm}`$ the interaction turns attractive. This observation is in line with Michael and Pennanen (1999) where strong attraction at $`r=0`$ due to gluonic effects is observed, and also with results in Mihály et al. (1997); Fiebig et al. (2000). A simple explanation is suggested by the $`SU(3)`$ color content of the two-body operators. Using standard nomenclature Lichtenberg (1978); Close (1979) we note that only the singlet from

$`\overline{\mathrm{𝟑}}\mathrm{𝟑}`$ $`=`$ $`\mathrm{𝟖}\mathrm{𝟏}`$ (67)
is used in the construction of (15). Thus, generically, the color-source structure of the operator $`\mathrm{\Phi }_1`$ is

$$\mathrm{\Phi }_1^{(\mathrm{𝟏})}(1,2)\varphi ^{(\mathrm{𝟏})}(1)\varphi ^{(\mathrm{𝟏})}(2),$$
(68)
where 1 and 2 denote the two color sources. The decomposition of the product of two color octets

$`\mathrm{𝟖}\mathrm{𝟖}`$ $`=`$ $`\mathrm{𝟐𝟕}\mathrm{𝟏𝟎}\mathrm{𝟖}\mathrm{𝟖}\overline{\mathrm{𝟏𝟎}}\mathrm{𝟏}`$ (69)
also contains a singlet and therefore mixes with (68). The construction of (16) involves two gauge link paths. Thus in some sense the color-source structure of $`\mathrm{\Phi }_2`$ is, schematically, described by the Clebsch-Gordon series

$$\mathrm{\Phi }_2^{(\mathrm{𝟏})}(1,2)\underset{i,j}{}\left(\begin{array}{c}\mathrm{𝟏}\end{array}|\begin{array}{cc}\mathrm{𝟖}& \mathrm{𝟖}\\ i& j\end{array}\right)\varphi _i^{(\mathrm{𝟖})}(1)\varphi _j^{(\mathrm{𝟖})}(2).$$
(70)
The interaction energies for one-gluon exchange in those states are proportional to the expectation values of $`2\stackrel{}{F}(1)\stackrel{}{F}(2)=(\stackrel{}{F}(1)+\stackrel{}{F}(2))^2\stackrel{}{F}(1)^2\stackrel{}{F}(2)^2`$ where $`\stackrel{}{F}^2`$ is an $`SU(3)`$ Casimir operator Close (1979). A simple calculation gives

$`2\stackrel{}{F}(1)\stackrel{}{F}(2)_{\mathrm{\Phi }_1^{(\mathrm{𝟏})}}`$ $`=`$ $`0`$ (71)
$`2\stackrel{}{F}(1)\stackrel{}{F}(2)_{\mathrm{\Phi }_2^{(\mathrm{𝟏})}}`$ $`=`$ $`6,`$ (72)
which indicates (possibly strong) attractive interaction at small $`r`$ such as described by excitations due to $`\mathrm{\Phi }_2`$. This is only a schematic picture, but is it appears to be consistent with the results of the lattice simulation.777Note that $`r=0`$ is a special case since there is no distinction in the color structure of (15) and (16), in fact $`\mathrm{\Phi }_1=\mathrm{\Phi }_2`$ at that point.

Finally, we turn to the transition matrix elements. Equations (43) and (48) relate the peak volumes $`Z_n`$ to (unnormalized) excitation probabilities of the states $`|n`$,

$`Z_n`$ $`=`$ $`\left|n|v_{n1}\mathrm{\Phi }_1(t_0)+v_{n2}\mathrm{\Phi }_2(t_0)|0\right|^2`$ (73a)
$``$ $`\left|n|\mathrm{\Phi }_1(t_0)|0\right|^2+\left|n|\mathrm{\Phi }_2(t_0)|0\right|^2.`$ (73b)
The coefficients $`v_{ni}`$ in (73a), being eigenvector components, ensure that $`Z_n`$ is maximal within the linear space of the available operators $`\mathrm{\Phi }_i(t_0),i=1,2`$, see Sect. IV. Accordingly we interpret $`Z_n`$, or

$$z_n=Z_n/(Z_1+Z_2),$$
(74)
as measures for the effectiveness of the set of operators to excite the state $`|n`$. The corresponding results, compiled from Tab. 3, are displayed in Fig. 10. The decrease of $`\overline{Z}_1`$ as $`r`$ becomes smaller means that the operators $`\mathrm{\Phi }_1`$ and $`\mathrm{\Phi }_2`$ start failing to capture some physics of the two-meson ground state at small relative distance. This is a sign of an emerging interaction mechanism that is is not well represented by $`\mathrm{\Phi }_{1,2}`$. Other operators should eventually be added to the set. From Fig. 10 we also see that the operator set $`\mathrm{\Phi }_{1,2}`$ is about 60–70% effective in creating the two-meson ground state, and the rest is left to excited state creation. It is possible that gauge link fuzzing and quark field smearing is responsible for the enhanced ground state presence. The normalized excitation rates are essentially independent of $`r`$, as Fig. 10 shows.

## VI Summary and conclusion

In an attempt to learn about the physics of hadronic interaction we have studied the spectrum of a meson-meson system consisting of two light and two static quarks. The static quarks define the relative distance $`r`$ between the mesons. We have supplemented the commonly used local two-meson operator $`\mathrm{\Phi }_1(t)`$ with a nonlocal one $`\mathrm{\Phi }_2(t)`$ which is built from spatially extended color singlets. The eigenvalues of the corresponding $`2\times 2`$ time correlation matrix were used in the spectral analysis thus allowing operator mixing. We have employed the maximum entropy method (MEM), a form of Bayesian inference, to extract the spectral densities of the ground and the excited states from the time correlation matrix.

The spectral analysis yields both energies of the ground and excited states of the meson-meson system as a function of the relative distance $`r`$. While, at large $`r`$, the ground state energy is weakly repulsive and flat, the excited state level is strongly $`r`$ dependent and decreases substantially as $`r`$ becomes small. By way of extrapolation, the salient feature of the spectrum is a level crossing of the ground and excited states at about $`r0.2\mathrm{fm}`$. There, the adiabatic potential changes from weak repulsion to strong attraction. Analysis of the eigenvalues of the correlation matrix, at asymptotic times, reveals that the interaction mechanism gradually changes from being dominated by quark degrees of freedom (quark-antiquark exchange) at large $`r`$ to gluon degrees of freedom( gluon pair exchange) as $`r`$ becomes smaller. This view is based on monitoring $`\mathrm{\Phi }_{1,2}`$ operator mixing as a function of $`r`$.

A side aspect of the current work relates to the various types of errors emerging in the analysis procedure. Between the gauge configuration statistical errors (which are small), the uncertainties native to the MEM analysis (which are comparable), and the widths of the spectral peaks (which are typically large) it remains a matter of judgment to decide which type of uncertainty is physically relevant. We have here taken the point of view that the spectral peak width should be considered the error of the simulation results because it is a measure of the information content in the spirit the Shannon-Jaynes entropy contained in the lattice data. By this measure we accept that the uncertainties generally exceed the gauge configuration statistical errors.

From the peak volumes of the spectral density functions we learn about the efficacy of the operators $`\mathrm{\Phi }_{1,2}`$ of coupling to the low-lying physical excitations of the meson-meson system. Although most of the important excitation mechanisms appear to be covered by $`\mathrm{\Phi }_{1,2}`$ more operators should be employed in more detailed studies of hadronic interaction mechanisms.

Technically, investigations into hadronic interaction are very demanding for both the need to extract small energy differences and having to dealing with (noisy) nonlocal operators involving loops that are a combination of gauge link products and light quark propagators. Anisotropic lattices and advanced analysis techniques seem essential tools in future studies.

###### Acknowledgements.

This material is based upon work supported by the National Science Foundation under Grant No. 0073362. Resources made available through the Lattice Hadron Physics Collaboration (LHPC) were used in this project.