Investigating efficient methods for computing four-quark correlation functions

TL;DR

This paper compares various computational methods for efficiently calculating four-quark correlation functions in lattice QCD, providing guidelines for future studies of multi-quark systems and meson interactions.

Contribution

It systematically evaluates and compares different techniques for computing four-quark correlation functions, offering practical recommendations for lattice QCD calculations.

Findings

Certain method combinations outperform others in efficiency.

The methods are applicable to various multi-quark systems.

Guidelines are provided for future lattice QCD studies.

Abstract

We discuss and compare the efficiency of various methods, combinations of point-to-all propagators, stochastic timeslice-to-all propagators, the one-end trick and sequential propagators, to compute two-point correlation functions of two-quark and four-quark interpolating operators of different structure including quark-antiquark type, mesonic molecule type, diquark-antidiquark type and two-meson type. Although we illustrate our methods in the context of the $a_0(980)$, they can be applied for other multi-quark systems, where similar diagrams appear. Thus our results could provide helpful guidelines on the choice of methods for correlation function computation for future lattice QCD studies of meson-meson scattering and possibly existing tetraquark states.

Figures

Tables

Table 1. Table 1: PACS-CS gauge link ensemble used (cf. also [ 25 ] ).

$\beta$	${\left(L/a\right)}^3\times T/a$	${\kappa }_{u,d}$	${\kappa }_s$	$c_{\text{SW}}$	$a$ [fm]	$m_{\pi }$ [MeV]	# configurations
1.90	${32}^3\times 64$	$0.1370$	$0.1364$	$1.715$	$0.091$	$300$	$500$

Equations

$${\mathcal{O}}^1$$

$${\mathcal{O}}^2$$

$${\mathcal{O}}^3$$

$${\mathcal{O}}^4$$

$${\mathcal{O}}^5$$

$${\mathcal{O}}^6$$

$$\displaystyle C_{jk}(t)=\Big{\langle}\mathcal{O}^{j}(t_{2})\mathcal{O}^{k{\dagger}}(t_{1})\Big{\rangle}\quad,\quad t=t_{2}-t_{1}>0,$$

$$S^{(q)}=\sum _{x,y}{\bar{q}}_{a,A}(x)D_{a,A;b,B}^{(q)}(x;y)q_{b,B}(y),$$

$$\sum _y{D}_{a,A;b,B}^{(q)}(x;y)G_{b,B;c,C}^{(q)}(y;z)={\delta }_{a;c}{\delta }_{A;C}\delta (x;z).$$

$$\displaystyle C_{11}(t)=-\frac{1}{V_{s}}\sum_{{\bf x},{\bf y}}\Big{\langle}\textrm{Tr}\Big{(}G^{(d)}({\bf x},t_{1};{\bf y},t_{2})G^{(u)}({\bf y},t_{2};{\bf x},t_{1})\Big{)}\Big{\rangle}_{U}\ \ =$$

$$\displaystyle\hskip 19.91684pt=-\frac{1}{V_{s}}\sum_{{\bf x},{\bf y}}\Big{\langle}\textrm{Tr}\Big{(}\gamma_{5}\Big{(}G^{(d)}({\bf y},t_{2};{\bf x},t_{1})\Big{)}^{\dagger}\gamma_{5}G^{(u)}({\bf y},t_{2};{\bf x},t_{1})\Big{)}\Big{\rangle}_{U},$$

$$\displaystyle C_{11}(t)=-\bigg{\langle}\sum_{\mathbf{y}}\textrm{Tr}\Big{(}\gamma_{5}\left(G^{(d)}({\bf y},t_{2};{\bf x},t_{1})\right)^{\dagger}\gamma_{5}G^{(u)}({\bf y},t_{2};{\bf x},t_{1})\Big{)}\bigg{\rangle}_{U},$$

$$\sum _y{D}_{a,A;b,B}^{(q)}(x;y){\varphi }_{b,B}^{(q)}(y)[c,C,z]={\xi }_{a,A}(x)[c,C,z],{\xi }_{a,A}(x)[c,C,z]={\delta }_{a;c}{\delta }_{A;C}\delta (x;z),$$

$$G_{b,B;a,A}^{(q)}(y;x)={\varphi }_{b,B}^{(q)}(y)[a,A,x],$$

$$\displaystyle C_{11}(t)=-\bigg{\langle}(\gamma_{5})_{A;B}\bigg{(}\sum_{\bf y}\phi^{(d)}(\mathbf{y},t_{2})[a,B,{\bf x},t_{1}]^{\dagger}\gamma_{5}\phi^{(u)}(\mathbf{y},t_{2})[a,A,{\bf x},t_{1}]\bigg{)}\bigg{\rangle}_{U}.$$

$$\sum _y{D}_{a,A;b,B}^{(q)}(x;y){\varphi }_{b,B}^{(q)}(y)[t_0,n]={\xi }_{a,A}(x)[t_0,n],{\xi }_{a,A}(x)[t_0,n]=\delta (x_0,t_0){\Xi }_{a,A}(\mathbf{x})[n],$$

$$\frac{1}{N}\sum _{n=1}^N{\Xi }_{a,A}^{\ast }(\mathbf{x})[n]{\Xi }_{b,B}(\mathbf{y})[n]={\delta }_{a;b}{\delta }_{A;B}\delta (\mathbf{x};\mathbf{y})+\text{unbiased noise}.$$

$$G^{(q)}(y;\mathbf{x},t_0)=\frac{1}{N}\sum _{n=1}^N{\varphi }^{(q)}(y)[t_0,n]\xi (\mathbf{x},t_0)[t_0,n{]}^{†}+\text{unbiased noise}.$$

$$C_{11}(t)=-\frac{1}{N(N-1)}\sum _{n\ne \tilde{n}}\frac{1}{V_s}$$

$$\displaystyle\hskip 19.91684pt\bigg{\langle}\bigg{(}\sum_{\mathbf{y}}\phi^{(d)}(\mathbf{y},t_{2})[t_{1},\tilde{n}]^{\dagger}\gamma_{5}\phi^{(u)}(\mathbf{y},t_{2})[t_{1},n]\Bigg{)}\bigg{(}\sum_{\mathbf{x}}\xi(\mathbf{x},t_{1})[t_{1},n]^{\dagger}\gamma_{5}\xi(\mathbf{x},t_{1})[t_{1},\tilde{n}]\bigg{)}\bigg{\rangle}_{U}.$$

$$\sum _{\mathbf{y}}{G}^{(q_1)}(x;\mathbf{y},t)\Gamma G^{(q_2)}(\mathbf{y},t;z),$$

$$\sum _y{D}_{a,A;b,B}^{(q_1)}(x;y){\varphi }_{b,B}^{(q_1)}(y)[t_0,n]={\xi }_{a,A}(x)[t_0,n]$$

$$\sum _y{D}_{a,A;b,B}^{(q_2)}(x;y){\tilde{\varphi }}_{b,B}^{(q_2)}(y)[t_0,\Gamma ,n]=({\gamma }_5{\Gamma }^{†}\xi {)}_{a,A}(x)[t_0,n],$$

$$\sum _{\mathbf{y}}{G}^{(q_1)}(x;\mathbf{y},t)\Gamma G^{(q_2)}(\mathbf{y},t;z)=\frac{1}{N}\sum _{n=1}^N{\varphi }^{(q_1)}(x)[t,n]{\tilde{\varphi }}^{(q_2)}(z)[t,\Gamma ,n{]}^{†}{\gamma }_5+\text{unbiased noise}.$$

$$\displaystyle C_{11}(t)=-\frac{1}{N}\sum_{n=1}^{N}\frac{1}{V_{s}}\bigg{\langle}\sum_{\mathbf{y}}\tilde{\phi}^{(d)}(\mathbf{y},t_{2})[t_{1},\mathbbm{1},n]^{\dagger}\gamma_{5}\phi^{(u)}(\mathbf{y},t_{2})[t_{1},n]\bigg{\rangle}_{U}.$$

$$\sum _{\mathbf{y}}{G}^{(q_1)}(x;\mathbf{y},t)\Gamma G^{(q_2)}(\mathbf{y},t;z).$$

$$\sum _y{D}_{a,A;b,B}^{(q_1)}(\mathbf{x},x_0;y){\psi }_{b,B}^{(q_1;q_2)}(y)[t,\Gamma ;c,C,z]=(\Gamma \varphi {)}_{a,A}^{(q_2)}(\mathbf{x},x_0)[c,C,z]\delta (x_0;t)$$

$$\displaystyle\bigg{(}\sum_{\mathbf{y}}G^{(q_{1})}(x;\mathbf{y},t)\Gamma G^{(q_{2})}(\mathbf{y},t;z)\bigg{)}_{A,B}^{a,b}=\psi_{a,A}^{(q_{1};q_{2})}(x)[t,\Gamma;b,B,z].$$

$$R^{(a),(b)}(t)=\frac{\Delta C^{(a)}(t)\cdot \sqrt{{\tau }^{(a)}}}{\Delta C^{(b)}(t)\cdot \sqrt{{\tau }^{(b)}}},$$

