Hadron-Hadron Interactions from $N_f=2+1+1$ Lattice QCD: isospin-1 $KK$ scattering length

TL;DR

This study computes the isospin-1 kaon-kaon scattering length using lattice QCD with $N_f=2+1+1$ flavors, providing a precise extrapolation to the physical point and continuum limit.

Contribution

It presents the first lattice QCD calculation of the $I=1$ $KK$ scattering length with $N_f=2+1+1$ flavors, including systematic error analysis and physical extrapolation.

Findings

Scattering length $a_0^{I=1} = -0.154(6)_{stat}(^{+0}_{-5})_{m_s}(^{+0}_{-2})_{Z_P}(2)_{r_f} ext{ fm}

Results are consistent with chiral perturbation theory predictions

Provides a benchmark for future lattice and experimental studies

Abstract

We present results for the interaction of two kaons at maximal isospin. The calculation is based on $N_f=2+1+1$ flavour gauge configurations generated by the European Twisted Mass Collaboration with pion masses ranging from about $230$ to $450\,\textrm{MeV}$ at three values of the lattice spacing. The elastic scattering length $a_0^{I=1}$ is calculated at several values of the bare strange and light quark masses. We find $M_K a_0 = -0.385(16)_{\textrm{stat}} (^{+0}_{-12})_{m_s}(^{+0}_{-5})_{Z_P}(4)_{r_f}$ as the result of a combined extrapolation to the continuum and to the physical point, where the first error is statistical, and the three following are systematical. This translates to $a_0 = -0.154(6)_{\textrm{stat}}(^{+0}_{-5})_{m_s} (^{+0}_{-2})_{Z_P}(2)_{r_f}\,\textrm{fm}$.

Tables

Table 1. Table 1 : The gauge ensembles used in this study. For the labelling of the ensembles we adopted the notation in Ref.  [ 12 ] . In addition to the relevant input parameters we give the lattice volume and the number of evaluated configurations, N conf subscript 𝑁 conf N_{\mathrm{conf}} .

ensemble	$\beta$	$a{\mu }_{\ell }$	$a{\mu }_{\sigma }$	$a{\mu }_{\delta }$	${\left(L/a\right)}^3\times T/a$	$N_{\mathrm{conf}}$
$A30.32$	$1.90$	$0.0030$	$0.150$	$0.190$	${32}^3\times 64$	$263$
$A40.20$	$1.90$	$0.0040$	$0.150$	$0.190$	${20}^3\times 48$	$268$
$A40.24$	$1.90$	$0.0040$	$0.150$	$0.190$	${24}^3\times 48$	$386$
$A40.32$	$1.90$	$0.0040$	$0.150$	$0.190$	${32}^3\times 64$	$244$
$A60.24$	$1.90$	$0.0060$	$0.150$	$0.190$	${24}^3\times 48$	$314$
$A80.24$	$1.90$	$0.0080$	$0.150$	$0.190$	${24}^3\times 48$	$305$
$A100.24$	$1.90$	$0.0100$	$0.150$	$0.190$	${24}^3\times 48$	$308$
$B35.32$	$1.95$	$0.0035$	$0.135$	$0.170$	${32}^3\times 64$	$235$
$B55.32$	$1.95$	$0.0055$	$0.135$	$0.170$	${32}^3\times 64$	$293$
$B85.24$	$1.95$	$0.0085$	$0.135$	$0.170$	${32}^3\times 64$	$290$
$D30.48$	$2.10$	$0.0030$	$0.120$	$0.1385$	${48}^3\times 96$	$369$
$D45.32sc$	$2.10$	$0.0045$	$0.0937$	$0.1077$	${32}^3\times 64$	$283$

Table 2. Table 2 : compilation of values for the Sommer parameter r 0 / a subscript 𝑟 0 𝑎 r_{0}/a , the lattice spacing a 𝑎 a and Z P subscript 𝑍 P Z_{\mathrm{P}} at 2 ​ GeV 2 GeV 2\ \mathrm{GeV} in the MS ¯ ¯ MS \overline{\mathrm{MS}} scheme determined with methods M1,M2 the three values of the lattice spacing. See Ref.  [ 17 ] for details.

$\beta$	$Z_P$ (M1)	$Z_P$ (M2)	$a\left[\mathrm{fm}\right]$	$r_0/a$
$1.90$	$0.529\left(07\right)$	$0.574\left(04\right)$	$0.0885\left(36\right)$	$5.31\left(8\right)$
$1.95$	$0.509\left(04\right)$	$0.546\left(02\right)$	$0.0815\left(30\right)$	$5.77\left(6\right)$
$2.10$	$0.516\left(02\right)$	$0.545\left(02\right)$	$0.0619\left(18\right)$	$7.60\left(8\right)$

Table 3. Table 3 : Values of the bare strange quark mass a ​ μ s 𝑎 subscript 𝜇 𝑠 a\mu_{s} used for the three β 𝛽 \beta -values. The lightest strange quark mass on the ensemble D30.48 is a ​ μ s = 0.0115 𝑎 subscript 𝜇 𝑠 0.0115 a\mu_{s}=0.0115 instead of a ​ μ s = 0.013 𝑎 subscript 𝜇 𝑠 0.013 a\mu_{s}=0.013 .

$\beta$	$1.90$	$1.95$	$2.10$
$a{\mu }_s$	0.0185	0.0160	0.013/0.0115
	0.0225	0.0186	0.015
	0.0246	0.0210	0.018

Table 4. Table 4 : Summary of the number of dilution vectors, N D subscript 𝑁 D N_{\mathrm{D}} , used in each index. We use a block scheme in time and an interlace scheme in eigenvector space.

${\left(L/a\right)}^3\times T/a$	$N_D\left(\text{time}\right)$	$N_D\left(\text{Dirac}\right)$	$N_D\left(\text{LapH}\right)$
${24}^3\times 48$	24	4	6
${32}^3\times 64$	32	4	4
${48}^3\times 96$	32	4	4

Table 5. Table 5 : Physical values of M s 2 superscript subscript 𝑀 𝑠 2 M_{s}^{2} for the three β 𝛽 \beta values. The stated values correspond to the continuum values of ( r 0 ​ M s phys ) 2 superscript subscript 𝑟 0 superscript subscript 𝑀 𝑠 phys 2 (r_{0}M_{s}^{\mathrm{phys}})^{2} equal to 1.33 ​ ( 7 ) 1.33 7 1.33(7) and 1.34 ​ ( 6 ) 1.34 6 1.34(6) for Z P subscript 𝑍 𝑃 Z_{P} from M1 and M2 , respectively.

$\beta$	${\left(a{M}_s^{\mathrm{phys}}\right)}^2$ (M1)	${\left(a{M}_s^{\mathrm{phys}}\right)}^2$ (M2)
1.90	$0.0473\left(28\right)$	$0.0475\left(26\right)$
1.95	$0.0400\left(22\right)$	$0.0402\left(20\right)$
2.10	$0.0231\left(12\right)$	$0.0232\left(11\right)$

Table 6. Table 6 : Physical values for M K ​ a 0 subscript 𝑀 𝐾 subscript 𝑎 0 M_{K}a_{0} obtained from the global fit of Eq. 24 to the data from the different approaches. We also give the χ 2 superscript 𝜒 2 \chi^{2} - and p 𝑝 p -values of the fit together with the best fit parameters Q 0 subscript 𝑄 0 Q_{0} - Q 2 subscript 𝑄 2 Q_{2} .

	M1A	M2A	M1B	M2B
${\left(M_K{a}_0\right)}^{\text{phys}}$	$-0.398\left(18\right)$	$-0.397\left(18\right)$	$-0.389\left(18\right)$	$-0.384\left(16\right)$
${\chi }^2$/dof	$2.23$/7	$2.43$/7	$3.07$/7	$4.94$/7
$p$-value	$0.95$	$0.93$	$0.88$	$0.67$
$Q_0$	$-0.69\left(12\right)$	$-0.74\left(12\right)$	$-0.67\left(12\right)$	$-0.70\left(12\right)$
$Q_1$	$2.4\left(6\right)$	$2.3\left(7\right)$	$2.0\left(6\right)$	$1.7\left(6\right)$
$Q_2$	$-0.39\left(2\right)$	$-0.39\left(2\right)$	$-0.38\left(2\right)$	$-0.38\left(2\right)$

Table 7. Table 7 : Values of a ​ μ s 𝑎 subscript 𝜇 𝑠 a\mu_{s} corresponding to the renormalized physical strange quark mass in lattice units for the three values of β 𝛽 \beta calculated from Eq.  21 .

$\beta$	$a{\mu }_s^{\mathrm{ref}}$ (M1)	$a{\mu }_s^{\mathrm{ref}}$ (M2)
1.90	$0.0202\left(12\right)$	$0.0204\left(11\right)$
1.95	$0.0182\left(10\right)$	$0.0181\left(9\right)$
2.10	$0.0150\left(8\right)$	$0.0151\left(8\right)$

