Bound states on the lattice with partially twisted boundary conditions

We propose a method to study the nature of exotic hadrons by determining the wave function renormalization constant $Z$ from lattice simulations. It is shown that, instead of studying the volume-dependence of the spectrum, one may investigate the dependence of the spectrum on the twisting angle, imposing twisted boundary conditions on the fermion fields on the lattice. In certain cases, e.g., the case of the $DK$ bound state which is addressed in detail, it is demonstrated that the partial twisting is equivalent to the full twisting up to exponentially small corrections.


I. INTRODUCTION
The search for the exotic states (tetraquarks, hybrids, hadronic molecules, etc) in the observed hadron spectrum has been a subject of both theoretical and experimental investigations for decades.
The exact pattern, how these states emerge, should be strictly determined by the underlying theory and should therefore contain important information about the behavior of QCD at low energies. In practice, however, extracting such information from the data encounters certain challenges, which are in part of a conceptual nature. In the present paper we wish to focus exactly on this issue.
In general, a state is called "exotic" if its quark content does not correspond to the "standard" constellation given by the non-relativistic quark model (qq for mesons and qqq for baryons). Consequently, one needs to use a particular model as a reference point to define how the exotic states are meant (note that the very notion of constituent quarks is, strictly speaking, model-dependent). Putting it differently, one has to agree on certain criteria formulated in terms of certain hadronic observables: if these observables are measured, or calculated on the lattice, and the results do not follow the pattern predicted by the quark model, this then should be interpreted as a signature for exotica.
A standard example for the exotic state candidates is given by the scalar nonet with the masses around 1 GeV. As it is well known, the observed mass hierarchy in this nonet is reversed as compared to, e.g., the pseudoscalar or vector multiplets. Such a mass ordering is counter-intuitive from the point of view of the naive quark model, but can be easily understood, if the scalar mesons were interpreted as tetraquark states (see, e.g., [1][2][3][4]). This is, however, not the only possible interpretation. In Refs. [5][6][7], the a 0 (980) and f 0 (980) were considered as hadronic molecules, whereas in Refs. [8] these states were described as a combination of a bare pole and the rescattering contribution. In the Jülich mesonexchange model, the f 0 (980) appears to be a bound KK state, whereas the a 0 (980) is a dynamically generated threshold effect [9]. Similar conclusions were inferred in Ref. [10] from the calculations in the unitarized ChPT with explicit resonance states. Finally, the investigations carried out within the framework of QCD sum rules are also indicative of the non-qq nature of a 0 (980) [11]. Given these multiple interpretations, it is natural to look for the clear-cut criteria based on the observables in order to minimize the model-dependence of the statements about the nature of the hadronic states in question.
In fact, such criteria are known for quite some time already. The "pole counting" method, considered in Refs. [12,13], relates the number of the S-matrix poles near threshold to the molecular nature of the states corresponding to these poles. Namely, it has been argued that the loosely bound states of hadrons (hadronic molecules) correspond to a single pole, whereas the poles corresponding to the tightly bound quark states (of standard or exotic nature) always come in pairs. A closely related criterion goes under the name of Weinberg's compositeness condition [14], which uses the quantity called the wave function renormalization constant Z, where 0 ≤ Z ≤ 1, to differentiate between the loosely bound states and tight QCD composites, the values Z ≃ 0 corresponding to the molecular states and vice versa. The application of these methods for the analysis of the data on scalar mesons are considered in Refs. [15][16][17][18][19][20], and the recent review on the subject may be found in Ref. [21]. Moreover, theoretically, one may study the dependence of the pole positions on the number of the colors N c (see Refs. [10,22,23]) or the quark masses (Refs. [24][25][26][27]). From the above studies, one can judge about the precise structure of these states beyond the simple alternative between a molecule and a tight quark composite.
Recent years have seen a renewed interest in the field, which is partly related to the progress in the lattice calculations of the QCD spectrum at the quark masses close to the physical values. It should be realized that the lattice studies have powerful tools at their disposal to analyze the nature of the states that emerge in QCD. Apart from the information about the dependence of the spectrum on quark masses, a valuable information comes from the volume dependence of the calculated spectrum as well as its dependence on the twisting angle in case of twisted boundary conditions, see Refs. [28][29][30][31][32]. Note that all this information is obtained from the first-principle calculations on the lattice and is thus in principle devoid of any model-dependent input.
In this paper we investigate the nature of the scalar states in the sector with one charm quark that is a natural generalization of our treatment of the light scalar mesons. We mainly focus on the case of the D * s0 (2317) meson [33], albeit the formalism, which we develop here, can be straightforwardly applied to the other cases where a bound state close to the elastic threshold emerges (note that, in this paper, we do not consider the generalization of the approach to the inelastic case. This forms a subject of a separate investigation.). The D * s0 (2317) does not fit very nicely to the quark-model picture, and its structure is still debated, see, e.g., Ref. [34] for a recent review. The molecular picture, due to the closeness of the DK threshold and a large coupling to the DK channel looks most promising among other alternatives. It would be highly desirable to verify this conjecture in a model-independent manner, on the basis of the lattice calculations. To this end, one may use the fact that the dependence of the bound-state energy on the kaon mass is very different for a molecule and a standard quarkmodel state, see Ref. [35]. Another possible method to address this issue has been described, e.g., in Refs. [29,30], where the authors propose to study the volume-dependence of the spectrum in order to apply the Weinberg's compositeness criterion on the lattice.
The exploratory study of light pseudoscalar mesons (π, K) of (D, D s ) in full lattice QCD has been carried out in Refs. [36][37][38]. In some isospin channels the study is plagued by the presence of disconnected contributions. The implementation of the method from Refs. [29,30], which implies carrying out calculations at different volumes, could be therefore quite expensive. In this paper we propose an alternative, which requires calculations at one volume, albeit with twisted boundary conditions. Moreover, we show that, in the study of D * s0 (2317), one may use partially twisted boundary conditions, despite the fact that the quark annihilation diagrams are present. The method used in the proof is the same as in Ref. [39]. Generally, one may expect that the simulations with partially twisted boundary conditions could be less expensive than working at different volumes, while they provide us the same information about the nature of the bound states in question.
This article is organized as follows. In Sect. II, we briefly review Weinberg's argument for the compositeness of particles. In Sect. III we describe the procedure of extraction of the parameter Z from the data with twisted boundary condition. Further, in Sect. IV we use some models and produce synthetic lattice data in order to check the procedure of the extraction in practice. The error analysis has also been carried out. Separately, in Sect. V, we discuss the use of the partially twisted boundary conditions and show that they are equivalent to the full twisting in our case. Sect. VI contains our conclusions.

