Lattice QCD investigation of the structure of the $a_0(980)$ meson

We investigate the quark content of the scalar meson $a_0(980)$ using lattice QCD. To this end we consider correlation functions of six different two- and four-quark interpolating fields. We evaluate all diagrams, including diagrams, where quarks propagate within a timeslice, e.g. with closed quark loops. We demonstrate that diagrams containing such closed quark loops have a drastic effect on the final results and, thus, may not be neglected. Our analysis shows that in addition to the expected spectrum of two-meson scattering states there is an additional energy level around the two-particle thresholds of $K + \bar{K}$ and $\eta + \pi$. This additional state, which is a candidate for the $a_0(980)$ meson, couples to a quark-antiquark as well as to a diquark-antidiquark interpolating field, indicating that it is a superposition of an ordinary $\bar{q} q$ and a tetraquark structure. The analysis is performed using AMIAS, a novel statistical method based on the sampling of all possible spectral decompositions of the considered correlation functions, as well as solving standard generalized eigenvalue problems.


I. INTRODUCTION
The mass ordering of the light scalar mesons σ, κ, f 0 (980) and a 0 (980), as observed in experiments, is inverted compared to what is expected based on the conventional quark-antiquark structure. A possible explanation of this mass ordering is to interpret these states as tetraquarks. Assuming such a four-quark structure the expected mass ordering is consistent with experimental results and the degeneracy of the f 0 (980) meson and the a 0 (980) meson is straightforward to understand [1].
Several lattice QCD studies of light scalar mesons have been published in the last couple of years (cf. e.g. [2][3][4][5][6][7][8][9][10]). Regarding the a 0 (980) meson the elaborate study presented in Ref. [10] using Lüscher's finite volume approach is of particular interest. The authors find a resonance close to the K +K threshold, which they interpret as the a 0 (980) meson. Moreover, cf. Ref. [11], where the same lattice data has been analyzed using chiral effective field theory.
In this work we continue our lattice QCD investigation of the a 0 (980) meson [12][13][14][15][16][17][18][19][20] with particular focus on clarifying its quark structure, i.e. whether it is of quarkantiquark type or rather a four-quark state. We include the computation of all diagrams, where quarks propagate within a timeslice, e.g. diagrams with closed quark loops, which were neglected in many previous investigations of scalar mesons. We show that these contributions are important and lead to the appearance of an additional qqlike state around the two-particle thresholds of η s +π and K +K, which could correspond to the a 0 (980) meson. On a technical level we apply the Athens Model Inde-pendent Analysis Scheme (AMIAS), an analysis method based on statistical concepts for extracting excited states from correlation functions. AMIAS is a novel analysis method, which has recently been used in a study of the nucleon spectrum [21]. AMIAS utilizes all the information encoded in the correlation function with the particular advantage of exploiting also data at small temporal separations, where statistical errors are typically small. In addition to AMIAS, we also use the standard generalized eigenvalue problem (GEVP) method, i.e. we solve a generalized eigenvalue problem and extract the spectrum from effective mass plateaus. Both GEVP and AMIAS provide information on the relative importance of the considered interpolating fields, which include quarkantiquark, four-quark and two-meson structures. Combining both methods allows to check the robustness of our results. This paper is organized as follows: In section II we describe the lattice techniques used and briefly discuss the inclusion and importance of diagrams, where quarks propagate within a timeslice. We demonstrate that these diagrams are essential both on a qualitative and quantitative level and, hence, may not be neglected. In section III we present the spectral decomposition of two-point correlation functions and explain possible complications arising due to the presence of multi-particle states. A short description of AMIAS is provided in section IV. In section V we present the analysis of a 6×6 correlation matrix containing various interpolating fileds using both GEVP and AMIAS. In section VI we summarize our findings and give our conclusions.