Table 8. Table 8 : Comparison of the scattering lengths of the L = 24 𝐿 24 L=24 ensembles determined from Lüscher’s formula to 𝒪 ​ ( L − 6 ) 𝒪 superscript 𝐿 6 \mathcal{O}(L^{-6}) , 𝒪 ​ ( L − 6 ) 𝒪 superscript 𝐿 6 \mathcal{O}(L^{-6}) without the effective range term and 𝒪 ​ ( L − 7 ) 𝒪 superscript 𝐿 7 \mathcal{O}(L^{-7}) , respectively at the lowest value of a ​ μ s 𝑎 subscript 𝜇 𝑠 a\mu_{s} .

Ensemble	$a_0$ at $𝒪\left(L^{-6}\right)$	$a_0$ at $𝒪\left(L^{-5}\right)$
A60.24	$-1.405\left(18\right)$	$-1.393\left(18\right)$
A80.24	$-1.412\left(14\right)$	$-1.400\left(14\right)$
A100.24	$-1.390\left(12\right)$	$-1.379\left(12\right)$
B85.24	$-1.592\left(20\right)$	$-1.572\left(19\right)$
A40.32	$-1.350\left(46\right)$	$-1.346\left(46\right)$
D30.48	$-2.143\left(13\right)$	$-2.130\left(12\right)$

Table 9. Table 9 : Comparison of M K ​ a 0 subscript 𝑀 𝐾 subscript 𝑎 0 M_{K}a_{0} with different orders of L 𝐿 L taken into account for determining a 0 subscript 𝑎 0 a_{0} . The data shown are the p 𝑝 p -value weighted medians over all fitranges for δ ​ E 𝛿 𝐸 \delta E at the lowest value of a ​ μ s 𝑎 subscript 𝜇 𝑠 a\mu_{s} .

Ensemble	$M_K{a}_0$ to $𝒪\left(L^{-6}\right)$	$M_K{a}_0$ to $𝒪\left(L^{-5}\right)$
A60.24	$-0.344\left(5\right)$	$-0.341\left(4\right)$
A80.24	$-0.360\left(4\right)$	$-0.357\left(3\right)$
A100.24	$-0.367\left(3\right)$	$-0.364\left(3\right)$
B85.24	$-0.368\left(5\right)$	$-0.363\left(4\right)$
A40.32	$-0.316\left(11\right)$	$-0.315\left(11\right)$
D30.48	$-0.322\left(19\right)$	$-0.320\left(18\right)$

Table 10. Table 10 : Single pion energy levels from Ref.  [ 12 ] , [ 13 ] and the finite size correction factors K M π subscript 𝐾 subscript 𝑀 𝜋 K_{M_{\pi}} and K M K subscript 𝐾 subscript 𝑀 𝐾 K_{M_{K}} computed in Ref.  [ 17 ] for M π subscript 𝑀 𝜋 M_{\pi} and M K subscript 𝑀 𝐾 M_{K} , respectively. The statistical uncertainty of K M K subscript 𝐾 subscript 𝑀 𝐾 K_{M_{K}} is only estimated. Where not given K M K subscript 𝐾 subscript 𝑀 𝐾 K_{M_{K}} is set to 1.

Ensemble	$a{M}_{\pi }$	$K_{M_{\pi }}$	$K_{M_K}$
A30.32	$0.12395\left(36\right)\left(14\right)$	$1.0081\left(52\right)$	$0.9954\left(1\right)$
A40.32	$0.14142\left(27\right)\left(42\right)$	$1.0039\left(28\right)$	$0.9974\left(1\right)$
A60.24	$0.17275\left(45\right)\left(23\right)$	$1.0099\left(49\right)$	$0.9907\left(1\right)$
A80.24	$0.19875\left(41\right)\left(35\right)$	$1.0057\left(29\right)$	$0.9950\left(1\right)$
A100.24	$0.22293\left(35\right)\left(38\right)$	$1.0037\left(19\right)$	$0.9970\left(1\right)$
B35.32	$0.12602\left(30\right)\left(30\right)$	$1.0069\left(32\right)$	$0.9951\left(1\right)$
B55.32	$0.15518\left(21\right)\left(33\right)$	$1.0027\left(14\right)$	$0.9982\left(1\right)$
B85.24	$0.19396\left(38\right)\left(54\right)$	$1.0083\left(28\right)$	$0.9937\left(1\right)$
D30.48	$0.09780\left(16\right)\left(32\right)$	$1.0021\left(7\right)$	$0.9986\left(1\right)$
D45.32	$0.12070\left(30\right)\left(10\right)$	$1.0047\left(14\right)$	$1.0000\left(1\right)$

Table 11. Table 11 : Fit ranges for the lowest value of a ​ μ s 𝑎 subscript 𝜇 𝑠 a\mu_{s} for the kaon correlation function C K subscript 𝐶 𝐾 C_{K} and the Ratio R 𝑅 R . The interval [ t i , t f ] subscript 𝑡 𝑖 subscript 𝑡 𝑓 [t_{i},t_{f}] denotes the lowest and largest timeslice considered in the fits, t min subscript 𝑡 min t_{\mathrm{min}} is the minimal extend of each fitrange.

		$C_K$		$R$
Ensemble	$a{\mu }_s$	$[t_i,t_f]$	$t_{\min }$	$[t_i,t_f]$	$t_{\min }$
A30.32	$0.0185$	$[12,32]$	$7$	$[12,32]$	$10$
A40.24	$0.0185$	$[12,24]$	$5$	$[11,24]$	$7$
A40.32	$0.0185$	$[15,32]$	$7$	$[12,32]$	$10$
A60.24	$0.0185$	$[12,24]$	$5$	$[11,24]$	$7$
A80.24	$0.0185$	$[12,24]$	$5$	$[11,24]$	$7$
A100.24	$0.0185$	$[12,24]$	$5$	$[11,24]$	$7$
B35.32	$0.0160$	$[15,32]$	$7$	$[12,32]$	$10$
B55.32	$0.0160$	$[15,32]$	$7$	$[14,32]$	$10$
B85.24	$0.0160$	$[12,24]$	$5$	$[11,24]$	$7$
D30.48	$0.0115$	$[16,48]$	$15$	$[8,41]$	$15$
D45.32	$0.0130$	$[15,32]$	$7$	$[12,32]$	$10$

Table 12. Table 12 : Same as Table  11 but for the medium value of a ​ μ s 𝑎 subscript 𝜇 𝑠 a\mu_{s} .

		$C_K$		$R$
Ensemble	$a{\mu }_s$	$[t_i,t_f]$	$t_{\min }$	$[t_i,t_f]$	$t_{\min }$
A30.32	$0.0225$	$[12,32]$	$5$	$[12,29]$	$10$
A40.24	$0.0225$	$[12,24]$	$5$	$[11,24]$	$7$
A40.32	$0.0225$	$[15,32]$	$7$	$[12,32]$	$10$
A60.24	$0.0225$	$[12,24]$	$5$	$[11,24]$	$7$
A80.24	$0.0225$	$[12,24]$	$5$	$[11,24]$	$7$
A100.24	$0.0225$	$[12,24]$	$5$	$[11,24]$	$7$
B35.32	$0.0186$	$[15,32]$	$7$	$[12,32]$	$10$
B55.32	$0.0186$	$[15,32]$	$7$	$[14,32]$	$10$
B85.24	$0.0186$	$[12,24]$	$5$	$[11,24]$	$7$
D30.48	$0.0150$	$[16,48]$	$15$	$[8,41]$	$15$
D45.32	$0.0150$	$[15,32]$	$7$	$[12,32]$	$10$

Table 13. Table 13 : Same as Table  11 but for the highest value of a ​ μ s 𝑎 subscript 𝜇 𝑠 a\mu_{s} .

		$C_K$		$R$
Ensemble	$a{\mu }_s$	$[t_i,t_f]$	$t_{\min }$	$[t_i,t_f]$	$t_{\min }$
A30.32	$0.02464$	$[12,32]$	$5$	$[12,32]$	$10$
A40.24	$0.02464$	$[12,24]$	$5$	$[11,24]$	$7$
A40.32	$0.02464$	$[15,32]$	$7$	$[12,32]$	$10$
A60.24	$0.02464$	$[12,24]$	$5$	$[11,24]$	$7$
A80.24	$0.02464$	$[12,24]$	$5$	$[11,24]$	$7$
A100.24	$0.02464$	$[12,24]$	$5$	$[11,24]$	$7$
B35.32	$0.0210$	$[15,32]$	$7$	$[12,32]$	$10$
B55.32	$0.0210$	$[15,32]$	$7$	$[14,32]$	$10$
B85.24	$0.0210$	$[12,24]$	$5$	$[11,24]$	$7$
D30.48	$0.0180$	$[16,48]$	$15$	$[8,41]$	$15$
D45.32	$0.0180$	$[15,32]$	$7$	$[12,32]$	$10$

Table 14. Table 14 : Parameters from Global fit of Eq.  19 to a ​ M K 2 𝑎 superscript subscript 𝑀 𝐾 2 aM_{K}^{2} with parameters from M1

$\beta$	$P_Z$	$P_r$	${\overline{P}}_0$	$P_1$	$P_2$	${\chi }^2$/dof
$1.90$	$0.524\left(7\right)$	$5.22\left(6\right)$	$5.53\left(20\right)$	$0.14\left(3\right)$	$5.21\left(1.61\right)$	$6.82$
$1.95$	$0.512\left(4\right)$	$5.84\left(5\right)$
$2.10$	$0.516\left(2\right)$	$7.57\left(8\right)$