$$\displaystyle\Big{(}\Delta C_{22}^{4\times\textrm{con}}(t)\Big{)}^{2}\propto\Big{\langle}\mathcal{O}[s^{(1)},s^{(2)},s^{(3)},s^{(4)}](t_{2})\tilde{\mathcal{O}}[s^{(1)},s^{(2)},s^{(3)},s^{(4)}]^{\dagger}(t_{1})\Big{\rangle}$$

$$\displaystyle\Big{(}\Delta C_{22}^{2\times\textrm{con}}(t)\Big{)}^{2}\propto\Big{\langle}\mathcal{O}[s^{(1)},s^{(1)},s^{(2)},s^{(2)}](t_{2})\tilde{\mathcal{O}}[s^{(3)},s^{(3)},s^{(4)},s^{(4)}]^{\dagger}(t_{1})\Big{\rangle}$$

$$\displaystyle\mathcal{O}[s_{1},s_{2},s_{3},s_{4}]=\bigg{(}\sum_{\bf{x}}\Big{(}{\bar{s}}_{1}({\bf x})\gamma_{5}u({\bf x})\Big{)}\Big{(}{\bar{d}}({\bf x})\gamma_{5}s_{2}({\bf x})\Big{)}\bigg{)}\bigg{(}\sum_{\bf{y}}\Big{(}{\bar{s}}_{3}({\bf y})\gamma_{5}u({\bf y})\Big{)}\Big{(}{\bar{d}}({\bf y})\gamma_{5}s_{4}({\bf y})\Big{)}\bigg{)}$$

$$\displaystyle\tilde{\mathcal{O}}[s_{1},s_{2},s_{3},s_{4}]=\sum_{\bf{x}}\Big{(}{\bar{s}}_{1}({\bf x})\gamma_{5}u({\bf x})\Big{)}\Big{(}{\bar{d}}({\bf x})\gamma_{5}s_{2}({\bf x})\Big{)}\Big{(}{\bar{s}}_{3}({\bf x})\gamma_{5}u({\bf x})\Big{)}\Big{(}{\bar{d}}({\bf x})\gamma_{5}s_{4}({\bf x})\Big{)}.$$

$$\frac{\Delta C_{22}^{4\times \text{connected}}(t)}{\Delta C_{22}^{2\times \text{connected}}(t)}\propto e^{-(2m_K-m_{\pi })t}.$$