II. LATTICE SETUP AND TECHNIQUES
To investigate the a 0 (980) meson, which has quantum numbers I(J P ) = 1(0 + ) and a mass around 980 MeV [22], we consider a 6 × 6 correlation matrix The interpolating fields O j , j = 1, . . . , 6 have either a two-quarkdu or a four-quarkduss structure. The fourquarks can be arranged as a meson-meson interpolating field (either a bound pair of mesons or two separated mesons) or in a diquark-antidiquark combination. We consider the interpolating fields where C is the charge conjugation matrix. The normalization factors N j are chosen such that C jj (t = a) = 1 (no sum over j), i.e. in a way that the six interpolating fields generate trial states with similar norm. All six interpolating fields couple to the a 0 (980) and other states with the same quantum numbers. For example, the interpolating fields O 5 and O 6 mostly generate the two-meson states K +K and η + π, respectively, which are expected to have masses close to the a 0 (980) state. In contrast to O 5 and O 6 , where the two mesons both have zero momentum, the interpolating fields O 2 and O 3 represent two mesons at the same point in space, i.e. resemble mesonic molecules. Previous results using these interpolating fields and Wilson clover fermions can be found in [18][19][20]. In this work several significant improvements have been carried out: (i) we have improved the statistical accuracy of the correlation matrix C jk (t), (ii) we include the propagation of strange quarks within a timeslice, (iii) we consider all 36 elements of the correlation matrix, (iv) we analyze the correlation matrix with both the standard GEVP and the AMIAS method.
We use an ensemble of around 500 gauge link configurations generated with N f = 2 + 1 dynamical Wilson clover quarks and the Iwasaki gauge action produced by the PACS-CS Collaboration [23]. The lattice size is 64 × 32 3 with lattice spacing a = 0.0907 (14) fm and pion mass m π ≈ 300 MeV. To optimize the coupling of the interpolating fields to the low-lying energy eigenstates, quark fields are Gaussian smeared with APE smeared spatial gauge links (cf. e.g. [24] for detailed equations).
For each diagram of the correlation matrix C jk (t) we have implemented and compared various combinations of techniques including point-to-all propagators, stochastic timeslice-to-all propagators, the one-end trick and sequential propagators. Based on the results of these comparisons we have chosen for each diagram individually the most efficient combination of techniques for our computation of C jk (t) used in the physics analysis in section V. Two example diagrams, which form the matrix element C 44 (t), are shown in Fig. 1. For the diagram on the left without closed quark loops we use four point-to-all propagators, i.e. its computation is technically simple and the statistical errors are quite small. Significantly more difficult is the diagram on the right with closed quark loops. In this case the use of four point-to-all propagators is not possible, because of the two closed quark loops. We found that using three point-to-all propagators and a stochastic timeslice-to-all propagator for one of the closed quark loops is the most efficient way of computing this particular diagram. For an extensive discussion we refer to [25], where each diagram of the 6×6 correlation matrix is studied in detail. Finding efficient methods is particularly important for diagrams, where quarks propagate within a timeslice, e.g. diagrams containing closed quark loops. These diagrams are significantly more noisy than their counterparts, where quarks do not propagate within a timeslice: their noise-to-signal ratio grows exponentially with increasing temporal separation as discussed in [25].
Considering diagrams, where quarks propagte within a timeslice, is vital for any solid study of a 0 (980), because they lead to non-vanishing correlations between two-quark and four-quark interpolating fields, i.e. allow forss creation and annihilation. Moreover, for correlation functions of two four-quark interpolating fields their contribution is sizable and should not be neglected as it has been done in the past, e.g. in [6,13]. This is demonstrated by Fig. 2, where we show C 44 (t) both with (blue dots) and without (red dots) closed quark loops. Similar findings have been reported in Refs. [9,26].