Table 15. Table 15 : Parameters from Global fit of Eq.  19 to a ​ M K 2 𝑎 superscript subscript 𝑀 𝐾 2 aM_{K}^{2} with parameters from M2

$\beta$	$P_Z$	$P_r$	${\overline{P}}_0$	$P_1$	$P_2$	${\chi }^2$/dof
$1.90$	$0.572\left(4\right)$	$5.19\left(6\right)$	$5.65\left(20\right)$	$0.16\left(3\right)$	$7.27\left(1.67\right)$	$6.85$
$1.95$	$0.547\left(2\right)$	$5.87\left(5\right)$
$2.10$	$0.545\left(2\right)$	$7.56\left(8\right)$

Table 16. Table 16 : Lattice results for a ​ M K 𝑎 subscript 𝑀 𝐾 aM_{K} , δ ​ E 𝛿 𝐸 \delta E , a 0 subscript 𝑎 0 a_{0} and M K ​ a 0 subscript 𝑀 𝐾 subscript 𝑎 0 M_{K}a_{0} for the smallest value of a ​ μ s 𝑎 subscript 𝜇 𝑠 a\mu_{s} on all ensembles used in this study. The first parentheses states the statistical uncertainty estimated from the bootstrap samples of the quantity, the second one states the systematic uncertainty estimated from the different fit ranges used for the correlation functions.

Ens	$a{\mu }_s$	$a{M}_K$	$a\delta E$	$a_0$	$\left(M_K{a}_0\right)$
A30.32	0.0185	$0.2292\left(2\right){(}_{-0}^{+0})$	$0.0025\left(1\right){(}_{-0}^{+1})$	$-1.306\left(82\right){(}_{-15}^{+29})$	$-0.299\left(19\right){(}_{-3}^{+7})$
A40.20	0.0185	$0.2385\left(5\right){(}_{-0}^{+0})$	$0.0126\left(2\right){(}_{-1}^{+2})$	$-1.523\left(19\right){(}_{-15}^{+19})$	$-0.363\left(5\right){(}_{-3}^{+5})$
A40.24	0.0185	$0.2364\left(3\right){(}_{-0}^{+0})$	$0.0065\left(1\right){(}_{-0}^{+1})$	$-1.423\left(17\right){(}_{-3}^{+20})$	$-0.336\left(4\right){(}_{-1}^{+5})$
A40.32	0.0185	$0.2342\left(2\right){(}_{-0}^{+0})$	$0.0025\left(1\right){(}_{-0}^{+0})$	$-1.346\left(46\right){(}_{-14}^{+9})$	$-0.315\left(11\right){(}_{-3}^{+2})$
A60.24	0.0185	$0.2449\left(3\right){(}_{-0}^{+0})$	$0.0061\left(1\right){(}_{-0}^{+1})$	$-1.393\left(18\right){(}_{-7}^{+18})$	$-0.341\left(4\right){(}_{-2}^{+4})$
A80.24	0.0185	$0.2548\left(2\right){(}_{-1}^{+1})$	$0.0059\left(1\right){(}_{-0}^{+1})$	$-1.400\left(14\right){(}_{-7}^{+16})$	$-0.357\left(3\right){(}_{-2}^{+4})$
A100.24	0.0185	$0.2642\left(2\right){(}_{-1}^{+1})$	$0.0056\left(1\right){(}_{-0}^{+0})$	$-1.379\left(12\right){(}_{-2}^{+2})$	$-0.364\left(3\right){(}_{-1}^{+1})$
B35.32	0.0160	$0.2053\left(2\right){(}_{-0}^{+0})$	$0.0035\left(1\right){(}_{-0}^{+0})$	$-1.606\left(58\right){(}_{-16}^{+20})$	$-0.330\left(12\right){(}_{-3}^{+4})$
B55.32	0.0160	$0.2153\left(2\right){(}_{-0}^{+0})$	$0.0030\left(1\right){(}_{-0}^{+0})$	$-1.491\left(47\right){(}_{-21}^{+18})$	$-0.321\left(10\right){(}_{-4}^{+4})$
B85.24	0.0160	$0.2312\left(3\right){(}_{-0}^{+0})$	$0.0075\left(1\right){(}_{-1}^{+1})$	$-1.572\left(19\right){(}_{-17}^{+16})$	$-0.363\left(4\right){(}_{-4}^{+4})$
D30.48	0.0115	$0.1504\left(1\right){(}_{-0}^{+0})$	$0.0018\left(1\right){(}_{-1}^{+1})$	$-2.130\left(123\right){(}_{-97}^{+94})$	$-0.320\left(18\right){(}_{-15}^{+14})$
D45.32	0.0130	$0.1657\left(3\right){(}_{-0}^{+0})$	$0.0066\left(2\right){(}_{-2}^{+0})$	$-2.307\left(51\right){(}_{-51}^{+12})$	$-0.382\left(9\right){(}_{-9}^{+2})$

Table 17. Table 17 : Same as Table  16 but for the medium value of a ​ μ s 𝑎 subscript 𝜇 𝑠 a\mu_{s}

Ens	$a{\mu }_s$	$a{M}_K$	$a\delta E$	$a_0$	$\left(M_K{a}_0\right)$
A30.32	0.0225	$0.2491\left(3\right){(}_{-2}^{+2})$	$0.0025\left(2\right){(}_{-2}^{+1})$	$-1.436\left(97\right){(}_{-104}^{+45})$	$-0.358\left(24\right){(}_{-26}^{+11})$
A40.20	0.0225	$0.2577\left(5\right){(}_{-0}^{+0})$	$0.0120\left(2\right){(}_{-1}^{+0})$	$-1.565\left(19\right){(}_{-11}^{+4})$	$-0.403\left(5\right){(}_{-3}^{+1})$
A40.24	0.0225	$0.2560\left(3\right){(}_{-0}^{+0})$	$0.0062\left(1\right){(}_{-1}^{+0})$	$-1.454\left(18\right){(}_{-21}^{+9})$	$-0.372\left(5\right){(}_{-5}^{+2})$
A40.32	0.0225	$0.2538\left(2\right){(}_{-0}^{+0})$	$0.0024\left(2\right){(}_{-0}^{+1})$	$-1.391\left(86\right){(}_{-23}^{+46})$	$-0.353\left(22\right){(}_{-6}^{+12})$
A60.24	0.0225	$0.2638\left(3\right){(}_{-1}^{+1})$	$0.0060\left(1\right){(}_{-1}^{+1})$	$-1.459\left(16\right){(}_{-28}^{+14})$	$-0.385\left(4\right){(}_{-7}^{+4})$
A80.24	0.0225	$0.2732\left(2\right){(}_{-1}^{+1})$	$0.0057\left(1\right){(}_{-0}^{+1})$	$-1.435\left(14\right){(}_{-7}^{+13})$	$-0.392\left(4\right){(}_{-2}^{+4})$
A100.24	0.0225	$0.2823\left(2\right){(}_{-0}^{+0})$	$0.0054\left(1\right){(}_{-0}^{+0})$	$-1.412\left(12\right){(}_{-10}^{+8})$	$-0.398\left(3\right){(}_{-3}^{+2})$
B35.32	0.0186	$0.2185\left(2\right){(}_{-1}^{+1})$	$0.0032\left(1\right){(}_{-1}^{+1})$	$-1.589\left(62\right){(}_{-47}^{+37})$	$-0.347\left(14\right){(}_{-10}^{+8})$
B55.32	0.0186	$0.2281\left(2\right){(}_{-0}^{+0})$	$0.0031\left(1\right){(}_{-0}^{+0})$	$-1.578\left(51\right){(}_{-9}^{+12})$	$-0.360\left(12\right){(}_{-2}^{+3})$
B85.24	0.0186	$0.2429\left(3\right){(}_{-1}^{+1})$	$0.0074\left(1\right){(}_{-0}^{+1})$	$-1.619\left(14\right){(}_{-9}^{+13})$	$-0.393\left(3\right){(}_{-2}^{+3})$
D30.48	0.0150	$0.1674\left(1\right){(}_{-0}^{+0})$	$0.0018\left(1\right){(}_{-1}^{+1})$	$-2.318\left(135\right){(}_{-59}^{+126})$	$-0.388\left(23\right){(}_{-10}^{+21})$
D45.32	0.0150	$0.1748\left(2\right){(}_{-3}^{+3})$	$0.0061\left(2\right){(}_{-3}^{+3})$	$-2.250\left(49\right){(}_{-77}^{+80})$	$-0.393\left(9\right){(}_{-13}^{+14})$

Table 18. Table 18 : Same as Table  16 but for the highest value of a ​ μ s 𝑎 subscript 𝜇 𝑠 a\mu_{s}

