Tetraquark interpolating fields in a lattice QCD investigation of the D∗ s0(2317) meson

Constantia Alexandrou, 2 Joshua Berlin, Jacob Finkenrath, Theodoros Leontiou, and Marc Wagner Department of Physics, University of Cyprus, P.O. Box 20537, 1678 Nicosia, Cyprus Computation-based Science and Technology Research Center, The Cyprus Institute, 20 Kavafi Street, 2121 Nicosia, Cyprus Goethe-Universität Frankfurt am Main, Institut für Theoretische Physik, Max-von-Laue-Straße 1, D-60438 Frankfurt am Main, Germany Computation-based Science and Technology Research Center, The Cyprus Institute, 20 Kavafi Street, 2121 Nicosia, Cyprus Department of Mechanical Engineering, Frederick University, 1036 Nicosia, Cyprus We investigate the D∗ s0(2317) meson using lattice QCD and considering correlation functions of several c̄s two-quark and c̄s(ūu+d̄d) four-quark interpolating fields. These interpolating fields generate different structures in color, spin and position space including quark-antiquark pairs, tetraquarks and two-meson scattering states. For our computation we use an ensemble simulated with pion mass mπ ≈ 0.296 GeV and spatial volume of extent 2.90 fm. We find in addition to the expected spectrum of two-meson scattering states another state around 60 MeV below the DK threshold, which we interpret as the D∗ s0(2317) meson. This state couples predominantly to a quark-antiquark interpolating field and only weakly to a DK two-meson interpolating field. The coupling to the tetraquark interpolating fields is essentially zero, rendering a tetraquark interpretation of the D∗ s0(2317) meson rather unlikely. Moreover, we perform a scattering analysis using Lüscher’s method and the effective range approximation to determine the D∗ s0(2317) mass for infinite spatial volume. We find this mass 51 MeV below the DK threshold, rather close to both our finite volume result and the experimentally observed value.


