The Optical Potential on the Lattice

The extraction of hadron-hadron scattering parameters from lattice data by using the L\"uscher approach becomes increasingly complicated in the presence of inelastic channels. We propose a method for the direct extraction of the complex hadron-hadron optical potential on the lattice, which does not require the use of the multi-channel L\"uscher formalism. Moreover, this method is applicable without modifications if some inelastic channels contain three or more particles.


Introduction
The Lüscher approach [1] has become a standard tool to study hadron-hadron scattering processes on the lattice. The use of this approach in case of elastic scattering is conceptually straightforward: besides technical complications, caused by partial-wave mixing, each measured energy level at a given volume uniquely determines the value of the elastic phase shift at the same energy.
In the presence of multiple channels, the extraction of the scattering phase becomes more involved. In case when only two-particle coupled channels appear, one can make use of the coupledchannel Lüscher equation [2][3][4][5][6][7][8] and fit a simple pole parameterization for the multi-channel Kmatrix elements to the measured energy spectrum in the finite volume [9]. A more sophisticated parameterization of the K-matrix elements, which is applicable in a wider range of the energies, can be obtained using unitarized chiral perturbation theory (ChPT) [10][11][12]. This approach has been successfully applied, e.g., in Ref. [13] to analyze coupled-channel ππ − KK P -wave scattering and to study the properties of the ρ-meson. Note the difference to the one-channel case: here, one has to determine several K-matrix elements (unknowns) from a single measurement of a finitevolume energy level. Hence, using some kind of (phenomenology-inspired) parameterizations of the multi-channel K-matrix elements becomes inevitable in practical applications.
In case when some of the inelastic channels contain three or more particles, the situation is far more complicated. Despite the recent progress in the formulation of the theoretical framework [14][15][16][17], it is still too cumbersome to be directly used in the analysis of the data. Moreover, the problem of the choice of the parameterization for three-particle scattering might become more difficult (and lead to even larger theoretical uncertainties) than in two-particle scattering.
From the above discussion it is clear that a straightforward extension of the Lüscher approach through the inclusion of more channels has its limits that are reached rather quickly. On the other hand, many interesting systems, which are already studied on the lattice, may decay into multiple channels. In our opinion, the present situation warrants a rethinking of the paradigm. One may for example explore the possibility to analyze the lattice data without explicitly resolving the scattering into each coupled channel separately. Such a detailed information is usually not needed in practice. Instead, in the continuum scattering problem, the effect of inelastic channels could be included in the so-called optical potential [18,19], whose imaginary part is non-zero due to the presence of the open inelastic channels. In many cases, it would be sufficient to learn how one extracts the real and imaginary parts of the optical potential from the lattice data, without resorting to the multi-channel Lüscher approach. In the present paper, we propose such a method, which heavily relies on the use of twisted boundary conditions [20][21][22][23][24]. Due to this, the method has its own limitations, but there exist certain systems, where it could in principle be applied. In particular, we have the following systems in mind: • The scattering in the coupled-channel πη − KK system in the vicinity of the KK threshold and the a 0 (980) resonance.
• The spectrum and decays of the XY Z states; namely, Z c (3900) ± that couples to the channels J/ψπ ± , h c π ± and (DD * ) ± (this system was recently studied in Ref. [25]) or the Z c (4025) that couples to the D * D * and h c π channels (see, e.g., Ref. [26]).
There certainly exist other systems where this method can be used. It should also be stressed that the systems, where the partial twisting (i.e., twisting only the valence quarks) can be carried out, are interesting in the first place -for an obvious reason. All examples listed above belong to this class. In general, the partial twisting can always be carried out when the annihilation diagrams are absent. In the presence of annihilation diagrams, each particular case should be analyzed separately, invoking the methods of effective field theories in a finite volume [24]. The present paper contains an example of such an analysis.
Further, note that there exists an alternative method for the extraction of hadron-hadron interaction potentials from the measured Bethe-Salpeter wave functions on the Euclidean lattice. This method goes under the name of the HAL QCD approach and its essentials are explained in Refs. [27,28]. Most interesting in the present context is the claim that the HAL QCD approach can be extended to the multi-channel systems, including the channels that contain three and moreparticles [29]. It should also be pointed out that this approach has already been used to study various systems on the lattice, including the analysis of coupled-channel baryon-baryon scattering (see, e.g., Ref. [30]). It would be interesting to compare our method with the HAL QCD approach.
The layout of the present paper is as follows. In Sect. 2 we discuss the theoretical framework for the extraction of the real and imaginary parts of the optical potential and provide an illustration of the method with synthetic data, generated by using unitarized ChPT. Further, in Sect. 3, we discuss the role of twisted boundary conditions for measuring the optical potential. Namely, the possibility of imposing partially twisted boundary conditions is explored in Sect. 3.1. Here, we also discuss the possibility of imposing the different boundary conditions on the quark and antiquark fields. The analysis of synthetic data, including an error analysis, is presented in Sect. 3.2. Finally, Sect. 4 contains our conclusions.
2 Optical potential in the Lüscher approach