Ens	$a{\mu }_s$	$a{M}_K$	$a\delta E$	$a_0$	$\left(M_K{a}_0\right)$
A30.32	0.0246	$0.2590\left(3\right){(}_{-2}^{+2})$	$0.0025\left(2\right){(}_{-1}^{+2})$	$-1.535\left(126\right){(}_{-26}^{+54})$	$-0.398\left(33\right){(}_{-7}^{+14})$
A40.20	0.0246	$0.2679\left(5\right){(}_{-0}^{+0})$	$0.0117\left(2\right){(}_{-0}^{+3})$	$-1.584\left(19\right){(}_{-5}^{+33})$	$-0.424\left(5\right){(}_{-1}^{+9})$
A40.24	0.0246	$0.2660\left(3\right){(}_{-1}^{+1})$	$0.0060\left(1\right){(}_{-1}^{+1})$	$-1.476\left(18\right){(}_{-23}^{+18})$	$-0.393\left(5\right){(}_{-6}^{+5})$
A40.32	0.0246	$0.2638\left(2\right){(}_{-0}^{+0})$	$0.0025\left(1\right){(}_{-0}^{+1})$	$-1.484\left(49\right){(}_{-20}^{+40})$	$-0.391\left(13\right){(}_{-5}^{+10})$
A60.24	0.0246	$0.2736\left(3\right){(}_{-0}^{+0})$	$0.0056\left(1\right){(}_{-0}^{+1})$	$-1.425\left(20\right){(}_{-3}^{+27})$	$-0.390\left(5\right){(}_{-1}^{+7})$
A80.24	0.0246	$0.2825\left(2\right){(}_{-2}^{+2})$	$0.0056\left(1\right){(}_{-1}^{+1})$	$-1.451\left(15\right){(}_{-18}^{+18})$	$-0.410\left(4\right){(}_{-5}^{+5})$
A100.24	0.0246	$0.2914\left(2\right){(}_{-1}^{+1})$	$0.0053\left(1\right){(}_{-0}^{+1})$	$-1.421\left(12\right){(}_{-8}^{+14})$	$-0.414\left(4\right){(}_{-2}^{+4})$
B35.32	0.0210	$0.2298\left(2\right){(}_{-0}^{+0})$	$0.0035\left(2\right){(}_{-1}^{+1})$	$-1.778\left(87\right){(}_{-27}^{+35})$	$-0.409\left(20\right){(}_{-6}^{+8})$
B55.32	0.0210	$0.2388\left(2\right){(}_{-0}^{+0})$	$0.0029\left(1\right){(}_{-1}^{+0})$	$-1.547\left(55\right){(}_{-25}^{+21})$	$-0.369\left(13\right){(}_{-6}^{+5})$
B85.24	0.0210	$0.2535\left(3\right){(}_{-0}^{+0})$	$0.0071\left(1\right){(}_{-1}^{+0})$	$-1.624\left(16\right){(}_{-12}^{+2})$	$-0.412\left(4\right){(}_{-3}^{+0})$
D30.48	0.0180	$0.1807\left(1\right){(}_{-0}^{+0})$	$0.0018\left(1\right){(}_{-0}^{+1})$	$-2.449\left(147\right){(}_{-43}^{+122})$	$-0.443\left(27\right){(}_{-8}^{+22})$
D45.32	0.0180	$0.1875\left(2\right){(}_{-1}^{+1})$	$0.0057\left(2\right){(}_{-1}^{+2})$	$-2.245\left(54\right){(}_{-36}^{+79})$	$-0.421\left(10\right){(}_{-7}^{+15})$

Equations

$$\begin{array}{rl}{r}_0& =0.470(12)\mathrm{fm}\text{(M1)},\\ r_0& =0.471(11)\mathrm{fm}\text{(M2)}.\end{array}$$

$$D_{\ell }=D_{\mathrm{W}}+m_0+i{\mu }_{\ell }{\gamma }_5{\tau }^3,$$

$$D_{\mathrm{h}}=D_{\mathrm{W}}+m_0+i{\mu }_{\sigma }{\gamma }_5{\tau }^1+{\mu }_{\delta }{\tau }^3.$$

$$D_s^{\pm }=D_{\mathrm{W}}+m_0\pm i{\mu }_s{\gamma }_5,$$

$${\mathcal{O}}_{\pi }(t)=\sum _{\mathbf{x}}\bar{u}(\mathbf{x},t)i{\gamma }_{5}d(\mathbf{x},t)$$

$${\mathcal{O}}_K(t)=\sum _{\mathbf{x}}{\bar{s}}^{+}(\mathbf{x},t)i{\gamma }_{5}d(\mathbf{x},t)$$

$$C_K(t-t^{\prime })=⟨{\mathcal{O}}_K(t){\mathcal{O}}_K^{†}(t^{\prime })⟩$$

$$C_K\propto \frac{1}{2}\left(e^{-M_{K}t}+e^{-M_K(T-t)}\right),$$

$${\mathcal{O}}_{KK}(t)=\sum _{\mathbf{x},{\mathbf{x}}^{\prime }}{\bar{s}}^{+}(\mathbf{x},t)i{\gamma }_{5}d(\mathbf{x},t){\bar{s}}^{+}({\mathbf{x}}^{\prime },t)i{\gamma }_{5}d({\mathbf{x}}^{\prime },t)$$

$$C_{KK}(t-t^{\prime })=⟨{\mathcal{O}}_{KK}(t){\mathcal{O}}_{KK}^{†}(t^{\prime })⟩.$$

$$C_{KK}\propto \frac{1}{2}\left(e^{-E_{KK}t}+e^{-E_{KK}(T-t)}\right)+\mathrm{const}.$$

$$R(t+1/2)=\frac{C_{KK}(t)-C_{KK}(t+1)}{C_K^2(t)-C_K^2(t+1)}$$

$$R(t+1/2)=A(\cosh (\delta E{t}^{\prime })+\sinh (\delta E{t}^{\prime })\coth (2E_K{t}^{\prime })),$$

$$a{M}_H^{\ast }:=\frac{a{M}_H}{K_{M_H}},$$

$$S=V_S{V}_S^{†}.$$

$$\widetilde{\Delta}^{ab}(x,y;U)=\sum_{k=1}^{3}\Bigl{\{}\widetilde{U}_{k}^{ab}(x)\delta(y,x+\hat{k})+\widetilde{U}^{ba}_{k}(y)^{\dagger}\delta(y,x-\hat{k})-2\delta(x,y)\delta^{ab}\Bigr{\}}\,,$$

$$\mathcal{Q}=\mathcal{S}{\Omega }^{-1}\mathcal{S}=V_{\mathrm{s}}(V_{\mathrm{s}}^{†}{\Omega }^{-1}{V}_{\mathrm{s}})V_{\mathrm{s}}^{†}.$$

$${\Omega }^{-1}\approx \frac{1}{N_{\mathrm{R}}}\sum _{r=1}^{N_{\mathrm{R}}}\sum _b{X}^{r[b]}{\rho }^{r[b]†},$$

$$\delta E=-\frac{4\pi a_0}{M_K{L}^3}\left(1+c_1\frac{a_0}{L}+c_2\frac{a_0^2}{L^2}+c_3\frac{a_0^3}{L^3}\right)-\frac{8{\pi }^2{a}_0^3}{M_K{L}^6}{r}_f+\mathcal{O}(L^{-7}),$$

$$c_1=-2.837297,c_2=6.375183,c_3=-8.311951.$$

$$M_s^2=M_K^2-M_{\pi }^2/2$$

$$(a{M}_K{)}^2=\frac{P_0}{P_r{P}_Z}(a{\mu }_l+a{\mu }_s)\left[1+P_1\frac{P_r}{P_Z}a{\mu }_l+\frac{P_2}{P_r^2}\right]$$

$${\chi }_{\mathrm{aug}}^2={\chi }^2+\sum _{\beta }\left[{\left(\frac{(r_0/a)(\beta )-P_r(\beta )}{\Delta r_0/a(\beta )}\right)}^2+{\left(\frac{Z_P(\beta )-P_Z(\beta )}{\Delta Z_P(\beta )}\right)}^2\right].$$

$$a{\mu }_s^{\mathrm{ref}}=\frac{(r_0{M}_K^{\mathrm{phys}}{)}^2{P}_Z}{P_r{P}_0[1+P_1{r}_0{m}_{\ell }^{\mathrm{phys}}+P_2{P}_r^{-2}]}-\frac{P_Z}{P_r}(r_0{m}_{\ell }^{\mathrm{phys}})$$

$$r_0{m}_s^{\mathrm{phys}}=\frac{(r_0{M}_K^{\mathrm{phys}}{)}^2}{P_0[1+P_1{r}_0{m}_{\ell }]}-(r_0{m}_{\ell }^{\mathrm{phys}}).$$

$$M_K{a}_0=\frac{M_K^2}{8\pi f_K^2}\left[-1+\frac{16}{f_K^2}\left(M_K^2{L}^{\prime }-\frac{M_K^2}{2}{L}_5+\zeta \right)\right].$$

$$M_K{a}_0=Q_0\frac{P_r}{P_Z}a{\mu }_{\ell }+Q_1\frac{1}{P_r^2}+Q_2,$$

$$(M_K{a}_0{)}^{\mathrm{phys}}=Q_0{r}_0{m}_{\ell }^{\mathrm{phys}}+Q_2.$$

$$w_X={\left[(1-2|p_X-0.5{|}^2)\cdot \min (\Delta X)/\Delta X\right]}^2,$$

$$\begin{split}m_{s}^{\textrm{phys.}}&=101.3(4.7)$\mathrm{MeV}$\qquad\textbf{(M1B)}\,\\ m_{s}^{\textrm{phys.}}&=99.4(4.4)$\mathrm{MeV}$\qquad\textbf{(M2B)}\,.\end{split}$$