I. INTRODUCTION
The D * s0 (2317) meson with quantum numbers I(J P ) = 0(0 + ), strangeness S = ±1 and charm C = S has mass m D * s0 = 2.3178(5) GeV, around 45 MeV below the DK threshold [1][2][3][4]. This experimental result is in contrast to theoretical predictions from quark models (see e.g. Refs. [5][6][7]), where the D * s0 (2317) meson is treated as acs quark-antiquark pair, which leads to a significantly larger mass in the range of 100 MeV to 200 MeV above the experimental value. Because of that discrepancy, there is an ongoing debate about the quark composition of the D * s0 (2317) meson. Besides a standard quark-antiquark structure it could also have a four-quark structure. For example Refs. [8][9][10] propose a tetraquark structure, while Ref. [11] provides arguments against such a scenario. Another possibility is a DK mesonic molecule structure as e.g. suggested by Refs. [12,13]. This picture is also supported by recent papers [14][15][16][17], where DK molecular components from around 60% to 75% are found. Other interesting approaches, which are able to explain the surprisingly low mass of the D * s0 (2317) meson, are e.g. presented in Ref. [18], where the D * s0 (2317) meson is a standardcs configuration with the coupling to the nearby DK threshold taken into account, and in Refs. [19][20][21], which are based on an SU(3) chiral Lagrangian. For a more detailed discussion of the properties of the D * s0 (2317) meson and existing literature we refer to the review articles [22,23].
Early quenched lattice QCD studies [24][25][26][27][28][29][30] of the D * s0 (2317) meson found masses significantly larger than the experimental result, similar to quark model predictions. There are also more recent lattice QCD studies [31][32][33][34][35][36][37][38][39], where only quark-antiquark interpolating fields of flavor structurecs were taken into account. The majority of these studies also finds masses for the D * s0 (2317) meson, which are larger than the experimental value, in particular if extrapolations to physical quark masses and to the continuum were performed (see e.g. Ref. [37]). If, however, in addition to quark-antiquark interpolating fields also two-meson DK interpolating fields are included, as done in the recent precision lattice QCD computations presented in Refs. [40][41][42], D * s0 (2317) masses below the DK threshold and close or consistent with the experimental result are found. In these investigations almost physical u and d quark masses were used, corresponding to pion masses m π ≈ 0.156 GeV and m π ≈ 0.150 GeV, respectively, and Lüscher's method was employed, to obtain the meson mass at infinite spatial volume. In this context it is also interesting to mention two closely related lattice QCD investigations. In Ref. [43] scattering of charmed and light pseudoscalar mesons was studied, including DK scattering, and by using SU(3) flavor symmetry the D * s0 (2317) mass was obtained in agreement with experiment and support for the interpretation as a DK molecule was found. In Ref. [44] scattering in the D * 0 sector was studied, however, with rather heavy quark masses, somewhere between the physical light and strange quark, corresponding to a pion mass m π ≈ 0.391 GeV, which led to some insights on the qualitative difference between the D * 0 (2400) meson and the D * s0 (2317) meson.
In this work we also study the D * s0 (2317) meson using lattice QCD with particular focus on tetraquark interpolating fields. As in the aforementioned lattice QCD investigations [40][41][42] we consider both quark-antiquark and two-meson interpolating fields. In addition, we include for the first time also tetraquark interpolating fields, where the four quarks are centered at the same point in space. We implemented color and spin contractions, where two standard meson interpolating fields of quarkantiquark type are put on top of each other (resembling DK and D s η mesonic molecules), as well as contractions, which have a diquark-antidiquark structure. Including such tetraquark interpolating fields might be essential, as it has recently been reported in Ref. [45] for the positive parity mesons a 0 (980) and K * 0 (700). In both cases a low-lying energy level is missed, if they are not taken into account. Moreover, we compute the couplings of the low lying states to different types of two-quark and fourquark interpolating fields and compare the spectra obtained from different subsets of interpolating fields. This might shed additional light on the question, whether the D * s0 (2317) meson is predominantly acs state or rather has a large tetraquark component.
We perform our computations in a single spatial volume of extent 2.90 fm and at unphysically heavy u and d quark mass corresponding to m π ≈ 296 MeV. For the analysis of correlation functions we apply the Athens Model Independent 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 previously been used in a study of the nucleon spectrum and the a 0 (980) meson [46,47]. 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 generalized eigenvalue problems and extract the spectrum from effective energy plateaus (cf. e.g. [48] and references therein). Note that both the GEVP and AMIAS provide information on the relative importance of the considered interpolating fields. Combining both methods allows to check the robustness of our results. This paper is organized as follows: In section II we describe the lattice setup and techniques with particular focus on the implemented interpolating fields. In section III we discuss the spectral decomposition of the corresponding two-point correlation functions. A short description of our two analysis methods, the GEVP and AMIAS, is provided in section IV. Section V is the main section of this work, where our numerical results are presented. First, in section V A, we show several finite volume spectra for the sector with D * s0 (2317) quantum numbers corresponding to different sets of interpolating fields. Based on these results the importance of each interpolating field is discussed. Then, in section V B, we perform a scattering analysis using Lüscher's method and the effective range expansion to determine the D * s0 (2317) mass at infinite volume. In section VI we summarize our findings and give our conclusions.
C denotes the charge conjugation matrix and the normalization factors N j are chosen such that C jj (t = a) = 1 (no sum over j; a is the lattice spacing), i.e. in a way that the interpolating fields generate trial states with similar norm. All interpolating fields couple to the D * s0 (2317) meson and to other states with the same quantum numbers. As in previous lattice QCD computations [40][41][42] we consider quark-antiquark interpolating fields, O qq, 1 and O qq, γ0 , as well as two-meson interpolating fields, O DK, 2part and O Dsη, 2part . In Refs. [40][41][42] it was shown that the latter interpolating fields are essential to determine the energy of the ground state and the first excitation reliably. were not considered in previous lattice QCD studies of the D * s0 (2317) meson. Thus, the main goal of this work is to explore, whether the inclusion of these tetraquark interpolating fields has an effect on the lattice QCD determination of the low-lying spectrum. Similar recent investigations of systems not including the D * s0 (2317) meson have led to different findings regarding the importance of tetraquark interpolating fields. While in Refs. [47,49] only marginal differences in the resulting spectra of the I = 1 hidden-charm and doubly-charmed sectors and the a 0 (980) sector were found, Ref. [45] observed additional energy levels both with K * 0 (700) and a 0 (980) quantum numbers. In this work, we also compute and compare the overlaps of the corresponding trial states O j |Ω (|Ω denotes the vacuum) to the lowest energy eigenstate, to obtain certain information about the quark composition of the D * s0 (2317) meson. This might contribute to the ongoing debate, whether the D * s0 (2317) meson is predominantly a quark-antiquark pair or a tetraquark (see the discussion in section I).
Note that the interpolating fields O 3 to O 7 do not generate orthogonal trial states. For example the terms with x = y in Eqs. (7) and (8) also appear in Eqs. (4) and (5). Similarly, one can relate two-meson combinations to diquark-antidiquark combinations via a Fierz identity, i.e. some of the terms present in Eqs. (4) and (5) are also part of Eq. (6) and vice versa. Even though the seven interpolating fields do not generate orthogonal trial states, they are not linearly dependent either, because each of them contains terms not present in any of the other six. Their non-orthogonality does not cause any particular problems during our analyses, because the two methods we use, the GEVP and AMIAS, are both able to deal with correlation matrices based on non-orthogonal trial states. We remark that on a technical level this work is similar to our lattice QCD investigation of the a 0 (980) meson [47], because to a large extent the same interpolating fields are used, just with different quark flavors.
To compute the correlation functions, 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 by the PACS-CS Collaboration [53]. The lattice size is 64 × 32 3 with lattice spacing a = 0.0907 (14) fm, i.e. the spatial lattice extent L is around 2.90 fm. The u and d quark mass and the s quark mass correspond to the pion mass m π ≈ 0.296 GeV and the kaon mass m K ≈ 0.597 GeV, i.e. are both heavier than in the real world, while the c quark mass corresponds to the D meson mass m D ≈ 1.845 GeV, i.e. is slightly lighter (see the detailed discussion in section V B). Note that the c quark only appears as a valence quark.
In a recent publication [54], we implemented and compared various combinations of techniques for the computation of propagators and correlation functions including point-to-all propagators, stochastic timeslice-to-all propagators, the one-end trick and sequential propagators. For each diagram of a similar 6 × 6 correlation matrix, which we used to study the a 0 (980) meson [47], we determined the most efficient combination of techniques. We have applied the same combinations of techniques in this work to compute the 7 × 7 correlation matrix (1) with the interpolating fields (2) to (8). 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 Ref. [54].

