Hadron-Hadron Interactions from $N_f=2+1+1$ Lattice QCD: $I=3/2$ $\pi K$ Scattering Length

In this paper we report on results for the s-wave scattering length of the $\pi$-$K$ system in the $I=3/2$ channel from $N_f=2+1+1$ Lattice QCD. The calculation is based on gauge configurations generated by the European Twisted Mass Collaboration with pion masses ranging from about $230$ to $450\,\text{MeV}$ at three values of the lattice spacing. Our main result reads $M_{\pi}\,a_0^{3/2,\text{phys}} = -0.059(2)$. Using chiral perturbation theory we are also able to estimate $M_{\pi}\,a_0^{1/2,\text{phys}} = 0.163(3)$. The error includes statistical and systematic uncertainties, and for the latter in particular errors from the chiral and continuum extrapolations.

= 0.163 (3). The error includes statistical and systematic uncertainties, and for the latter in particular errors from the chiral and continuum extrapolations.

I. INTRODUCTION
For understanding the strong interaction sector of the standard model (SM) it is not sufficient to compute masses of stable particles. Gaining insight into interactions of two or more hadrons and resonances is a must. Due to the non-perturbative nature of low energy quantum chromodynamics (QCD), computations of interaction properties from lattice QCD are highly desirable. While ultimately the phase shift in a given partial wave is to be computed, also the scattering length is in many cases a useful quantity, in particular when the two-particle interaction is weak.
Due to the importance of chiral symmetry in QCD the investigation of systems with two pseudoscalar mesons is of particular interest. Here, chiral perturbation theory (ChPT) is able to provide a description of the pion mass dependence. And any non-perturbative computation in turn allows to check this dependence. Naturally, ChPT works best for two pion systems, while convergence is unclear for pion-kaon or two kaon systems. The two pion system is studied well experimentally, also in the different isospin channels. However, as soon as one or both pions are replaced by kaons, experimental results become sparse. On the other hand, this gap starts to be filled by lattice QCD calculations. For the pion-kaon system with isospin I = 3/2 there are by now a few lattice results available focusing on the scattering length [1][2][3][4][5]. The most recent computation in Ref. [4] uses one lattice at physical pion and kaon masses and lattice spacing a ≈ 0.114 fm. For the sea and valence sector they use N f = 2+1 Möbius domain wall fermions and an Iwasaki gauge action. In Ref. [2] a systematic study of the elastic scattering lengths for the light pseudoscalar mesons was carried out with N f = 2 + 1 order O (a)-improved Wilson quarks at pion masses ranging from 170 MeV to 710 MeV and a lattice spacing a ≈ 0.09 fm. Furthermore Refs. [1,3] use N f = 2 + 1 flavors on the MILC configurations with a rooted staggered sea quark action. Whereas Ref. [3] calculates the scattering length at a lattice spacing a ≈ 0.15 fm, a slightly smaller lattice spacing a ≈ 0.125 fm has been used in Ref. [1]. The pion masses in Ref. [1] range from 290 MeV to 600 MeV using domain wall valence quarks with a chiral extrapolation done in mixed-action chiral perturbation theory (MAChPT) [6,7]. The range of pion masses, 330 MeV to 466 MeV, for the Asqtad improved staggered fermions of Ref. [3] is a bit smaller compared to Ref. [1]. In Ref. [5] the phaseshifts and scattering lengths for π-K-scattering in I = 3/2 and I = 1/2 in the s-wave and the p-wave has been determined. The gauge action is a N f = 2 tree level improved Wilson-Clover action. The authors include the strange quark as a valence quark only which then corresponds to pion and kaon masses of M π = 266 MeV and M K = 522 MeV, respectively.
In this paper we are going to present results for the s-wave scattering length of the pion-kaon system in the elastic region with isospin I = 3/2. The investigation is based on gauge configurations produced by the European Twisted Mass Collaboration (ETMC) with N f = 2 + 1 + 1 dynamical quark flavors [8]. In contrast to previous computations, we are able to perform a continuum extrapolation owing to 11 ensembles with M π ranging from 230 MeV to 450 MeV distributed over 3 different lattice spacing values. We employ in total 4 different extrapolation methods to also estimate systematic uncertainties associated with our computation.
Finally, since this paper is the fourth in a series of publications [9][10][11] concerning elastic scattering of two pions in different channels and kaon-kaon with I = 1, we are able to compare results of two pseudoscalar mesons at maximal isospin involving different amounts of strangeness. The leading order ChPT predictions for the dependence on the reduced mass ensemble β aµ aµ σ aµ δ (L/a) 3 Table I: The gauge ensembles used in this study. For the labeling of the ensembles we adopted the notation in Ref. [13]. In addition to the relevant input parameters we give the lattice volume and the number of evaluated configurations, N conf .
divided by the relevant decay constant are identical for the three systems and differences appear only at NLO. This paper is organized as follows: We first introduce the lattice details of our calculation. After the discussion of the analysis methods we present the main result, followed by a detailed discussion of the analysis details. We close with a discussion and summary. Technical details can be found in the appendix.