$$\displaystyle S_{\textrm{quark}}[q,\bar{q},U]=\sum_{q=u,d,s}\bigg{(}\sum_{x}\bar{q}_{x}q_{x}-\kappa_{q}c_{\textrm{SW}}\sum_{x,\mu,\nu}\frac{i}{2}\bar{q}_{x}\sigma_{\mu\nu}F_{x,\mu\nu}q_{x}$$

Full text

Investigating efficient methods for computing four-quark correlation functions

**Abdou Abdel-Rehim1,2, Constantia Alexandrou1,3, Joshua Berlin4,

Mattia Dalla Brida5, Jacob Finkenrath1, Marc Wagner4**

1Computation-based Science and Technology Research Center, The Cyprus Institute,

20 Kavafi Street, 2121 Nicosia, Cyprus

2Department of Engineering, Science and Mathematics SUNY Polytechnic Institute, Utica, New York 13502, USA

3Department of Physics, University of Cyprus, P.O. Box 20537, 1678 Nicosia, Cyprus

4Goethe-Universität Frankfurt am Main, Institut für Theoretische Physik, Max-von-Laue-Straße 1, D-60438 Frankfurt am Main, Germany

5 Dipartimento di Fisica, Università di Milano-Bicocca & INFN, sezione di Milano-Bicocca, Piazza della Scienza 3, I-20126 Milano, Italy

January 25, 2017

Abstract

We discuss and compare the efficiency of various methods, combinations of point-to-all propagators, stochastic timeslice-to-all propagators, the one-end trick and sequential propagators, to compute two-point correlation functions of two-quark and four-quark interpolating operators of different structure including quark-antiquark type, mesonic molecule type, diquark-antidiquark type and two-meson type. Although we illustrate our methods in the context of the $a_{0}(980)$, they can be applied for other multi-quark systems, where similar diagrams appear. Thus our results could provide helpful guidelines on the choice of methods for correlation function computation for future lattice QCD studies of meson-meson scattering and possibly existing tetraquark states.

1 Introduction

In recent years the study of four-quark systems has revived a lot of interest within the field of lattice QCD. These systems allow to study meson-meson scattering, as well as to investigate the existence and structure of possibly existing four-quark states, e.g. of mesonic molecule type or of diquark-antidiquark type. In particular the latter is motivated by experimental results providing strong indications for the existence of tetraquarks (e.g. the recently observed charged $Z_{c}$ and $Z_{b}$ states [1, 2, 3] or the mass ordering of the nonet of light scalar mesons, which is inverted compared to the expectation from a standard quark-antiquark picture [4]).

In order to perform this study within lattice QCD we need to compute two-point correlation functions involving both two-quark and four-quark interpolating operators. Depending on the details of these operators diagrams involving up to four quark propagators with a non-trivial spacetime structure are present. In contrast to simple quark-antiquark correlation functions, which can be computed in a straightforward way, e.g. by using standard point-to-all propagators, these four-quark correlation functions require in many cases all-to-all propagators (propagators from any point in space on a timeslice to any other point in space on another timeslice, what makes them considerably expensive) and, hence, more advanced techniques.

One possibility to compute them is the distillation method [5], which provides all-to-all propagators between specifically smeared quark field operators (Laplacian Heaviside smearing). Distillation has recently been applied to such four-quark correlation functions (cf. e.g. [6, 7, 8]). The distillation method however comes with an expensive overhead and only pays if a much larger number of interpolating fields is used. Moreover, for large volume lattices it is impractical and one has to apply its stochastic variant, which introduces additional stochastic noise terms for each quark propagator.

In this work we explore a different strategy, namely combining several traditional techniques to compute quark propagators and four-quark correlation functions: (A) point-to-all propagators (cf. e.g. [9, 10]); (B) stochastic timeslice-to-all propagators (cf. e.g. [11, 12]); (C)  the one-end trick (cf. e.g. [13, 14]); (D) sequential propagators (cf. e.g. [15]). There are many different types of diagrams and for each type we discuss several methods (combinations of the above mentioned techniques (A) to (D)) and determine numerically the most efficient method. We study the $a_{0}(980)$ channel (quantum numbers $I(J^{P})=1(0^{+})$) at a lattice spacing of $a\approx 0.09$fm and spacetime volume of $(32a)^{3}\times 64a$ with a variety of interpolating operators including quark-antiquark type, mesonic molecule type, diquark-antidiquark type and two-meson type (the latter describes two independent mesons with total zero momentum). In other words our work is intended to provide guidelines for future lattice QCD work concerned with arbitrary four-quark correlation functions, in particular guidelines for a quick and hence time-saving decision, which diagrams of a given correlation matrix to compute with which combinations of the above listed techniques.

Parts of this work have already been presented at recent conferences [16, 17, 18, 19, 20, 21].

The paper is structured as follows. In section 2 we introduce the six interpolating operators we investigate and the resulting $6\times 6$ correlation matrix. In section 3 we recapitulate the above mentioned four techniques for propagator and correlation function computation. Section 4 is the main section of this work, where various methods are discussed for each element and diagram of the correlation matrix, with a numerical comparison of their efficiency. In section 5 we conclude, in particular we summarize a few general rules regarding the choice of an efficient method, which seem to hold for most correlation matrix elements investigated. Section A outlines briefly the used lattice setup and gives detailed examples for the calculation of certain correlation functions.

2 Interpolating operators and correlation matrix

To study the $a_{0}(980)$ meson, we consider several interpolating operators $\mathcal{O}^{j}$, which create states with quantum numbers $I(J^{P})=1(0^{+})$, when applied to the vaccuum:

[TABLE]

($V_{s}$ is the spatial volume, $C$ is the charge conjugation matrix). These interpolating operators are of different structure. $\mathcal{O}^{q\bar{q}}$ generates a quark-antiquark pair, while the other operators (2) to (6) generate two quarks and two antiquarks. \mathcal{O}^{\text{K\bar{K}, point}} and \mathcal{O}^{\text{\eta_{s}\pi, point}} are of mesonic molecule type, i.e. resemble a $K\bar{K}$ pair or $\eta_{s}\pi$ pair111${\bar{s}}_{\bf x}\gamma_{5}{s}_{\bf x}$ excites a meson-like structure composed of an $s\bar{s}$ pair, which is expected to have significant overlap to both $\eta$ and $\eta^{\prime}$ and, hence, is denoted by $\eta_{s}$. centered around the same spatial point $\mathbf{x}$. \mathcal{O}^{\text{Q\bar{Q}}} corresponds to a diquark-antidiquark pair222We consider only the lightest (anti-)diquarks, which have spin structure $C\gamma_{5}$ [22, 23, 24].. These three operators are candidates to model the structure of a possibly existing bound four-quark state, i.e. of a tetraquark. The remaining two operators \mathcal{O}^{\text{K\bar{K}, 2part}} and \mathcal{O}^{\text{\eta_{s}\pi, 2part}} also generate meson pairs ($K\bar{K}$ and $\eta_{s}\pi$), this time, however, at independent spatial points $\mathbf{x}$ and $\mathbf{y}$. They should be suited to resolve low-lying two-meson states within the $I(J^{P})=1(0^{+})$ sector.

The interpolating operators (1) to (6) enter a $6\times 6$ correlation matrix,

[TABLE]

where $\langle\ldots\rangle$ denotes the lattice QCD path integral expectation value. Computing such expectations values at affordable numerical costs and with small statistical errors is a highly non-trivial task, in particular, when four-quark operators like (2) to (6) are involved. One possible strategy, which we follow throughout this work, is to suitably combine standard techniques for propagator and correlator computation (point-to-all propagators, stochastic propagators, the one-end trick, sequential propagators). One of the main goals of this work is to discuss, which combinations of techniques are possible in principle for each of the matrix elements $C_{jk}$ and, in a second step, to determine numerically which combination is most efficient, i.e. results in the smallest statistical error at comparable computational cost.

Even though we focus on the $a_{0}(980)$ meson with quantum numbers $I(J^{P})=1(0^{+})$, our findings are of general interest regarding the study of tetraquark systems and unstable mesonic resonances. The study of such systems usually requires the computation of correlation matrices of identical or similar structure as (7) with interpolating operators (1) to (6). For example, after replacing the quark flavors according to $(u,d,s)\rightarrow(c,s,u/d)$, one would obtain a matrix suited to study the $J^{P}=0^{+}$ $D_{s0}^{\ast}(2317)$ meson, and probe whether two- or for-quark structures are dominating. Similarly, $(u,d,s)\rightarrow(u,d,c)$ would allow to study certain charmonium states, for which are tetraquark structure is frequently discussed, e.g. the $Z(4430)$ meson.

To ease notation, we will often picture correlation matrix elements in a diagrammatic way, where quark propagators are represented by arrows. These diagrams do not exhibit information about color and spin, but clearly display the spacetime structure. In particular one can read off, which combinations of methods is suited to compute a diagram. In our case, each correlation matrix element (7) corresponds to either a single diagram or a sum of two diagrams. Due to the flavor structure of the interpolating operators (either $u\bar{d}$ or $u\bar{d}s\bar{s}$) there are within each diagram either two or four quark propagators connecting the timeslices at $t_{1}$ and $t_{2}$. The diagrams are, hence, denoted as “$2\times$ connected” and “$4\times$ connected”, respectively (for an example cf. Figure 1). The diagrammatic representation of the full $6\times 6$ correlation matrix is shown in Figure 2. Note that certain correlation matrix elements have identical diagrammatic representations, since they only differ in their color and spin structure.

3 Techniques for propagator computation

The lattice action for quarks is bilinear in the quark fields $q$ and $\bar{q}$,

[TABLE]

where $D^{(q)}$ denotes the Dirac operator for quark flavor $q$ and indices $a,b,\ldots$, $A,B,\ldots$ and $x,y,\ldots$ label color, spin and spacetime, respectively. The propagator $G^{(q)}$ is the inverse of the Dirac operator, i.e. the solution of the linear system

[TABLE]

In the following subsections we discuss several standard techniques for propagator computation which are well-know in the literature, namely: (1) point-to-all propagators, (2) stochastic propagators, (3) the one-end trick and (4) sequential propagators. We illustrate some of these techniques in the context of a simple example, the correlation function of the interpolating operator \mathcal{O}^{1}=\mathcal{O}^{\text{q\bar{q}}}=(1/\sqrt{V_{s}})\sum_{\bf x}\bar{d}({\bf x})u({\bf x}) (eq. (1)), i.e. correlation matrix element $C_{11}$ in (7) and Figure 2. Integrating over the Grassmann valued quark fields allows us to express the correlation function in terms of quark propagators,

[TABLE]

where $\gamma_{5}$ hermiticity has been used, $\textrm{Tr}(\ldots)$ denotes the trace in spin and color space and $\langle\ldots\rangle_{U}$ is the average over gauge link configurations distributed proportionally to $e^{-S_{\textrm{eff}}}=e^{-(S_{\textrm{gauge}}-\ln(\det(Q)))}$.

3.1 Point-to-all propagators

The exact computation of a propagator $G^{(q)}(x;y)$ from any point in spacetime $x$ to any other point $y$ is numerically not feasible, because for typical lattices with e.g. $32^{4}$ lattice sites both $D^{(q)}$ and $G^{(q)}$ are matrices with $\mathcal{O}(10^{7})$ entries. However, using translation invariance it is often possible to simply compute propagators from a single spacetime point $x$ to any other point $y$. For example in (10), $\sum_{\textbf{x}}$ can be replaced by the spatial volume $V_{s}$,

[TABLE]

where $\mathbf{x}$ denotes an arbitrary but fixed point in space. The spacetime column index of both propagators is now $(\mathbf{x},t_{1})$, i.e. it is not anymore necessary to compute the full matrix $G^{(u/d)}$, but only 12 columns of this matrix ($\textrm{3 color}\times\textrm{4 spin}$), which is numerically feasible.

Explicitly, one has to solve 12 linear systems,

[TABLE]

where $c=1,2,3$ and $C=1,2,3,4$ label 12 different point sources. (Each solution $\phi$ corresponds to a single column of the inverse of the Dirac matrix $D$; $\phi$ are, therefore, also commonly called “inversions”). The propagator ending at spacetime point $x$, a so-called point-to-all propagator (cf. e.g. [9, 10]), is then

[TABLE]

The example correlation function (10) expressed in terms of such point-to-all propagators is

[TABLE]

Of course, for each diagram, translation invariance allows one to replace only a single spatial sum $\sum_{\mathbf{x}}$ by an arbitrary fixed $\mathbf{x}$. Diagrams, where all propagators either start or end at the same spacetime point can, hence, be expressed exclusively in terms of point-to-all propagators (e.g. $C_{11}$, $C_{12}$, left diagram of $C_{22}$). However, since this is not the case for the majority of diagrams, additional methods to compute propagators are necessary.

After replacing a spatial sum $\sum_{\mathbf{x}}$ by a fixed $\mathbf{x}$, spatial averaging to reduce statistical errors is not implemented anymore. For each set of 12 inversions (12) and each diagram only a single sample is computed. Of course, it is possible to compute additional samples by choosing different x values. This however requires additional sets of 12 inversions for each different $\mathbf{x}$, and quickly becomes expensive.

3.2 Stochastic timeslice-to-all propagators

While it is not possible in practice to compute the propagator from any spacetime point $x$ to any other spacetime point $y$ exactly, one can at least estimate it stochastically. Quite common are so-called stochastic timeslice-to-all propagators (cf. e.g. [11, 12]): stochastically estimated propagators from any space point $\mathbf{x}$ in a given time-slice $t_{0}$, to any other spacetime point $y$.

Again linear systems have to be solved labeled now by $n=1,\ldots,N$,

[TABLE]

where $\Xi_{a,A}(\mathbf{x})[n]$ are random numbers satisfying

[TABLE]

A convenient choice is $\Xi_{a,A}(\mathbf{x})[n]\in\mathbb{Z}(2)\times\mathbb{Z}(2)$ which results in an unbiased noise proportional to {\rm O}\Big{(}\frac{1}{\sqrt{N}}\Big{)}. As usual since the noise average and the average over the gauge field commute, in practice one can take a fairly small number $N$ of noise sources per gauge configuration, but not smaller than the number of propagators in the diagram.

Using (15) and (16) it is straightforward to show

[TABLE]

The example correlation function (10) expressed in terms of stochastic timeslice-to-all propagators is

[TABLE]

Note that each propagator needs to be estimated by a different pair of stochastic sources $\xi[n]$ and corresponding inversions $\phi[n]$ (guaranteed here by $\sum_{n\neq\tilde{n}}$).

Stochastic timeslice-to-all propagators are most flexible, i.e. replacing a spatial sum $\sum_{\mathbf{x}}$ by a fixed $\mathbf{x}$ as in the case of point-to-all propagators is not necessary. In principle all diagrams of the correlation matrix (7) can be computed using exclusively stochastic timeslice-to-all propagators. A severe drawback of these propagators is, however, that they introduce additional stochastic noise. The number of stochastic noise terms is $\approx V_{s}^{M}\times(\textrm{number of signal terms})$, where $M$ is the number of stochastic timeslice-to-all propagators in a diagram333For example the number of signal terms in (18) is $\propto V_{s}^{2}$ (all $(\mathbf{x},t_{1})$ connected by a pair of propagators with all $(\mathbf{y},t_{2})$), while the number of stochastic noise terms is $\propto V_{s}^{4}$.. While using a single stochastic propagator, i.e. $M=1$, typically leads to acceptable signal-to-noise ratios, the noise grows quite rapidly with the number of stochastic propagators. Already for $M>2$, the signal can easily be lost in the noise if these techniques are applied naively. Therefore, a promising strategy might be to combine a single stochastic timeslice-to-all propagator with several point-to-all propagators, as we shall see later on.

3.3 The one-end trick

The one-end trick is an efficient technique to estimate the product of two propagators stochastically (cf. e.g. [13, 14]). The product must be of the form

[TABLE]

i.e. the propagators are connected at spacetime point $(\mathbf{y},t)$ with a sum over $\mathbf{y}$, but no further propagators starting or ending at $(\mathbf{y},t)$. The one-end trick is thus particularly suited to compute correlation matrix elements, where at least one of the two interpolating operators is either a $q\bar{q}$ or a two-meson operator, i.e. \mathcal{O}^{\text{q\bar{q}}}, \mathcal{O}^{\text{K\bar{K}, 2part}} or \mathcal{O}^{\text{\eta_{s}\pi, 2part}} defined in section 2.

One has to solve $2N$ linear systems labeled by $n=1,\ldots,N$,

[TABLE]

where $\xi$ is a stochastic timeslice source defined in (15). With the resulting $\phi$ and $\tilde{\phi}$ the product of propagators can be estimated as,

[TABLE]

Applying the one-end trick to the example correlation function (10) results in

[TABLE]

Note that there is a big advantage in efficiency when using the one-end trick, compared to using two ordinary stochastic timeslice-to-all propagators discussed in the previous subsection. The one-end trick introduces $\approx V_{s}\times(\textrm{number of signal terms})$ stochastic noise terms, i.e. for the example correlation function (23) $\propto V_{s}^{3}$ stochastic noise terms. In contrast to that using two stochastic timeslice-to-all propagators (18) will generate $\propto V_{s}^{4}$ stochastic noise terms.

3.4 Sequential propagators

Quite similar to the one-end trick the technique of sequential propagators (cf. e.g. [15]) is applicable when two propagators are connected at spacetime point $(\mathbf{y},t)$ with a sum over $\mathbf{y}$, but no further propagators start or end at $(\mathbf{y},t)$, i.e. as in (19),

[TABLE]

In practice, sequential propagators need to be combined with the methods discussed in the previous subsections, point-to-all propagators, stochastic timeslice-to-all propagators and/or the one-end trick. In the following we explain the basic idea of using sequential propagators in combination with point-to-all propagators. Corresponding equations for stochastic timeslice-to-all propagators or the one-end trick are straightforward to derive and, therefore, not presented.

In a first step a point-to-all propagator for flavor $q_{2}$ from a fixed spacetime point $z$ to any other spacetime point $y$ is computed by solving 12 linear systems as explained in section 3.1. This propagator $\phi_{a,A}^{(q_{2})}(x)[b,B,y]=G_{a,A;b,B}^{(q_{2})}(x;y)$ is then used as the right hand side of further sets of 12 linear systems,

[TABLE]

(labeled by indices $c=1,2,3$, $C=1,2,3,4$ and $t$), which have to be solved with respect to $\psi$. Then

[TABLE]

The use of sequential propagators is quite efficient when $t=z_{0}$, i.e. when $q_{2}$ propagates within timeslice $t=z_{0}$. Then the total number of inversions is limited to 24 (12 to obtain $\phi$, then 12 more to obtain the final result $\psi$). The result is then exact. If one is, however, interested in several $t\neq z_{0}$, one has to solve 12 linear systems for each value of $t$ to obtain the corresponding $\psi$, which can be a rather computer time consuming task.

Even though sequential propagators have a similar application as the one-end trick (cf. (19) and (24)) with the additional practical limitation to $t=z_{0}$, it is essential to use this technique, when computing the correlation matrix (7). In contrast to the one-end trick, which is a self contained technique, sequential propagators always need to be combined with other techniques and, hence, offer a lot of flexibility. For example the triangular diagram appearing e.g. in $C_{15}$ can be computed in a very efficient way, when combining sequential propagators and the one-end trick (for details cf. section 4.3.1).

3.5 Diagrammatic representation of propagators and correlation functions

In the following we introduce a diagrammatic representation of propagators and correlation functions, which will be particularly useful in section 4. There we discuss how to compute the elements of the correlation matrix 7 using the previously introduced numerical techniques.

Propagators are represented by arrows, which point from the sources (denoted by $\xi$ in sections 3.1 to 3.4) to the corresponding solutions of the linear systems (denoted by $\phi$). Point sources (12) (which are restricted to a specific spacetime point) are represented by black boxes, while stochastic timeslice sources (15) and the solutions of the linear systems (which have non-zero values for any space or spacetime point, respectively) are represented by black circles. Point-to-all propagators are colored in blue, stochastic timeslice-to-all propagators in red and a one-end trick combination of two propagators in green. Since sequential propagators are always used in combination with one of the three techniques mentioned in the previous subsection, they are colored accordingly in blue, red or green. Two propagators expressed in terms of a sequential propagator can easily be identified by a big circle at their junction.

After correlation functions are expressed in terms of sources $\xi$ and solutions $\phi$ using the techniques discussed in sections 3.1 to 3.4, there are typically one or more spatial sums. These sums are indicated in the diagrams by $\sum\longleftrightarrow$ in black. When such sums are implicitly realized by either the one-end trick or by sequential propagators, they are indicated by $\sum\longleftrightarrow$ in gray.

To illustrate this diagrammatic representation, we show a couple of examples in Figure 3: the correlation function $C_{11}$ expressed in terms of point-to-all propagators (14), stochastic timeslice-to-all propagators (18) and the one-end trick (23) as well as the sequential propagator (26).

4 Computation of the correlation matrix

Each diagram of the correlation matrix (7) (cf. also Figure 2) can be computed in several ways (denoted as methods in the following) by combining the techniques for propagator computation discussed in section 3. An efficient method satisfies the following three criteria,

• (A)

requires only a small number of inversions,

• (B)

averages the diagram over space (to reduce statistical fluctuations originating from the gauge links),

• (C)

introduces no or only a moderate number of stochastic noise terms (due to stochastic timeslice-to-all propagators or the one-end trick).

In practice ideal methods do not exist. In particular (B) and (C) exclude each other: to avoid additional stochastic noise terms, one has to use exclusively point-to-all-propagators, which do not average the diagram over space. Therefore, the best compromise has to be found for each diagram. We do this in this section by comparing the efficiency of different methods numerically.

To compare two methods $(a)$ and $(b)$ quantitatively, we define the quality ratio

[TABLE]

where $\Delta C^{(x)}(t)$ denotes the statistical error of the diagram at temporal separation $t$ obtained with method $(x)$ and $\tau^{(x)}$ the corresponding computing time. $R^{(a),(b)}(t)<1$ indicates that method $(b)$ is inferior to method $(a)$, i.e. that the statistical error obtained with method $(a)$ at comparable computing time is smaller, while $R^{(a),(b)}(t)>1$ indicates the opposite. In the following we label the investigated combinations of techniques for each diagram by $(a)$, $(b)$, $(c)$, … according to decreasing quality, i.e. method $(a)$ is always the best. We compare the remaining methods $(b)$, $(c)$, … exclusively to method $(a)$, hence, it is convenient to use the notation $R^{(b)}(t)=R^{(a),(b)}(t)$.

The numerical comparisons presented in the following are based on gauge link configurations generated by the PACS-CS collaboration [25]. Details regarding these configurations are collected in appendix A. When using point-to-all propagators, we did a single set of 12 inversions per gauge link configuration corresponding to a set of 12 point sources at the same randomly chosen point in space. For stochastic timeslice-to-all propagators, we used 15 independently chosen stochastic sources per gauge link configuration, i.e. $N=15$ in section 3.2. When using the one-end trick, we inverted 3 independently chosen stochastic sources per gauge link configuration, i.e. $N=3$ in section 3.3. The resulting quality ratios $R^{(x)}$ are only weakly dependent on $N$, since larger $N$ lead to smaller statistical errors $\Delta C^{(x)}$ (suppressed by approximately $1/\sqrt{N}$)), but also to larger computing times $\tau^{(x)}$ (increased by approximately $N$), which has been taken into account in the definition (27).