III. CORRELATION FUNCTIONS FOR PERIODIC TEMPORAL DIRECTION
A correlation function computed on a lattice with periodic temporal direction of extension T can be expanded according to with energy eigenstates |m , corresponding energy eigenvalues E m , possibly complex c j m,n = m|O j |n and Z = m e −EmT . Using the QCD symmetries charge conjugation and time reversal one can show that all elements of the correlation matrix (1) with interpolating fields (2) to (8) are real. Moreover, one can rewrite Eq. (9) in more convenient form, with real c j m,n and H jk (x) = − sinh(x) for j = 2, k = 2 and j = 2, k = 2 + cosh(x) otherwise .
Since the elements of the correlation matrix are either symmetric with respect to the reversal of time, , it is sufficient to restrict the following discussion to temporal separations 0 ≤ t ≤ T /2. For sufficiently large T , where Z ≈ e −EΩT (Ω denotes the vacuum), and for sufficiently large t, Eq. (10) reduces to if the correlation function is not contaminated by effects related to multi-hadron states as discussed below.
denotes the sum over a finite number of low-lying energy eigenstates in the sector with D * s0 (2317) quantum numbers, which is probed by the interpolating fields (2) to (8) (in the following we assume the ordering For temporal separations t around T /2 the correlation functions C jk (t) have a more complicated expansion than Eq. (12), if there are low-lying multi-hadron states with the same quantum numbers. In our case, i.e. for D * s0 (2317) quantum numbers, the lowest multi-hadron state is a DK scattering state, which has an energy only slightly above the mass of the D * s0 (2317) meson. Clearly, the interpolating fields (2) to (8) do not only excite such a DK scattering state, when applied to the vacuum |Ω , but also yield non-vanishing matrix elements D|O j |K and K|O j |D , i.e. annihilate a kaon and create a D meson and vice versa. For example a significant contribution to C jj (t) is as can be seen from Eq. (10). Assuming coefficients |c j D,K | ≈ |c j DK,Ω |, where DK denotes a low-lying DK scattering state, one can see that in the region of t ≈ T /2 the corresponding terms in Eq. (10) are comparable in magnitude. Therefore, terms as in Eq. (13) have to be taken into account, when extracting energy levels from correlation functions at large temporal separations t ≈ T /2. For smaller temporal separations t such contributions may be neglected, since they are exponentially suppressed ∝ e −2m K (t−T /2) with decreasing t. Analytical estimates as well as numerical experiments have shown, that within our setup this is the case for t < ∼ 15 a, which is an upper bound for all t fitting ranges used in the following.

IV. ANALYSIS METHODS
To analyze the 7 × 7 correlation matrix discussed in section III and various submatrices, we use both the GEVP method and the AMIAS method. While the GEVP method is quite common and very well known, the AMIAS method has proven to be particularly suited to study excited states [46] and was succesfully used in our related previous lattice QCD study of the a 0 (980) meson [47]. In section V we will show that both methods yield consistent results, which we consider to be an important cross-check, in particular due to the fact that the signalto-noise ratios of the elements of the correlation matrix grow rapidly with increasing temporal separations. In the following we summarize both methods and discuss the details of our analyses.

A. GEVP method
A commonly used method to extract several energy levels from an N × N correlation matrix is to solve the generalized eigenvalue problem (see e.g. Ref. [48] 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 ) (m = 0, . . . , N − 1) and t 0 ≥ a an input parameter. We use t 0 = a, which is a typical choice. A number of N effective energies E eff,m (t) can be obtained by solving for each eigenvalue λ m (t, t 0 ). At sufficiently large, but not too large temporal separations, i.e. in a t-region, where Eq. (12) is a valid parameterization of the correlation matrix, the effective energies E eff,m (t) exhibit plateaus. The values of these plateaus correspond to the N lowest energy levels in the sector probed by the interpolating fields, i.e. to E m . We determine each energy level E m by first fitting to the eigenvalue λ m (t, t 0 ) in the region t min ≤ t ≤ t max , where A 0 , A 1 and E 0 < E 1 are fitting parameters. t min and t max are chosen as follows: • t min is the smallest temporal separation t, where and ∆E eff,m (t) is the statistical error of E eff,m (t)).
• t max is the largest temporal separation t, where as well as This definition of t min and t max guarantees that the effective energy is consistent with a plateau within statistical errors for t ≥ t min and that its statistical errors are still reasonably small at t = t max . The energy level E m is then determined by averaging E f eff,m (t) over the fitting region, (see also Ref. [55], where a similar procedure was used).
The components of the eigenvectors v m (t, t 0 ) obtained by solving the GEVP (14) provide information about the structure of the corresponding energy eigenstates: for sufficiently large t, where the ≈ sign denotes the expansion of the energy eigenstate |m within the subspace spanned by the trial states O j † |Ω . We found that for t ≥ t min the eigenvector components are constant within statistical errors. Thus, we average the eigenvector com- and normalize via v j m → v j m /|v m |.