A. Action
The lattice details for the investigation presented here are very similar to the ones we used to study the kaon-kaon scattering length [11]. We use N f = 2 + 1 + 1 flavor lattice QCD ensembles generated by the ETM Collaboration, for which details can be found in Refs. [8,12,13]. The parameters relevant for this paper are compiled in Table I: we give for each ensemble the inverse gauge coupling β = 6/g 2 0 , the bare quark mass parameters µ , µ σ and µ δ , the lattice volume and the number of configurations on which we estimated the relevant quantities.
The ensembles were generated using the Iwasaki gauge action and employ the N f = 2+1+1 twisted mass fermion action [14][15][16]. For orientation, the β-values 1.90, 1.95 and 2.10 correspond to lattice spacing values of a ∼ 0.089 fm, 0.082 fm and 0.062 fm, respectively, see also Table III. The ensembles were generated at so-called maximal twist, which guarantees automatic order O (a) improvement for almost all physical quantities [14]. The renormalized light quark mass m is directly proportional to the light twisted quark mass via with Z S the non-singlet scalar renormalization constant.
As noted in Refs. [13,17], the renormalized sea strange quark masses across the "A", "B" and "D" ensembles vary by up to about 20% and in a few cases differ from the physical strange quark mass to the same extent. For D30.48 and D45.32sc at the finest lattice spacing, the sea strange quark mass on the former ensemble overshoots the physical strange quark mass while it is consistent on the latter ensemble. In order to correct for these mistunings and to avoid the complicated flavor-parity mixing in the unitary non-degenerate strange-charm sector [8], we adopt a mixed action ansatz with so-called Osterwalder-Seiler (OS) [16] valence quarks, while keeping order O (a) improvement intact. We denote the OS bare strange quark parameter with µ s . It is related to the renormalized strange quark mass by For each ensemble we investigate three values of µ s which are compiled in Table II. More details on the mixed action approach can be found in Ref. [11]. As a smearing and contraction scheme we employ the stochastic Laplacian-Heaviside approach, described in Ref. [18]. Details of our parameter choices can be found in Refs. [9,11].

B. Lattice Operators and Correlation Functions
For reasons which will become clear later we need to estimate the masses of the pion, the kaon and the η meson on our ensembles. The masses for the pion and kaon are obtained from the large Euclidean time dependence of two point functions of the form where X ∈ {π, K}. The operators for the charged pion and kaon projected to zero momen- with For the η (and η ) meson we use the two operators From these we build a two-by-two correlator matrix by taking the disconnected diagrams into account. The η (principal) correlator is determined by solving a generalized eigenvalue problem as described in detail in Ref. [19]. A complete discussion of the analysis of the η (and η ) meson is beyond the scope of this paper and the full analysis will be presented in a future publication [20]. In addition to the aforementioned meson masses, we also need to estimate the energy E πK of the interacting pion-kaon two particle system. For the case of maximal isospin, i.e. I = 3/2, the corresponding two particle operator reads It is used to construct the two-particle correlation function E πK can then be determined from the large Euclidean time dependence of C πK .