For some techniques the number of inversions and, hence, the computing time is proportional to the number of temporal separations $t$, for which the correlation matrix element is computed. In such cases we limited our computations to the range $0\leq t/a\leq 14$, which are typical temporal separations where we have signal in our correlators. Since these larger numbers of inversions are always done to compute the strange quark propagator, and are thus relatively cheap, the quality ratios $R^{(x)}$ are only weakly dependent on variations of the computed temporal range.

We do not compare all possible methods, but focus on promising methods, i.e. methods, which fulfill (A) to (C), the criteria listed above, at least to some extent. In particular, we do not consider methods making excessive use of stochastic techniques and, hence, introduce a large number of stochastic noise terms. The maximum we allow is either a stochastic timeslice-to-all propagator in combination with the one-end trick or twice the one-end trick444The number of stochastic noise terms is $\approx V_{s}^{M}\times(\textrm{number of signal terms})$, where $M$ counts the number of stochastic timeslice-to-all propagators and applications of the one-end trick (cf. sections 3.2 and 3.3). We set the maximum to $M=2$.. Moreover, whenever possible, we combine the techniques for propagator computation in such a way that additional inversions for each considered temporal separation $t$ are not necessary (cf. the previous paragraph). In particular, to avoid a very large number of inversions, we exclude methods where additional inversions for each considered temporal separation $t$ are done for point-to-all propagators (point-to-all propagators require 12 inversions per separation).