B. AMIAS method
In practice, 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 energy levels E m difficult. Therefore, in addition to the GEVP method we employ an alternative analysis method called AMIAS [46,56,57].
In section III we have discussed that lattice QCD results for correlation functions C jk (t) (see Eq. (1)) with interpolationg fields O j , j = 1, . . . , 7 (see Eqs. (2) to (8)) can be parameterized according to Eq. (12). In the t range we are going to consider, a ≤ t ≤ 15 a, and for the energy levels E m expected, cosh and sinh can be approximated by exponential functions, resulting in fit functions The fit parameters E m and c j m,Ω are real. In the following they 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 in uncorrelated χ 2 minimizing fits. C jk (t) denotes the correlation functions computed using lattice QCD with corresponding statistical errors ∆C jk (t), while C fit j,k (t) is given by Eq. (24). In principle one can also use a correlated χ 2 . Then, however, one has to estimate a covariance matrix, which requires rather precise data and computations on a large number of gauge link configurations (cf. e.g. Ref. [58] for a detailed discussion).
To obtain the PDF Π(A r ) for a specific fit parameter A r , one has to integrate Eq. (27) 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 use a parallel tempering scheme combined with the Metropolis algorithm as described in detail in Ref. [46]. The parallel tempering scheme prevents the algorithm from getting stuck in a region around a local minimum of χ 2 and guarantees ergodicity of the algorithm.
While we use t max = 15 a in Eq. (28), we vary in our analyses both t min and the number of terms in the truncated sum in Eq. (24), until we find a stable region with no observable change in the PDFs for the low-lying energy eigenstates of interest. For a detailed example see Ref. [46].
The coefficients c j Ω,m = Ω|O j |m = m|O j † |Ω in the fit function (24) 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 inverting Eq. (30) 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. jṽ j m c j Ω,n ≈ δ m,n .
Note that the coefficientsṽ j m are equivalent to the eigenvector components v j m obtained by solving a GEVP (see Eq. (23)) and, thus, the resuls from the two methods can be compared in a meaningful way, after choosing the same normalization (ṽ m ) 2 = 1. Our goal in this section is to determine the two lowest energy levels in the sector with D * s0 (2317) quantum numbers in the finite spatial volume L 3 of the lattice. From previous results [40][41][42] we expect that one of the corresponding energy eigenstates is the lowest DK scattering state, while the other has a somewhat smaller energy and represents the D * s0 (2317) meson. A precise determination of these two energy levels is necessary to study the infinite volume limit using Lüscher's finite volume method in section V B.
Moreover, in this section we will also investigate the quark content and arrangement of the D * s0 (2317) state by studying the eigenvector components v j m and the coefficientsṽ j m introduced in section IV.