II. COMPOSITENESS OF BOUND STATES
As mentioned before, in view of the plethora of candidates of exotic hadrons, it is very important to make model-independent statements on the nature of these states. Model-independence requires that we can only study the physical observables which can be defined in terms of the matrix elements between asymptotic states. In particular, we would like to ask a question, whether a given particle, corresponding to the S-matrix pole, can be regarded as "elementary" or rather as a bound state (molecule) of other hadrons. The central place in this identification belongs to the so-called wave function renormalization constant Z, which has been used to distinguish composite particles from elementary ones since the early 1960's [14,[40][41][42][43][44]. To see its role, we will first discuss a non-relativistic quantum mechanical system, following the discussion of Ref. [14].
In this section, we will restrict our discussion to the infinite volume. Let us consider a two-body system with a Hamiltonian H = H 0 + V , where H 0 is the free Hamiltonian, and V specifies the interaction. Both H and H 0 have a continuum spectrum. Let us assume that there is a bound state solution of the Schrödinger equation with a binding energy E B , and H 0 also has a discrete spectrum which are the bare elementary particles. For simplicity, we will assume that there is only one such state, denoted by |B 0 . In the Hilbert space spanned by the eigenstates of the free Hamiltonian, the completeness relation is thus given by where µ = m 1 m 2 /(m 1 + m 2 ) is the reduced mass. Thus, the probability for the physical state |B overlapping with the elementary state |B 0 which, by definition, equals to Z, is given by where Eq. (1) is used. The quantity 1 − B 0 |B 2 then describes the probability of the physical state not being the elementary state or finding the physical state in the two-particle state. In other words, Z ≃ 1 corresponds to a mostly elementary state whereas a state with Z ≃ 0 can be interpreted as a predominately molecular one.
In general, the above integral depends on the matrix element q |V |B , which is not directly measurable. However, for loosely bound states, the quantity Z can be related to the observables. Consider, for instance, an S-wave bound state with a small binding energy. The binding energy should be much smaller than the inverse of the range of forces so that the matrix element q |V |B can be approximated by a constant g NR . We get from Eq. (3) Note that, in the past, this equation has been often applied to distinguish composite particles from elementary ones, see e.g. [14,16,21,45]. The non-relativistic coupling constant g 2 NR coincides with the residue of the non-relativistic scattering matrix at the bound state pole. This can be immediately seen, considering the Low equation in the vicinity of the pole [14,44]. Here, E q = q 2 /(2µ).
Finally, we would like to relate the quantity Z to the physical observables, namely, to the scattering length a and effective range r. Here, we are closely following the path of Ref. [14]. It is important to note that these relations can be derived when the binding energy is much smaller than the inverse of the range of forces. We start with the twice-subtracted dispersion relation for the inverse of t(E) where the two subtraction constants have been determined from Eq. (5). The S-wave transition matrix element is related to the non-relativistic S-wave scattering amplitude f (k) = 1/[k cot δ(k) − i k] as f (k) = −µ t(E)/(2π) with k = √ 2µE and δ(k) being the S-wave phase shift. Thus, one gets Im t −1 (w) = µ √ 2µw/(2π). Inserting this into Eq. (6), we obtain where R = 1/ √ 2µE B denotes the characteristic distance between the constituents in the twobody bound system. Comparing the above expression with the effective range expansion t −1 (E) = −µ/2π −1/a + r k 2 /2 − i k , and using Eq. (4), one can express the scattering length and effective range in terms of the binding energy and compositeness [14] Therefore, for an S-wave shallow two-body bound state, the compositeness can be measured by measuring the low-energy scattering parameters.
Next, we turn to the compositeness condition within the framework of the quantum field theory. For simplicity, let us first consider the situation when a scalar particle described by a field Φ(x) with the bare mass M 0 couples with two scalars φ 1,2 (x) with the masses m 1,2 . The interaction Lagrangian takes the form L int = g 0 Φφ 1 φ 2 .
Consider now the two-point function of the field Φ(x) Summing up one-loop bubble diagrams to the two-point function, one arrives at the expression (see Fig. 1) where the one-loop self-energy is given by The relativistic scattering amplitude for the process φ 1 φ 2 → φ 1 φ 2 in the same approximation is given by (see Fig. 1 1 Here, in order to be consistent with the non-relativistic formalism, the sign convention S = 1−iT is used in the definition of the T -matrix. FIG. 1. The scattering matrix for the process φ 1 φ 2 → φ 1 φ 2 (a) and the two-point function of the field Φ (b).
The relativistic and the non-relativistic scattering matrices are the same up to an overall normalization.
In the rest frame of the bound system, the relation takes the form where w i (k) = m 2 i + k 2 . Now, let us consider the behavior of the scattering amplitude in the vicinity of the bound-state pole. The two-point function has the following behavior where M is the physical mass.
The residue of the propagator determines the wave function renormalization constant for the particle Φ: where g 2 = Z g 2 0 is the renormalized coupling constant, and G ′ (M 2 ) = d ds G(s) s=M 2 . In order to establish the relation of the quantity Z, defined by Eq. (15), with its non-relativistic counterpart, we perform the contour integration over q 0 of the loop integral in Eq. (11): where In the rest frame of the bound state, one has P = 0. Taking derivative with respect to s, and then taking the non-relativistic approximation which amounts to ω 1 ≃ m 1 + q 2 /(2m 1 ) and ω 2 ≃ m 2 + q 2 /(2m 2 ), we get where we have used E B = m 1 +m 2 −M . Taking into account the difference between relativistic and nonrelativistic normalizations, we finally arrive at the relation Comparing now this relation with Eq. (3), one immediately sees that the wave function renormalization constant Z is the same as its non-relativistic counterpart and thus the compositeness condition for an S-wave bound state can be written as One might treat the above argumentation with a grain of salt, since it is based on certain approximations.
Namely, the amplitude is given as a sum of one-loop diagrams only. It is, however, clear that the result is valid beyond this approximation, if bound states close to an elastic threshold are considered. The justification is provided by the statement that such bound states can be consistently described within a non-relativistic effective field theory, which is perturbatively matched to the underlying relativistic theory (see, e.g., Ref. [46] for a review on the subject). Such an effective theory is equivalent to the non-relativistic quantum mechanics (the number of particles is conserved) and hence the compositeness can be rigorously defined along the lines discussed above. Finally, we would like to mention that the quantity Z, which is defined in Eq. (15), is ultraviolet finite, since the quantity g is defined through the residue of the renormalized scattering amplitude.

III. COMPOSITENESS FROM LATTICE DATA
As stated above, the wave function renormalization constant, Z, gives an overlap of the physical state with the elementary state and hence could be used as a parameter that describes the compositeness of a given state. Lattice calculations provide a model-independent way to determine Z from the volume dependence of the spectrum [30,[47][48][49][50][51], or -as we propose in this paper -from the dependence on the twisting angle. In this section we set up a finite-volume formalism, which describes the dependence of the bound-state mass on the volume or twisting angle.

A. Finite volume formalism
We consider elastic scattering of particles with the masses m 1 and m 2 in the S-wave 2 . Then, generally, a unitary partial-wave amplitude in infinite volume is given by where is the relative momentum squared in the center of mass (c.m.) frame. Further, the function V −1 (s) ("the inverse potential") is a regular function in the vicinity of the threshold. The notation used here is reminiscent of that of unitarized Chiral Perturbation Theory, but Eq. (19) may in fact describe any elastic unitary amplitude, with the particular dynamics encoded in the function V (s). The loop function G(s) is given by Eqs. (11) and (16). This function contains a unitarity cut. Across this cut, we have Im G(s) = −k/(8π √ s). Other (distant) cuts that may be also present are included in V (s). The loop function G(s) is divergent and has to be renormalized. Here we do the renormalization with a subtraction constant. As it will be seen below, the extension to the finite volume is independent of any regulator.
When the particles are put in a finite box of size L, their momenta become discretized due to boundary conditions. So, the continuum spectrum, which gives rise to the cut in the infinite volume, becomes a discrete set of two-particle levels. In order to obtain the spectrum in a finite volume, one should replace the momentum integrals by the sums over the discretized momenta in the expression of the scattering amplitude. Then, the "finite volume scattering amplitude"T contains poles on the real axis that correspond to the discrete two-particle levels. It should be noted that the finite-volume effects in V (s) are exponentially suppressed (see, e.g., [52]), so the the finite volume scattering amplitude can be obtained just by changing the loop function by its finite volume counterpartG θ Here I( q ) denotes the integrand in Eq. (16), and q n the allowed momenta in a finite volume, whose value depends on the box size L and the boundary conditions used. For the periodic boundary conditions we have q n = 2π L n, n ∈ Z 3 . In case of twisted boundary conditions, the momenta also depend on the twisting angle θ according to q n = 2π L n + θ L , 0 ≤ θ i < 2π. Using the methods of Ref. [53], it can be shown that ∆G θ L can be related to the modified Lüscher function Z θ 00 , see Appendix A, wherek = kL/(2π) and the dots stand for terms that are exponentially suppressed with the volume size L [53].
In this paper, we are going to apply Lüscher formalism to study shallow bound states, where the finite-volume effects are exponentially suppressed. Since, for such states, the binding momentum κ is presumed to be much smaller than the lightest mass in the system, the exponentially suppressed corrections emerging, e.g., from the potential V (s) could be consistently neglected as compared to the corrections ∼ e −κL that arise from Z θ 00 (1,k 2 ). Note however that, if masses of the constituents increase for a fixed binding energy, then the magnitude of the binding momentum also increases and, for the bound states of heavy mesons, may become comparable to the pion mass. In this case, further study of the problem is necessary. A recent example of such a study (albeit in the light quark sector) is given in Ref. [54]. In the present paper this issue is not addressed.
Finally, note that the divergences arising at Λ → ∞ in Eq. (20) cancel between the sum and the integral, so we can safely send the cutoff to infinity. Thus, ∆G θ L does not depend on any regulator. In Appendix A we show in detail, how ∆G θ L could be calculated below threshold for different types of boundary conditions.
where ψ(k 2 ) is the analytic continuation of k cot δ(k) for arbitrary complex values of k 2 , which is needed since the bound state is located below threshold, k 2 B < 0. On the other hand, the discrete levels in a finite volume are obtained as the poles of the finite-volume scattering amplitudeT and, in particular, the bound state pole gets shifted to M L , with binding momentum k L ≡ iκ L , given bỹ Note that, below threshold, both T −1 and ∆G θ L are real, so the pole position is real. The discrete scattering levels above threshold are real as well (as they should be), since the imaginary part of ∆G θ L cancels exactly with that of T −1 .
Next, we relate the finite-volume pole position with the infinite-volume quantities as the bound state mass, M , and the coupling, g 2 (defined as the residue of the scattering amplitude at the pole s = M 2 ).
To this end, we expand ψ(k 2 L ) around the infinite-volume pole position, k B = iκ, where the prime denotes a derivative respect to k 2 . Then, evaluating the residue at M 2 in Eq. (19) we where the derivative dk 2 /ds is to be evaluated at s = M 2 . Finally, using Eqs. (23) and (24), we obtain for the pole position shift This equation gives the bound state pole position, κ L (or, equivalently, M L = m 2 1 − κ 2 L + m 2 2 − κ 2 L ) as a function of the infinite-volume parameters g 2 and κ. It is worth noting that, within the approximation (24), the position of the bound state pole in a finite volume depends only on these two parameters.
This approximation works remarkably well in all cases considered in this paper.
If the difference κ L − κ is small enough, Eq. (26) can be solved iteratively. For periodic boundary conditions, with the use of Eq. (A3), it can be shown that the lowest-order iterative solution reads which coincides with the result given in Refs. [30,48,50]. However, it will be shown below that, for shallow bound states, where κ is very small, one should take more than just the first term in the sum (A3). Moreover, in some cases, the iterations converge very slowly, if at all. Therefore, in our opinion, it is safer to consider solving Eq. (26) numerically, without further approximations, in order to obtain the finite volume pole position κ L . This is the way we proceed.
Using Eq. (26), it is possible to fit the infinite-volume parameters M and g 2 from the bound state levels κ L , obtained through lattice simulations at different L or θ. This, in turn, allows one to determine the compositeness parameter from Eq. (18). However, in actual lattice simulations, the measured energy levels have some uncertainty, and the number of different volumes or different twisting angles might be not very large. Therefore, it is important to know in advance, at which accuracy should be the lattice measurements carried out, in order to render the extraction of the parameter Z reliable. We address this question in some exactly solvable models with a given V (s), producing "synthetic lattice data," adding random errors and trying to extract back the infinite volume parameters M, g 2 and Z from data.

A. A toy model
The potential in this model is given by a "bare state pole", which depends on two parameters: a bare pole position s 0 and a bare coupling constant g 0 . By appropriately choosing the value of the bare parameters, we can reproduce a bound state with any given mass M and coupling g.
If our model describes the interaction of two particles, where a bound state with the mass M is present, the scattering partial wave amplitude (19) should have a pole at s = M 2 , The physical coupling of the bound state, g, is given by the residue of the scattering partial-wave amplitude at the bound state pole One can use above equations to trade the bare parameters for the physical ones in the expression of the scattering amplitude and write the latter in terms of M and Z: Note that the above amplitude does not depend on the subtraction constant that renders G(s) finite.
This model can describe a bound state with any given value of the wave function renormalization constant.
Next, we study the finite volume effects in the bound-state mass. In the actual calculations, we take can be seen in Fig. 2, Eq. (26) is able to reproduce the synthetic lattice results very accurately. On the other hand, note that for shallow bound states the binding momentum κ is small, so no wonder that the expansion in ∆G θ L converges rather slowly. Consequently, retaining only the leading-order term and constructing iterative solution, see Eq. (27), might not be sufficient in all cases.
In the right panel of the same figure we show the dependence of the bound-state mass on the twisting angle θ = (θ, θ, θ) for the fixed value of Lm π = 3. We see that, for such a choice of twisting, the size of the effect of twisting for a fixed L is almost the double of the maximal effect caused by the variation of L from the same value to infinity (periodic boundary conditions). Thus, using (partially) twisted boundary conditions to determine Z, besides being cheaper, could give more accurate results than a method based on the study of the volume-dependence of the energy level. Note also that, for the above choice of the twisting angle, the twisting effect is maximal. Other choices, e.g., θ = (0, 0, θ) lead to a smaller effect. Now we turn our attention to the realistic case of the hadronic bound state D * s0 (2317) in the DK scattering channel with isospin I = 0 and strangeness S = 1. When isospin symmetry is exact, this state is stable under strong interactions, since it does not couple to the lighter hadronic channels (the observed decay D * s0 (2317) → D s π breaks isospin symmetry). Thus, the formalism above, tailored for stable bound states, does apply in this case. The case of quasi-bound states, which are coupled to inelastic channels, requires special treatment and is not addressed here.
A popular view on the D * s0 (2317) meson is that this state is dynamically generated as a pole through the S-wave interactions between the D-meson and the kaon in the isoscalar channel [36,[55][56][57][58][59]. We shall study this system, using the model used from Ref. [56], which is based on the leading-order heavy flavor chiral Lagrangian [60][61][62] and unitarizes the amplitude [6,7,63,64]. Namely, the infinite-volume amplitude is obtained from Eq. (19) with the S-wave-projected potential where x = cos θ is the cosine of the scattering angle, f π ≃ 92.4 MeV is the pion decay constant, and s and u are usual Mandelstam variables. We regularize the loop function with a subtraction constant a(µ), as done in Refs. [56,65]. Its value at the scale µ = m D is taken to be a(m D ) = −0.71. With this value of the subtraction constant, we find a bound state pole, associated with the D * s0 (2317), at M = 2316.9 MeV, and the coupling to DK, which is given by the residue of the pole, takes the value g = 10.7 GeV. One can easily calculate the compositeness parameter of the bound state as well, using Eq. (18). The calculation yields Z = 0.29. Hence, in this model, the D * s (2317) is predominately a molecular state.
Next, we study this model in a finite volume and consider twisting of different quarks, from which the D and K mesons consist. The net effect is that these mesons get different momenta as a result of such twisting, so the expression for G θ L changes. Note that this issue is important in view of the fact that partial twisting is allowed only for certain quarks (see Section V for more details).
In Fig. 3, we display the volume dependence of the bound state mass for different twisting angles which are again chosen as θ = (θ, θ, θ). In the left panel, we plot the L-dependence for three different values of the twisting angle, when twisted boundary conditions are applied to the u-quark. In the right panel, twisted boundary conditions are applied to the s-quark. As we shall see later, in the latter case the use of partial twisting gives the same results as using fully twisted boundary conditions. The size of the finite volume effects, using twisted boundary conditions for the c-quark, is very small, so we do not discuss this case. In this model, we test again that the predictions obtained from Eq. and M . On the other hand, we see that retaining only the leading exponential in the expansion of G θ L will have a large impact on the accuracy. Consequently, the first few terms should be retained. We see that the convergence is satisfactory: e.g., taking n max ≥ 3, where n max denotes the number of terms retained in the expansion, we see that the largest difference between the synthetic data and the prediction from Eq. (26) is less than 0.1 MeV.
Analyzing Fig. 3, we again come to the conclusion that the use of (partially) twisted boundary conditions can provide a better way to extract the compositeness parameter Z from lattice results. This can already be seen by comparing the curves for θ = 0 and θ = π. One namely observes that the size of the effect due to twisting at a fixed volume is almost twice as big as due to changing the volume for periodic boundary conditions. In Fig. 4, for three different volumes, we show the dependence of the bound-state mass on the twisting angle both for u-and s-quark twisting. On the other hand, taking the results of the θ-dependence (like in Fig. 4) at a fixed volume for granted, one could fit the value of the infinite-volume mass and coupling constant to these data, using Eq. (26). After this, it is straightforward to obtain the value of Z. In fact, producing four synthetic lattice data points at a fixed Lm π = 2.5 and θ = 0, π/3, 2π/3, π (either for u-or s-quark twisting), we were able to obtain values for M and g that differ less than 1% from those calculated from the infinite volume model by fitting the solution to Eq. (26) (with n max = 5) to the synthetic data.
Real lattice simulations, however, produce results which carry uncertainties. Hence, the question arises, how big these errors could be in order to be still able to determine Z with a desired accuracy. Since, as seen from the figures, the finite volume effects (for reasonable volume sizes, say, above Lm π = 2.5) are at most around 10 MeV, one expects that a relatively high accuracy will be needed in the measurement of the bound-state energy. In order to determine, how high this accuracy should actually be, we assign an uncertainty to the synthetic data that we generate from our model. In particular, using the von Neumann rejection method, from the "exact" data points we generate a new, "randomized" data set, where the central values of each data point are shifted randomly, following the Gaussian distribution centered at exact data values and with a standard deviation, given by the lattice data error. Repeating this process several times, we obtain several sets of synthetic lattice data with errors and central values shifted accordingly. We then fit each of the randomized data sets and obtain a corresponding value for M and g (and therefore, for Z), one for each set, ending up with as many values for the parameters, as many randomized data sets we have generated. We can obtain then the mean and standard deviation of the distributions for M , g and Z. Thus, for a given data error, we can estimate the accuracy of the parameter extraction.
For the case of the s-quark twisting, we construct 5000 sets of randomized data at a fixed volume, for different input errors ∆M L and different number of data points per set. Fitting the parameters to each set, we obtain the corresponding distributions of 5000 points for each parameter M , g and Z.
In Index Channel Quark content the input lattice data is even bigger. If we increase the number of lattice data points, we get slightly better results but, in general, the dependence on the increase of the size of the data set is very mild.
For example, we need to use around 20 data points to achieve an accuracy of order 0.1 in Z, given an input error ∆M L = 2 MeV and volume Lm π = 2.5.

V. PARTIALLY TWISTED BOUNDARY CONDITIONS IN THE DK SYSTEM
The partial twisting, unlike the full twisting, is more affordable in terms of computational cost in lattice simulations, because one does not need to generate new gauge configurations. Thus, it is very interesting to study whether it is possible to extract any physically relevant information from simulations using this kind of boundary conditions. Problems may arise when there are annihilation channels present, as is the case in the DK scattering in the isoscalar channel, where light quarks may annihilate. An analysis of Lüscher approach with partial twisting for scattering problem in the presence of annihilation channels was recently addressed in [39]. Namely, a modified partially twisted Lüscher equation was derived for the πη − KK coupled channel scattering in the framework of non-relativistic EFT.
Here, we address the same problem in the context of the DK scattering. The method is described in Ref. [39], to which the reader is referred for further details. Consider first the scattering in the infinite volume. We start from building the channel space by tracking the quarks of different species following where T , V and G are given by 3 × 3 matrices.
The free Green function is given by where G(s) is defined in Eqs. (11) and (16), supplemented by the prescription that the integral is performed in dimensional regularization after expanding the integrand in powers of 3-momenta (see Refs. [39,66] for details). The minus sign on the diagonal of the matrix G arises due to fermionic nature of D and K mesons composed of valence and (commuting) ghost quarks.
The crucial point now is that there exist linear symmetry relations between various elements of T due to equal valence, sea and ghost quark masses. Note that scattering matrix elements are given by residues of the 4-point Green functions Γ ij of the bilinear quark operators at the poles, corresponding to the external mesonic legs. Decomposing Γ ij into connected t c and disconnected t d pieces through Wick contractions (see Fig. 5) and noting that quark propagators are the same for all light quark species, we Since in our case there are no neutral states and thus no mixing occurs, following the argumentation given in Ref. [39], it is easy to show that T -matrix obeys the same symmetry relations as Γ Here T 11 = t corresponds to the physical elastic DK scattering amplitude, i.e scattering in the sector with valence quarks only. Other diagonal entries are unphysical in the sense that they correspond to scattering of particles, composed of sea and ghost light quarks. Non-diagonal elements of T -matrix describe coupling between valence and sea/ghost sectors through disconnected diagrams. Furthermore, it is straightforward to check from Eq. (34) that the elements of potential matrix V satisfy the same symmetry relations as T and can be expressed in the following form Let us now turn to the case of a finite volume and derive the Lüscher equation for a couple of particular choices of partially twisted boundary conditions. Note that the potential V remains the same (up to exponentially suppressed in terms L ), while in the loop functions the integration is substituted by summation over lattice momenta.
1. Twist the s/c-quark, leaving u and d-quarks to obey periodic boundary condition. In this case, the matrix of the Green functions is diag G θ L ,G 0 L , −G 0 L . The solution of the Lippmann-Schwinger equation in a finite volume for the physical amplitude t is given by whereG θ L is the loop function G(s) in a finite volume. We see that the finite-volume spectrum in case of the partial twisting is determined from the Lüscher equation in the same way as in the full-twisting case. Thus, the results obtained by using of the partially twisted boundary conditions on the c-or s-quark are equivalent to those using full twisting.
2. Twist the valence u-and d-quarks simultaneously, leaving s-and c-quarks obey periodic boundary condition. In this case, the ghost light quarks also need to be twisted, and the matrix of the Green The Lüscher equation determining the finite volume spectrum now takes the form Vanishing of the first bracket on the r.h.s gives the Lüscher equation with no twisting. Note also that the quantity τ − υ is in fact the connected part of the scattering potential for the isoscalar DK system, which is identical to the DK scattering potential in the isovector channel. Hence, vanishing of the second bracket is equivalent to the fully twisted Lüscher equation for the isovector DK scattering 4 .

VI. SUMMARY AND CONCLUSIONS
i) Lattice QCD does not only determine the hadron spectrum. Under certain circumstances, it may provide information about the nature of hadrons, which renders lattice simulations extremely useful for the search and the identification of exotic states. Note that the lattice QCD possesses unique tools at its disposal (e.g., the study of the volume and quark mass dependence of the measured quantities), which are not available to experiment.
ii) In the present paper, we concentrate on the identification of hadronic molecules on the lattice.
Experimentally, one may apply Weinberg's compositeness condition to the near-threshold bound states, in order to distinguish the molecular states from the elementary ones. To this end, one may use the value of the wave function renormalization constant Z which obeys the inequalities 0 ≤ Z ≤ 1. The vanishing value of the parameter Z corresponds to the purely molecular state.
In this paper we consider the lattice version of the Weinberg's condition.
iii) It is known that the quantity Z can be extracted from lattice data by studying the volume dependence of the measured energy spectrum. We have shown that the same result can be achieved by measuring the dependence of the spectrum on the twisting angle in case of twisted boundary conditions. Moreover, within the method proposed, the expected effect is approximately twice as large in magnitude and comes at a lower computational cost. Further, we have analyzed synthetic data to estimate the accuracy of the energy level measurement which is required for a reliable extraction of the value of Z on the lattice. iv) As an illustration of the method, we consider the D * s0 (2317) meson, which is a candidate of a DK molecular state. It is proven that, despite the presence of the so-called annihilation diagrams, one may still use the partially twisted boundary conditions for the extraction of Z from data if the charm or strange quark is twisted. The effects which emerge due to partial twisting, are suppressed at large volumes. We compute the scattering amplitude in a finite volume by replacing the loop function G by its finite volume counterpartG θ L = G + ∆G θ L and obtain synthetic data from the poles of the finite volume scattering amplitude. In particular, the pole below threshold gives the mass of the bound state in a finite volume. For the case of a level below threshold, there exists a fairly simple way to calculate ∆G θ L defined by Eq. (20), so that the equation (26) for κ L can be easily solved. Here, we consider three different cases, one with periodic boundary conditions, and two with twisted boundary conditions. Depending on which quarks are twisted, the momenta of the mesons are modified accordingly.