Multichannel potential, projection operators
In the continuum scattering theory, the inelastic channels can be effectively included in the so-called optical potential by using the Feshbach projection operator technique [18]. Namely, let us start from the multi-channel T -matrix which obeys the Lippmann-Schwinger equation Here, V is the potential and G 0 = (E − H 0 ) −1 denotes the free Green's function with E the total energy in the center-of-mass system. The quantities T, V, G 0 are all N × N matrices in channel space.
In case when only two-particle intermediate states are present, using dimensional regularization together with the threshold expansion, it can be shown that the Lippmann-Schwinger equation (1) after partial-wave expansion reduces to an algebraic matrix equation (see, e.g., Ref. [31]). With the proper choice of normalization, the matrix G 0 (E) in this case takes the form G 0 (E) = diag (ip 1 (E), · · · , ip n (E)) , where p k (E) denotes the magnitude of the center-of-mass three-momentum, i.e., and m (k) 1,2 are the masses of particles in the k th scattering channel. Hence, if dimensional regularization is used in case of two-particle channels, the potential V coincides with the multi-channel K matrix. The latter quantity can always be defined, irrespectively of the used regularization. Our final results are of course independent of the use of a particular regularization.
Suppose further that we focus on the scattering in a given two-particle channel. Let us introduce the projection operators P and Q = 1 − P , which project on this channel and on the rest, respectively. In the following, we refer to them as the primary (index P ) and the secondary (index Q) channels. The secondary channels may contain an arbitrary number of particles. It is then straightforward to show that the quantity T P (E) = P T (E)P obeys the following single-channel Lippmann-Schwinger equation where W (E) = P V + V Q 1 E − H 0 − QV Q QV P and G P (E) = P G 0 (E)P .
It is easily seen that, while V is Hermitean, W (E) above the secondary threshold(s) is not. The imaginary part of W (E) is expressed through the transition amplitudes into the secondary channels where For illustration, let us consider scattering in the πη − KK coupled channels. Let KK and πη be the primary and secondary channels, respectively. Then, the formulae for the S-wave scattering take the following form (we suppress the partial-wave indices for brevity): Here, p KK , p πη denote the magnitude of the relative three-momenta in the center-of-mass frame in the respective channel, as given in Eq. (3). It is often useful to introduce the so-called M -matrix M = V −1 . In terms of this quantity, the above formula can be rewritten in the following form: Using the latter form can be justified, when a resonance near the elastic threshold exists that shows up as a pole on the real axis in V . In contrast, the quantity M is smooth in this case and can be Taylor-expanded near threshold. In a finite volume, one may define a counterpart of the scattering amplitude T KK→KK (E). Imposing, e.g., periodic boundary conditions leads to the modification of the loop functions (for simplicity, we restrict ourselves to the S-waves from here on and neglect partial wave mixing) whereas the potential V remains unchanged up to exponentially suppressed corrections. In the above expressions, L is the size of the cubic box and Z 00 denotes the Lüscher zeta-function. The energy levels of a system in a finite volume coincide with the poles of the modified scattering amplitude. The position of these poles is determined from the secular equation The positions of these poles on the real axis are the quantities that are measured on the lattice.