The following subsections focus on the results of the comparisons of different methods. Technical details regarding the implementation of some of the diagrams can be found in appendix B.

4.1 Two-quark – two-quark correlation function

4.1.1 $C_{11}$ ($1\equiv q\bar{q}$)

Promising methods for numerical computation are the following (cf. also section 3, where this diagram has extensively been discussed, to illustrate various techniques for propagator computation).

• $(a)$

one-end trick:

2 $u/d$ inversions

• $(b)$

2 point-to-all propagators:

12 $u/d$ inversions

The one-end trick averages the diagram over space (reflected by $\sum\longleftrightarrow$ both at timeslice $t_{1}$ and at timeslice $t_{2}$ in $(a)$), but at the same time introduces stochastic noise terms. A numerical comparison using the lattice setup described in appendix A results in $R^{(b)}(t)\approx 0.3\ldots 0.4$ and $\bar{R}^{(b)}=0.35$, where $\bar{R}^{(b)}$ is the average of $R^{(b)}(t)$ in the range $5\leq t/a\leq 14$, which are typical temporal separations to extract energy levels from our correlation functions (cf. Figure 4). Consequently, the one-end trick leads to statistical errors smaller by a factor around $2$ compared to point-to-all propagators, when investing similar amounts of computing time. Hence, this is our preferred method to compute $C_{11}$. It should, however, be noted that the efficiency might depend to some extent on the lattice setup, in particular could be somewhat different, when another $u/d$ quark mass or spatial volume is used555For example in [26] it has been found that the efficiency of the two methods is quite similar, when the $u/d$ quark mass is replaced by the much heavier charm quark mass. Similarly, counting the number of signal terms and the number of stochastic noise terms in (23) ($\propto V_{s}^{2}$ and $\propto V_{s}^{3}$, respectively) suggests that the one-end trick becomes even more efficient for larger spatial volumes $V_{s}$..

4.2 Two-quark – four-quark correlation functions

4.2.1 $C_{12}$, $C_{13}$, $C_{14}$ ($1\equiv q\bar{q}$, $2\equiv K\bar{K}\textrm{, point}$, $3\equiv\eta_{s}\pi\textrm{, point}$, $4\equiv Q\bar{Q}$)

Promising methods for numerical computation are the following (cf. also appendix B.1, where certain technical aspects are discussed).

• $(a)$

3 point-to-all propagators:

12 $u/d$ inversions, 12 $s$ inversions

• $(b)$

one-end trick:

2 $u/d$ inversions

stochastic timeslice-to-all propagator:

$\#t$ $s$ inversions ($\#t$ denotes the number of temporal separations computed)

• $(c)$

2 point-to-all propagators:

12 $u/d$ inversions

stochastic timeslice-to-all propagator:

$\#t$ $s$ inversions

Method $(b)$ averages the diagram over space (reflected by $\sum\longleftrightarrow$ both at timeslice $t_{1}$ and at timeslice $t_{2}$), but at the same time introduces a rather large number of stochastic noise terms. Methods $(a)$ and $(c)$ on the other hand compute only a single sample. For all three diagrams $C_{12}$, $C_{13}$ and $C_{14}$ method $(a)$, which does not introduce any stochastic noise terms, is clearly more efficient, as shown in Figure 5), with corresponding quality ratios $\bar{R}^{(b)},\bar{R}^{(c)}\approx 0.35\ldots 0.63$.

4.3 Two-quark – two-meson correlation functions

4.3.1 $C_{15}$ ($1\equiv q\bar{q}$, $5\equiv K\bar{K}\textrm{, 2part}$)

Since the $s$ quark propagates within timeslice $t_{1}$, $C_{15}$ can be computed efficiently using a sequential propagator (the $s$ quark propagator together with either the $u$ or the $d$ quark propagator). Promising methods are the combination of a sequential propagator with the one-end trick or with point-to-all propagators as shown in the figure below (cf. also appendix B.2, where certain technical aspects are discussed).

• $(a)$

one-end trick (with 1 sequential propagator):

2 $u/d$ inversions, 1 $s$ inversion

• $(b)$

2 point-to-all propagators (1 is a sequential propagator):

24 $u/d$ inversions, 12 $s$ inversions

As for the matrix element $C_{11}$ discussed in section 4.1.1, the one-end trick is more efficient than point-to-all propagators (cf. Figure 6).

4.3.2 $C_{16}$ ($1\equiv q\bar{q}$, $6\equiv\eta_{s}\pi\textrm{, 2part}$)

$C_{16}$ is a product of two disconnected parts. The quark loop can either be computed using a point-to-all or a stochastic timeslice-to-all propagator. Promising methods to compute this are the following.

• $(a)$

one-end trick:

1 $u/d$ inversion

stochastic timeslice-to-all propagator:

1 $s$ inversion

• $(b)$

one-end trick:

1 $u/d$ inversion

point-to-all propagator:

12 $s$ inversions

• $(c)$

2 point-to-all-propagators:

12 $u/d$ inversions

stochastic timeslice-to-all propagator:

1 $s$ inversions

Method $(a)$ introduces a rather large number of stochastic noise terms. However, in contrast to method $(b)$ and method $(c)$, which also introduce stochastic noise terms, method $(a)$ averages the diagram over space and results in being more efficient (cf. Figure 7).

4.4 Four-quark – four-quark correlation functions

4.4.1 $4\times$ connected $C_{22}$, $C_{23}$, $C_{24}$, $C_{33}$, $C_{34}$, $C_{44}$ ($2\equiv K\bar{K}\textrm{, point}$, $3\equiv\eta_{s}\pi\textrm{, point}$, $4\equiv Q\bar{Q}$)

The correlation matrix elements $C_{22}$, $C_{23}$, $C_{24}$, $C_{33}$, $C_{34}$ and $C_{44}$ have different spin and color structures, but are identical with respect to spacetime. Therefore, it is appropriate to discuss them together. Since both at timeslice $t_{1}$ and at timeslice $t_{2}$ four quarks are located at the same point in space, there is only a single possibility to compute this type of diagram efficiently.

• $(a)$

4 point-to-all-propagators:

12 $u/d$ inversions, 12 $s$ inversions

4.4.2 $2\times$ connected $C_{22}$, $C_{23}$, $C_{24}$, $C_{33}$, $C_{34}$, $C_{44}$ ($2\equiv K\bar{K}\textrm{, point}$, $3\equiv\eta_{s}\pi\textrm{, point}$, $4\equiv Q\bar{Q}$)

Again there is only a single possibility to compute this type of diagram efficiently.

• $(a)$

3 point-to-all-propagators:

12 $u/d$ inversions, 12 $s$ inversions

stochastic timeslice-to-all propagator:

$\#t$ $s$ inversions

4.4.3 Comparison of $4\times$ connected and $2\times$ connected diagrams

The statistical error of each of the six correlation matrix elements $C_{22}$, $C_{23}$, $C_{24}$, $C_{33}$, $C_{34}$ and $C_{44}$ is a combination of the statistical errors of the corresponding $4\times$ connected and $2\times$ connected diagrams and will be dominated by the larger of the two errors. Therefore, it is interesting to compare the statistical errors of the $4\times$ connected and the $2\times$ connected diagrams. For that purpose we use again the quality ratio defined in (27), $R^{2\times\textrm{con}}(t)\equiv R^{4\times\textrm{con},2\times\textrm{con}}(t)$.

In the top row of Figure 8 these quality ratios are shown for all six correlation matrix elements $C_{22}$, $C_{23}$, $C_{24}$, $C_{33}$, $C_{34}$ and $C_{44}$. They rapidly decrease with the temporal separation $t$ and are mostly below $0.1$ or even significantly smaller for $5\leq t\leq 14$, where effective mass plateaus are typically read off. Consequently, the statistical errors of the correlation matrix elements will be dominated by the $2\times$ connected diagrams. Thus, it is much more important to determine the optimal method of computation for the $2\times$ connected than for the $4\times$ connected diagrams.

To investigate this drastic difference in the statistical errors in more detail, we plot in the bottom row of Figure 8 a related quantity, $\log(R^{2\times\textrm{con}}(t)/R^{2\times\textrm{con}}(t+a))$. Assuming an exponential behavior $R^{2\times\textrm{con}}(t)\propto e^{-\alpha t}$, the plotted quantity $\log(\ldots)$ corresponds to the exponent $\alpha$. For all six cases $\log(\ldots)$ fluctuates around the same constant, $\alpha\approx 0.4$, which indeed shows that the statistical errors of the $2\times$ connected diagrams increase exponentially in comparison to the $4\times$ connected diagrams, i.e. proportionally to $e^{\alpha t}$.