Periodic boundary conditions
In the case of periodic boundary conditions, the meson momenta in a box are given by q n = 2π n L , n ∈ Z 3 .
We can evaluate the sum in Eq. (20), using the Poisson summation formula n δ(n − x) = n e 2πinx .
Transforming the sum into the integral gives Next, we note that the integrand I( q ) can be approximated by 1 2 √ s 1 k 2 − q 2 , since the difference is exponentially suppressed [53]. Here, k 2 is the three-momentum squared of the particles in the center of mass (c.m.) frame. Then, for k 2 < 0, ∆G θ L reads The function ∆G 0 L can be expressed in terms of the Lüscher zeta-function Z 00 (1,k 2 ), as follows [53]: wherek = kL/(2π).

Twisted boundary conditions: both momenta shifted
In the case of twisted boundary conditions, when the momenta of both particles are shifted but the particles still are in the c.m. frame, the allowed momenta in a box are: where θ is the twisting angle. Now, acting in the same way, we can evaluate the sum in Eq.
Again, we can express ∆G θ L in terms of the Lüscher zeta-function with twisted boundary conditions, Z θ 00 (1,k 2 ), as follows,

Twisted boundary conditions: only one momentum shifted
Finally, in the case of twisted boundary conditions, when only the momentum of one of the particles (say, particle 1) is shifted, the allowed momenta in a box are The particles are not in the c.m. frame any more: the c.m. momentum is equal to P = θ/L. Hence, we have to evaluate ∆G θ L in a moving frame with momentum P , Again, we can approximate the integrand by [67] I( q ) = − 1 2P 0 1 ( q ′ ) 2 − ( q ′ · P ) 2 /P 2 0 − k 2 + · · · , q ′ = q − µ P , where µ = 1 2 1 − where n and n ⊥ are the components parallel and perpendicular to P of n, and γ = P 0 / √ s is the relativistic gamma-factor. Once again, we can relate ∆G θ L in this case with the Lüscher zeta function in the moving frame Z d 00 (1; (q * ) 2 ) [68], see also Refs. [67,69,70]: where d = P L/2π = θ/2π. For the case of θ = (θ, θ, θ), the first few terms in the above expansion are ∆G (θ,θ,θ) L (M ) = − 1 8πM L 6 √ 3 cos(µθ) In the case of shallow bound states, the exponential factor κ will be usually quite small, so in order to reproduce accurately the full function, one should take several terms in the expansion for ∆G θ L above.