Continuation to the complex energy plane
The main question, which we are trying to answer, can now be formulated as follows: Is it possible to extract the real and imaginary parts of W (E) from the measurements performed on the lattice? We expect that the answer exists and is positive, for the following reason. Let us imagine that all scattering experiments in Nature are performed in a very large hall with certain boundary conditions imposed on its walls. It is a priori clear that nothing could change in the interpretation of the results of this experiment, if the walls are moved to infinity. Consequently, there should exist a consistent definition of the infinite-volume limit in a finite-volume theory that yields all quantities defined within the scattering theory in the continuum. Since the optical potential is one of these, there should exist a quantity defined in a finite volume, which coincides with the optical potential in the infinite-volume limit.
In order to find out, which quantity corresponds to the optical potential in a finite volume and how the infinite-volume limit should be performed, let us follow the same pattern as in the infinite volume. Namely, we apply the one-channel Lüscher equation for the analysis of data, instead of the two-channel one. As a result, we get: .
The left-hand side of this equation is measured on the lattice at fixed values of p KK , corresponding to the discrete energy levels in a finite volume. Methods to measure W −1 L are discussed in Sect. 3. The quantity on the right-hand side is proportional to the cotangent of the so-called pseudophase, defined as the phase extracted with the one-channel Lüscher equation [2,3]. It coincides with the usual scattering phase in the absence of secondary channels. Fig. 1 shows the real and imaginary parts of the quantity W −1 (E) that is constructed by using a simple parameterization of the two-channel T -matrix, based on unitarized ChPT (see Ref. [32]). For comparison, the finite-volume counterpart W −1 L (E), which is defined by Eq. (13), is also shown. If the secondary channels were absent, W −1 (E) would be real and equal to W −1 L (E), up to exponentially suppressed contributions. Fig. 1 clearly demonstrates the effect of neglecting the secondary channels. While the "true" function W −1 (E) is a smooth (and complex) function of energy, the (real) function W −1 L (E) has a tower of poles and zeros. The (simple) zeros of W −1 L (E) (poles of W L (E)) emerge, when E coincides with one of the energy levels in the interacting πη system. The background, obtained by subtracting all simple poles, is a smooth function of E. It should be stressed that this statement stays valid even in the presence of multiple secondary channels, some of which containing three or more particles. The only singularities that emerge in general are the simple poles that can be traced back to the eigenvalues of the total Hamiltonian  Figure 3: The real and imaginary parts of the quantityŴ −1 L (E + iε) for ε = 0.02 GeV (solid black lines) and ε = 0.05 GeV (dashed blue lines) versus the real and imaginary parts of the infinite-volume counterpart W −1 (E) (dotted red lines). All quantities are given in units of GeV.
restricted to the subspace of the secondary states 1 .
It is important to note that, if L tends to infinity, the optical potential does not have a welldefined limit at a given energy. As the energy levels in the secondary channel(s) condense towards the threshold, the quantity W −1 L (E) at a fixed E oscillates from −∞ to +∞. Thus, the question arises, how the quantity W −1 (E) can be obtained in the infinite-volume limit.
It should be pointed out that this question has been already addressed in the literature in the past. In this respect, we find Ref. [33] most useful. In this paper it is pointed out that, in order to give a correct causal description of the scattering process, one should consider adiabatic switching of the interaction. This is equivalent to attaching an infinitesimal imaginary part E → E + iε to the energy. Further, as argued in Ref. [33], the limits L → ∞ and ε → 0 are not interchangeable. A correct infinite-volume limit is obtained, when L → ∞ is performed first (see Ref. [34] for a more detailed discussion of this issue). Physically, this statement is clear. The quantity ε defines the available energy resolution, and the distance between the neighboring energy levels tends to zero for L → ∞. If this distance becomes smaller than the energy resolution, the discrete levels merge into a cut and the infinite-volume limit is achieved. It is also clear, why the infinite-volume limit does not exist on the real axis: ε = 0 corresponds to an infinitely sharp resolution and the cut is never observed.
The above qualitative discussion can be related to Lüscher's regular summation theorem [35]. On the real axis above threshold, the zeta-function Z 00 (1; q 2 πη ) in Eq. (13) does not have a welldefined limit. Assume, however, that the energy E gets a small positive imaginary part, E → E +iε. The variable q 2 πη also becomes imaginary: It is immediately seen that above threshold, E > M η + M π , the quantity ε is strictly positive. Now, for real energies E, the nearest singularity is located at the distance ε from the real axis, so the regular summation theorem can be applied. It can be straightforwardly verified that the remainder term in this theorem vanishes as exp(−ε L) (modulo powers of L), when L → ∞. The above argumentation can be readily extended to the cases when intermediate states contain any number of particles. Consider a generic loop diagram in the effective field theory where these particles appear as internal lines. It is most convenient to use old-fashioned time-ordered perturbation theory, where the integrand contains the energy denominator (E + iε − w 1 (p 1 ) . . . − w n (p n )) −1 . Here, w i (p 1 ) , i = 1, . . . , n stand for the (real) energies of the individual particles in the intermediate state. It is clear that, if ε = 0, the denominator never vanishes, and the regular summation theorem can be applied. The remainder, as in the two-particle case, vanishes exponentially when ε = 0.
The analytic continuation into the complex plane can be done as follows. Suppose one can measure the quantity W −1 L (E) on the real axis. Bearing in mind the above discussion, one may fit this function by a sum of simple poles plus a regular background. Fig. 2 shows the result of such a fit which was performed by using the function to fit a sample of the exact W −1 L without errors. The exact values of the fit parameters are not listed here since Fig. 2 is given for the illustrative purposes only. In the actual numerical simulation of Sect. 3.2, the order of the polynomial is varied.
The continuation into the complex plane is trivial: one uses Eq. (15) with fixed values of Z i , Y i , D i and makes the substitution E → E + iε. The real and imaginary parts of the quantitŷ W −1 L (E + iε) for ε = 0.02 GeV and ε = 0.05 GeV are shown in Fig. 3. For comparison, the real and imaginary parts of the infinite-volume counterpart W −1 (E) are also given. As seen, the finitevolume "optical potential" oscillates around the true one and the magnitude of such oscillation grows larger, when ε becomes smaller. On the other hand, the artifacts caused by a finite ε grow, when ε becomes large.

Infinite-Volume Extrapolation
From the above discussion it is clear that, performing the limit L → ∞ for a fixed ε, and then taking ε → 0, the infinite-volume limit is restored fromŴ −1 L (E + iε). For the actual extraction on the lattice, however, taking the large volume limit could be barely feasible. An alternative to this procedure is to "smooth" the oscillations arising from Eq. (15) if evaluated at complex energies at a finite L and ε. This allows one to perform the extraction of the optical potential at a reasonable accuracy even at sufficiently small values of L. As in the present study the true optical potential is known, the validity of this procedure can be tested. We would like to stress that LM π = 5 used in this study is rather small and thus not completely beyond reach.
In the present section we test two different algorithms for smoothing the quantityŴ −1 L (E + iε). In both cases, the result is calledŴ −1 , i.e., the estimate of the true infinite-volume potential W −1 . The final results of the numerical studies, presented in Sect. 3.2 are evaluated with both methods.