Full text

1 Introduction

Shortly after the Big Bang the universe is believed to have been in a quark gluon plasma state of matter. Apart from the inside of neutron stars the only places where this state of matter appears and can be studied are detectors investigating heavy ion or proton-proton collisions like the STAR detector at the Relativistic Heavy Ion Collider (RHIC) at BNL [1] or the ALICE experiment at the LHC at CERN [2]. The collisions taking place at such sites yield in their final states numerous light hadrons like pions and kaons. Due to the mass difference between kaons and pions the produced kaons carry much lower momenta than the pions, therefore being much more likely to interact elastically. The interaction of two kaons is determined by Quantum Chromodynamics (QCD), which is non-perturbative at low energies. The understanding and interpretation of the results of the aforementioned experiments make a non-perturbative investigation of kaon-kaon interactions highly desirable. While this can be formulated in chiral perturbation theory (ChPT), it is theoretically interesting to check if the effective approach is able to properly describe kaon-kaon scattering. Lattice QCD provides a non-perturbative ab initio method to perform such a study.

Hadron-hadron scattering has become more and more accessible to lattice QCD simulations over the last years. This is on the one hand due to Lüschers finite volume formalism, and on the other hand due to lattice QCD ensembles becoming ever more realistic. For kaon-kaon scattering in the isospin-1 channel only a few lattice QCD calculations have been performed [3, 4] where the result of the former calculation has been used in Ref. [2] for the ALICE results. In the maximal isospin channel kaon-kaon scattering resembles the well studied pion-pion case [5, 6, 7, 8, 9, 10, 4]: there are no fermionic disconnected diagrams and only one light quark is replaced by a strange quark. Since we already investigated pion-pion scattering in the isospin-2 channel [11] a lot of our analysis tools can be carried over to the present investigation.

In this paper we present the first study of $K^{+}K^{+}$ scattering from lattice QCD based on $N_{f}=2+1+1$ ensembles of the European Twisted Mass Collaboration (ETMC) [12, 13] covering three values of the lattice spacing. These ensembles, which employ up to five values of the light quark mass per lattice spacing value allow us to perform reliable chiral and continuum extrapolations of our results.

For the strange quark we employ a mixed action approach with so-called Osterwalder-Seiler valence quarks on the Wilson twisted mass sea [14]. This allows us to tune the valence strange quark mass value to its physical value without spoiling the automatic $\mathcal{O}(a)$-improvement guaranteed by Wilson twisted mass lattice QCD at maximal twist [15]. However, while unitarity breaking effects vanish in the continuum limit, this ansatz also introduces partial-quenching effects, which we cannot control in the present calculation. However, in previous calculations with this setup, no sizable effects were found, see e.g. [16, 17]. The mixed-action approach for the strange quark also allows us to avoid the parity-flavour mixing present in the $1+1$ (strange-charm) sea sector of Wilson twisted mass lattice QCD at maximal twist with $N_{f}=2+1+1$ flavours.

Our final result differs by about $2\sigma$ from the determinations by NPLQCD [3] and about $4\sigma$ from the determination of PACS-CS [4]. This deviation can likely be attributed to lattice artefacts: NPLQCD works mainly at a single lattice spacing with the exception of one ensemble at a second lattice spacing value. PACS-CS works at a single lattice spacing only. However, we can also not exclude residual unitarity breaking effects in our calculation. Interestingly, our result is actually equal to the leading order ChPT prediction for $M_{K}a_{0}$.

2 Lattice action

We use gauge configurations generated by the ETM collaboration with $N_{f}=2+1+1$ dynamical quark flavours [12]. The Iwasaki gauge action [18] is used in combination with the Wilson twisted mass fermion discretisation. There are three values of the lattice spacing available, with $\beta=1.90$, $\beta=1.95$ and $\beta=2.10$ corresponding to $a\sim 0.089\ \mathrm{fm}$, $a\sim 0.082\ \mathrm{fm}$ and $a=0.062\ \mathrm{fm}$, respectively. The ensembles we used are compiled in table 1. The lattice scale for the ensembles has been determined in Ref. [17] using $f_{\pi}$. Also in Ref. [17] the pseudoscalar renormalisation constant $Z_{P}$, the inverse of which is the quark mass renormalisation constant in the twisted-mass approach, has been determined for each lattice spacing and then converted to the $\overline{\textrm{MS}}$ scheme at a scale of $2\ \mathrm{GeV}$.

The computation of $Z_{P}$ employs the RI-MOM renormalisation scheme and further makes use of two different methods which are labelled M1 and M2 by the authors. The two methods, M1 and M2, give results which differ by lattice artefacts. As an intermediate length scale we use the Sommer parameter $r_{0}/a$ determined in Ref. [12] for each value of the light quark mass $m_{l}$ and extrapolated to the chiral limit in Ref. [17], assuming either a linear or quadratic dependence on the light quark mass. The value of $r_{0}$ in fm was determined in Ref. [17] using chiral perturbation theory (ChPT) employing methods M1 and M2 for $Z_{P}$, reading

[TABLE]

We keep the two values separate here, because we will use them to estimate systematic uncertainties. The values for $Z_{P}$, the lattice spacing $a$ and $r_{0}/a$ are summarised in table 2 for the three $\beta$-values. For details we refer to Ref. [17]. Note that $\mu_{\sigma}$ and $\mu_{\delta}$ are kept fixed for all $\mu_{\ell}$ values at $\beta=1.90$ and $\beta=1.95$. Between the two ensembles D30.48 and D45.32sc they differ slightly.

In order to set the strange quark mass, we use $M_{K}$ in physical units as input. We use $M^{\textrm{phys}}_{K}=$494.2(3)\text{,}\mathrm{MeV}$$ corrected for electromagnetic and isospin breaking effects [19].