III. CORRELATION FUNCTIONS ON A PERIODIC LATTICE
A correlation function computed on a lattice with periodic temporal direction of extension T can be expanded according to with E m,n = E m − E n . For the interpolating fields defined by Eqs.
Moreover, c j m,n (c k m,n ) * is real and without restriction one can choose real c j m,n , which we do in the following. Consequently, Eq. (9) can be simplified to Since the correlators are symmetric with respect to the reversal of time, C jk (t) = C jk (T − t), it is sufficient to restrict the follwing discussion to temporal separations 0 ≤ t ≤ T /2. For j = k, for sufficiently large T , where Z ≈ e −EΩT , and for sufficiently large t, Eq. (11) reduces to if the correlation function is not contaminated by effects related to multi-particle states as discussed below (|Ω and |0 denote the vacuum and the lowest state in the I(J P ) = 1(0 + ) sector probed by the interpolating fields defined by Eqs.
(2) to (7), respectively). Consequently, the energy difference E0 ,Ω can be extracted by fitting to the lattice QCD results for the correlation function C jj (t) at sufficiently large t with fitting parameters E0 ,Ω and A. Alternatively, one can solve the equation with respect to E eff (t), where E eff (t) ≈ E0 ,Ω = const for large t. In other words, a plateau-like behavior of E eff (t) indicates the mass E0 ,Ω . A common method to extract several energy levels from an N × N correlation matrix is to solve the GEVP (cf. e.g. [27] and references therein) where C(t) is the correlation matrix with entries C jk (t), j, k = 1, . . . , N , v m (t, t 0 ) the eigenvector corresponding to the eigenvalue λ m (t, t 0 ) and t 0 ≥ a a parameter, where a typical choice is t 0 = a. N effective energies E eff,m (t) can be obtained by solving for each eigenvalue λ m (t, t 0 ), m = 0, . . . , N − 1. In practice, however, effective energies E eff,m (t) often exhibit strong statistical fluctuations, in particular for large t and m > 0, rendering a reliable identification of plateaus and extraction of masses Em ,Ω difficult (|0 , |1 , . . . denote the lowest states in the I(J P ) = 1(0 + ) sector probed by the interpolating fields defined by Eqs.
(2) to (7)). Therefore, besides using the GEVP, we also employ an alternative analysis approach, AMIAS, which is discussed in detail in section IV. When low-lying multi-particle states are present in the investigated sector, the determination of masses becomes even more difficult. For example in the I(J P ) = 1(0 + ) sector the lowest state is expected to be a two-meson state composed of an η meson and a π meson. Clearly, the interpolating fields defined by Eqs. (2) to (7) do not only generate an η + π state, when applied to the vacuum |Ω , but also yield non-vanishing matrix elements π|O j |η and η|O j |π(u ↔ d) , i.e. anihilate an η meson and generate a π meson and vice versa. Consequently, a significant contribution to C jk (t) is (cf. Eq. (11)). Similarly, a term providing information about the mass of the η + π state is Assuming coefficients |c j π,η | ≈ |c j η+π,Ω | it is easy to see that in the region of t ≈ T /2 the two terms are comparable in magnitude. Therefore, Eq. (17) needs to be taken into account, when extracting masses from the correlation matrix at large temporal separations. Only for sufficiently small temporal separations t the contribution from Eq. (17) may be neglected, since the ratio of Eq. (18) and Eq. (17) grows exponentially ∝ e −2mπ(t−T /2) with decreasing t. This is illustrated in Fig. 3, where we show the effective mass E eff (t) defined in Eq. (14) for lattice QCD results for the correlation function C 66 (t) corresponding to the interpolating field O ηsπ, 2part . Note that we have neglected closed quark loops and quark propagation within a timeslice for this computation, to obtain sufficiently precise results at large temporal separations. We also compare with the effective mass (14), where C jj (t) is the sum of Eqs. (17) and (18), the analytically expected dominating terms at large temporal separations t, i.e.
with m η a = 0.364, m π a = 0.138 (the η meson and π meson masses in our lattice setup) and m η+π = m η +m π , and find excellent agreement as can be seen in Fig. 3.