Parametric method
The basic idea of this method is to fit the optical potentialŴ −1 L (E + iε) from Eq. (15) at complex energies in the whole energy range with a suitable Ansatz. Model selection is performed with LASSO regularization (as explained in detail later) in combination with cross validation. Such methods have the advantage that basic properties of the optical potential, like Schwartz's reflection  Figure 4: Left: The χ 2 as a function of the degree of the fit polynomial, n max . While the χ 2 of the unconstrained fits (gray squares) monotonically decreases, a finite penalty factor of λ =λ opt = 0.2 for P 2 stabilizes the result (red triangles). Right: Cross validation. The χ 2 of the fits to the training set according to Eq. (17) are shown with gray squares; the χ 2 V of these fits, evaluated for the test/validation set, are indicated with red triangles; the χ 2 t of these fits evaluated for the (unknown) true optical potential according to Eq. (20) are displayed with blue circles. The minimum of the χ 2 V (λ) of the test/validation set (red) estimates the penalty factorλ opt ∼ 0.15 − 0.2 which is very close to the truly optimal λ opt ∼ 0.2 − 0.3 (blue). The absolute and relative scales of the different χ 2 's are irrelevant.
principle and threshold behavior, can be built in explicitly. In our problem, this is particularly simple because the only non-analyticity is given by the branch point at the πη threshold. In more complex problems, additional non-analyticities like resonance poles or complex branch points from multi-channel states [36,37] have to be included in the parameterization. Yet, all these non-analyticities are situated on other than the first Riemann sheet. The parametric and nonparametric methods proposed here use an extrapolation from finite, but positive ε to ε → 0, i.e., an extrapolation performed on the first Riemann sheet that is analytic by causality.
A suitable yet sufficiently general parameterization of the optical potential is given bŷ with real parameters a j , b j . The only non-analyticity ofŴ −1 is given by the cusp function i p πη , evaluated at the complex energy E (see Eq. (10)), that is therefore explicitly included in the Ansatz; the rest is then analytic and can be expanded in a power series around a real E 0 chosen in the center of the considered energy region, in order to reduce correlations among fit parameters (the actual value of E 0 is irrelevant).
To perform the effective infinite-volume extrapolation through smoothing, we consider the minimization of the χ 2 , where P i are penalty functions specified below. The absolute scale of the χ 2 is irrelevant. The quantityŴ −1 L is fitted by sampling at the complex energies E k = E min + k δE + iε (ε = 0.05 GeV)  over the considered energy range E min ≤ E ≤ E max with a step δE = 10 MeV, and assigning an arbitrary error of σ k = σ = 1 GeV. Note that in cross validation (to be specified later), the position of the minimal χ 2 determines the size of the penalty, i.e., the size of σ is irrelevant. The infinite-volume optical potential is then obtained by simply evaluatingŴ −1 at real energies, i.e., setting ε = 0. If we assume for the moment that the penalty function P i in Eq. (17) is zero, then it is clear that the minimized χ 2 is a monotonically decreasing function of the degree of the fit polynomial n max . This is demonstrated by the gray squares in Fig. 4 (left panel). Apparently, the fit stabilizes first for n max = 3 − 6, which might lead to the wrong conclusion that an optimal smoothing had been obtained. Then, for higher n max , another plateau is reached at n max = 7 − 9 and then another one for n max = 10 − 14. Thus, without an additional criterion, one cannot decide which n max is optimal.
In general, for a small n max , the smoothing will be too aggressive (large χ 2 ), while for too large values of n max the fit will start following the oscillations (Fig. 3), resulting in a low χ 2 but missing the point of smearing the optical potential. These two extreme cases are illustrated in Fig. 5 with the thin dashed and thin solid lines, respectively 2 . There is obviously a sweet spot for n max . Model selection refers to the process of determining this spot as outlined in the following.
Model selection for the fit (16) is formally introduced through a penalty P (a j , b j ) imposed on the fit parameters. The penalty is formulated using the LASSO method developed by Tibshirani in 1996 [38]. See also Refs. [39,40] for an introduction into the topic. The LASSO method has been recently applied in hadronic physics for the purpose of amplitude selection [41].
A natural choice to suppress oscillations is to penalize the modulus of the second derivative [39], where the integral is performed along a straight line in the complex plane. Another choice is to penalize only the polynomial part of the Ansatz (16), i.e., removing the p πη factor that has an inherently large second derivative at the πη threshold, Including λ to the fourth power is done in order to have a clearer graphical representation of the penalty factor in subsequent plots. Imposing a penalty, the decrease of χ 2 with n max is eventually stabilized, as shown by the red triangles in Fig. 4 (left panel) for some yet to be determined value of λ. Clearly, the minimized χ 2 from Eq. (17) is a monotonically increasing function of λ as demonstrated by the gray squares in Fig. 4 (right panel) for the penalty function P 2 .
The fitted data (ε = 0.05 GeV) form the so-called training set [38]. The main idea of cross validation to determine the sweet spot of λ is as follows (for more formal definitions and k-fold cross validation, see Refs. [38][39][40]): after a random division of a given data set into training and test/validation sets, the fit obtained from the training set is used to evaluate its χ 2 with respect to the test/validation set, called χ 2 V in the following (without changing fit parameters and setting P i = 0). For too large values of λ, both values of χ 2 from training and from test/validation sets will be large. For too small λ, the fit to the training set is too unconstrained and sensitive to unwanted random properties such as fluctuations in the training data. However, those unwanted random properties are different in the validation set, leading to a worse χ 2 V for too small λ. It is then clear that χ 2 V (λ) exhibits a minimum at the sweet spot λ =λ opt . Here, we cannot meaningfully divide the data set randomly. Instead, we have to look for data, for which the physical property (infinite-volume optical potential) is unchanged, but the unphysical property (oscillations from finite-volume poles) is changed. This is naturally given byŴ −1 L but evaluated for a substantially different value of ε (we choose ε = 0.15 GeV). The analytic form of Eq. (16) ensures that the infinite-volume optical potential can be analytically continued to different values of ε, and only the unwanted finite-volume oscillations are different for different ε. Indeed, as indicated with the red triangles in Fig. 4 (right panel), χ 2 V exhibits a clear minimum at λ =λ opt ∼ 0.2. The potential dependence of the this value on the chosen ε is discussed below. Furthermore, in this example, we know the underlying optical potential and can simply determine the (generally unknown) truly optimal value for λ, λ opt by evaluating the χ 2 of the estimate of the optical potential,Ŵ −1 , with respect to the true optical potential on the real axis, W −1 , Note that the quantity χ 2 t (λ) (implicitly) depends on λ, because the quantityŴ −1 (Re E k ) was determined at a fixed value of λ. The quantity χ 2 t is shown with the blue filled circles in Fig. 4 (right panel). Its minimum at λ = λ opt is very close to the minimum of the validation χ 2 V at λ =λ opt , demonstrating that cross validation [39] is indeed capable of estimating the optimal penalty in our case.
Instead of using the penalty function P 2 , one can also choose P 1 , see Eqs. (18) and (19). The estimatedλ opt given by the minimum of χ 2 V will, of course, change. But, again, it was checked that the newλ opt is very close to the new λ opt given by the minimum of χ 2 t . Similarly, we have checked other forms of penalization, with the same findings: imposing penalty on the third derivative, variation of the value of ε for the training set, and variation of the value of ε for the test/validation set. The only restriction is that the ε of the test/validation set has to be chosen sufficiently larger than ε of the training set for a minimum in χ 2 V to emerge -if the two ε's are too similar, the oscillations are too similar and no minimum in χ 2 V is obtained. Also, n max has to be chosen high enough so that, at a given ε for the training set, the fit is capable of fitting oscillations (for small λ) which is a prerequisite for a minimum in χ 2 V to appear. In all simulations we have chosen n max = 18 although n max ∼ 7 would suffice as the left panel of Fig. 4 shows.
For the initially considered case, using P 2 for the penalty, ε = 0.05 GeV for the training set, and ε = 0.15 GeV for the test/validation set, the resulting optical potential is shown in Fig. 5 with the thick black solid lines. For comparison, the true optical potential is shown with the thick red (dashed) lines. The optical potential is well reconstructed over the entire energy range. At the πη threshold, the reconstructed potential reproduces the square-root behavior due to the explicit factor p πη in the parameterization (16). The reconstructed potential explicitly fulfills Schwartz's reflection principle and its imaginary part is zero below threshold. At the highest energies, small oscillations become visible originating from the upper limit of the fitted region at E max = 1.7 GeV. Here, the smoothing algorithm, that is an averaging in energy, has simply no information onŴ −1 L beyond E max . Note that in the numerical simulation of the next section, that uses re-sampling techniques and realistic error bars, these small oscillations themselves average out over the Monte-Carlo ensemble, simply resulting in a widened, but smooth, error band at the highest energies.
For illustration, we also show in Fig. 5 a largely under-constrained result (too small λ, thin solid lines) in which the oscillations from the finite-volume poles inŴ −1 L survive. The opposite case, i.e., an over-constrained fit with too large λ, is shown with the thin dashed lines exhibiting too large of a penalization on the second derivative.