This behavior can be understood in the following way. The squared statistical error of a diagram is proportional to the squared diagram and, hence, can be expressed as the correlation function of two eight-quark operators with appropriately chosen spin, color and spacetime structure. As an example consider,

[TABLE]

with

[TABLE]

Note that four different degenerate strange quarks $s^{(1)},\ldots,s^{(4)}$ have been introduced in such a way that the correlation functions (28) and (29) reproduce exactly the squared diagrams of interest shown in Figure 9. Since both $\mathcal{O}$ and $\tilde{\mathcal{O}}$ with flavor structure $[s^{(1)},s^{(2)},s^{(3)},s^{(4)}]$ generate quantum numbers $I(J^{P})=2(0^{+})$ and strangeness $S^{(1,3)}=+1$ and $S^{(2,4)}=-1$ for the four strange flavors, the correlation function (28) will decay asymptotically according to $e^{-4m_{K}t}$. Similarly, both $\mathcal{O}$ and $\tilde{\mathcal{O}}$ with flavor structures $[s^{(1)},s^{(1)},s^{(2)},s^{(2)}]$ and $[s^{(3)},s^{(3)},s^{(4)},s^{(4)}]$ will generate quantum numbers $I(J^{P})=2(0^{+})$, but strangeness $S^{(1)}=S^{(2)}=S^{(3)}=S^{(4)}=0$ for all four strange flavors and, hence, the correlation function (29) will decay asymptotically according to $e^{-2m_{\pi}t}$. Consequently, one has

[TABLE]

Inserting the masses $am_{K}\approx 0.274$ and $am_{\pi}\approx 0.137$ yields $\alpha=0.411$ in agreement with the numerical findings from Figure 8.

In conclusion the relative exponential increase of the statistical errors of the $2\times$ connected diagrams is not associated with the method of computation, but rather it is an intrinsic property of these diagrams.

4.5 Four-quark – two-meson correlation functions

4.5.1 $4\times$ connected $C_{25}$, $C_{26}$, $C_{35}$, $C_{36}$, $C_{45}$, $C_{46}$ ($2\equiv K\bar{K}\textrm{, point}$, $3\equiv\eta_{s}\pi\textrm{, point}$, $4\equiv Q\bar{Q}$, $5\equiv K\bar{K}\textrm{, 2part}$, $6\equiv\eta_{s}\pi\textrm{, 2part}$)

The correlation matrix elements $C_{25}$, $C_{26}$, $C_{35}$, $C_{36}$, $C_{45}$ and $C_{46}$ have different spin and color structures, but their $4\times$ connected diagrams are identical with respect to spacetime and, hence, can be discussed together. Promising methods for numerical computation are the following.

• •

$C_{25}$, $C_{35}$, $C_{45}$:

• $(a)$

2$\times$ one-end trick:

2 $u/d$ inversions, 2 $s$ inversions

• $(b)$

2 point-to-all-propagators:

12 $u/d$ inversions, 12 $s$ inversions

one-end trick:

1 $u/d$ inversion, 1 $s$ inversion

• $(c)$

4 point-to-all-propagators:

12 $u/d$ inversions, 12 $s$ inversions

Note that there are two variants of method $(b)$, when the operator $\mathcal{O}^{6}$ is considered, either applying the one-end trick to the $\pi$ meson (method $(b_{1})$) or to the $\eta_{s}$ meson (method $(b_{2})$).

• •

$C_{26}$, $C_{36}$, $C_{46}$:

• $(a)$

2$\times$ one-end trick:

1 $u/d$ inversions, 1 $s$ inversions

• $(b_{1})$

one-end trick:

1 $u/d$ inversion

2 point-to-all-propagators:

12 $s$ inversions

• $(b_{2})$

2 point-to-all-propagators:

12 $u/d$ inversions

one-end trick:

1 $s$ inversion

• $(c)$

4 point-to-all-propagators:

12 $u/d$ inversions, 12 $s$ inversions

A numerical comparison of these methods for all six cases is shown in Figure 10. The conclusions are similar as before. Method $(a)$ is most efficient, despite a relatively large number of additional stochastic noise terms is introduced, because the diagrams are averaged over space. In contrast to that method $(b)$ and method $(c)$ consider only a single sample.

It is interesting to note that even though in method $(b)$ stochastic noise terms are present, it performs better than method $(c)$. It seems that omitting a sum over space where only two quark lines end (as in method $(b)$), is more cost effective than omitting a sum over space where four quark lines end (as in method $(c)$).

Finally, our findings are consistent with the expectation that method $(b_{1})$ is superior to method $(b_{2})$. One reason is that the larger number of 12 inversions needed for the point-to-all propagators requires much less computing time for the heavier $s$ quark than for the lighter $u/d$ quarks. Secondly, the one-end trick has been found to be more efficient than point-to-all-propagators when applied to light quarks [26].

4.5.2 $2\times$ connected $C_{25}$, $C_{35}$, $C_{45}$, ($2\equiv K\bar{K}\textrm{, point}$, $3\equiv\eta_{s}\pi\textrm{, point}$, $4\equiv Q\bar{Q}$, $5\equiv K\bar{K}\textrm{, 2part}$)

Similarly to $C_{15}$, one of the $s$ quarks propagates within a timeslice and therefore using a sequential propagator is expected to be quite efficient. One can combine the sequential propagator either with the one-end trick or with point-to-all propagators.

• $(a)$

one-end trick (with 1 sequential propagator):

2 $u/d$ inversions, 1 $s$ inversion

stochastic timeslice-to-all propagator:

$\#t$ $s$ inversions

• $(b)$

2 point-to-all propagators (1 is a sequential propagator):

24 $u/d$ inversions, 12 $s$ inversions

stochastic timeslice-to-all propagator:

$\#t$ $s$ inversions

Treating the light quarks with the one-end trick, which only introduces a moderate number of additional stochastic noise terms but averages the diagram over space, is more efficient than point-to-all propagators (cf. Figure 11). This is consistent with the results for matrix element $C_{15}$ discussed in section 4.3.1.

4.5.3 $2\times$ connected $C_{26}$, $C_{36}$, $C_{46}$, ($2\equiv K\bar{K}\textrm{, point}$, $3\equiv\eta_{s}\pi\textrm{, point}$, $4\equiv Q\bar{Q}$, $6\equiv\eta_{s}\pi\textrm{, 2part}$)

In this case, at least one of the closed quark loops should be computed using a stochastic timeslice-to-all propagator. Therefore, promising methods are the following.

• $(a)$

one-end trick:

1 $u/d$ inversion

point-to-all propagator:

12 $s$ inversions

stochastic timeslice-to-all propagator:

$\#t$ $s$ inversions

• $(b)$

3 point-to-all-propagators:

12 $u/d$ inversions, 12 $s$ inversions

stochastic timeslice-to-all propagator:

$\#t$ $s$ inversions

A numerical comparison of these methods is shown in Figure 12. It is interesting to note that method $(a)$ is superior to method $(b)$ for $C_{26}$ and $C_{36}$, while the opposite is true for $C_{46}$. This is the only case we have investigated, where the spin and color structure of a diagram has a significant impact on the efficiency of the methods used for its computation.

Similarly to what has been observed for the $4\times$ connected $C_{25}$, $C_{26}$, $C_{35}$, $C_{36}$, $C_{45}$ and $C_{46}$ (cf. section 4.5.1) method (a) and method (b) perform on a similar level, even though a significantly larger amount of noise terms is introduced by method (a). This supports the conclusion that omitting a sum over space where only two quarks are present (as in method $(a)$), is more cost effective than omitting a sum over space where four quarks are present (as in method $(b)$).

4.5.4 Comparison of $4\times$ connected and $2\times$ connected diagrams

A comparison of $4\times$ connected and $2\times$ connected diagrams for $C_{25}$, $C_{35}$, $C_{45}$, $C_{26}$, $C_{36}$ and $C_{46}$ is shown in Figure 17. Results are qualitatively very similar to those shown in section 4.4.3 and the conclusions are essentially the same. Compared to their $4\times$ connected counterparts, the statistical errors of the $2\times$ connected diagrams increase exponentially as $e^{(2m_{K}-m_{\pi})t}$, when the temporal separation $t$ is increased. This independently of the methods used for the diagrams computation.

In this respect, note that the quality ratio $R^{2\times\textrm{con}}(t)$ is particularly small for $C_{26}$, $C_{36}$, $C_{46}$. This is so because even the most efficient method to compute the $2\times$ connected diagrams (method (a) in section 4.5.3) requires extensive use of stochastic techniques, while it still does not average the diagrams over space.

4.6 Two-meson – two-meson correlation functions