IV. BASICS OF THE AMIAS METHOD
In this section we summarize the basics of the AMIAS method [21,28,29]. A more detailed description and application to the excited nucleon spectrum can be found in Ref. [21].
Lattice QCD results for correlation functions C jk (t) with O j , j = 1, . . . , 6 defined in Eq. (2) to Eq. (7) can be parameterized according to Eq. (11). After approximating Z ≈ e −EΩT and truncating the sums m,n to a limited number of important terms, typically terms corresponding to small energy differences |E m,n |, one obtains fit functions of the form which are appropriate for sufficiently large t and T (the factor 2 is included, because of u ↔ d flavor symmetry; states m, n and states n(u ↔ d), m(u ↔ d) contribute with identical terms and, thus, should be combined in the fit functions). In practice it is convenient to use the equivalent fit functions instead of Eq. (20), where X is a superindex replacing {m, n}, a X ≡ √ 2e −(Em,Ω+En,Ω)T /4 c j m,n and E X ≡ E m,n . The fit parameters E X and a j X are real. In the following these fit parameters are collectively denoted by A r .
AMIAS determines a probability distribution function (PDF) Π(A r ) for each fit parameter A r . The estimates for the values of the fit parameters and their uncertainties are the expectation values and the standard deviations of the corresponding PDFs, AMIAS is able to handle a rather large number of parameters using Monte Carlo techniques, i.e. it is suited to study several energy eigenstates, if the lattice QCD results for correlation functions are sufficiently precise. The PDF for the complete set of fit parameters is defined by with appropriate normalization N and which is the well-known χ 2 used e.g. in χ 2 minimizing fits. C jk (t) denotes the correlation functions computed using lattice QCD with corresponding statistical errors σ jk (t), while C fit j,k (t) is given by Eq. (20). To obtain the PDF Π(A r ) for a specific fit parameter A r , one has to integrate Eq. (24) over all other parameters. In particular, the probability for the parameter A r to be inside the interval This multi-dimensional integral can be computed with standard Monte Carlo methods. We implemented a parallel-tempering scheme combined with the Metropolis algorithm as described in detail in Ref. [21]. The paralleltempering scheme prevents the algorithm from getting stuck in a region around a local minimum of χ 2 and guarantees ergodicity of the algorithm.
A common choice for t min and t max in Eq. (25) is t min = a and t max = T /2. With AMIAS one can determine the maximum number of parameters, to which the lattice QCD results for the correlation functions are sensitive on, i.e. the number of terms considered in the truncated sum in Eq. (21). The strategy is to increase the number of fit parameters, until there is no observable change in the PDFs for the low-lying energy eigenstates of interest. For a detailed example cf. Ref. [21]. For correlation functions with large statistical errors it might be helpful to also vary t min and t max and check the stability of the results, as e.g. done in section V C.
As an example we show in Fig. 4 the AMIAS analysis of the correlation function C 66 (t) corresponding to the interpolating field O ηsπ, 2part , which we already discussed at the end of section III. We have found that using t min = a, t max = T /2 = 32a and three terms in the truncated sum of the fit function (Eq. (20)) allows to determine two energy differences rather precisely, E 0 corresponding to E η,π a = (E η − E π )a ≈ 0.23 and E 1 corresponding to E η+π,Ω = (E η+π − E Ω )a ≈ m η + m π ≈ 0.50. E 2 should not be interpreted as a specific energy difference, since it is unstable under a variation of the number of terms in the truncated sum. In this section we analyse the 6 × 6 correlation matrix discussed in section III and various submatrices using both the GEVP method and the AMIAS method. The latter has proven to be particularly suited to study excited states [21]. Both methods yield consistent results, which we consider to be an important cross-check, in particular due to the fact that the signal-to-noise ratios of the elements of the correlation matrix grow rapidly with increasing temporal separations.
To extract energy differences E m,n with the GEVP method, effective energies are computed as defined in Eq. (16) (we always use t 0 = a in Eq. (15)). The plateau values at sufficiently large temporal separation correspond to E m,n as determined by fitting constants. The same energy differences E m,n are computed with the AMIAS method using the fit function (21) as explained in detail in section IV.
The components of the eigenvectors v m (t, t 0 ) obtained by solving a GEVP (15) provide information about the structure of the corresponding energy eigenstates: for sufficiently large t, where the approximate equality sign denotes the expansion of the energy eigenstate |m within the subspace generated by the interpolating fields via O j † |Ω . We always normalize the eigenvectors according to (v m (t, t 0 )) 2 = 1, before plotting them. One can easily convert the amplitudes a j X extracted with AMIAS to c j Ω,m introduced in section III via Since c j Ω,m = m|O j † |Ω , these are the coefficients of the expansions of the trial states O j † |Ω in terms of the energy eigenstates |m , i.e.
More interesting, however, is invertig Eq. (29) and writing the extracted energy eigenstates in terms of the trial states, One can show that the matrix formed by the coefficients v j m is the inverse of the matrix formed by the coefficients c j Ω,m up to exponentially small corrections, i.e.
Note that the coefficientsṽ j m are equivalent to the eigenvector components v j m (t, t 0 ) obtained by solving a GEVP (15) and, thus, the resuls from the two methods can be compared in a meaningful way, after choosing the same normalization (ṽ m ) 2 = 1.