Non-parametric method
The advantage of non-parametric methods lies in its blindness of analytic structures, which, however, also leads to the fact that threshold behavior and Schwartz' reflection principle cannot be implemented easily. As a particular method, we utilize an approach, commonly used in image processing applications. This approach goes under the name of Gaussian smearing. The basic idea of the Gaussian smearing is quite simple: for a given set of uniformly distributed data, any data point is replaced by a linear combination of its neighboring data points (within a given radius r), with individual weights, w(x) given by Here, x and σ 0 denote the distance of the individual points from the central one and the standard deviation. Typically, the latter value is chosen to match the radius of the smearing by σ 0 = r/2. Therefore, the only undetermined quantity is given by the smearing radius r.
The general prescription to determine the smearing radius should rely on the properties of the original data only. Recall that the latter is determined by the functionŴ −1 L in Eq. (15), which splits up into a real and an imaginary set, when evaluated at the energy E + iε for a fixed ε > 0 and uniformly distributed values of E. Therefore, after the fits to the (synthetic) lattice data are performed, the scale of the structures to be smeared is determined by the distance between two poles, see Fig. 8. Of course, since the poles are not distributed uniformly over the whole energy range, one could argue in favor of using different values of r for different energies. It is also clear that constraint on the standard deviation σ 0 = r/2 affects the result of the smoothing. However, in order not to over-complicate the procedure, in the following we choose the smearing radius to be twice as large as the typical (average) distance between two poles. If the radius is much larger than this, physical information (i.e. the functional form of the optical potential) will be smeared out. If, however, the radius is much smaller than this value, then the (unphysical) oscillations will remain, preventing the reconstruction of the underlying optical potential. The situation is in fact very similar to the under-and over-constrained results, discussed in the previous section for the too small and too large values of λ.
After the parameters of the smearing kernel (21) are fixed, the method is applied to the sets of real and imaginary parts ofŴ −1 L at a fixed ε > 0. Then the procedure is repeated, each time assuming slightly smaller value of ε than before. In the final step, a simple (polynomial) extrapolation is performed to real energies, i.e. ε = 0, to obtain the final result of this procedure, namelyŴ −1 (E).
In this section, we have demonstrated that the real and imaginary parts of the optical potential can be reconstructed from the pseudophase measured on the lattice for real energies, W −1 L , if the analytic continuation into the complex plane is performed. Two distinct methods are presented to smear the oscillations which emerge from the analytic continuation, and to recover the optical potential for real energies. It remains to be seen, how the pseudophase can be measured in practice. This issue will be considered in the Sect. 3 where a realistic numerical simulation will be carried out as well.

Reconstruction of the optical potential
The quantity W −1 L (E), which is used to extract the optical potential, along with the energy E, depends on other external parameters, say, on the box size L, boundary conditions, etc. In the fit to W −1 L (E), the values of these parameters have to be fixed. Otherwise, for example, the position of the poles in W −1 L (E) will be volume-dependent and a fit is not possible. Hence, we are quite restricted in the ability to scan the variable E: the knob, which tunes E, must leave all other parameters in the pseudophase intact.