D and K meson masses and DK threshold
As a preparatory step we computed the masses of the pseudoscalar mesons D and K within our lattice setup. It is rather straightforward to obtain precise values for m K and m D from correlation functions of standard interpolating fields We find m D = 1.8445(9) GeV (35) m K = 0.5965(4) GeV (36) using the AMIAS analysis method (extracting the masses from the corresponding effective energies, as explained in the context of the GEVP method in section IV A, leads to compatible results, however, with somewhat larger statistical errors) [59]. Consequently, which is the lowest two-meson threshold in the sector with D * s0 (2317) quantum numbers and, thus, plays an important role in the interpretation of further results. Also of interest is the energy of a non-interacting DK pair with one quantum of relative momentum p min = 2π/L,

Preselection of interpolating fields
To reduce the seven interpolating fields (2) to (8) to a somewhat smaller set of interpolating fields, which are most important to resolve the lowest energy eigenstates with D * s0 (2317) quantum numbers, we performed GEVP as well as AMIAS analyses using individual correlation functions or 2 × 2 correlation matrices. From these analyses it became clear, that three of the interpolating fields (2) to (8) are less relevant.
• O qq, γ0 (Eq. (3)): One can determine the energy of a low lying state, which we will identify below as the D * s0 (2317) meson, using the correlation function of one of the two quark-antiquark interpolating fields, i.e. either of O qq, 1 or of O qq, γ0 . For the latter, however, the effective energy plateau is reached at larger temporal separation and statistical errors are larger as well. An analysis of the corresponding 2 × 2 correlation matrices gives the same low lying state and a second rather noisy effective energy significantly above, which most likely receives contributions from several excited states. The eigenvector components v j m indicate a strong dominance of the interpolating field O qq, 1 for the ground state. In view of these findings we consider O qq, 1 superior to O qq, γ0 and do not use the latter interpolating field in any of the following analyses.
• O Dsη, point and O Dsη, 2part (Eqs. (5) and (8)): Correlation functions containing either O Dsη, point or O Dsη, 2part exhibit large statistical errors. This seems to be a consequence of the "η interpolator" uγ 5 u +dγ 5 d, which is part of O Dsη, point as well as of O Dsη, 2part (note that a lattice QCD study of the η meson is quite challenging by itself, partly because of strong statistical fluctuations; see e.g. Refs. [60,61] for a detailed discussion and a recent computation). Moreover, the mass of the D * s0 (2317) meson is close to the DK threshold, while the D s η threshold is around 155 MeV above [4]. This suggests that including the interpolating fields O Dsη, point and O Dsη, 2part is not essential to resolve the two lowest energy eigenstates. This is supported by Refs. [14][15][16][17], where the molecular components for the D * s0 (2317) were found to be from around 60% to 75% for DK and below 15% for D s η. Thus, we do not use O Dsη, point and O Dsη, 2part in any of the following analyses.
and (7)) are used in the following to determine the two lowest energy levels in the sector with D * s0 (2317) quantum numbers.