B. Neglecting quark propagation within a timeslice
At first we neglect diagrams, where quarks propagate within a timeslice. Thus,ss pair creation and annihilation is excluded. Consequently, the quark-antiquark interpolating field O 1 probes a different sector than the four-quark interpolating fields O 2 to O 6 . Thus, within this subsection we restrict the analysis to the 5 × 5 correlation matrix formed by the interpolating fields O 2 to O 6 . Neglecting quark propagation within a timeslice leads to results with rather small statistical uncertainties and, thus, allows to cross-check our analysis methods and to compare with our previous study of the a 0 (980) meson [13], where we used a different lattice discretization and setup. Note, however, that contributions from diagrams, where quarks propagate within a timeslice, are sizeable (cf. e.g. Fig. 1 and [ 9,25]) and have to be taken into account to arrive at full QCD results, which can be compared to experimental data in a meaningful way.
In Fig. 5 we show effective energies obtained by solving the GEVP for the 5×5 correlation matrix. In the absence of quark propagation within a timeslice all five effective energies exhibit convincing plateaus and the corresponding energy differences can be determined in a straightforward and precise way by fitting constants at large temporal separations, e.g. at t > ∼ 10 a . . . 15 a. The plateaus are consistent with the two-particle thresholds • 2m K ≈ 1192 MeV as shown in Fig. 5 (note that the light quarks are unphysically heavy, corresponding to m π ≈ 300 MeV, m K ≈ 596 MeV and m η ≈ 792 MeV). The lowest momentum excitations are given by • 2(m 2 K + p 2 min ) 1/2 ≈ 1466 MeV, where p min = 2π/L denotes one quantum of momentum and L is the spatial lattice extent. We do not observe any sign of an additional energy level near the expected mass of the a 0 (980) meson, i.e. in the region around 1100 MeV to 1200 MeV, which could be interpreted as the a 0 (980) meson [30].
The same energy differences are found, when using the AMIAS method. The corresponding PDFs, generated with the fit function given in Eq. (21) and eight terms in the truncated sum [31], is shown in Fig. 6. To generate quantitative results including uncertainties we compute the mean values and widths of the PDFs. We find perfect agreement with the expected energies of states with two particles at rest, (expectation: m π + m η ≈ 1092 MeV), • E 3 = 1194(9) MeV (expectation: 2m K ≈ 1192 MeV),  (7)).
and fair agreement, when these particles have one quantum of relative momentum, (expectation: 2(m 2 K + p 2 min ) 1/2 ≈ 1466 MeV). When including correlation matrix data for temporal separations around t = T /2, which has small statistical uncertainties, when quark propagation within a timeslice is neglected, AMIAS also finds two significantly smaller energy differences E 0 and E 1 in the region of m K − mK = 0 and m η − m π . This is expected and has been discussed in detail in section III.
To support our interpretation of the states corresponding to E m , m = 2, . . . , 5, we show the corresponding coefficientsṽ j m in Fig. 6. The states corresponding to E 2 and E 3 are clearly two-particle states η + π and K +K with both mesons at rest, since the coefficients v j m indicate a strong domination of interpolating fields O 5 = O KK, 2part and O 6 = O ηsπ, 2part . The states corresponding to E 4 and E 5 exhibit significant contributions from interpolating fields O 2 = O KK, point and O 3 = O ηsπ, point , which is consistent with their interpretation as two-particle states with one quantum of relative momentum. A similar GEVP analysis [20] provides consistent results, i.e. eigenvector components v j m (t, t 0 ), which are in agreement with the coefficientsṽ j m .