Partially twisted boundary conditions
In certain systems, there indeed exists a possibility to scan the energy within a given range in this manner. It is provided by the use of twisted boundary conditions and can be realized, e.g., in the coupled-channel πη − KK scattering. Namely, as was discussed in Refs. [3,24], in this system it is possible to apply (partially) twisted boundary conditions so that, when the twisting angle is changed continuously, the KK threshold moves, whereas the πη threshold stays intact. This can be achieved, for example, by twisting the light u, d quarks by the same angle and leaving the s-quark with periodic boundary conditions. This will lead to the modification of the secular equation (12), The expression for W −1 L (E) remains the same and does not contain the twisting angle θ. The method can be used to study the isospin I = 1 scattering in the πη −KK system. As shown in Ref. [24], despite the presence of the annihilation diagrams, the partial twisting in this case is equivalent to the full twisting, if the light quarks are twisted, whereas twisting of the s-quark does not lead to an observable effect. As a rule of thumb, one expects that the partial twisting of a given quark will be equivalent to full twisting, only if this quark line goes through the diagram without being annihilated (of course, a rigorous proof of this statement should follow by using effective field theory methods [24]). In our case, we could choose to work with the state with maximal projection  of the isospin, say I = 1, I 3 = 1. This state contains one u-quark and oned-quark, which cannot be annihilated. Choosing the same twisting angle for both quarks, the system stays in the center-ofmass frame and the pseudophase becomes independent from the twisting angle, as required. From the above discussion it is also clear that using our method for the extraction of the optical potential in the channel with isospin I = 0 implies the use of full twisting instead of partial twisting. The same trick can be used to study the Z c (3900) and Z c (4025) states, which both have isospin I = 1. Twisting u-and d-quarks by the same angle, the D-and D * -mesons will get additional momenta proportional to the twisting angle, whereas the J/ψ, h c and π-mesons will not. Consequently, one may choose the channels containing the D and D * mesons as the primary ones (in our nomenclature) and regard every other channel as secondary. For this choice, the pseudophase will not depend on the twisting angle.
Last but not least, an unconventional twisting procedure was used in the study of the J/ψφ scattering from Y (4140) decays [42]. Namely, in that work the c-and s-quarks were twisted by the angles θ and −θ, respectively, whereas their Hermitean conjugatesc,s were subject to periodic boundary conditions. Albeit in the particular case of J/ψφ scattering the twisting cannot be used for the extraction of the optical potential, one could not exclude a possibility that this kind of twisting could be applied in other systems for this purpose. For this reason, we consider in detail this case of (unconventional) twisting in App. B.  Figure 7: Comparison of different scenarios with respect to the number of poles reconstructed below the primary threshold. The curves were produced by using the parameters of the perfect fit from the Sect. 2, but neglecting a certain number of poles below the KK threshold.