GEVP analysis
The results of a GEVP analysis of the 4 × 4 correlation matrix containing the four interpolating fields identified in the previous subsection are collected in FIG. 1. The upper plot shows effective energies as functions of the temporal separation. The four plots below contain the squared eigenvector components (v j m ) 2 . There are two convincing effective energy plateaus with small statistical errors close to the DK threshold. One is around 60 MeV below, while the other is somewhat above, but significantly closer to the DK threshold than to the energy of a non-interacting DK pair with one quantum of relative momentum as indicated by the horizontal gray lines (see also Eqs. (37) and (38)). Thus, there is an additional low-lying state compared to the non-interacting DK spectrum. The eigenvector components clearly indicate that the lowest state resembles a quark-antiquark pair ((v qq, 1 0 ) 2 > ∼ 0.90; black bar chart), while the first excitation is a DK scattering state similar to a non-interacting two-meson state with both mesons at rest ((v DK, 2part 1 ) 2 ≈ 0.85; red bar chart). These eigenvector components suggest to identify the lowest state as the D * s0 (2317) meson. Energy levels E m are determined from effective energies as discussed in detail in section IV A. For m = 0, 1 they are listed in Table I.
The third effective energy E eff,2 (t) has large statistical errors and is somewhat above the estimated energy of a non-interacting DK pair with one quantum of relative momentum. It seems likely that it corresponds to a linear superposition of several DK scattering states with  non-vanishing relative momenta. This interpretation is supported by the eigenvector components, which indicate a dominance of the O DK, point interpolating field ((v DK, point 2 ) 2 ≈ 0.90; green bar chart), which by construction excites DK states with many different relative momenta. The fourth effective energy E eff,2 (t) has even larger statistical errors and is around 1 GeV above the DK threshold. Most likely it represents a superposition of a larger number of highly excited states.
From (v QQ, γ5 m ) 2 < ∼ 0.05 for m = 0, 1, 2 one can conclude that the diquark-antidiquark interpolating field O QQ, γ5 is not important to resolve any of the three lowest energy eigenstates. In particular the ground state seems to be predominantly a standard quark-antiquark pair ((v qq, terpret this as indication that the D * s0 (2317) meson has no sizable tetraquark component.

AMIAS analysis
The results of an AMIAS analysis of the 4 × 4 correlation matrix are collected in FIG. 2. The upper plot shows the PDFs generated with the fit function given in Eq. (24) and six terms in the truncated sum [62]. The four plots below show the squared coefficients (ṽ j m ) 2 . When comparing the PDFs to the effective energies in FIG. 1 one can see, that the energy levels obtained with AMIAS are consistent with those from the GEVP analysis. Statistical errors for the AMIAS results are somewhat smaller than for the GEVP results (see Table I). The coefficients (ṽ j m ) 2 are also in reasonable agreement with the GEVP eigenvector components (v j m ) 2 from FIG. 1, supporting that the ground state is mostly a quark antiquark-pair.
To cross-check the obtained results, in particular to confirm our findings regarding the quark composition and interpretation of the low-lying energy eigenstates, it is useful to compare the above 4 × 4 AMIAS analysis to analogous analyses using the the four possible 3 × 3 submatrices as input (for the latter five terms in the truncated sum in Eq. There is essentially no difference between the 3 × 3 and 4 × 4 PDFs for the three lowest energy levels. This confirms that the diquark-antidiquark interpolating field O QQ, γ5 is not important to resolve the low-lying energy eigenstates. This in turn supports our conclusion from section V A 3 that the D * s0 (2317) meson does not have a sizable tetraquark component. level of the second excitation is, however, significantly larger. This indicates that the interpolating field O DK, point is useful to resolve higher momentum excitations, while it is not essential for a determinaton of the lowest two energy levels.
The lowest two energy levels are slightly larger compared to the 4 × 4 result, but they are still compatible, because of their drastically larger statistical errors (see Table I). Thus, it is possible to excite the D * s0 (2317) meson with only four-quark interpolating fields, i.e. it seems to have a non-vanishing, but small DK component, which is in agreement with our findings from section V A 3. The lowest energy level is consistent with the 4 × 4 result, even though it has a much larger statistical error. The energy level of the first excitation, however, cannot be determined reliably anymore. The PDF has a large width and its peak is localized at energies significantly above the DK threshold. This confirms that the interpolating field O DK, 2part is of central importance for a determination the energy of the lowest DK scattering state.

Summary of finite volume results and conclusions
A summary plot of the obtained energy levels with the 4 × 4 GEVP as well as the 4 × 4 and 3 × 3 AMIAS analyses is shown in FIG. 4. Again it can be seen that the most important interpolating fields to determine the two lowest energy levels are O qq, 1 and O DK, 2part . Analyses using these two interpolating fields (4 × 4 GEVP, 4 × 4 AMIAS, 3 × 3 AMIAS (A) and (B)) yield consistent energy levels with small statistical errors. From the GEVP eigenvector components v j m and the AMIAS coefficientsṽ j m we conclude that the lowest energy level mostly corresponds to a quark-antiquark bound state, possibly similar to the D * s0 (2317) meson in the infinite volume. There seems to be a small DK component, but no sign of any sizable tetraquark component. Moreover, the components v j m and the coefficients v j m clearly indicate that the first excitation is a DK scattering state. These two energy levels will be used for the finite volume analysis in section V B.

B. Scattering analysis and infinite volume limit
The energy levels collected in Table I were computed at finite spatial volume L 3 with periodic boundary conditions. One can determine the mass of the D * s0 (2317) meson, which is the infinite volume limit of the ground state energy, from the lowest two energy levels E 0 and E 1 at finite volume by performing a scattering analysis continued to imaginary momenta, i.e. using Lüscher's finite volume method [63]. This approach has been used in a lattice QCD study of the D * s0 (2317) meson for the first time in Refs. [40,41] and later also in Ref. [42]. It was also used to study other systems (see e.g. Refs. [44,64]).
For a recent review on scattering in lattice QCD see Ref. [65].
The first step is to determine the squared scattering momenta k 2 0 and k 2 1 via where E 0 and E 1 can be taken from Table I and m D and m K are the D meson and K meson masses obtained within the same lattice setup (see Eqs. (35) and (36)). With Lüschers finite volume method one can then compute the corresponding two phase shifts δ 0 (k 0 ) and δ 0 (k 1 ), Here Z 00 denotes the generalized zeta function and L ≈ 2.90 fm the spatial lattice extent (see section II). k cot(δ 0 (k)) can be written as a Taylor series in k 2 , where a 0 is the S wave scattering length and r 0 the S wave effective range. For sufficiently small k 2 one can neglect terms of order k 4 and parameterize k cot(δ 0 (k)) by the first two terms on the right hand side of Eq. (41). a 0 and r 0 are then fixed by the two data points cot(δ 0 (k 0 )) and cot(δ 0 (k 1 )) obtained via Eq. (40). This parameterization is called effective range expansion. A stable D * s0 (2317) meson manifests itself as a pole in the scattering amplitude, i.e. corresponds to cot(δ 0 (k D * s0 )) = i, where k D * s0 denotes the position of the pole. Combining this condition with the parameterization (41) leads to which can easily be solved with respect to k 2 (note that for our data one of the two solutions has to be discarded, because it is far outside the region of validity of the effective range expansion (41), where O(k 4 ) terms cannot be neglected). The mass of the D * s0 (2317) meson is given by the right hand side of Eq.
In Table II we show the results obtained for the lowest two energy levels E 0 and E 1 , for the squared scattering momenta k 2 0 and k 2 1 , for the phase shifts, for the S wave scattering length a 0 and effective range r 0 as well as for the position of the pole. To verify that the effective range expansion (41) is a reasonable approximation, we also provide s(k 2 0 ), s(k 2 1 ) and s(k 2 D * s0 ), where s(k 2 ) = |a 0 r 0 k 2 /2| corresponds to the ratio of the O(k 2 ) term and the O(k 0 ) term in Eq. (41). We find values 1 for the two scattering momenta as well as for the position of the pole, which gives certain indication that higher order terms are suppressed, i.e. that O(k 4 ) terms in Eq. (41) are indeed negligible. In Table II we also list m D * s0 , the resulting mass of the D * s0 meson, and m D + m K − m D * s0 , the binding energy with respect to the DK threshold. All results are provided both for the 4 × 4 GEVP and the 4 × 4 AMIAS determination of the lowest two energy levels E 0 and E 1 discussed in sections V A 3 and V A 4.
In FIG. 5 we show the parameterization of k cot(δ 0 (k)) with the effective range expansion (41) together with the two data points k 0 cot(δ 0 (k 0 )) and k 1 cot(δ 0 (k 1 )). Note that the effective range expansion is equivalent to the right hand side of Eq. (43). We also show the left hand side of that equation, ik = − √ −k 2 . The intersection of the two curves corresponds to k 2 D * s0 , the binding momentum of the D * s0 (2317) meson. In FIG. 6 we illustrate the pole in the scattering amplitude by plotting in the complex k plane, i.e. Eq. (42) with the parameterization (41) inserted. The color reflects the quality of the effective range expansion (41) and indicates that the pole is in a region, where O(k 4 ) terms should be negligible. It is important to note that a direct comparison of our result for m D * s0 to the corresponding experimental result m D * s0,exp = 2317.8(5) MeV [4] is not meaningful, because the quark masses in our simulation differ from their experimental counterparts: • The light u and d quark mass is unphysically heavy, reflected by the pion mass m π ≈ 0.296 GeV.
• The c quark is unphysically light, because m D ≈ 1.845 GeV, i.e. below m D,exp ≈ 1.867 GeV.
Also for the binding energy of the D * s0 (2317) meson with respect to the DK threshold, m D + m K − m D * s0 , it is not clear a priori, whether quark masses, which differ from their physical values, will result in a value similar to the corresponding experimental value m D,exp + m K,exp − m * Ds0,exp ≈ 45 MeV. One reason for this is that the threshold m D + m K clearly depends on the light quark mass, while m D * s0 , which according to section V A E0/GeV E1/GeV k 2 0 /GeV 2 k 2 1 /GeV 2 k0 cot(δ0(k0))/GeV k1 cot(δ0(k1))/GeV GEVP, 4   is mostly acs state, should be almost independent of the light quark mass (see also the discussion in Ref. [41]). Note, however, that we find m D + m K − m D * s0 ≈ 51 MeV rather close to the experimentally observed 45 MeV, which indicates that with respect to the D * s0 meson we The pole corresponds to the D * s0 (2317) meson. The color of the plotted surface is related to the value of s(k 2 ) (green: s(k 2 ) < 0.1; yellow: 0.1 ≤ s(k 2 ) < 0.2; red: 0.2 ≤ s(k 2 )). Thus, it reflects the quality of the effective range expansion (41) and indicates that the pole at k D * s0 is in a region, where O(k 4 ) terms should be negligible. might be in a similar situation as in real world QCD, even though we are not precisely at physical quark masses. Because of this and since m D * s0 is close to the lowest energy level E 0 obtained at finite lattice volume (around 10 MeV difference as can be seen from Table II), we expect that our findings and statements from section V about the importance of the two-quark and the four quark interpolating fields also apply for physical quark masses and the infinite volume limit. This is further supported by the qualitative agreement of our results for a 0 and r 0 and corresponding results obtained in lattice QCD computations at almost physical quark masses [41,42].

VI. SUMMARY AND CONCLUSIONS
We studied the D * s0 (2317) meson with lattice QCD using interpolating fields of different structure. In addition to quark-antiquark interpolating fields and two-meson interpolating fields, which were already considered in previous lattice QCD studies, we implemented and explored the importance of tetraquark interpolating fields. For these tetraquark interpolating fields the four quark operators are centered at the same point in space and their color and spin structure corresponds to either a mesonmeson pair or a diquark-antidiquark pair.
In the finite spatial volume of our lattice with extent L ≈ 2.90 fm we find two low-lying energy eigenstates, one around 60 MeV below the DK threshold, the other slightly above the DK threshold. The GEVP eigenvector components and the AMIAS coefficients and PDFs clearly indicate that the state below threshold, which corresponds to the D * s0 (2317) meson, is mostly of quark-antiquark type with only a small DK component, while the state above threshold is a DK scattering state. The tetraquark interpolating fields explored in this work turned out to be essentially irrelevant, when extracting the corresponding two energy levels, i.e. the couplings of the state below threshold to these interpolating fields is close to zero. We interpret this as indication that the D * s0 (2317) meson is mainly a quark-antiquark state and not a tetraquark, as discussed or proposed by various existing papers.
It is important to keep in mind that our computation was carried out for a single spatial volume and at quark masses different from those in the real world, in particular a u and d quark mass corresponding to a heavier pion, m π ≈ 0.296 GeV. We performed a scattering anal-ysis using Lüscher's method to determine the mass of the D * s0 (2317) in the infinite volume limit. We find this mass 51 MeV below the DK threshold, rather close to our finite volume result as well as to the experimental value 45 MeV. Thus we expect that our findings concerning the importance of various interpolating fields as well as the quark composition of the D * s0 (2317) meson will also apply to infinite volume and physical quark masses at least on a qualitative level. Of course, it would be worthwhile and interesting to perform similar computations at physical quark masses and for several volumes in the future, in particular to check the approximate independence of the GEVP eigenvector components or the AMIAS coefficients from the spatial volume. 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.
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 [66] with a multigrid solver [67].