4.6.1 $4\times$ connected $C_{55}$, $C_{66}$ ($5\equiv K\bar{K}\textrm{, 2part}$, $6\equiv\eta_{s}\pi\textrm{, 2part}$)

The two disconnected parts of these diagrams are in fact identical to $C_{11}$. Therefore, the methods discussed in section 4.1.1 are applied to both parts.

• •

$C_{55}$:

• $(a)$

2$\times$ one-end trick:

2 $u/d$ inversions, 2 $s$ inversions

• $(b)$

2 point-to-all-propagators:

12 $u/d$ inversions, 12 $s$ inversions

one-end trick:

1 $u/d$ inversion, 1 $s$ inversion

Similar to $C_{26}$, $C_{36}$ and $C_{46}$ there are two variants of method $(b)$ for $C_{66}$, either applying the one-end trick to the $\pi$ meson (method (b1)) or to the $\eta_{s}$ meson (method (b2)).

• •

$C_{66}$:

• $(a)$

2$\times$ one-end trick:

1 $u/d$ inversions, 1 $s$ inversions

• $(b_{1})$

one-end trick:

1 $u/d$ inversion

2 point-to-all-propagators:

12 $s$ inversions

• $(b_{2})$

2 point-to-all-propagators:

12 $u/d$ inversions

one-end trick:

1 $s$ inversion

A numerical comparison of these methods is shown in Figure 14. The conclusions are essentially the same as in section 4.5.1 for $C_{25}$, $C_{26}$, $C_{35}$, $C_{36}$, $C_{45}$ and $C_{46}$. Method $(a)$ is more efficient, because the diagrams are averaged over space, while method $(b)$ considers only a single sample. Method $(b_{1})$ is superior to method $(b_{2})$, because it is cheaper to perform the larger number of 12 inversions for the heavier $s$ quark than for light $u/d$ quarks and the one-end trick has been found to be more efficient than point-to-all-propagators when applied to light quarks [26].

4.6.2 $2\times$ connected $C_{55}$ ($5\equiv K\bar{K}\textrm{, 2part}$)

Since both the $s$ quark at $t_{1}$ and the $s$ quark at $t_{2}$ propagate within the timeslice, the diagram can be computed rather efficiently using two sequential propagators combined with the one-end trick.

• $(a)$

2$\times$ one-end trick (with 2 sequential propagators):