Analysis of synthetic data
In the following, we shall reconstruct the optical potential from a synthetic lattice data set generated by the chiral unitary approach of Ref. [32]. Twisted boundary conditions are applied as described above, and the box size is taken to be L = 5M −1 π . In the first stage of our analysis we have observed that more than 100 energy eigenvalues are required to extract the potential in the considered, and quite wide, energy range from E = 2M K to E = 1.7 GeV. To produce the synthetic data, we consider the following set of six different twisting angles For these values, Z θ 00 (1; q 2 KK ) has the smallest number of poles. This requirement is important, when the energy eigenvalues are measured with a finite accuracy. Then, in proximity of its poles, the function Z θ 00 (1; q 2 KK ) will exhibit a very large uncertainty. Solving Eq. (12) with Z 00 (1; q 2 KK ) replaced by Z θ 00 (1; q 2 KK ) for each of the aforementioned angles we were able to extract 186 energy eigenvalues above and 3 below the KK threshold. Further, in any realistic lattice simulation, the eigenvalues will be known only up to a finite precision. To check the feasibility of the proposed method, it is important to account for this error, ∆E, and to see how this uncertainty 3 is reflected in the final result as studied with re-sampling techniques in the following. Therefore, we start from a sufficiently large number (∼ 1000) of re-sampled lattice data sets, normally distributed around the (189) synthetic eigenvalues with a standard deviation of ∆E. An example of 75 synthetic lattice data sets with ∆E = 1 MeV is presented in Fig. 6.
In the next step, we determine the parameters of Eq. (15) for each of these sets. Prior to doing so, we have to clarify several questions: • Range of applicability. Below the KK threshold, the function Z θ 00 (1; q 2 KK ) does not depend on θ up to exponentially suppressed contributions. Therefore, only a limited number of energy  13) to the synthetic lattice data as described in the main text. Different curves represent fits to different sets of re-sampled synthetic lattice data corresponding to the notation of Fig. 6. The gray dashed line shows the actual amplitude W −1 L (E) to guide the eye. eigenvalues can be determined. A reliable extraction of positions and residua of all four lowest poles is not possible because the twisting cannot generate the necessary scan of W −1 L in this energy region. This means that, on the one hand, this approach does not allow one to extract the optical potential below the primary (KK) threshold. On the other hand, it is crucial to recall that, due to smearing applied in the complex energy plane, this failure will yield the wrong real and especially imaginary parts of the reconstructedŴ −1 (E). This is demonstrated in Fig. 7, which was produced by using the test parameters of the perfect fit from the last section, but neglecting a certain number of poles below the KK threshold. It is seen that the imaginary part ofŴ −1 at the primary threshold deviates by about 50%, if no poles are considered below this threshold. However, already the inclusion of the first pole below the primary threshold improves the description drastically. Therefore, all poles above as well as the one below the primary threshold should be considered in the fit to the (synthetic) lattice data. Note also that if the secondary channels open above the primary channel, none of these complications arise.
• Number of poles -starting values. We found that, for sufficiently many eigenvalues and ∆E of the order of several MeV, the number of poles above the primary threshold to be fitted can be determined, searching for a rapid sign change of Z θ 00 (1; q 2 KK ). The corresponding energy eigenvalues serve us as limits on the pole positions, while the residua are allowed to vary freely.
• Highest order of the polynomial part. In principle, the order of the polynomial part of Eq. (15) is not restricted a priori. We have tested explicitly that adding terms of fourth  or fifth order in energy to the fit function yields only a small change of the reconstructed potential. This part may be further formalized by conducting combined χ 2 -and F -tests on the χ 2 defined below.
• Definition of χ 2 . The uncertainty of the (synthetic) lattice data is given by ∆E only. Therefore, a proper definition of χ 2 d.o.f. should account for the difference between the measured {E i |i = 1, ..., N } and fitted eigenvalues {E f i |i = 1, ..., N } compared to ∆E for all N data points. The E f i eigenvalues are defined as the solutions of the following equation which is technically very intricate. The problem of finding such solutions can be circumvented by expanding both sides of the latter equation in powers of (E f i − E i ) around E i for each i = 1, ..., N . Up to next-to-leading order in this expansion, the correct quantity to minimize  reads where n is the number of free parameters and θ i is the twisting angle corresponding to the energy eigenvalue E i . Note that the χ 2 in Eq. (25) differs from the usual definition by a correction factor in the denominator, given by the difference of the derivatives of the Lüscher and the fit function.
For every member of the data sets, each consisting of 188 energy eigenvalues (186 above and 2 below threshold), we perform a fit, minimizing χ 2 d.o.f. given in Eq. (25). Note that the two closest energy eigenvalues below the KK threshold, which are included in the fit, are assigned a weight factor of 6, because they are measured at every value of θ of Eq. (23) and do not depend on its value up to exponentially suppressed contributions. Further, the number of free parameters n is set to 32, consisting of 13(1) pole positions and 13(1) residua above(below) KK threshold, as well as 4 parameters in the polynomial part. The minimization is performed by using the Minuit2 (5.34.14) library from Ref. [43]. A representative subset (75 synthetic lattice data sets) of the results of the fits is shown in Fig. 8. It is seen that the data are described fairly well by all fits in a large energy region starting above the KK threshold. At and below this threshold, there is much larger spread of the fit curves describing the data. Especially the pole at ∼ 0.9 GeV is not fixed very precisely which is quite natural, keeping in mind the small number of synthetic data points in this energy region.
For each of the above fits we proceed as described in Sect. 2. First, the functionŴ −1 L (E) is evaluated at the complex energies. Second, using the Gaussian smearing as well as the parametric method discussed in Sect. 2.3, the real and imaginary parts of the potential are smoothened. The penalty factor λ = 0.28 (see App. A) and the smearing radius r = 0.2 GeV are used in these methods, respectively. Finally, for every energy, we calculate the average and the standard deviation σ. The result of this procedure is presented in Fig. 9. It is seen that both smearing methods yield very similar results. Overall, the exact solution (the dashed line) in the considered energy region lies within 1 or 2 sigma bands around the reconstructed potential. The error band appears to be comfortably narrow, but becomes broader around the KK threshold and E max = 1.7 GeV. This effect is a natural consequence of the missing information outside the energy region, which influences the prediction within the energy region via smearing during the intermediate steps of the potential reconstruction.
Furthermore, we have repeated the whole procedure of synthetic lattice data generation, fitting and recovering of the optical potential for higher uncertainty of the energy eigenvalues, ∆E = 2 MeV and ∆E = 3 MeV. The results are presented in Fig. 10 and show that the error bars grow roughly linearly with ∆E and that the real part of the reconstructed amplitude remains quite stable. The imaginary part is more sensitive to the value of ∆E. Further, at even higher values of ∆E ∼ 10 MeV, the fit is not reliable anymore and the imaginary part becomes very small.