III. ANALYSIS METHODS
We focus in this work on pion-kaon scattering in the elastic region. For small enough squared scattering momentum p 2 one can perform the effective range expansion for partial wave with phase shift δ and scattering length a . For the pion-kaon system it is, to a very good approximation, sufficient to study the s-wave, i.e. = 0. In lattice QCD the phase shift or the scattering length can only be computed from finite volume induced energy shifts. The relevant energy shift here is given by Using again the effective range expansion, one arrives at the Lüscher formula [21] relating δE directly to the scattering length a 0 , the reduced mass of the pion-kaon system and the spatial extent of the finite volume L. The coefficients read [21] c 1 = −2.837297 , c 2 = 6.375183 .   Given δE, µ πK and L, Lüscher's formula allows one to determine the scattering length a 0 by solving Equation (14) for a 0 . In what follows, we will describe how we extract δE and the other relevant bare quantities from correlation functions. Then we will give details on our approach to inter-or extrapolate the results to physical conditions and to the continuum limit.
In order to gain some understanding of systematic uncertainties, we perform the analysis in two different ways once the bare data has been extracted. Combined chiral and continuum extrapolations are performed at fixed strange quark mass using next to leading order ChPT (NLO ChPT) and a variant thereof referred to as the Γ-method, as described in Ref. [1].

A. Physical Inputs
For the analysis presented below, we require physical inputs for the pion, the kaon and η-meson masses as well as the pion decay constant. To this end, we employ the values in the isospin symmetric limit, M π and M K , as determined in chiral perturbation theory [22] and given in Ref. [23] as For the η meson mass we use the average obtained by the Particle Data Group [24]: For the decay constant, we use the phenomenological average determined by the Particle Data Group given in Ref. [25] as As an intermediate lattice scale, we employ the Sommer parameter r 0 [26]. It was determined in Ref. [17] from the ensembles we use here to be In the parts of the analysis which require r 0 , we use parametric bootstrap samples with central value and width given in Equation (19). Where r 0 /a values enter as fit parameters, we constrain the corresponding fit parameters using Gaussian priors in the augmented χ 2 function given as

B. Energy Values from Correlation Functions
The energies of the two point correlation functions as given in Equation (4) are extracted from fits of the form to the data. While for M K and M π the signal extends up to T /2, for the η we have to face more noise . We deal with this by applying the excited state subtraction method used and described in Refs. [19,27].
In the determination of the energy shift δE, the total energy E πK of the interacting π-K system must be computed. However, in the spectral decomposition of the two-particle correlation function, unwanted time dependent contributions, so-called thermal pollutions, appear. Taking into account that our π-K correlation function is symmetric around the T /2 point, the leading contributions in the spectral decomposition can be cast into the form where only the first line corresponds to the energy level we are interested in. However, at finite T -values, the second contribution might be sizable, in particular at times close to T /2. Moreover, the thermal pollution cannot be separated easily from the signal we are interested in. We have studied two different methods, labeled E1 and E2, to extract E πK from C πK (t), where E1 has already been discussed in Ref. [28].
• E1: weighting and shifting: To render one of the polluting terms in Equation (22) time independent, the correlation function first gets weighted by a factor exp((E K − E π )t). We chose this factor, because exp(−E π T ) is significantly larger than exp(−E K T ). The resulting constant term can then be removed by the shifting procedure, which thus replaces C πK (t) by where δt is a fixed number of time slices.
Subsequently, we multiply C w πK (t) by exp(−(E K − E π )t), which (mostly) recovers the original time dependence in the contribution of interest We then extract the total energy of the π-K system, E πK , from a fit of the form, Note that in contrast to Ref. [28], where correlator matrices with various thermal pollutions are considered, we are able to takeÃ 1 as an additional fit parameter in order to account for this sub-leading term.
• E2: dividing out the pollution: To improve on method E1, we assume that the decomposition given in Equation (22) allows one to neglect any further thermal pollutions. This leads to dividing out the time dependent part explicitly. With we then proceed to calculate from which E πK can be extracted through a fit of the form We remark that for both methods E1 and E2 the energies E π and E K , i.e. M π and M K for zero momentum, are required as an input. They are determined from the corresponding two-point correlation functions. For the error analysis bootstrap samples are used to fully preserve all correlations. After solving Equation (14) for a 0 up to order O (L −5 ) on every ensemble for each strange quark mass of Table II, we have three parameters in which we want to extra-or interpolate: the lattice spacing a, the strange quark mass m s and the light quark mass m . To evaluate a 0 at the physical point we follow a two step procedure. We first fix the strange quark mass to its physical value and do a combined chiral and continuum extrapolation afterwards. For the extrapolation in the light quark mass, we use the continuum ChPT expression in terms of meson masses and decay constants, as detailed in Refs. [3,7,29].