C. Taking into account quark propagation within a timeslice
When including quark propagation within a timeslice, correlation functions of the quark-antiquark interpolat-  (7)).
ing field (Eq. (2)) and the four-quark interpolating fields (Eqs. (3) to (7)) are non-zero. Thus, all interpolating fields probe the same sector and the quark-antiquark interpolating field can be included in the analysis.
In Fig. 7 we show effective energies obtained by solving the GEVP for two different correlation matrices: In comparison to the effective energies from Fig. 5, where quark propagation within a timeslice has been neglected, statistical uncertainties drastically increase. For a detailed discussion concerning the reason for this increase cf. e.g. [25], section 4.4.3. The signal is essentially lost for temporal separations t > ∼ 7 a. For temporal separations t < ∼ 6 a effects related to multiparticle states, as discussed in section III, are negligible and, thus, can be ignored throughout this subsection. Due to the large statistical uncertainties, the identification of plateaus and energy differences is rather difficult. Nevertheless, there is a clear qualitative difference between the results from the 4 × 4 and the 6 × 6 correlation matrix. In the 4 × 4 plot there seem to be only two low-lying states around 1100 MeV to 1200 MeV, while the next two states are significantly above, somewhere in the energy region of momentum excitations. Thus, the observed spectrum is consistent with the expected spectrum of two-meson states.
In contrast to that, in the 6×6 plot, i.e. when taking also the quark-antiquark and the diquark-antidiquark interpolating fields into account, there is an additional third state around 1200 MeV, which could correspond to the a 0 (980) meson. Similar plots have been obtained for 5×5 correlation matrices, where either the quark-antiquark interpolating field is considered, but not the diquarkantidiquark interpolating field, or vice versa. However, the effective energy of the additional state is somewhat larger and has larger statistical errors, when O 1 = O qq is excluded. In summary, these results indicate that in addition to the expected two-meson states there is another low-lying state in the region of 1100 MeV to 1200 MeV, which is predominantly of quark-antiquark type with a smaller, but non-negligible diquark-antidiquark component. This interpretation is confirmed by the squared eigenvector components (v j m (t, t 0 )) 2 , which are plotted for the three lowest states in Fig. 8 for t = 5a (within statistical errors these eigenvector components are independent of t). The lowest state is mostly an η + π two-meson state, the second-lowest state mostly a quark-antiquark state and the third-lowest state a K +K two-meson state. The lowest two states also contain a small diquarkantidiquark contribution. To support these findings we have also determined energy differences from the full 6 × 6 correlation matrix with the AMIAS method. The corresponding PDFs for t min = 2a, t max = 8a and eight terms in the truncated sum of the fit function (21) is shown in Fig. 9. There are two clear peaks consistent with the expected energies of two-particle η + π and K +K states with both mesons at rest: (expectation: m π + m η ≈ 1092 MeV), • E 2 = 1192(11) MeV (expectation: 2m K ≈ 1192 MeV).
Moreover, there is another peak signaling the existence of an additional third state in the same energy region with i.e. significantly below the expectation for the lowest twoparticle states with one quantum of relative momentum, (m 2 η + p 2 min ) 1/2 + (m 2 π + p 2 min ) 1/2 ≈ 1422 MeV.  To clarify the structure of the extracted energy eigenstates, we also show the corresponding squared coefficients (c j Ω,m ) 2 , m = 0, 1, 2 in Fig. 9 [32]. The two lowest states corresponding to the energy differences E 0 and E 1 contribute to the two-meson trial state O 6 |Ω = O ηsπ, 2part |Ω as well as to the quark-antiquark and the diquark-antidiquark trial states O 1 |Ω = O qq |Ω and O 4 |Ω = O QQ |Ω . The second excitation corresponding to the energy difference E 2 contributes almost exclusively to the two-meson trial state O 5 |Ω = O KK, 2part |Ω . These AMIAS results are in agreement with the GEVP results discussed above and, thus, confirm our previous interpretation that there is an additional state in the energy region of 1100 MeV to 1200 MeV, which could correspond to the a 0 (980) meson. This additional state has both a quark-antiquark and a diquark-antidiquark component.
One of the advantages of the AMIAS method, compared to e.g. the GEVP method, is that one can use an arbitrary selection of correlation matrix elements for an analysis. To check the correctness and stability of our results, in particular the existence of an additional low-lying state with significant quark-antiquark component, we compare the PDFs for energy differences based on two different analyses and sets of correlation matrix elements in Fig. 10: (a) the full 6 × 6 correlation matrix (same data as in Fig. 9), (b) as in (a), but the diagonal element C 11 (t) is excluded from the analysis, i.e. j = k = 1 is omitted in the sum j,k in Eq. (25); this implies that the correlation of the quark-antiquark interpolating field O 1 = O qq with itself is excluded, while correlations with the other four-quark interpolating fields are still included. Fig. 10 10. PDFs for energy differences obtained with the AMIAS method (quark propagation within a timeslice taken into account; full 6×6 correlation matrix and the same matrix without C11(t)).