As further inputs we use the average up/down quark mass, $m_{l}^{\textrm{phys}}=$3.70(17)\text{,}\mathrm{MeV}, from Ref. [[17](#bib.bib17)] as well as the neutral pion mass, $M^{\textrm{phys}}_{\pi^{0}}=$134.98\text{\,}\mathrm{MeV} [20].

In more detail, for the sea quarks we use the Wilson twisted mass action with $N_{f}=2+1+1$ dynamical quark flavours. The Dirac operator for the light quark doublet reads [21]

[TABLE]

where $D_{\mathrm{W}}$ denotes the standard Wilson Dirac operator and $\mu_{\ell}$ the bare light twisted mass parameter. $\tau^{3}$ and in general $\tau^{i},i=1,2,3$ represent the Pauli matrices acting in flavour space. $D_{\ell}$ acts on a spinor $\chi_{\ell}=(u,d)^{T}$ and, hence, the $u$ ($d$) quark has twisted mass $+\mu_{\ell}$ ($-\mu_{\ell}$).

For the heavy doublet of $c$ and $s$ quarks [14] the Dirac operator is given by

[TABLE]

The bare Wilson quark mass $m_{0}$ has been tuned to its critical value $m_{\mathrm{crit}}$ [22, 12]. This guarantees automatic order $\mathcal{O}\left(a\right)$-improvement [15], which is one of the main advantages of the Wilson twisted mass formulation of lattice QCD. For a discussion on how to tune to $m_{\mathrm{crit}}$ we refer to Refs. [22, 12].

The splitting term in the heavy doublet Eq. 3 introduces parity and flavour mixing between strange and charm quarks which would render the present analysis very complicated. For this reason we rely in this paper on a mixed-action approach for the strange quark: in the valence sector we use the so-called Osterwalder-Seiler (OS) discretisation [23] with Dirac operator

[TABLE]

with bare strange quark mass $\mu_{s}$. Formally, this introduces two valence strange quarks with $\pm\mu_{s}$ as bare quark mass. We will denote these two as $s^{\pm}$ and they will coincide in the continuum limit. Hence, observables computed using the one or the other will differ by $\mathcal{O}(a^{2})$ lattice artefacts. It was shown in Ref. [23] that $\mathcal{O}(a)$-improvement stays intact when $m_{0}$ is set to the same value $m_{\mathrm{crit}}$ as used in the unitary sector. For each $\beta$-value, we choose a set of three bare strange quark masses $a\mu_{s}$ as listed in Table 3. The mass values are chosen such as to bracket the physical strange quark mass independently of the light quark mass.

We remark here that in twisted mass lattice QCD the quark masses renormalise multiplicatively with $1/Z_{P}$ [21]. Since OS and unitary actions agree in the chiral limit, also the OS strange quark mass renormalises multiplicatively with $1/Z_{P}$.

2.1 Lattice Operators and Correlation Functions

For the charged pion we use the interpolating operator

[TABLE]

projected to zero momentum. As interpolating operator with the quantum-numbers of the kaon we use

[TABLE]

projected to zero momentum. We use the combination of a strange quark with $+|\mu_{s}|$ and the down quark with $-|\mu_{\ell}|$, because it is known that observables employing this combination are subject to milder lattice artefacts compared to the combination with same signs. The corresponding two-point function reads

[TABLE]

and likewise the pseudo-scalar two point function $C_{\pi}$ with $\mathcal{O}_{K}$ replaced by $\mathcal{O}_{\pi}$. From the behaviour of $C_{K}$ ($C_{\pi}$) at large Euclidean time

[TABLE]

the kaon mass $aM_{K}$ ($aM_{\pi}$) can be extracted. In order to compute the finite volume energy shift $\delta E=E_{KK}-2M_{K}$, needed in Lüschers formula to obtain the scattering length $a_{0}$, we have to determine the energy of the two kaon system in the interacting case. Using the isospin $I=1$ operator

[TABLE]

one defines the correlation function

[TABLE]

It shows a dependence on Euclidean time similar to $C_{K}$ with the addition of a time independent piece, the so-called thermal pollution

[TABLE]

To determine $\delta E$ from $C_{KK}$ we use a method which was devised in Ref. [8] for the $\pi\pi$ system with $I=2$. In this method, we consider the ratio

[TABLE]

which can be shown to have the large Euclidean time dependence

[TABLE]

with $t^{\prime}=t+1/2-T/2$ and amplitude $A$.

The kaon and pion masses are affected by (exponentially suppressed) finite size effects. The corresponding ChPT corrections $K_{M_{\pi}}=M_{\pi}(L)/M_{\pi}(L=\infty)$ and $K_{M_{K}}=M_{K}(L)/M_{K}(L=\infty)$ were determined from the data in Ref. [17] and we reuse these values, which are collected in Table 10. From here on we only work with finite size corrected hadron masses:

[TABLE]

for $H=\pi,K$ and drop the asterisk to ease the notation.

2.2 Stochastic LapH

As a smearing scheme we employ the so-called stochastic Laplacian-Heaviside (sLapH) method [24, 25]. In this approach the quark field under consideration is smeared with the so-called smearing matrix

[TABLE]

The matrices $V_{S}$ are matrices obtained by stacking the eigenvectors of the lattice Laplacian,

[TABLE]

columnwise. The complete set of eigenvectors spans the so called LapH-space. The indices $a$,$b$ denote different colours, the variables $x$,$y$ space-timepoints and $\widetilde{U}$ (possibly smeared) $SU(3)$-gauge link matrices. The index $S$ on $V_{S}$ denotes a truncation of the eigenspectrum of $\widetilde{\Delta}$ such that excited state contaminations of the quark field are maximally suppressed. In addition we smear the gauge fields appearing in Eq. 14 with 3 iterations of 2 level HYP smearing [26], with parameters $\alpha_{1}=\alpha_{2}=0.62$. To build correlation functions we denote quark lines connecting source and sink timeslices with

[TABLE]

where $\Omega^{-1}$ denotes the quark propagator and $\mathcal{P}=(V_{\mathrm{s}}^{\dagger}\Omega^{-1}V_{\mathrm{s}})$ is called perambulator. We use all-to-all propagators to calculate the correlation functions which can get prohibitively expensive when done exactly. Therefore, we employ a stochastic method with random vectors diluted in time, Dirac-space and LapH-subspace. The all-to-all propagator then reads

[TABLE]

with the number of random vectors $N_{\mathrm{R}}$ and the compound index $r[b]$ counting the total number of random vectors and the total number of dilution vectors $N_{\textrm{D}}$. For the kaon correlation functions we reused the light quark propagators already calculated for the $\pi\pi$ paper, Ref. [11]. The number of dilution vectors for the light quark propagators, therefore, is the same. An exception is ensemble D30.48 which was not included in the $\pi\pi$ paper. For this volume of $L/a=48$ the values for the several $N_{\textrm{D}}$ are collected in Table 4 together with the values of $N_{\textrm{D}}$ for the other lattice sizes. Concerning the newly calculated strange quark propagators we adopted the same dilution scheme.

An investigation of the number of random vectors $N_{\mathrm{R}}$ yielded no further error reduction for the energy shift $\delta E$ when increasing $N_{\mathrm{R}}$ from 4 to 5 random vectors for each strange quark perambulator. Thus we decided to take 4 random vectors per strange quark perambulator into account for the current analysis.

3 Analysis Methods

3.1 Lüscher Method

We are interested in the limit of small scattering momenta for the kaon-kaon system with $I=1$ below inelastic threshold. Very much like in the case of $\pi\pi$ scattering with $I=2$, the scattering length $a_{0}$ can be related in the finite range expansion to the energy shift $\delta E$ by an expansion in $1/L$ as follows [27]

[TABLE]

with coefficients [27, 28]

[TABLE]

Here, $r_{f}$ is the effective range parameter. Eq. 17 can be solved for the scattering length $a_{0}/a$ given $L/a$, $a\delta E$ and $aM_{K}$ if the terms up to $\mathcal{O}(1/L^{5})$ are taken into account. This approach is valid only if the residual exponentially suppressed finite volume effects are negligible compared to the ones related for $\delta E$. Moreover, by truncating Eq. 17 at $\mathcal{O}(1/L^{5})$, one assumes that the effective range has no sizeable contribution. We estimate the effect of this truncation in Appendix A and find it to be negligible.

3.2 Chiral and Continuum Extrapolations

The values of $\delta E$ and $a_{0}$ are calculated for each combination of $a\mu_{s}$ and $a\mu_{\ell}$. In order to arrive at our final values for the scattering length, we need to perform interpolations in the strange quark mass, extrapolations in the light quark mass and the continuum extrapolation. We adopt the following strategy: we will first tune the renormalised strange quark to its physical value for all $\beta$-values and light quark masses. Next we interpolate $M_{k}a_{0}$ in the strange quark mass for all ensembles to this value. The value for $M_{K}a_{0}$ obtained from this interpolation are finally extrapolated to the physical point and the continuum limit in a combined fit.

We use two different strategies, from here on denoted by A and B, to tune the renormalised strange quark mass to its physical value:

1. A:

as a strange quark mass proxy we use

[TABLE]

which is directly proportional to the strange quark mass at leading order in ChPT. We interpolate $M_{K}a_{0}$ linearly in $(aM_{s})^{2}$ to the value where $M_{s}^{2}$ assumes its physical value for each ensemble separately. This requires the physical value of $M_{K}$ and $M_{\pi}$ and the lattice spacing as an input. The bare strange quark mass is not explicitly used in this case. 2. B:

here we are going to use the bare strange quark mass parameter $\mu_{s}$ explicitly. To determine the renormalised, physical value of the strange quark mass, we first perform a global fit of the NLO SU$(2)$ ChPT prediction for $M_{K}^{2}$

[TABLE]

to all our data for $aM_{K}$ simultaneously. Note that in SU$(2)$ ChPT there are no chiral logarithms in $M_{K}^{2}$ predicted at NLO. Here we have three global fit parameters $P_{0}$, $P_{1}$ and $P_{2}$. In addition, we have $\beta$-dependent fit parameters $P_{r}(\beta)$ and $P_{Z}(\beta)$ for $r_{0}/a$ and $Z_{P}$, respectively, which we constrain using Gaussian priors based on the determinations of these from Ref. [17].

Hence, we have in total nine fit parameters for which we define the augmented $\chi^{2}$ function:

[TABLE]

Using the best fit parameters, $a\mu_{s}^{\mathrm{ref}}$ can be determined from

[TABLE]

using the input values specified before.

This allows us to interpolate $M_{K}a_{0}$ in $a\mu_{s}$ to the reference value $a\mu_{s}^{\mathrm{ref}}$ for each ensemble separately. In the continuum limit, the physical value of the renormalised strange quark mass, $r_{0}m_{s}^{\mathrm{phys}}$, is then given by

[TABLE]

In the following we will denote the combination of M1 with strategy A as M1A and likewise M1B, M2A and M2B.

The values of $M_{K}a_{0}$ interpolated as explained above are now to be understood at fixed renormalised strange quark mass. The quark mass dependence of $M_{K}a_{0}$ is known from ChPT and is given at NLO [29, 30, 31] by

[TABLE]

Here, $L_{5}$ is a low energy constant (LEC) and $L^{\prime}$ a combination of LECs. $\zeta$ is a known function with chiral logarithms, which can be found in the references above. We can rewrite Eq. 23 in terms of the quark masses by replacing $M_{K}^{2}$ and $f_{K}$ by their corresponding LO ChPT expressions. Note that we use the convention with $f_{\pi}=130\ \mathrm{MeV}$.

As we will see later, our data for $M_{K}a_{0}$ is not sufficiently precise to resolve terms beyond leading order, in contrast to $M_{K}^{2}$. Including lattice artefacts of order $a^{2}$, we therefore resort to the following effective fit ansatz for $M_{K}a_{0}$ linear in $\mu_{\ell}$ and $a^{2}$

[TABLE]

with three free fit parameters $Q_{0}$, $Q_{1}$ and $Q_{2}$. The continuum and chiral limit for $M_{K}a_{0}$ is then given by

[TABLE]

For the fit we use again an augmented $\chi^{2}$ like in Eq. 20 to take the errors on $r_{0}/a$ and $Z_{P}$ into account.

All errors are computed using the (chained) bootstrap with $1500$ bootstrap samples. Values not determined by ourselves, e.g. for $r_{0}/a$ or $Z_{P}$ are included in the bootstrap analysis using the parametric bootstrap. Where relevant, fits are fully correlated. The configurations used are well separated in HMC trajectories and we have checked explicitly for autocorrelation using a blocked bootstrap.

4 Results

In this section we present the results for the energy shift $\delta E$, the scattering length $a_{0}$ and the chiral and continuum extrapolations of $M_{K}a_{0}$. From the four approaches M1A, M1B, M2A and M1B we obtain four values for $M_{K}a_{0}$, which we combine into our final result. The spread between the four values is used to estimate the systematic uncertainty.

4.1 Energy Shift $\delta E$

The energy shift is calculated from fitting Eq.13 to the data of the ratio defined in Eq.12. Because of the $\cosh$-like behaviour of $C_{K}$ and $C_{KK}$, we symmetrize the correlation functions. For the kaon masses we use the results of fully correlated fits to the two-point correlation function Eq. 8. We repeat our fits for multiple fit ranges for each correlation function. The systematic uncertainties of the fitting procedures are then estimated using the approach introduced in Ref. [11]. The energy value is determined as the median of the weighted distribution over the fit ranges. The weight assigned to each fit reads

[TABLE]

where $X=E_{K},\delta E$. $p_{X}$ is the $p$-value of the fit and $\Delta X$ denotes the statistical uncertainty of the considered quantity $\braket{X}$. An estimate of the systematic uncertainty is then calculated from the 68.54% confidence interval of the weighted distribution of $X$. The statistical error comes from bootstrapping this procedure.

In order to choose the fit ranges for obtaining $M_{K}$ from $C_{K}$ and $\delta E$ from $R$, we require several criteria to be fulfilled. Concerning the initial timeslice $t_{i}$, we demand that excited states, both in $C_{K}$ and $R$ have decayed away sufficiently. For $C_{K}$, we visually inspect the effective mass. Since $C_{K}$ does not suffer from exponential error growth at late times we set $t_{f}=T/2$. Thus we vary $t_{i}$ and $t_{f}$ within the constraints above. In the case of the ratio, $t_{f}$ is set to the timeslice where $R$ starts to deviate significantly from the behaviour suggested by Eq. 13. The minimal number of timeslices for a fit range is chosen with the same criterion as for $C_{K}$. The values for $t_{i}$, $t_{f}$ and $t_{\textrm{min}}$ for $C_{K}$ and $R$ are compiled in Tables 11-13 for each value of $a\mu_{s}$ in the Appendix B.

In Fig. 1 we show exemplary fits of the ratio Eq. 13 to the data for several ensembles and selected fit ranges. At least for the Ensembles with $L=24$ the tendency of an upward bend of the data at late times can be seen clearly.

As mentioned before, for Eq. 17 to be valid residual exponentially suppressed finite volume effects must be negligible. Moreover, the terms in Eq. 17 of order $1/L^{6}$ and higher must be negligible. We can test the latter for ensembles A40.20, A40.24 and A40.32, which differ only in the volume. In Fig. 2 we plot $\delta E$ as a function of $1/L$ for these three ensembles and $a\mu_{s}=0.0185$. The other two $\mu_{s}$-values give similar results. We have solved Eq. 17 including all terms up to order $1/L^{6}$ for $a_{0}$ and $r_{f}$ using A40.24 and A40.32 only, the result of which is shown as the solid line with errorband in the left panel of the figure. It leads to $M_{K}a_{0}=-0.292(20)$. Including also A40.20, we perform a two parameter fit with three data points finding $M_{K}a_{0}=-0.318(9)$. The corresponding fit is shown in the right panel of the figure. Leaving out the effective range term at order $1/L^{6}$ results in unreasonably large $\chi^{2}$-values.

Noting that solving Eq. 17 up to order $1/L^{5}$ for $a_{0}$ for ensemble A40.32 gives $M_{K}a_{0}=-0.315(11)$, which agrees within error with the two estimates from above, we conclude that $L/a=32$ is sufficiently large, while $L/a=24$ is at the border. $L/a=20$ is certainly too small to extract $M_{K}a_{0}$ from a single volume neglecting the effective range term.

We checked the impact of the inclusion of $r_{f}$ on the extraction of $M_{K}a_{0}$ in Appendix A. With a LO ChPT estimation of $r_{f}$ included in the extraction of $a_{0}$ the values for $M_{K}a_{0}$ vary by about one standard deviation. The central values for the $L/a=24$ lattices change by about 1% on the inclusion of the order $1/L^{6}$ terms (cf. Table 9). Thus we attribute a conservatively estimated systematic uncertainty of 1% to our chiral and continuum extrapolated value of $M_{K}a_{0}$.

4.2 Scattering Length

Given the values of $a\delta E$ and $aM_{K}$, the scattering length $a_{0}$ is determined using Eq. 17.

The number of fit ranges for extracting $a\delta E$ is low, compared to the $\pi\pi$-case of Ref. [11]. Thus an estimate of the systematic effects stemming from the fitting procedure is likely to be incorrect. Therefore, instead of estimating the systematic uncertainty introduced by the fitting procedure after the chiral extrapolations we consider the $p$-value weighted median over the fitranges. This procedure is further supported by the fact that the statistical uncertainties of $M_{K}a_{0}$ do essentially not differ from the uncertainties obtained by adding statistical and systematic uncertainties in quadrature. The final results for $aM_{K}$, $a\delta E$, $a_{0}/a$ and $M_{K}a_{0}$ are compiled in Tables 16–18 for all ensembles.

4.3 Strategies M1A and M2A: $M_{K}a_{0}$ from fixed $M_{s}^{2}$

To evaluate $M_{K}a_{0}$ at the physical strange quark mass, we convert $M_{s}^{2}$ to lattice units using $r_{0}/a$ listed in Table 2. First, we express $M_{s}^{2}$ in units of $r_{0}$ using the estimates in Eq. 1, which gives $(r^{\mathrm{M1}}_{0}M_{s}^{\mathrm{phys}})^{2}=1.33(7)$ with $Z_{P}$ from M1 and $(r^{\mathrm{M2}}_{0}M_{s}^{\mathrm{phys}})^{2}=1.34(6)$ with $Z_{P}$ from M2. In lattice units at our three lattice spacings, these correspond to the values given in Table 5.

For each ensemble, we then interpolate $M_{K}a_{0}$ by performing a correlated linear fit to the data at the three values of $a\mu_{s}$ (the independent variable being $a^{2}M_{s}^{2}$). An example of this is given in Fig. 6 in Appendix B.

Having interpolated $M_{K}a_{0}$ on all ensembles, the data is extrapolated to the physical point and to the continuum in a global fit using Eq. 24. In Fig. 3 the dimensionless product $M_{K}a_{0}$ is shown as a function of $r_{0}m_{l}$ together with the global fit for each value of $\beta$ for M1A in the left and for M2B in the right panel, respectively. Note that we take into account all correlation between data which enters through the procedure for fixing the strange quark mass at each value of the lattice spacing. The results of the fits can be found in Table 6.

4.4 Strategies M1B and M2B: $M_{K}a_{0}$ from fixed $m_{s}^{R}$

Analysis B involves as a first step a global fit of Eq. 19 to the values of $aM_{K}$. As an example the fits to the data of the A ensembles are shown in Appendix B in Figure 7 for $Z_{P}$ from M1 (left panel) and M2 (right panel), respectively.

The fit takes into account the correlation between data at different values of $a\mu_{s}$ but the same $a\mu_{\ell}$-value. The results of the global fits are compiled in Tables 14 and 15. The fitted parameters allow us to calculate the renormalised strange quark mass, $m_{s}^{\mathrm{phys}}$, from Eq. 22. As input, we use $r_{0}$ from Eq. 1, $Z_{P}$ from Table 2, $m_{l}^{\mathrm{phys}}$ and $M_{K}^{\mathrm{phys}}$.

For the physical values of the strange quark mass at $2\ \mathrm{GeV}$ in the $\overline{\textrm{MS}}$ scheme, we find

[TABLE]

These values compare well to the corresponding results from Ref. [17]:

[TABLE]

We can convert the values from Eq. 26 to lattice units for the three $\beta$-values, which we compiled in Table 7. Next we interpolate $M_{K}a_{0}$ in $a\mu_{s}$ to these values for all ensembles. As an example we show the linear correlated fit for ensemble B55.32 in Appendix B in Fig. 8 for M1B in the left and M2B in the right panel.

Interpolated to the reference strange quark mass, the values of $M_{K}a_{0}$ are shown as a function of the renormalised light quark mass in Fig. 4 in units of $r_{0}$. We also show the best fit function for each $\beta$-value and the continuum extrapolation. The continuum extrapolated values at the physical point $(M_{k}a_{0})^{\mathrm{phys}}$ are indicated by the diamonds. Note again that due to the strange quark mass fixing procedure all points for a single lattice spacing are correlated.

In Table 6 we give our final results for $(M_{k}a_{0})^{\mathrm{phys}}$ for the four different approaches M1A, M2A, M1B and M2B together with the best fit parameters $Q_{1,2,3}$, the $\chi^{2}/\mathrm{dof}$- and the $p$-value of the fit.

As the final result we quote the $p$-value weighted median over the four determinations

[TABLE]

The statistical uncertainty comes from the bootstrap procedure. The systematic uncertainty coming from the two methods to estimate $Z_{P}$ is estimated as follows: we first compute the weighted average of only M1A and M1B and also of only M2A and M2B. The systematic uncertainty is then taken as the deviation between these two weighted averages and the final result, Eq. 27. For the systematic uncertainty from setting the strange quark mass we proceed in the same way, just that we first compute the weighted average of only M1A and M2A and also of only M1B and M2B. As the last error we quote the systematic uncertainty from neglecting higher order terms in the calculation of the scattering length. Using $M_{K}^{\textrm{phys}}$ we obtain for the scattering length

[TABLE]

5 Discussion

We have used four methods to determine $M_{K}a_{0}$ at the physical light and strange quark mass value in the continuum limit. The differences between these methods are lattice artefacts. From Table 6 it becomes clear that all four methods give results which are well compatible within statistical uncertainties. This gives us confidence in our procedure and in our final result Eq. 27. The four different estimates can still serve as an estimate of systematic effects, which are, however, smaller than the statistical uncertainty of about 4%. The largest fraction of this statistical uncertainty stems from the uncertainty in the scale.

It turns out that lattice artefacts are not negligible in $M_{K}a_{0}$: from $\beta=1.90$ to the continuum a roughly 20% relative change in the result is observed. From our finest lattice spacing we still see a change of about 8%. It is interesting to note that our central value equals within errors to the LO ChPT estimate

[TABLE]

A possibly still uncontrolled systematic uncertainty could come from our chiral and continuum extrapolation. In lattice ChPT, usually the $a^{2}$ term is taken to be of higher order than the term linear in $\mu_{\ell}$. For this we would need to include higher orders in the quark mass as well. However, the precision in our data is not sufficient to resolve such terms. But the need for the $a^{2}$ term is evident. Therefore, we decided to stick to a power counting with $a^{2}\propto\mu_{\ell}$. An alternative and probably better chiral representation of $M_{k}a_{0}$ in terms of $M_{K}/f_{K}$ was used in Ref. [3]. This representation turned out to be not feasible for us, because we have only very little spread in $M_{K}/f_{K}$. Smaller uncertainties on $M_{K}a_{0}$ might enable the investigation of the light and strange quark mass dependence using mixed action ChPT at NLO.

We also cannot estimate the effects from partial quenching of the strange quark. However, it should be noted that the kaon masses that we obtain in the OS valence sector at the physical strange quark mass, as set via either method A or B, deviate from those of the unitary kaon mass published in Ref. [32] by a few percent at most. Partial quenching effects in analyses using OS valence fermions on a $N_{f}=2+1+1$ twisted mass sea have been shown to be small for other observables in the past. Moreover, we would like to remark that the dependence of $M_{K}a_{0}$ on $\mu_{s}$ is not very pronounced. Finally, our estimate in Appendix A indicates that the $\mathcal{O}(L^{-6})$-terms in the Lüscher formula Eq. 17 are indeed negligible for our case. Nevertheless we do not have a sufficient number of volumes available to determine it from the data.

Two other lattice calculations of $M_{K}a_{0}$ are available. The NPLQCD collaboration used three-flavour mixed action ChPT to obtain $M_{K}a_{0}=-0.352(16)$, with statistical and systematic uncertainties combined in quadrature [3]. They worked with domain wall valence quarks on a sea of $N_{f}=2+1$ asqtad-improved rooted staggered quarks. A second calculation was performed by the authors of Ref. [4] with $N_{f}=2+1$ dynamical flavours of non-perturbatively $O(a)$-improved Wilson quarks. Their result reads $M_{K}a_{0}=-0.310(17)(32)$. The discrepancy between these determinations and our final result, Eq. 27, is quite substantial. In the NPLQCD determination predominantly one lattice spacing of $a=0.125\ \mathrm{fm}$ was considered in the chiral extrapolation. One ensemble with a finer lattice spacing was included in the analysis to attempt a quantification of discretisation errors, but it should be noted that the uncertainty on this point was about a factor of three larger than on all other points in the analysis. The PACS-CS collaboration used only one lattice spacing value with $a\sim 0.09\ \mathrm{fm}$, very close to our coarsest lattice spacing value. PACS-CS included one ensemble with $M_{\pi}=170\ \mathrm{MeV}$ in their analysis, which is, however, giving very noisy results. Both collaborations use one strange quark mass value which was tuned to be close to physical.

In Fig. 5 we compare our result at the coarsest lattice spacing, i.e. the A-ensembles, interpolated to the physical strange quark mass with method M1A to the results of the other two collaborations. There is no obvious conclusion from this comparison. But the errors of the PACS-CS results appear to be large enough to explain the observed differences, given the fact that the PACS-CS result is at one lattice spacing value only. The comparison to the NPLQCD data points is more difficult, in particular since the one NPLQCD point with a finer lattice spacing points towards an even smaller absolute value for $M_{K}a_{0}$, though with a large statistical uncertainty. This can only be resolved with continuum extrapolations for the other formulations. However, the NPLQCD and our result agree within two standard deviations.

6 Summary

We investigated the scattering length of the $K^{+}$-$K^{+}$ system by means of finite volume methods for lattice QCD devised by M. Lüscher. The lattice formulation is Wilson twisted mass lattice QCD at maximal twist and $N_{f}=2+1+1$ dynamical quark flavours. The gauge configurations, involving 11 pion masses at 3 different lattice spacings, were generated by the ETMC. To the author’s knowledge our result represents the first study of the $K^{+}$-$K^{+}$ system controlling lattice artefacts using three lattice spacing values and up/down, strange and charm dynamical quarks. For the strange quark we used a mixed action approach with so-called Osterwalder-Seiler valence strange quarks to be able to correct for a slight mistuning of the sea strange quark mass value.

In total, we followed four different strategies to arrive at the continuum extrapolated value for $M_{K}a_{0}$ at physical light and strange quark masses. All four show very good agreement indicating that the corresponding extrapolations are well controlled. Our final result for the scattering length is

[TABLE]

from the weighted median over the four strategies. In our calculation we find that the continuum extrapolation is vital in obtaining the final number: from the coarsest to the continuum result we observe a roughly 20% difference. We think that this is also the reason for the discrepancy we observe when comparing to the two previous lattice calculations of $M_{K}a_{0}$, because for the other two results a continuum extrapolation could not be performed.

In the near future we will extend the analysis performed here to the pion-kaon case.

Acknowledgements

We thank the members of ETMC for the most enjoyable collaboration. The computer time for this project was made available to us by the John von Neumann-Institute for Computing (NIC) on the Jureca and Juqueen systems in Jülich. We thank A. Rusetsky for very useful discussions. We thank S. Simula for the estimates of the finite size corrections to $M_{\pi}$ and $M_{K}$. This project was funded by the DFG as a project in the Sino-German CRC110. The open source software packages tmLQCD [33, 34, 35], Lemon [36], Eigen [37], Boost [38], SciPy [39] and R [40] have been used. In addition we employed QUDA [41, 42] for calculating propagators on GPUs.

Appendix A Effective Range from ChPT

We start from the partial wave expansion for the scattering amplitude $T^{I}(s,t,u)$ [43]

[TABLE]

which depends on the Legendre polynomials $P_{\ell}(\cos\vartheta)$, and the partial wave amplitudes $t^{I}_{\ell}(s)$. The amplitudes $t^{I}_{\ell}(s)$ can be expanded in terms of the scattering momentum $q$ and the slope parameters:

[TABLE]

Since we are interested in maximal isospin and the $s$-wave, we take $I=1$ and $\ell=0$. This yields

[TABLE]

In Ref. [43] $T(s,t,u)$ for $K^{+}K^{-}\rightarrow K^{+}K^{-}$ is given to leading order by:

[TABLE]

To turn this into an amplitude valid for $K^{+}K^{+}$-scattering we employ crossing symmetry which interchanges the Mandelstam variables $s$ and $u$. With that the partial wave amplitude becomes

[TABLE]

where we expressed $s$ with the momentum transfer $q$: $s=4(M_{K}^{2}+q^{2})$. Expanding Eq. 33 in a Taylor series gives:

[TABLE]

Comparing Eq. 34 with Eq. 30 we can extract $b^{1}_{0}$ and use $r_{f}=-2M_{K}b^{1}_{0}$ to get

[TABLE]

To estimate the effective range we use the physical value of the kaon mass $M_{K}=$494.2\text{,}\mathrm{MeV} and the ChPT value $f_{\pi}=$94.2\text{\,}\mathrm{MeV}. Converting to a length unit with $\hbar c=$197.37\text{,}\mathrm{MeV}\text{,}\mathrm{fm}$$ gives

[TABLE]

We can use this to estimate the influence of the $\mathcal{O}(L^{-6})$-terms on the determination of the scattering length $a_{0}$ from Lüscher’s formula. To this end we compare the results for the scattering length up to order $\mathcal{O}(L^{-5})$ to the ones of up to order $\mathcal{O}(L^{-6})$ with and without the term involving the scattering length. Table 8 gives an overview over these differences .

For converting $r_{f}$ back to lattice units we use parametric bootstrap-samples of the lattice spacing $a$. In Table 9 the results for $M_{K}a_{0}$ for $a_{0}$ up to $\mathcal{O}(L^{-6})$ and $a_{0}$ truncated at $\mathcal{O}(L^{-5})$ are compared.

As visible from the table the inclusion of the terms to order $\mathcal{O}(L^{-6})$ in the determination of the scattering length does not change the values of $M_{K}a_{0}$ beyond one standard deviation.

Appendix B Data Tables and Plots