Conclusions
i) In the present paper, we formulate a framework for the extraction of the complex-valued optical potential, which describes hadron-hadron scattering in the presence of the inelastic channels, from the energy spectrum of lattice QCD. An optical potential, defined in the present article, is obtained by using causal prescription E → E + iε for the continuation into the complex energy plane. It converges to the "true" optical potential in the limit L → ∞, ε → 0. A demonstration of the effectiveness of the method has been carried out by utilizing synthetic data.
ii) The approach requires the precise measurement of the whole tower of the energy levels in a given interval. The optical potential is then obtained through averaging over all these levels.
iii) Moreover, the availability of this approach critically depends on our ability to take the lattice data at neighboring energies without changing the interaction parameters in the secondary channels. This can be achieved, e.g., by using (partially) twisted boundary conditions that affects the pripary channel only. In the paper, we consider several systems, where the method can be applied. It is remarkable that some candidates for the QCD exotica are also among these systems.
We would like to emphasize that the use of twisted boundary conditions is only a tool, which is used to perform a continuous energy scan of a certain interval. Whatever method is used to measure the dependence of the pseudophase on energy (all other parameters fixed), our approach, based on the analytic continuation into the complex plane, could be immediately applied.
iv) The approach could be most useful to analyze systems, in which the inelastic channels contain three or more particles. Whereas direct methods based on the use of multi-particle scattering equations in a finite volume will be necessarily cumbersome and hard to use, nothing changes, if our approach is applied. The reason for this is that, in case of an intermediate state with any number of particles, the single poles are the only singularities in any Green's function in a finite volume.
A Penalty factor for a realistic set of the synthetic data In Sect. 2.3, where the parametric method for the smearing was introduced, we assumed that the quantity W −1 L can be measured with no uncertainties and at all energies from E min = M π + M η to E max = 1.7 GeV. We now turn to a more realistic case, studied in the numerical simulation in Sect. 3. For this, the search forλ opt is adapted to the interval from E min = 2M K to E max = 1.7 GeV, using severalŴ −1 L 's from the Monte-Carlo ensemble (see description there). Fig. 11 shows the χ 2 behavior for the training set, the test/validation set χ 2 V , and the true χ 2 t for one arbitrarily chosen fit of the Monte-Carlo ensemble of differentŴ −1 L 's. Both variants of the penalty, P 1 and P 2 from Eqs. (18,19), are shown in the left and right panels, respectively.
As Fig. 11 shows, the minima of χ 2 V (red triangles) are even more pronounced than in the previously discussed, idealized case, leading toλ opt = 0.34 for P 1 andλ opt = 0.28 for P 2 (minima of the curves shown with red triangles). The minima of the χ 2 t occur almost at the same respective values of λ (blue filled circles) which again demonstrates the applicability of the method. For both penalties, we also show the moduli of the Fourier coefficients |c n |, n = 1, . . . , 4 in the respective right panels, where Here, the infinite-volume quantityŴ −1 (E) implicitly depends on λ. These coefficients indicate the weight of the available frequencies to built up the optical potential over a finite energy range. As long as the potential is smooth, we expect the lowest |c n | to dominate. For decreasing values of λ, eventually a point is reached at which the oscillations will become noticeable and coefficients |c n | with larger n will become more relevant. Indeed, the figure shows that, close to the respectivê λ opt 's, the coefficients |c 2 | to |c 4 | exhibit a very pronounced rise. In all simulations, which were • n=3 • n=3  Figure 11: Determination ofλ opt for a realistic numerical simulation. Notation as in Fig. 4. Left two graphs: Using the penalization P 1 of Eq. (18). Right two graphs: Using the penalization P 2 of Eq. (19). For each case, the χ 2 (training set), χ 2 V (test/validation set) and χ 2 t (true χ 2 ) are displayed. Additionally, the moduli of the Fourier coefficients |c n |, n = 1, . . . , 4 are shown for each case. For further explanations, see text.
carried out, we have observed this behavior. This suggests that the λ-dependence of the Fourier coefficients can be used as a tool to cross-check the results from cross validation.
As a final remark, the value ofλ opt itself carries uncertainty that can be estimated by k-fold cross validation [39,40]. Using this uncertainty, the simplest model is in principle obtained by the 1-σ rule, i.e., the maximal λ compatible with the uncertainty ofλ opt [39,40]. For the numerical simulations, we simply choose one value ofλ opt = 0.28 for the penalty P 2 , because uncertainties are dominated by the statistics of the lattice measurements. As mentioned above, the valueλ opt = 0.28 corresponds to one randomly chosen fit from the Monte-Carlo ensemble, but we have made sure that this value is representative.

B Partial twisting
In this section, we would like to examine in detail the unconventional twisting prescription, which was introduced in Ref. [42], in the context of studying J/ψφ scattering from Y (4140) decays. We remind the reader that, within this prescription, only quark fields are twisted, whereas the antiquark fields are subject to the periodic boundary conditions. One could ask whether such a prescription is rigorously justified.
We address this problem by using the same methods as in Ref. [24]. In order to simplify things, we restrict ourselves to the case of elastic J/ψφ scattering and neglect the coupling to the inelastic channels. In order to treat the partial twisting, we introduce valence (v), sea (s) and ghost (g) quarks for each quark flavor, subject to twisting. Only valence and ghost quarks are twisted, whereas the sea quarks are not. In total, 9 different J/ψφ states are possible 1) (c vcv ) (s vsv ) 2) (c vcv ) (s sss ) 3) (c vcv ) (s gsg ) 4) (c scs ) (s vsv ) 5) (c scs ) (s sss ) 6) (c scs ) (s gsg ) 7) (c gcg ) (s vsv ) 8) (c gcg ) (s sss ) 9) (c gcg ) (s gsg ) .
The free Green's function is given by a a diagonal 9 × 9 matrix. Taking into account the sign cc cc ss ss x y c y s b Figure 12: The fully connected piece (x), the partially connected pieces (y c and y s ) and the fully disconnected piece (b) of the J/ψφ scattering amplitude. where and G ± = 0 due to the conservation of the total momentum, if θ is not equal to a multiple of 2π. As seen from Eq. (32), the finite-volume scattering matrix at θ = 0 contains two towers of poles, determined by the equations 1 = 0 and 4 = 0, respectively, where the former depends on the parameter θ and the latter does not. The explicit expression of the scattering matrix element in the valence sector is given by It is also clear that the θ-dependent singularities are determined by the fully connected part of the scattering amplitude, whereas the θ-independent part contains the full amplitude. Consequently, the approach of Ref. [42] can be safely used if and only if the contribution of the disconnected diagrams is much smaller than the connected one (in fact, this was mentioned already in Ref. [42]). In this case, i.e., whenb =c =s = 0, the double and triple poles in Eq. (34) vanish and one arrives at the expression that was expected from the beginning For the particular problem, considered here, one expects that the disconnected contributions will be strongly suppressed, according to the OZI rule. Consequently, the justification of the method, proposed in Ref. [42], heavily rests on the effectiveness of the OZI suppression.