2 $u/d$ inversions, $(\#t+1)$ $s$ inversions

4.6.3 $2\times$ connected $C_{66}$ ($6\equiv\eta_{s}\pi\textrm{, 2part}$)

At least one of the closed quark loops needs to be computed using a stochastic timeslice-to-all propagator. Therefore, there is only a single possibility to compute this type of diagram efficiently.

• $(a)$

one-end trick:

1 $u/d$ inversions

point-to-all propagator:

12 $s$ inversions

stochastic timeslice-to-all propagator:

$\#t$ $s$ inversions

4.6.4 $4\times$ connected $C_{56}$ ($5\equiv K\bar{K}\textrm{, 2part}$, $6\equiv\eta_{s}\pi\textrm{, 2part}$)

Similar to $4\times$ connected $C_{55}$ and $C_{66}$ there are two promising methods.

• $(a)$

2$\times$ one-end trick:

1 $u/d$ inversions, 1 $s$ inversions

• $(b_{1})$

one-end trick:

1 $u/d$ inversion

2 point-to-all-propagators:

12 $s$ inversions

• $(b_{2})$

2 point-to-all-propagators:

12 $u/d$ inversions

one-end trick:

1 $s$ inversion

A numerical comparison of these methods is shown in Figure 15. The conclusions are essentially the same as for the diagrams discussed in section 4.5.1 and section 4.6.1. Method $(a)$ is more efficient, because the diagrams are averaged over space, while method $(b)$ considers only a single sample.

4.6.5 $2\times$ connected $C_{56}$ ($5\equiv K\bar{K}\textrm{, 2part}$, $6\equiv\eta_{s}\pi\textrm{, 2part}$)

Again the triangular loop can be computed efficiently with a sequential propagator. There are several promising methods, which combine that propagator with point-to-all propagators, stochastic timeslice-to-all propagators and the one-end trick.

• $(a)$

one-end trick (with 1 sequential propagator):

2 $u/d$ inversion, $\#t$ $s$ inversions

point-to-all propagator:

12 $s$ inversions

• $(b)$

one-end trick (with 1 sequential propagator):

2 $u/d$ inversion, 1 $s$ inversions

stochastic timeslice-to-all propagator:

$\#t$ $s$ inversion

• $(c)$

2 point-to-all propagators (1 is a sequential propagator):

24 $u/d$ inversions, 12 $s$ inversions

stochastic timeslice-to-all propagator:

$\#t$ $s$ inversions

Note that there is another possibility very similar to method $(b)$, which also uses the one-end trick (with 1 sequential propagator) and a stochastic timeslice-to-all propagator, and requires the same number of inversions. The only difference is that the stochastic noise introduced by the one end-trick is located on the timslice of the closed quark loop (i.e. also at $t_{1}$), while in method $(b)$ it is located on the opposite timeslice (i.e. at $t_{2}$). We did not explore this additional method because in numerical studies of related diagrams we found that distributing the noise on the two timeslices normally results in smaller statistical errors.

A numerical comparison of the methods is shown in Figure 16. Method $(a)$ and method $(b)$ perform on a similar level. At first glance this is a bit surprising, because in previous subsections we typically observed that a method averaging a diagram over space is superior when only a moderate number of additional stochastic noise terms is introduced (cf. e.g. the correlation matrix element $C_{16}$ discussed in section 4.3.2, which has a very similar structure, and which can be computed by essentially the same combinations of techniques). Note, however, that this time both method $(a)$ and method $(b)$ require a comparably large number of $s$ inversions ($>\#t$), which is a consequence of the sequential propagator. This is not the case in many of the previously discussed cases, e.g. $C_{16}$. There the computation of the closed quark loop with a stochastic timeslice-to-all propagator requires only a single $s$ inversion, while here a point-to-all propagator needs to be computed, which requires a significantly larger amount of $s$ inversions i.e. 12 inversions.

4.6.6 Comparison of $4\times$ connected and $2\times$ connected diagrams

A comparison of $4\times$ connected and $2\times$ connected diagrams for $C_{55}$, $C_{56}$ and $C_{66}$ is shown in Figure 17. Results and conclusions are again very similar to previous comparisons. The statistical errors of the $2\times$ connected diagrams increase exponentially as $e^{(2m_{K}-m_{\pi})t}$ with respect to the temporal separation $t$ compared to their $4\times$ connected counterparts independent of the methods employed.

On the other hand, it is interesting to note that the quality ratio $R^{2\times\textrm{con}}(t)$ for $C_{55}$ is orders of magnitude larger than for all other diagrams. This seems to be the case due to a combination of the following two reasons. (1) The $2\times$ connected $C_{55}$ diagram is a single connected loop, i.e. it is not a combination of several disconnected pieces as e.g. for $C_{56}$ or $C_{66}$. (2) The method of computation used for $C_{55}$ introduces only a moderate number of noise terms (only 1 application of the one-end trick), but still averages the diagram over space, i.e. it is very efficient. Consequently, up to temporal separation $t/a<10$ the statistical error of the $2\times$ connected diagram is smaller than that of the $4\times$ connected diagram.

5 Conclusions

We have explored various methods to compute correlation functions of two- and four-quark interpolating operators of different structures including quark-antiquark pairs, mesonic molecules, diquark-antidiquark pairs and two independent mesons. Computing such correlation functions in an efficient way and with small statistical errors is essential to study tetraquark candidates or meson-meson scattering with lattice QCD.

The investigations have been performed in the context of a long-term project focused on the $a_{0}(980)$ meson, which might have a significant tetraquark component. We expect, however, that our conclusions hold at least qualitatively also for other tetraquark candidates and four-quark systems, e.g. for the positive parity $D_{s}$ mesons $D_{s0}^{\ast}(2317)$ and $D_{s1}(2460)$ or for certain $Z_{c}$ systems. Of course, the quality ratios presented and discussed in the previous section, might be somewhat different, when changing quark masses or quark flavors, when studying different lattice spacings or spacetime volumes, or when using another lattice QCD quark action. The general tendencies we have observed should, however, be the same also for other lattice QCD setups and other four-quark systems. Therefore, the investigations reported in this paper might provide a comprehensive overview of existing methods and helpful guidelines, which methods to choose for which kind of diagrams.

We have considered combinations of several techniques for quark propagator and correlation function computation: (A) point-to-all propagators, (B) stochastic timslice-to-all propagators, (C) the one-end trick and (D) sequential propagators. For applications of the distillation method to the computation of four-quark correlation functions we refer to [6, 7, 8].

For each diagram of the $6\times 6$ correlation matrix (7) we have implemented up to four promising methods, i.e. combinations of the above mentioned techniques (A) to (D). While details can be found in section 4, there are a few general observations which seem to hold for most diagrams:

• •

Methods, which average a diagram over space (i.e. which are not using any point-to-all propagators), are usually quite efficient. Exceptions are methods which extensively use stochastic methods ($M\geq 2$, where $M$ counts the number of stochastic timeslice-to-all propagators, and applications of the one-end trick, i.e. the number of stochastic sources) for diagrams, which can also be computed without introducing any stochastic noise (and without averaging over space). An example of the latter is method $(a)$ in section 4.2.1. A rule of thumb, which is fulfilled by the majority of the diagrams and methods explored in this work, is the following.

Method $(a)$ averaging a diagram over space is more efficient than method $(b)$ not averaging a diagram over space, if $M^{(a)}\leq M^{(b)}+1$.

Exceptions have been found in section 4.5.1, method $(b)$ versus method $(c)$ and section 4.5.3, method $(a)$ versus method $(b)$ (cf. the next item) and section 4.6.5, method (a) versus method $(b)$ (specific reasons are discussed in the section).

• •

When using one or more point-to-all propagators a spatial sum has to be omitted in the correlation function, i.e. the diagram is not averaged over space.

Omitting spatial sums over only two quark field operators, seems to result in much smaller statistical errors than omitting a sum over four quark field operators.

Examples are presented in section 4.5.1, method $(b)$ versus method $(c)$ and section 4.5.3, method $(a)$ versus method $(b)$.

• •

Sequential propagators are particularly useful for $2\times$ connected diagrams, where at least one operator \mathcal{O}^{\text{K\bar{K}, 2part}} is involved, i.e. to compute the strange quark propagating between different points in space within a timeslice.

• •

Stochastic methods are more efficient for lighter quarks (where the propagators exhibit rather strong statistical fluctuations due to the gauge links) and less efficient for heavy quarks (where the opposite is the case). An example is presented in section 4.6.5, method $(b_{1})$ versus method $(b_{2})$ (cf. also [26]).

A severe problem when computing a correlation matrix like (7), is that statistical errors of $2\times$ connected diagrams are increasing exponentially with respect to the temporal separation compared to their $4\times$ connected counterparts. This is independent of the method of computation and inherent to these diagrams as explained in detail in section 4.4.3. Their computation with sufficient statistical precision is, therefore, extremely challenging and the statistical errors of the correlation matrix elements will be dominated by the $2\times$ connected contributions.

In the past several studies have neglected $2\times$ connected diagrams (cf. e.g. [27, 28]). While in specific cases this might be justified, it is certainly not a general rule that these diagrams are negligible compared to their $4\times$ connected counterparts. A recent study [29] has stressed this and we moreover have found that the contribution of 2$\times$ connected diagrams is in many cases sizable. For illustration we show a few cases in Figure 18, where the ratio $|C_{jk}^{2\times\textrm{con}}|/(|C_{jk}^{4\times\textrm{con}}|+|C_{jk}^{2\times\textrm{con}}|)$ is plotted, i.e. the percentage of contribution of the $2\times$ connected diagram. For $C_{44}$ the $2\times$ connected diagram contributes more than the corresponding $4\times$ connected diagram to the whole correlation function, but also for the other examples, $C_{22}$, $C_{33}$ and $C_{56}$ both contributions are at least of similar importance.

On the one hand $2\times$ connected diagrams may not be neglected, on the other hand they are the reason why certain elements of the correlation matrix are very noisy. Therefore, performing a standard analysis by solving a generalized eigenvalue problem (cf. [30] and references therein) might not be the most successful strategy. Clearly, the numerical solution of a generalized eigenvalue problem requires the knowledge of all elements of a correlation matrix and, hence, the quality of its result will be dictated by the statistical errors of the most imprecise elements. A more promising strategy might be to fit a sum of decaying exponential functions, where the statistical error is rather related to the most precise elements of the correlation matrix, and, where in contrast to the generalized eigenvalue problem, strongly fluctuating matrix elements can be excluded from the analysis or will not affect the fitting result in a significant way. An advanced related variant, which we have recently explored in the context of this work is the AMIAS method [31].

After significantly increasing the statistical accuracy using the best method for each diagram as detailed in the main section of this paper, we plan to carry out a physics analysis using exponential fitting and the AMIAS method in the near future. The corresponding results regarding $a_{0}(980)$ meson will be part of an upcoming publication.

Appendix A Lattice setup

The numerical investigations presented in this paper are based on gauge link configurations generated by the PACS-CS collaboration [25]. The gluonic action is the Iwasaki gauge action [32] and the quark action for 2+1 flavors is a non-perturbatively $\mathcal{O}(a)$ improved Wilson quark action,

[TABLE]

with the field strength

[TABLE]

where $U_{x,\mu\nu}(j)$, $j=1,\ldots,4$ denotes the four plaquettes in the $\mu\nu$-plane attached to $x$.

We consider a single ensemble with gauge coupling $\beta=1.90$, corresponding lattice spacing $a\approx 0.091\,\textrm{fm}$ and pion mass $m_{\pi}\approx 300\,\textrm{MeV}$. Further details are listed in Table 1.

To enhance the ground state overlap of trial states $\mathcal{O}^{j}|\Omega\rangle$ generated by our interpolating operators (1) to (6), we apply standard smearing techniques. For the quark field operators $u_{\mathbf{x}}$, $d_{\mathbf{x}}$ and $s_{\mathbf{x}}$ we use Gaussian smearing with APE smeared spatial links (parameters $N_{\textrm{Gauss}}=50$, $\kappa_{\textrm{Gauss}}=0.5$, $N_{\textrm{APE}}=20$, $\alpha_{\textrm{APE}}=0.45$; cf. [33] for detailed equations). The widths of these smeared quark field operators can be estimated, $D\approx\sqrt{2N_{\textrm{Gauss}}\kappa_{\textrm{Gauss}}/(1+6\kappa_{\textrm{Gauss}})}a\approx 0.3\,\textrm{fm}$ (eq. (27) in [33]), i.e. they excite the corresponding quark fields in spherically extended regions with radius size $D$ centered around $\mathbf{x}$.

Appendix B Computation of the correlation matrix: technical aspects

B.1 $C_{12}$, $C_{13}$, $C_{14}$ ($1\equiv q\bar{q}$, $2\equiv K\bar{K}\textrm{, point}$, $3\equiv\eta_{s}\pi\textrm{, point}$, $4\equiv Q\bar{Q}$)

Even though the diagrams for $C_{12}$, $C_{13}$ and $C_{14}$ are identical (cf. Figure 2), the corresponding mathematical expressions differ due to a different spin and color structure of the interpolating operators (2), (3) and (4),

[TABLE]

The three combinations of techniques investigated in section 4.2.1 are rather independent of the spin and color structure introduced by either \mathcal{O}^{\text{K\bar{K}, point}}, \mathcal{O}^{\text{\eta_{s}\pi, point}} or \mathcal{O}^{\text{Q\bar{Q}}} and, hence, can be applied in essentially the same way to $C_{12}$, $C_{13}$ and $C_{14}$. Therefore, we present equations only for $C_{12}$.

• •

Method $(a)$:

[TABLE]

• •

Method $(b)$:

[TABLE]

• •

Method $(c)$:

[TABLE]

B.2 $C_{15}$ ($1\equiv q\bar{q}$, $5\equiv K\bar{K}\textrm{, 2part}$)

With respect to spin and color $C_{15}$ is identical to $C_{12}$ (cf. eq. (35)). The crucial difference is that the $K$ and the $\bar{K}$ meson are not created at the same point in space, but at different points ${\bf x}$ and ${\bf z}$,

[TABLE]

Equations corresponding to the two combinations of techniques investigated in section 4.3.1 are the following.

• •

Method $(a)$:

[TABLE]

where

[TABLE]

are the analogs of (20) and (25) for combining the one-end trick with a sequential propagator, which are straightforward to derive.

• •

Method $(b)$:

[TABLE]

where $\psi^{(u;s)}$ is the sequential propagator $G^{(u)}\gamma_{5}G^{(s)}$ according to (26).

Acknowledgments

We thank Mario Gravina for his contributions in an early stage of this project. We acknowledge useful conversations with Stefan Meinel and Christian Wiese.

J.B. and M.W. acknowledge support by the Emmy Noether Programme of the DFG (German Research Foundation), grant WA 3000/1-1. M.D.B. is grateful to the Theoretical Physics Department at CERN for the hospitality and support.

This work was supported in part by the Helmholtz International Center for FAIR within the framework of the LOEWE program launched by the State of Hesse.

This work was cofunded by the European Regional Development Fund and the Republic of Cyprus through the Research Promotion Foundation (Project Cy-Tera ΝΕΑ $\Upsilon\Pi\text{O}\Delta\text{OMH}/\Sigma\text{TPATH}/0308/31$ ) by the grand cypro914.

Calculations on the LOEWE-CSC and on the on the FUCHS-CSC high-performance computer of the Frankfurt University were conducted for this research. We would like to thank HPC-Hessen, funded by the State Ministry of Higher Education, Research and the Arts, for programming advice.

Computations have been performed using the Chroma software library [34].