VI. SUMMARY AND CONCLUSIONS
We have computed the low-lying spectrum of the I(J P ) = 1(0 + ) sector using lattice QCD. To this end, we have considered a variety of interpolating fields representing quark-antiquark and four-quark bound states as well as two-meson scattering states. In addition to the expected two-meson scattering states we found another state in the energy region of 1100 MeV to 1200 MeV, which is a candidate for the a 0 (980) meson. This state is predominantly generated by the quark-antiquark interpolating field, but also receives sizable contributions from the diquark-antidiquark interpolating field, i.e. a likely interpretation is that it is mainly a quark-antiquark state with a minor tetraquark component. To some extent this is supported by a computation, where we have neglected quark propagation within a timeslice. Then the quark-antiquark interpolating field decouples from the four-quark interpolating fields and an analysis of the four-quark correlation matrix only yields the expected two-meson scattering states.
We have cross-checked and confirmed our results by applying two different analysis methods, the well-known GEVP method and the rather new AMIAS method. The latter is particularly suited to handle lattice QCD data with large uncertainties, because it also exploits correlation functions at small temporal separations.
A current source of systematic error, which we plan to eliminate in the future, is the unphysically large u and d quark masses corresponding to m π ≈ 300 MeV. A resonance study of a 0 (980) using Lüscher's finite volume approach was presented in Ref. [10] using, however, u and d quark masses corresponding to m π ≈ 391 MeV. In Ref. [11], which is an analysis of the same lattice QCD data using chiral effective field theory, it was pointed out that lighter u and d quark masses are essential to obtain more precise results. However, the extension of the method to lighter u and d quark masses is very demanding, since it requires the identification of all energy states below that of the a 0 (980). Therefore, the work presented here constitutes an important preparatory step.
Another obvious direction to extend this research is to investigate the D * s0 (2317) meson, which is also considered as a tetraquark candidate. Our techniques and code can be used with only minor changes, since both a 0 (980) and D * s0 (2317) have the same quantum numbers J P = 0 + and the same flavor structure, i.e. a quark-antiquark pair of different flavorq 1 q 2 and possibly another quarkantiquark pair of the same flavorq 3 q 3 . Such an investigation could be of particular interest, because lattice QCD studies like those presented in [33][34][35], which do not include four-quark interpolating fields, found masses significantly above the experimental result. Other lattice QCD studies, e.g. [36,37], which include four-quark interpolating fields found a state below the DK threshold much closer to the experimental result. Thus, it would be interesting to further explore the existence and mass of a D * s0 (2317) state below the DK threshold, by varying the set of interpolating fields considered in the analysis, and to investigate its internal structure, e.g. by employing a variety of tetraquark interpolating fields.