C. Fixing the strange quark mass
In order to fix the strange quark mass we adopt the following procedure: we match the quantity which is proportional to the strange quark mass at the leading order of ChPT, to its physical value (M phys using our determinations of M 2 K at three valence strange quark masses on a per-ensemble basis. For each ensemble, we then interpolate all valence strange quark mass dependent observables, i.e. µ πK a 3/2 0 , M K , M η and µ πK , in M 2 s to this reference value.

D. Chiral extrapolation
With the strange quark mass fixed, the extrapolation to the physical point can be carried out using ChPT. The first NLO calculation of the scattering amplitude and scattering lengths was done in Ref. [30]. From the formulae for the isospin even (odd) scattering lengths a + (a − ) in Ref. [29] the NLO ChPT formulae for µ πK a I 0 , I ∈ { 1/2, 3/2 }, can be derived as sketched in Appendix A, giving Equation (33) depends on the masses of the pion and the kaon, their reduced mass as defined in Equation (15), the η mass and the pion decay constant. In addition, the equation depends on the low energy constants (LECs) L 5 and L πK while χ 3/2 NLO is a known function, see Appendix A 2.
We express Equation (33) in terms of the meson masses and decay constants as they are determined on the lattice, which has the benefit that their ratios can be determined with high statistical precision without the need for explicit factors of the lattice scale. Formally we fix the scale-dependent LECs at the renormalization scale Λ χ = f Automatic order O (a) improvement of Wilson twisted mass fermions at maximal twist guarantees that the leading lattice artifacts are of order O (a 2 ) or better. For instance, for the I = 2 ππ s-wave scattering length, discretization effects start only at order O (a 2 M 2 π ) [31]. A corresponding theoretical result for πK is missing so far. However, our numerical data suggest that also for πK lattice artefacts are very small. Still, we include a term c · f (a 2 ) accounting for possible discretization effects, with fit parameter c and f (a 2 ) either equal to a 2 /r 2 0 or to a 2 M 2 X , with M 2 X one of the masses or mass combinations M 2 π , M 2 K , M 2 K + 0.5M 2 π , µ 2 πK . In the following analysis we will include the term c · f (a 2 ) into our fit for every choice of f (a 2 ) and thus investigate a possible dependence of our data on the lattice spacing.
To summarize, our fit parameters are the LECs L 5 and L πK , and c, where L πK is the combination of renormalized LECs Let us mention already here that the fits to the data described in the next section turn out to be not sensitive to L 5 . Therefore, we include it as a prior in the fit with the value taken from Ref. [23]. In slight abuse of language we will denote this extrapolation method as NLO ChPT.
E. Extrapolations Using the Γ Method at Fixed m s Next, we describe an alternative way to extrapolate our data. For this we fix the strange quark mass as described in the previous subsection and interpolate the other quantities. in the continuum limit at the physical point obtained with the different methods used in this paper. The fit ranges decrease with increasing index as described in Table IV. The inner error band represents the statistical error only, while the outer error band represents the statistical and systematic errors added in quadrature.
Using the interpolated data we can compute the following quantity [1] Γ with next to leading order ChPT functions χ ± NLO given in Appendix A 2. Using the ChPT expressions for the isospin even/odd scattering lengths a + and a − , Γ can also be expressed as [1] Γ For the derivation see Appendix A. Hence, a linear fit of Equation (36) to the data obtained via Equation (35) allows one to determine the LECs L 5 and L πK . Given L 5 and L πK from the fit one can compute µ πK a 3/2 0 at the physical point using Equation (33). Again, it turns out we are not sensitive to L 5 in our fits. Therefore, we use a prior as discussed before. This extrapolation method we denote as Γ method.

IV. RESULTS
In this section we present our main result for µ πK a 3/2 0 in the continuum limit and extrapolated to the physical point.
We use two thermal state pollution removal methods, E1 and E2, for E πK . Next we employ the two (related) ChPT extrapolations, Γ-method and NLO ChPT, as discussed  before. For each of the two ChPT extrapolation methods we use three fit ranges as compiled in Table IV. Hence, we have twelve estimates for each quantity at the physical point in the continuum limit available, which we use to estimate systematic uncertainties. We remark that the fit for the Γ method is in terms of M K /M π and for NLO ChPT in terms of µ πK /f π . Thus, we vary the fit range at the lower end for the Γ method and at the upper end for NLO ChPT.
For µ πK a 3/2 0 the twelve estimates are shown in Figure 1. The final result is obtained as the weighted average over all of these, as shown in the figure as the horizontal bold line. The weight is computed according to with p the p-value of the corresponding ChPT fit and ∆ the statistical uncertainty obtained from the fit. The statistical uncertainty of the final results is determined from the bootstrap procedure. For µ πK a 3/2 0 this is shown in Figure 1 as the inner error band. In addition, we determine three systematic uncertainties: The first is obtained from the difference between using only E1 or only E2 results. The second from the difference between using only the Γ method or only NLO ChPT. Finally, we use the maximal difference of the weighted average to the twelve estimates as a systematic uncertainty coming from the choice of fit ranges.
The results of all twelve fits can be found in Table VII for the Γ method and Table VIII for NLO ChPT fits. The fit range indices used in Figure 1 are resolved in Table IV. With this procedure and all errors added in quadrature we quote This translates to as our final results. The error budget is compiled in Table V. While the dominating contribution to the error for both µ πK a 3/2 0 and L πK summed in quadrature is coming from the fit range and the statistical uncertainty, also the choice of the thermal state removal method contributes significantly. The contribution from the different chiral extrapolation methods is negligible. If the errors were added (not in quadrature), the total error would become a factor ∼ 1.7 larger.
We remark that these results have been obtained with L 5 as an input, because the fits are not sufficiently sensitive to determine L 5 directly. We use the most recent determination from a N f = 2 + 1 + 1 lattice calculation by HPQCD [23], which is extrapolated to the continuum limit. At our renormalization scale it reads  and L πK .

A. Error Analysis, Thermal Pollution and Choice of Fit Ranges
The error analysis is performed using the stationary blocked bootstrap procedure [32]. In order to determine an appropriate average block length, we compute the integrated autocorrelation time τ int for the correlation functions C X (t) at all source-sink separations, with X being π, K, η or πK. In the case of πK, C X (t) is of course first suitably transformed for the extraction of the interaction energy as discussed in section III B. The computation of τ int is detailed in Ref. [33]. The average block length is then chosen to be the ceiling of the maximum integrated autocorrelation time observed over all correlation functions at all source-sink separations b = max on a per-ensemble basis. We have confirmed explicitly that this method produces a block length at which the estimate of the statistical error plateaus and are thus confident that we properly take into account the effect of autocorrelations on our quoted statistical errors.
Using the so-determined block length on a per-ensemble basis, we generate N = 1500 samples from which we estimate statistical errors throughout our analysis. As discussed in Section III B, we employ methods E1 and E2 to remove unwanted thermal pollutions from the πK two particle correlation function. Both methods allow us to describe the data rather well, but the choice of best fit range depends on the method used to remove the thermal pollution. This in turn affects the value of the extracted E πK and, subsequently, the value of µ πK a 3/2 0 obtained from the energy difference.
To demonstrate the quality of our fits, we look at the ratio where C [E1,E2] are defined in Equations (24) and (29), respectively, and the fit functions f [E1,E2] (t) are given in Equations (25) and (30), respectively. The ratio is shown in Figure 2 for the two ensembles A40.24 and A40.32 for two fit ranges for which both methods describe the data well. The choice of the fit ranges to determine energy levels is always difficult. In the past, we have used many fit ranges and weighted them according to their fit qualities [9,11]. However, this procedure relies on properly estimated variance-covariance matrices, which is notoriously difficult. For the pion-kaon correlation functions needed in this paper we have observed several cases where the fit including the variance-covariance matrix did not properly describe the data after visual inspection. Therefore, we use fits here assuming independent data points with the correlation still taken into account by the bootstrap procedure.
As a consequence, we cannot apply the weighting procedure used in Refs. [9,11] any longer and have to choose fit ranges. The procedure is as follows: Due to exponential error growth of C πK we fix t f = T /2−4a and vary t i , beginning from a region where excited states do not contribute significantly anymore. From these fits we choose one fit range where the ratio of Equation (41) is best compatible with 1. The statistical error is calculated from the bootstrap samples as discussed before. We then estimate the systematic uncertainty from the remaining fit ranges. To this end we determine the difference of the mean value to the upper and lower bound of values for E πK . This procedure results in an asymmetric estimate of the systematic uncertainty of E πK . The results are compiled in Table XI. Since C K and C π do not suffer from exponential error growth at late times we set t f = T /2. Table VI gives an overview of our chosen values of t i and t f for all ensembles.
It is always difficult to include systematic uncertainties in the analysis chain. Since we see systematic uncertainties on extracted energies on the same level as the statistical one, we adopt the following procedure to include this uncertainty: Because µ πK a  on each ensemble after the data have been interpolated in the strange quark mass. To this end we define a scaling factor s via the standard error ∆X and the average of the systematic uncertainties Q X over the 3 strange quark masses for each ensemble: where the average Q X is the simple mean over the six systematic errors. In order to extract δE we first determine M K and M π from fitting Equation (21) to our data for C π (t) and C K (t). We then calculate the reduced mass µ πK via Equation (15) for all combinations of fit ranges. M π , M K and µ πK are listed in Table X. The two methods E1 and E2 give us two estimates of E πK as outlined in Section III B, from which we determine δE and hence the scattering length using Equation (14). The values for E πK and δE are collected in Tables XI and XII. We introduce factors K FSE X for X ∈ {M K , M π , f π } to correct our lattice data for finite size effects. They have already been calculated in Ref. [17] and are listed in Table IX. We apply these factors for e.g. M π like We correct every quantity of the set named above and drop the asterisk in what follows to improve legibility. For M η statistical uncertainties are too big to resolve finite volume effects, see also Ref. [34].
For the two methods E1 and E2 we solve Equation (14) for a 0 up to order O(L −5 ) numerically. The values for a 0 and its product with the reduced mass, µ πK a 3/2 0 are collected in Table XIII. Since the finite size behavior of the scattering length is unknown, we do not apply finite size corrections to the reduced mass appearing in µ πK a 3/2 0 , either.

C. Strange quark mass fixing
Before extrapolating to the continuum limit and the physical point, we interpolate all data to reference strange quark masses as discussed before. The data for the three strange quark masses are strongly correlated because the same stochastic light perambulators were used for all light-strange observables. As a consequence, the variance-covariance matrix was sometimes not sufficiently well estimated such that its inverse was unreliable. As a result, we resort to performing these fits using uncorrelated χ 2 which results in best fit parameters which describe the data much better. It should be noted that all statistical covariance is still fully taken into account by the bootstrap procedure and our final statistical errors on all fit parameters are correctly estimated.
As an example, we show in Figure 3 the interpolation of µ πK a 3/2 0 in M 2 K − 0.5M 2 π for ensemble B55.32 comparing methods E1 and E2. The large uncertainty in the interpolation variable stems mainly from the uncertainty in the scaling quantity r 0 . Furthermore the errors of the three data points are highly correlated. The interpolation to the reference point is shown as a red diamond. In general, the strange quark mass dependence of µ πK a 3/2 0 is mild and stems mainly from the reduced mass µ πK . The values thus determined are compiled in Appendix B 1. They serve as input data for the subsequent chiral extrapolations.

D. Chiral and Continuum Extrapolations
Having interpolated all our lattice data to a fixed reference strange quark mass corresponding to the physical strange quark mass at leading order, we will describe below possible systematic errors in our chiral and continuum extrapolations.

Γ-Method
In this section we present results employing the determination of L πK using the linear fit introduced in Section III E. Figure 4 shows the chiral extrapolations in terms of M K /M π for pollution removal E1 and E2. Since we work at fixed strange quark mass, the light quark mass decreases from left to right in the figure. In order to check how our extraction of L πK is affected by the range of included pion masses we employ three different fit ranges In Table VII we compile the results for the fits corresponding to the data points of Figure 4 for all 3 fit ranges. As the fit range is restricted to our lightest ensembles, the value extracted for L πK tends up, while the absolute value of the extracted µ πK a 3/2 0 decreases. It is worth noting that this behavior is only observed for the pollution removal E2, whereas for E1 the values for L πK and µ πK a   Figure 4 correspond to the largest fit range in the table.

Chiral Perturbation Theory at NLO
We first fit the continuum ChPT formula Equation (33) with c = 0, i.e. assuming no lattice artefacts. In the right column of plots in figure Figure 5, we show the lattice data for µ πK a 3/2 0 interpolated to the reference strange quark mass as a function of µ πK /f π for the two thermal pollution removal methods E1 and E2. The solid line corresponds to the leading order, parameter free ChPT prediction. Plotting our best fit curve with NLO ChPT together with the data is difficult, because µ πK a 3/2 0 depends on meson masses and f π besides µ πK /f π . Therefore, in order to demonstrate that the fit is able to describe our data, we indicate the relative deviation δ r (µ πK a 3/2 0 ) between the fitted points and the original data in Figures 5b and 5d. The indicated error bars are statistical only and it is clear that within these uncertainties, our data is reasonably well described by the fit. As in Section V D 1, to investigate the validity of Equation (33) across our entire range of pion masses, we studied three different fit intervals for µ πK /f π , namely µ πK f π ∈ { (0; 1.35), (0; 1.41), (0; 1.60) } , where now the first range corresponds to only using our lightest pion masses and the third range includes all of our ensembles. The resulting trend in the extracted values of L πK and µ πK a 3/2 0 is shown in Table VIII. Just as in the study of the Γ-method, including heavier pion masses leads to smaller values of L πK and correspondingly smaller values of µ πK a 3/2 0 . To investigate possible discretization effects we now allow c = 0 in Equation (33). Fitting Equation (33) for the different choices of f (a 2 ), we are neither able to obtain a statistically significant result for the fit parameter c, nor do we see significant differences in the extracted values of L πK and µ πK a 3/2 0 . We conclude that we are not able to resolve lattice artifacts in this quantity given our statistical uncertainties.

VI. DISCUSSION
Let us first discuss the main systematics of our computation: In contrast to the pionpion or kaon-kaon systems, there are time dependent thermal pollutions in the correlation  functions relevant for the extraction of the pion-kaon s-wave scattering length. This very fact turns out to represent one of the major systematic uncertainties in the present computation.
We have investigated two methods to remove the leading thermal pollutions, denoted as E1 and E2. With both we are able to describe the data for the correlation functions. However, there is uncertainty left, because we remove only the leading pollution and the removal procedure requires input estimated from other two point functions. Thus, we eventually decided to use both methods E1 and E2 and include the differences in the systematic   Figure 5 correspond to the largest fit range in the table. uncertainty. Secondly, we perform a mixed action simulation for the strange quark. We use this to correct for small mis-tuning in the sea strange quark mass value used for the gauge configuration generation. This leads -at least in principle -to a small mismatch in the renormalization condition used for the continuum extrapolation. We cannot resolve the corresponding effect on our results quantitatively given our statistical uncertainties. But, since we study quantities which mainly depend on the valence quark properties we expect them to be small.
Thirdly, in the ChPT determination of L πK the remaining LEC, L 5 , entered as a prior to numerically stabilize our fits. The HPQCD value of Ref. [35] stems from an independent lattice simulation, but is extrapolated to the continuum limit. In Ref. [35] L 5 is given at scale M η , which we translated to our renormalization scale given by the pion decay constant.
In addition, the extrapolation from our data to the physical point is quite long. Here, a computation directly with physical pion mass would improve our confidence in the result. The final error on our determination is only as small as it is due to the highly constraining ChPT description of µ πK a 3/2 0 . Finally, although we are not able to resolve lattice artefacts in our determination of µ πK a 3/2 0 our statistical errors and limited set of gauge ensembles especially at the finest lattice spacing might make us unable to resolve possible lattice artefacts.

VII. SUMMARY
In this paper we have presented a first lattice computation of the pion-kaon s-wave scattering length for isospin I = 3/2 extrapolated to the continuum limit. By varying our methodology we estimate the systematic uncertainties in our results. Our errors cover statistical uncertainties, continuum and chiral extrapolations as well as the removal of thermal pollutions.
In the left panel of Figure 6 we compare the results presented in this paper with previous lattice determinations. The inner (darker) error bars show the purely statistical errors whereas the outer (lighter) ones correspond to the statistical and systematic errors added in quadrature. Even though the four other determinations lack the extrapolation to the continuum limit, overall agreement within errors is observed. However, concerning the final uncertainty, our determination improves significantly on the previous determinations by controlling more sources of uncertainty. from various lattice computations [1][2][3][4][5]. The unfilled point denotes the LO extrapolation to the physical point using the data of Ref. [5]. Right: s-wave scattering lengths for the pion-pion, pion-kaon and kaon-kaon maximum isospin channels as a function of the squared reduced mass µ 2 of the system divided by f 2 π for pion-pion and pion-kaon and by f 2 K for kaon-kaon.
As mentioned in the introduction, the three two particle systems pion-pion, pion-kaon and kaon-kaon are very similar. Therefore, it is interesting to compare the data for pionpion [9], kaon-kaon [11] and pion-kaon in a single plot. This is done in the right panel of Figure 6 where we show µ · a 0 as a function of (µ/f ) 2 . Here µ is the reduced mass of the corresponding two particle system and f is the pion decay constant f π for the pion-pion and pion-kaon and f K for the kaon-kaon system. The dashed line in the right panel of Figure 6 is the leading order, parameter-free ChPT prediction all three systems share. The three symbols (and colors) represent our data for the three different systems, respectively.
It can be seen that for all three systems, the deviations from LO ChPT are small. For the pion-kaon system, a parametrization in terms of f K · f π would bring the points even closer to the LO line, while increasing the deviation of the final result from the LO estimate. For the kaon-kaon system, instead, a parametrization in terms of f π rather than f K (which is perfectly valid at this order of ChPT) would render the deviation from the LO line more severe.
It is somewhat surprising that ChPT appears to work so well for all three systems, especially for the heavier points in our simulations and even more for the kaon-kaon system, where the expansion parameter becomes large. A possible reason for this finding might be the fact that all three systems are only weakly interacting.
providing us with the data for f π and S. Simula for the estimates of the finite size corrections to M π , M K and f π . This project was funded by the DFG as a project in the Sino-German CRC110. The open source software packages tmLQCD [36], Lemon [37], QUDA [38][39][40], R [41], python [42] and SciPy [43] have been used. depends on L πK . From Equation (A1) it follows that which immediately carries over to the scattering lengths a 1/2 and a 3/2 yielding Inserting the ChPT formulae for a + and a − into Equation (A5) one obtains the ChPT formulae Equation (33). Moreover, after the insertion of a + and a − into Equation (A5) the terms containing the LECs L 5 and L πK can be collected at one side in order to arrive at Equation (36). The expression Γ in Equation (36) is then given by Equation (35) The functions ν ( ) X are given by From these isospin even/odd functions the definite isospin functions χ 3/2 NLO and χ 1/2 NLO can be derived in the same way as the scattering lengths of Equation (A5) Here the functions κ ( ) X are given by   (12)