Quark mass dependence of the low-lying charmed mesons at one loop in HH$\chi$PT

TL;DR

This paper investigates how light and heavy quark masses influence low-lying charmed mesons using one-loop heavy-hadron chiral perturbation theory, fitting lattice QCD data to determine relevant constants.

Contribution

It introduces a systematic approach to extract low energy constants from lattice data and discusses discretization effects related to charm quark mass.

Findings

Determined low energy constants from lattice QCD data.

Analyzed discretization effects due to charm quark mass.

Enhanced precision in modeling charmed meson mass dependence.

Abstract

We study the light and heavy quark mass dependence of the low-lying charmed mesons in the framework of one-loop HH$\chi$PT. The low energy constants are determined by analyzing the available lattice data from different LQCD simulations. Model selection tools are implemented to determine the relevant parameters as required by data with a higher precision. Discretization and other effects due to the charm quark mass setting are discussed.

Tables

Table 1. Table 1: Comparison of the charm quark mass settings of the different collaborations through their value of m D s subscript 𝑚 subscript 𝐷 𝑠 m_{D_{s}} . Note that the physical value is m D s = 1968.35 ± 0.07 subscript 𝑚 subscript 𝐷 𝑠 plus-or-minus 1968.35 0.07 m_{D_{s}}=1968.35\pm 0.07  MeV Zyla et al. ( 2020 ) . L 𝐿 L and a 𝑎 a are given in fm, and the masses in units of MeV.

Col.	$a$	$L$	$m_{\pi }$	$m_{D_s}$
ETMC1	0.0619	3.0	224	1988
ETMC2	0.0619	3.0	224	2001
PACS-CS	0.0907	2.9	156	1809
HSC	0.12	3.8	239	1967
CLS	0.08636	2.1	280	1981
Prelovsek et al.	0.1239	2.0	266	1657
RQCD	0.071	4.5	150	1977

Table 2. Table 2: Collection of parameters from different fits. The subindex (I) means ”before LASSO” and (II) means ”after LASSO”.

Fit	$\frac{\sigma m_{\pi }}{B_0^{}{}_{}{}^2}\cdot {10}^{-4}$	$\frac{a{m}_{\pi }}{B_0^{}{}_{}{}^2}\cdot {10}^{-4}$	$\frac{b{m}_{\pi }^3}{B_0^{}{}_{}{}^2}\cdot {10}^{-10}$	$\frac{c{m}_{\pi }^3}{B_0^{}{}_{}{}^2}\cdot {10}^{-10}$	$\frac{d{m}_{\pi }^3}{B_0^{}{}_{}{}^2}\cdot {10}^{-10}$	$\frac{{\Delta }^{\left(\sigma \right)}{m}_{\pi }}{B_0}\cdot {10}^{-5}$	$\frac{{\Delta }^{\left(a\right)}{m}_{\pi }}{B_0}\cdot {10}^{-5}$
ETMC${}_{\text{(I)}}$	$0.89\pm 0.08$	$2.08\pm 0.13$	$-0.55\pm 0.21$	$6.51\pm 0.20$	$3.39\pm 0.18$	$5.78\pm 3.43$	$-2.16\pm 2.27$
ETMC${}_{\text{(II)}}$	$1.02\pm 0.01$	$2.02\pm 1.18$	-	$6.16\pm 0.24$	$3.24\pm 0.27$	$5.88\pm 3.85$	$-2.17\pm 2.24$
MILC${}_{\text{(I)}}$	$3\pm 22$	$3\pm 18$	$7\pm 108$	$2\pm 55$	$1.56\pm 7.08$	$4\pm 255$	$36\pm 478$
MILC${}_{\text{(II)}}$	$1.74\pm 0.40$	$2.40\pm 0.73$	-	$5.61\pm 1.66$	-	-	-
PACS${}_{\text{(I)}}$	$0.90\pm 0.24$	$2.62\pm 0.20$	$-2.38\pm 1.81$	$7.03\pm 1.24$	$3.53\pm 0.74$	$3.60\pm 0.16$	$-0.76\pm 0.44$
PACS${}_{\text{(II)}}$	$1.22\pm 0.02$	$2.36\pm 0.03$	-	$5.40\pm 0.03$	$4.50\pm 0.03$	$3.59\pm 0.16$	$-0.75\pm 0.44$
GLOBAL${}_{\text{(I)}}$	$1.84\pm 0.01$	$2.54\pm 0.02$	$0.57\pm 0.05$	$4.47\pm 0.04$	$2.34\pm 0.03$	$2.70\pm 0.21$	$-1.43\pm 0.47$
GLOBAL${}_{\text{(II)}}$	$1.17\pm 0.01$	$2.58\pm 0.02$	-	$4.89\pm 0.05$	$2.27\pm 0.03$	$2.56\pm 0.20$	$-1.47\pm 0.46$

Table 3. Table 3: Tree level parameters of the ETMC fits. The remaining parameters are shown in the first two rows of Table  2 .

Fit		$m_H$ (MeV)		$\Delta$ (MeV)
Col.	$a$ (fm)	(I)	(II)	(I)	(II)
ETMC1	$0.0619$	$1995\pm 12$	$1988\pm 16$	$122\pm 21$	$121\pm 24$
	$0.0815$	$2014\pm 13$	$2007\pm 17$	$129\pm 25$	$128\pm 27$
	$0.0885$	$2015\pm 14$	$2008\pm 18$	$132\pm 28$	$131\pm 29$
ETMC2	$0.0619$	$2008\pm 12$	$2001\pm 16$	$123\pm 21$	$123\pm 23$
	$0.0815$	$2038\pm 13$	$2031\pm 16$	$126\pm 26$	$126\pm 27$
	$0.0885$	$2047\pm 14$	$2041\pm 18$	$131\pm 25$	$130\pm 27$
ETMC	0	$1972$	$1965$	$114$	$114$

Table 4. Table 4: Tree level parameters of the MILC fit. The remaining parameters are shown in the first two rows of Table  2 .

Collaboration	$m_H$ (MeV)		$\Delta$ (MeV)
	(I)	(II)	(I)	(II)
MILC	$1863\pm 1271$	$1959\pm 40$	$110$	$110$

Table 5. Table 5: Tree level parameters of the PACS fit. The remaining parameters are shown in the first two rows of Table  2 .

Collaboration	$m_H$ (MeV)		$\Delta$ (MeV)
	(I)	(II)	(I)	(II)
PACS	$1767\pm 47$	$1706\pm 3$	$99\pm 2$	$99\pm 2$

Table 6. Table 6: Tree level parameters of the global fit. The lattice spacing with a ∗ , corresponding to HSC, are the temporal lattice spacing and they are shown in units of GeV -1 . The remaining parameters are shown in the first two rows of Table  2 .

Collaboration		$m_H$ (MeV)		$\Delta$ (MeV)
Label	$a$ (fm)	(I)	(II)	(I)	(II)
ETMC1	$0.0619$	$1937\pm 4$	$1945\pm 4$	$140\pm 7$	$141\pm 7$
	$0.0815$	$1949\pm 5$	$1956\pm 5$	$148\pm 14$	$148\pm 14$
	$0.0885$	$1955\pm 5$	$1963\pm 5$	$152\pm 12$	$153\pm 12$
ETMC2	$0.0619$	$1950\pm 3$	$1958\pm 3$	$139\pm 7$	$140\pm 7$
	$0.0815$	$1973\pm 5$	$1980\pm 5$	$145\pm 13$	$146\pm 13$
	$0.0885$	$1988\pm 4$	$1996\pm 4$	$148\pm 13$	$149\pm 13$
ETMC	0	$1916$	$1925$	$131$	$132$
HSC	$1/{6079}^{\ast }$	$1903\pm 5$	$1910\pm 5$	$99\pm 7$	$100\pm 7$
	$1/{5667}^{\ast }$	$1874\pm 5$	$1881\pm 5$	$91\pm 8$	$92\pm 8$
CLS	$0.08636$	$1940\pm 1$	$1945\pm 1$	$95\pm 2$	$96\pm 2$
Prelovsek	$0.1239$	$1574\pm 2$	$1582\pm 2$	$97\pm 2$	$98\pm 2$
RQCD	$0.0714$	$1910\pm 1$	$1917\pm 1$	$98\pm 1$	$98\pm 1$
MILC	$0$	$1913\pm 4$	$1920\pm 4$	$111\pm 4$	$112\pm 4$

Table 7. Table 7: Results of the analyses with the tree level formulas, Eqs. ( 33 ) and ( 34 ). The remaining parameters are shown in Table  8

Fit	$\frac{\sigma m_{\pi }}{B_0^{}{}_{}{}^2}$$\cdot {10}^{-4}$	$\frac{a{m}_{\pi }}{B_0^{}{}_{}{}^2}$$\cdot {10}^{-4}$	$\frac{{\Delta }^{\left(\sigma \right)}{m}_{\pi }}{B_0}$$\cdot {10}^{-5}$	$\frac{{\Delta }^{\left(a\right)}{m}_{\pi }}{B_0}$$\cdot {10}^{-5}$
ETMC	$0.16\pm 0.17$	$2.03\pm 0.07$	$6.22\pm 4.78$	$-2.02\pm 2.24$
MILC	$0.83\pm 0.26$	$2.45\pm 0.37$	$6.07\pm 1.75$	$3.75\pm 4.47$
GLOBAL	$0.69\pm 0.01$	$1.98\pm 0.03$	$0.48\pm 0.23$	$-1.88\pm 0.47$

Table 8. Table 8: Results of the analyses with the tree level formulas, Eqs. ( 33 ) and ( 34 ). The remaining parameters are shown in Table  7 .

Fit	Label	$a$ (fm)	$m_H$ (MeV)	$\Delta$ (MeV)
ETMC	ETMC1	$0.0619$	$1994\pm 12$	$137\pm 34$
		$0.0815$	$2010\pm 14$	$145\pm 39$
		$0.0885$	$2014\pm 14$	$147\pm 42$
	ETMC2	$0.0619$	$2006\pm 12$	$138\pm 33$
		$0.0815$	$2034\pm 13$	$143\pm 39$
		$0.0885$	$2047\pm 14$	$146\pm 39$
	ETMC	$0$	$1971$	$129$
MILC	MILC	$0$	$1952\pm 21$	$140$
GLOBAL	ETMC1	$0.0619$	$1959\pm 3$	$178\pm 8$
		$0.0815$	$1972\pm 5$	$188\pm 16$
		$0.0885$	$1974\pm 4$	$194\pm 14$
	ETMC2	$0.0619$	$1971\pm 3$	$177\pm 8$
		$0.0815$	$1996\pm 5$	$185\pm 14$
		$0.0885$	$2007\pm 4$	$190\pm 15$
	HSC	$1/6079$	$1929\pm 5$	$129\pm 8$
		$1/5667$	$1899\pm 5$	$124\pm 9$
	CLS	$0.08636$	$1966\pm 1$	$124\pm 2$
	Prelovsek	$0.1239$	$1596\pm 1$	$128\pm 3$
	RQCD	$0.0714$	$1934\pm 1$	$127\pm 2$
	MILC	$0$	$1938\pm 3$	$144\pm 5$

Table 9. Table 9: D ( s ) ( ∗ ) subscript superscript 𝐷 𝑠 D^{(*)}_{(s)} and J / Ψ 𝐽 Ψ J/\Psi meson masses. Data collected from Guo et al. ( 2018b ) provided by the authors of Kalinowski and Wagner ( 2015 ) . The charm quark mass is fixed by reproducing the physical J / ψ 𝐽 𝜓 J/\psi meson mass. Lattice spacing uncertainties are discused on Carrasco et al. ( 2014 ) . N F = 2 + 1 + 1 subscript 𝑁 𝐹 2 1 1 N_{F}=2+1+1 .

$a$ (fm)	${\left(L/a\right)}^3\times T/a$	$a{m}_{\pi }$	$a{m}_K$	$a{m}_D$	$a{m}_{D_s}$	$a{m}_{D^{\ast }}$	$a{m}_{D_s^{\ast }}$	$a{m}_{J/\Psi }$
$0.0619\left(18\right)$	${48}^3\times 96$	0.0703	0.1697	0.5905	0.6236	0.6466	0.6770	0.9715
	${48}^3\times 96$	0.0806	0.1738	0.5906	0.6234	0.6506	0.6763	0.9697
	${48}^3\times 96$	0.0975	0.1768	0.5913	0.6229	0.6486	0.6764	0.9703
$0.0815\left(30\right)$	${32}^3\times 64$	0.1074	0.2133	0.7840	0.8159	0.8568	0.8905	1.2791
	${32}^3\times 64$	0.1549	0.2279	0.7895	0.8183	0.8678	0.8950	1.2828
	${24}^3\times 48$	0.1935	0.2430	0.7934	0.8175	0.8745	0.8965	1.2818
$0.0885\left(36\right)$	${32}^3\times 64$	0.1240	0.2512	0.8514	0.8953	0.9356	0.9806	1.3890
	${32}^3\times 64$	0.1412	0.2569	0.8544	0.8972	0.9363	0.9802	1.3895
	${24}^3\times 48$	0.1440	0.2589	0.8552	0.8978	0.9403	0.9844	1.3906
	${24}^3\times 48$	0.1988	0.2764	0.8599	0.8950	0.9487	0.9841	1.3882

Table 10. Table 10: The same as Table 9 but with D ( s ) ( ∗ ) subscript superscript 𝐷 𝑠 D^{(*)}_{(s)} and η c subscript 𝜂 𝑐 \eta_{c} meson masses. The charm quark mass is fixed by reproducing the physical η c subscript 𝜂 𝑐 \eta_{c} meson mass.

$a$ (fm)	${\left(L/a\right)}^3\times T/a$	$a{m}_{\pi }$	$a{m}_K$	$a{m}_D$	$a{m}_{D_s}$	$a{m}_{D^{\ast }}$	$a{m}_{D_s^{\ast }}$	$a{m}_{{\eta }_c}$
$0.0619\left(18\right)$	${48}^3\times 96$	0.0703	0.1697	0.5947	0.6279	0.6506	0.6809	0.9351
	${48}^3\times 96$	0.0806	0.1738	0.5949	0.6277	0.6546	0.6803	0.9332
	${48}^3\times 96$	0.0975	0.1768	0.5955	0.6271	0.6526	0.6804	0.9335
$0.0815\left(30\right)$	${32}^3\times 64$	0.1074	0.2133	0.7946	0.8263	0.8664	0.9001	1.2312
	${32}^3\times 64$	0.1549	0.2279	0.8004	0.8291	0.8777	0.9049	1.2342
	${24}^3\times 48$	0.1935	0.2430	0.8039	0.8278	0.8840	0.9059	1.2314
$0.0885\left(36\right)$	${32}^3\times 64$	0.1240	0.2512	0.8677	0.9114	0.9506	0.9953	1.3370
	${32}^3\times 64$	0.1412	0.2569	0.8708	0.9132	0.9511	0.9949	1.3379
	${24}^3\times 48$	0.1440	0.2589	0.8714	0.9137	0.9545	0.9990	1.3382
	${24}^3\times 48$	0.1988	0.2764	0.8753	0.9102	0.9627	0.9980	1.3325

Table 11. Table 11: Data for heavy meson masses from Mohler and Woloshyn ( 2011 ) and light meson masses collected from Aoki et al. ( 2009 ) [PACS-CS]. The charm quark hopping pa rameter κ c subscript 𝜅 𝑐 \kappa_{c} has been tuned to the value where the spin-averaged kinetic mass ( M D s + 3 ​ M D s ∗ ) / 4 subscript 𝑀 subscript 𝐷 𝑠 3 subscript 𝑀 superscript subscript 𝐷 𝑠 4 (M_{D_{s}}+3M_{D_{s}^{*}})/4 takes its physical value. The lowest pion mass data is also named as Ensemble 2 in Lang et al. ( 2014 ) . Also, ( L / a ) 3 × T / a = 32 3 × 64 superscript 𝐿 𝑎 3 𝑇 𝑎 superscript 32 3 64 (L/a)^{3}\times T/a=32^{3}\times 64 and N F = 2 + 1 subscript 𝑁 𝐹 2 1 N_{F}=2+1 .

$a$ (fm)	$a{m}_{\pi }$	$a{m}_K$	$a{m}_D$	$a{m}_{D_s}$	$a{m}_{D^{\ast }}$	$a{m}_{D_s^{\ast }}$	$a{m}_{{\eta }_c}$	$a{m}_{J/\Psi }$
$a=0.0907\left(13\right)$	0.32242	0.36269	0.84022	0.86273	0.91237	0.93297	1.22924	1.27877
	0.26191	0.32785	-	0.85158	-	0.91864	1.22098	1.26860
	0.18903	0.29190	0.79580	0.84000	0.86327	0.90429	1.21369	1.25988
	0.13593	0.27282	0.78798	0.83929	0.85776	0.90429	1.21268	1.25904
	0.07162	0.25454	0.77646	0.83149	0.83656	0.89268	1.20765	1.25258

Table 12. Table 12: Data collected from Cheung et al. ( 2021 ) for N F = 2 + 1 subscript 𝑁 𝐹 2 1 N_{F}=2+1 [HSC]. The charm quark mass parameters have been determined using the η c subscript 𝜂 𝑐 \eta_{c} physical mass, obtaining the values, ( a t ​ m η c ) 2 = 0.2412 superscript subscript 𝑎 𝑡 subscript 𝑚 subscript 𝜂 𝑐 2 0.2412 (a_{t}m_{\eta_{c}})^{2}=0.2412 and 0.2735 0.2735 0.2735 for the light and heavy pion mass respectively (see Figs. 1 of Cheung et al. ( 2016 ); Liu et al. ( 2012 ) ).

$a_t^{-1}$ (MeV)	${\left(L/a_s\right)}^3\times \left(T/a_t\right)$	$m_{\pi }$	$m_K$	$m_D$	$m_{D_s}$	$m_{D^{\ast }}$	$m_{D_s^{\ast }}$
6079	${32}^3\times 256$	0.03928(18)	0.08344(7)	0.30923(11)	0.32356(12)	0.33058(24)	0.34448(15)
5667	${(16,20,24)}^3\times 128$	0.06906(13)	0.09698(9)	0.33303(31)	0.34441(29)	0.35494(46)	0.36587(35)

Table 13. Table 13: Data collected from Prelovsek et al. ( 2021 ) where the spin average M a ​ v ​ g = ( m η c + 3 ​ m J / Ψ ) / 4 subscript 𝑀 𝑎 𝑣 𝑔 subscript 𝑚 subscript 𝜂 𝑐 3 subscript 𝑚 𝐽 Ψ 4 M_{avg}=(m_{\eta_{c}}+3m_{J/\Psi})/4 has been tune to the physical point, for ( L / a ) 3 × T / a = 24 3 × 128 superscript 𝐿 𝑎 3 𝑇 𝑎 superscript 24 3 128 (L/a)^{3}\times T/a=24^{3}\times 128 and N F = 2 + 1 subscript 𝑁 𝐹 2 1 N_{F}=2+1 [CLS]. Masses given in MeV.

$a$ (fm)	$m_{\pi }$	$m_K$	$m_D$	$m_{D_s}$	$m_{D^{\ast }}$	$m_{D_s^{\ast }}$	$M_{avg}$
$a=0.08636\left(138\right)$	280	467	1927	1981	2048	2102	3103

Table 14. Table 14: Data collected from Ensemble 1 in Lang et al. ( 2014 ) for ( L / a ) 3 × T / a = 16 3 × 32 superscript 𝐿 𝑎 3 𝑇 𝑎 superscript 16 3 32 (L/a)^{3}\times T/a=16^{3}\times 32 and N F = 2 subscript 𝑁 𝐹 2 N_{F}=2 [Prelovsek et al.]. The spin-average ( m η c + 3 ​ m J / ψ ) / 4 subscript 𝑚 subscript 𝜂 𝑐 3 subscript 𝑚 𝐽 𝜓 4 (m_{\eta_{c}}+3m_{J/\psi})/4 is tuned to reproduce its physical value.

$a$ (fm)	$a{m}_{\pi }$	$a{m}_K$	$a{m}_D$	$a{m}_{D_s}$	$a{m}_{D^{\ast }}$	$a{m}_{D_s^{\ast }}$	$a{M}_{avg}$
$a=0.1239\left(13\right)$	0.1673	0.3467	0.9801	1.04075433	1.0629	1.11635248	1.52499

Table 15. Table 15: Data collected from table I in Bali et al. ( 2017 ) for N F = 2 subscript 𝑁 𝐹 2 N_{F}=2 [RQCD]. The charm quark mass is fixed through the experimental value of the spin-averaged 1 ​ S 1 𝑆 1S charmonium mass, m 1 ​ S = 3068.5 subscript 𝑚 1 𝑆 3068.5 m_{1S}=3068.5 MeV. Masses given in MeV.

$a$ (fm)	${\left(L/a\right)}^3\times T/a$	$m_{\pi }$	$m_K$	$m_D$	$m_{D_s}$	$m_{D^{\ast }}$	$m_{D_s^{\ast }}$
$a=0.071$	${24}^3\times 48$	306.9	540	1907	-	2038	-
	${32}^3\times 64$	294.6	528	1902	-	2030	-
	${40}^3\times 64$	288.8	527	1901	-	2030	-
	${64}^3\times 64$	289.5	526	1898	-	2030	-
	${48}^3\times 64$	159.7	500	1880	-	2007	-
	${64}^3\times 64$	149.7	497	1877	1976.9	1996	2094.9

Equations

$$H_l^{(Q)}=\frac{(1+/v)}{2}(P_{l\mu }^{\ast (Q)}{\gamma }^{\mu }-P_l^{(Q)}{\gamma }_5).$$

$${\bar{H}}_l^{(Q)}={\gamma }^0{H}_l^{(Q)†}{\gamma }^0=(P_{l\mu }^{\ast (Q)†}{\gamma }^{\mu }+P_l^{(Q)†}{\gamma }_5)\frac{(1+/v)}{2}.$$

$$H_l^{(Q)}\to S_Q(H^{(Q)}{U}^{†}{)}_l$$

$${\bar{H}}_l^{(Q)}\to (U{\bar{H}}^{(Q)}{)}_l{S}_Q^{†},$$

$$H_v^{(Q)}=\sqrt{m_H}{e}^{i{m}_H}/v{v}_{\mu }{x}^{\mu }{H}^{(Q)}(x),$$

$${\mathcal{L}}^0=$$

$$+2g\mathrm{Tr}{\bar{H}}_v{H}_v(S_{lv}\cdot A)$$

$$V^{\mu }=\frac{1}{2}(\xi {\partial }^{\mu }{\xi }^{†}+{\xi }^{†}{\partial }^{\mu }\xi ),A^{\mu }=\frac{i}{2}(\xi {\partial }^{\mu }{\xi }^{†}-{\xi }^{†}{\partial }^{\mu }\xi ),$$

$$\varphi =$$

$${\mathcal{L}}_v^{/s}=$$

$$-\frac{\Delta }{8}\mathrm{Tr}{\bar{H}}_v{\sigma }^{\mu \nu }{H}_v{\sigma }_{\mu \nu }-\frac{{\Delta }^{(a)}}{8}\mathrm{Tr}{\bar{H}}_v{\sigma }^{\mu \nu }{H}_v{\sigma }_{\mu \nu }{M}_{\xi }-\frac{{\Delta }^{(\sigma )}}{8}\mathrm{Tr}{M}_{\xi }\mathrm{Tr}{\bar{H}}_v{\sigma }^{\mu \nu }{H}_v{\sigma }_{\mu \nu },$$

$$\frac{1}{4}(M_{P_l}+3M_{P_l^{\ast }})$$

$$M_{P_l^{\ast }}-M_{P_l}$$

$$m_H=m_H^{\prime }+r_H\mathcal{O}(a^k)$$

$$\Delta ={\Delta }^{\prime }+r_{\Delta }\mathcal{O}(a^k),$$

$${\chi }^2(\vec{p},\vec{q})=\sum _i^N\frac{{\left(y_i-f_i(\vec{p},\vec{q})\right)}^2}{{\sigma }_{y_i}^2}+P,$$

$$P=\frac{\lambda }{10}\sum _{i=1}^7|p_i|.$$

$$\text{AIC}=2k-2\log (L)=2k+{\chi }^2,$$

$$\text{AICc}=\text{AIC}+\frac{2k(k+1)}{n-k-1},$$

$$\text{BIC}=k\log (n)-2\log (L)=k\log (n)+{\chi }^2,$$

$$m_H=m_H^{\prime }+r_H{a}^2$$

$$\Delta ={\Delta }^{\prime }+r_{\Delta }{a}^2,$$

$${\alpha }_l=2\sigma (2m+m_s)+2a(m{\delta }_{1,2l}+m_s{\delta }_{3l})$$

$$\displaystyle\beta^{(\pi)}=g^{2}\left({\begin{array}[]{c}\frac{3}{2}\\ \frac{3}{2}\\ 0\end{array}}\right);\,\beta^{(K)}=g^{2}\left({\begin{array}[]{c}1\\ 1\\ 2\end{array}}\right);\beta^{(\eta)}=g^{2}\left({\begin{array}[]{c}\frac{1}{6}\\ \frac{1}{6}\\ \frac{2}{3}\end{array}}\right)$$

$$({\gamma }_l^{(\pi )}-{\lambda }_l^{(\pi )}{\alpha }_l)=$$

$$({\gamma }_l^{(K)}-{\lambda }_l^{(K)}{\alpha }_l)=$$

$$+3g^{2}a(m_s-m)({\delta }_{1l}-{\delta }_{3l})$$

$$({\gamma }_l^{(\eta )}-{\lambda }_l^{(\eta )}{\alpha }_l)=$$

$$-\frac{\sigma }{3}(2m+4m_s)$$

$$({\gamma }_l^{(X)}-{\lambda }_l^{(X)}\Delta )=$$

$$c_l=$$

$$+4d(2m^2+m_s^2)$$

$$\delta c_l=$$

$$m=\frac{M_{\pi }^{02}}{2B_0},m_s=\frac{M_K^{02}}{B_0}-\frac{M_{\pi }^{02}}{2B_0}.$$

$$\frac{1}{4}(M_{P_l}$$

$$m_H+2\sigma (2m+m_s)+2a(m{\delta }_{1,2l}+m_s{\delta }_{3l}),$$

$$M_{P_l}^{\ast }$$

$$\Delta +2{\Delta }^{(\sigma )}(2m+m_s)+2{\Delta }^{(a)}(m{\delta }_{1,2l}+m_s{\delta }_{3l}).$$

Full text

I Introduction

The discovery of exotic hadrons in the heavy quark sector, which cannot be accomodated in terms of $q\bar{q}$ mesons or $qqq$ baryons, like those with a tetraquark or pentaquark structure, has manifested the relevance of hadronic loops in order to explain the masses and other properties of many states in the hadron spectrum Guo et al. (2018a).

Even there has been claims since long time of such effects for some low-lying states. Such is the case, for example, of the scalar mesons. An approach for (pseudoscalar) meson - (pseudoscalar) meson scattering based on chiral symmetry and unitarity, the so-called Unitarized Chiral Perturbation Theory (UChPT) Truong (1988); Dobado and Pelaez (1993, 1997); Nieves and Ruiz Arriola (2000); Garcia-Recio et al. (2004); Oller et al. (1999); Gomez Nicola and Pelaez (2002) has predicted the properties of the $\sigma$, $a_{0}$ and $f_{0}$, being in reasonably good agreement with experiment Oller and Oset (1997); Oller et al. (1998); Kaiser (1998); Locher et al. (1998); Nieves and Ruiz Arriola (2000). The study of the large $N_{c}$ behaviour concludes that, while the low-lying scalars are non-ordinary mesons, the vector mesons are dominated by the $q\bar{q}$ interaction Pelaez (2016). From the Lattice QCD (LQCD) side, however, both, $q\bar{q}$ and $\pi\pi$ operators are needed in order to obtain the properties of the $\rho$ meson Aoki et al. (2007); Gockeler et al. (2008); Feng et al. (2011); Lang et al. (2011); Pelissier and Alexandru (2013); Bali et al. (2016); Guo et al. (2016); Erben et al. (2020); Wilson et al. (2015); Dudek et al. (2013); Bulava et al. (2016); Feng et al. (2015); Alexandrou et al. (2017); Fu and Wang (2016); Metivet (2015); Fischer et al. (2020). The role of the coupling of the $K\bar{K}$ channel to $\pi\pi$ in the $\rho$-meson mass has been studied in Guo et al. (2016); Hu et al. (2016, 2017); Molina and Ruiz de Elvira (2020).

In the charm sector, it could be easier to dissentangle ordinary from non-ordinary mesons, since many of the observed particles are narrow states. In fact, a plethora of exotic states, for which hadron loops might play a relevant role, especially when these states are close to thresholds, have been discovered in the recent years Olsen et al. (2018); Chen et al. (2022). Unexpectedly, it has been found that, when hadron loops are incorporated in the quark model, mass shifts can be quite large even for the low-lying charmonium states Barnes and Swanson (2008).

In the heavy-light meson sector, heavy quark symmetry arises in the heavy-quark limit ($m_{Q}\to\infty$). This involves both Heavy Quark Spin Symmetry (HQSS) and Heavy Quark Flavor Symmetry (HQFS) Isgur and Wise (1989, 1990, 1991); Manohar and Wise (2000). Heavy quark spin symmetry emerges in this limit because the color magnetic moment of a heavy quark is proportional to the inverse of its mass, $\mu^{c}\propto 1/m_{Q}$, and thus the interaction becomes invariant under the SU(2) spin symmetry, i. e., for each $Q\bar{q}$ state there will be another one degenerated in mass which can be obtained by a flip of the spin of the heavy quark. The spin-parity of the heavy and light degrees of freedom are then separately conserved. Beyond that, the exact value of the heavy quark mass in the interactions with the light degrees of freedom plays no role and this gives rise to heavy flavor symmetry. Obviously, these symmetries are broken for finite quark masses and terms that go like $1/m_{Q}$ could be relevant for the charm quark.

Thus, in the heavy quark limit, one can use the light degrees of freedom (d.o.f) to classify the heavy-light meson states. The total angular momentum of the light d.o.f are $j_{l}=L\pm\frac{1}{2}$, being $L$ the angular momentum. In the ground state, two degenerate states with $J^{P}=\{0,1\}^{-}$ appear for $L=0$ ($j_{l}^{P}=\frac{1}{2}^{-}$), while for $L=1$ there are a couple of doublets, one with $J^{P}=\{0,1\}^{+}$ ($j_{l}^{P}=\frac{1}{2}^{+}$), and another for $J^{P}=\{1,2\}^{+}$ ($j_{l}^{P}=3/2^{+}$). Actually, one can find in the PDG the $D$ and $D^{*}$ mesons which correspond to the ground state doublet, while the $D^{*}_{0}(2300)$ and $D_{1}(2420)$ could be related to the $J^{P}=\{0,1\}^{+}$ doublet, and the $D_{1}(2430)$ and $D^{*}_{2}(2460)$ would correspond to the doublet with $J^{P}=\{1,2\}^{+}$.

Similarly, we find in the PDG their strange counterparts, the $(D_{s},D^{*}_{s})$, $(D_{s0}(2317)$, $D_{s1}(2460))$, and $(D_{s1}(2536)$, $D_{s2}(2573))$, which should correspond to the doublets with $J^{P}=\{0,1\}^{-}$, $\{0,1\}^{+}$ and $\{1,2\}^{+}$ respectively. Regarding these states, the constituent quark model predicts broad states decaying to $D^{(*)}K$ for the $J^{P}=\{0,1\}^{+}$ doublet with masses higher in about $100$ MeV than the experimental ones Godfrey and Isgur (1985); Godfrey and Kokoski (1991); Di Pierro and Eichten (2001). In contrast, the observed states, the $D_{s0}(2317)$ and $D_{s1}(2460)$, are very narrow resonances which are placed relatively close to the $DK$ and $DK^{*}$ thresholds respectively in the spectrum Aubert et al. (2003); Besson et al. (2003). For these reasons, molecular Barnes et al. (2003); Kolomeitsev and Lutz (2004); Faessler et al. (2007a, b); Lutz and Soyeur (2008); Albaladejo et al. (2018) and also tetraquark explanations Cheng and Hou (2003); Terasaki (2003); Maiani et al. (2005); Dmitrasinovic (2005) have been proposed. On the other hand, the predictions of the constituent quark model for the $J^{P}=\{1,2\}^{+}$ ($j_{l}^{P}=3/2^{+}$) doublet are in agreement with experiment Godfrey and Isgur (1985); Godfrey and Kokoski (1991); Di Pierro and Eichten (2001). Still, molecular components for the same quantum numbers cannot be discarded. In fact, in Molina et al. (2010), assuming that the $D^{*}_{2}(2573)$ is a $D^{*}K^{*}$ molecular state, the recenlty observed $T_{cs}$ and $T_{c\bar{s}}$ states Aaij et al. (2020a, b); Collaboration (2022) are predicted, see also Molina and Oset (2020, 2022).

In general, states with an additional strange quark are in principle expected to be $100$ MeV higher in mass, since $m_{s}/m_{d}\simeq 20$ MeV with $m_{d}\approx 5$ MeV, however, this is not the case for the charmed mesons. As reported by the PDG, the masses of the $D^{*}_{0}(2300)$ and $D_{1}(2420)$ are similar to those of their strange partners, $D_{s0}^{*}(2317)$ and $D_{s1}(2460)$. The situation for the $L=1$ charmed mesons is more complex. The measurements of the masses of the $0^{+}$ state lie between $2300$ and $2400$ MeV Abe et al. (2004); Aubert et al. (2009); Link et al. (2004). Theoretically, a two pole structure has been proposed in coupled-channel approaches based on chiral lagrangians, unitarity and heavy quark symmetry, resulting in two states with masses of around $2100$ MeV and $2400-2450$ MeV Kolomeitsev and Lutz (2004); Guo et al. (2006, 2009); Albaladejo et al. (2017); Du et al. (2018). The existence of a lower mass pole is supported by lattice QCD, which consistently finds the mass of the $D^{*}_{0}(2300)$ well below that of the $D^{*}_{s0}(2317)$, in contrast to the experiment Moir et al. (2016); Gayer et al. (2021).

A strongly interacting system at low energies composed of heavy and light mesons with the heavy meson being of $Q\bar{q}$ type can be studied by means of Heavy Hadron Chiral Perturbation Theory (HH$\chi$PT) Wise (1992); Burdman and Donoghue (1992); Yan et al. (1992). This is an effective field theory based on both chiral and heavy quark symmetries, where a systematic two double expansions in the parameters $Q/\Lambda_{\chi}$ and $\Lambda_{QCD}/m_{Q}$ can be made, being $Q$ the momentum transfer (soft scale), $Q\sim m_{\pi}\sim p_{\pi}$, $m_{Q}$, is the mass of the heavy quark, and $\Lambda_{\chi}=4\pi f\simeq 1$ GeV (hard scale). HH$\chi$PT has been applied to evaluate the ground state heavy meson masses in the heavy quark limit Isgur and Wise (1989, 1990). Leading order corrections to these masses due to the finite heavy quark mass have been studied in Rosner and Wise (1993); Randall and Sather (1993); Di Bartolomeo et al. (1995); Blok et al. (1997); Jenkins (1994); Yeh and Lee (1995) for non-zero light quark masses. The $j^{P}_{l}=\frac{1}{2}^{+}$ doublet has also been introduced as an explicit degree of freedom in Mehen and Springer (2005); Ananthanarayan et al. (2007); Alhakami (2016) where the masses of both doublets were studied including one-loop corrections.

Given the large amount of parameters in one-loop HH$\chi$PT it has not been possible yet to pin down precisely the parameters of the theory. Because of the recent progress in lattice QCD it is worth coming back to this issue. In this article we study the quark mass dependence of the ground state charmed mesons, $D$, $D^{*}$, $D_{s}$ and $D^{*}_{s}$ within one-loop HH$\chi$PT Jenkins (1994) by analyzing the available lattice data. There are several reasons to do this: (1) to test the predictive power of HH$\chi$PT with the available lattice data, (2) to better understand the significance of one-loop corrections involving light mesons for the ground state charmed mesons, (3) to provide an overview of the trends of charmed mesons in different lattice QCD simulations, and (4), since many exotic states have been observed in the recent years Guo et al. (2018a); Aubert et al. (2003); Besson et al. (2003); Aaij et al. (2020a, b, 2022), the analysis done here can be used in the future to investigate the quark mass dependence of these new states. This could help to discriminate between diferent pictures (tetraquark, molecule …).

We investigate both, the light and heavy quark (near the physical point) mass dependence of the ground state charmed mesons, $D_{(s)}$ and $D_{(s)}^{*}$. We investigate which parameters of one-loop HH$\chi$PT are more relevant as required by the LQCD data implementing model selection tools, as the Least Absolute Shrinkage and Selection Operator (LASSO) method in combination with cross validation and information criteria James and Tibshirani (2009); Tibshirani (1996); Landay et al. (2017). This can be possible now since there are sufficient lattice simulations to better constraint the parameters of HH$\chi$PT. For this analysis, we have taken as input the lattice data on the ground state charmed mesons of Ref. Kalinowski and Wagner (2015); Carrasco et al. (2014) (ETMC), the PACS-CS ensembles Mohler and Woloshyn (2011); Aoki et al. (2009), the ensembles from Hadron Spectrum Collaboration (HSC) Cheung et al. (2016, 2021), the data of Table 1 of Prelovsek et al. (2021), the data from Ref. Bruno et al. (2017) (CLS), and that of Refs. Bali et al. (2017, 2015) (RQCD). Most of these data are collected in Tables tables VIII-XIV of Guo et al. (2018b), where a previous lattice data analysis using the chiral lagrangians developed in Lutz and Soyeur (2008); Kolomeitsev and Lutz (2004) was performed. In addition, we incorporate the latest studies of charmed meson masses from MILC in $N_{f}=2+1+1$ Bazavov et al. (2018a, b), which have been done for different quark masses and lattice spacings.

II Formalism

Heavy meson masses to one-loop order in HH$\chi$PT

In HH$\chi$PT  Wise (1992); Burdman and Donoghue (1992); Yan et al. (1992) it is customary to introduce a single field, $H_{l}$, with $l=1,2,3$, to refer to $u,d,s$ respectively, including both the pseudoscalar and vector meson fields defined as, $P_{l}^{(*)}=(D^{(*)0},D^{(*)+},D_{s}^{(*)})$, for the charmed mesons (and similarly for the bottom mesons),

[TABLE]

with $Q=c,b$. Similarly, the conjugated meson field is defined as,

[TABLE]

These heavy meson fields transforms as a doublet under heavy quark spin symmetry and like an anti-triplet under $SU(3)_{V}$ flavour symmetry, i.e.,

[TABLE]

and

[TABLE]

for the heavy meson field and its conjugate, respectively. A velocity dependent field can also be defined as Georgi (1990),

[TABLE]

which satisfies $\not{v}H_{v}=H_{v}$ and $H_{v}\not{v}=-H_{v}$. The mass $m_{H}$ in Eq. (5) can be expanded in terms of powers of $1/m_{Q}$ starting in $m_{H}=m_{Q}$. The chiral lagrangian involving heavy meson fields defined as in Eq. (5) and pseudoscalar fields, has a systematic expansion in $1/m_{H}$, derivative and loop expansion. At leading order, the chiral lagrangian reads Jenkins (1994),

[TABLE]

where $S_{lv}$ stands for the velocity dependent spin operator $S^{\mu}_{lv}$ Jenkins and Manohar (1991), which satisfies $v\cdot S_{lv}=0$, and in the rest frame of the heavy meson, $v^{\mu}=(1,0)$, this is, $S^{\mu}_{l}=(0,\vec{\sigma}/2)$. The coupling $g$ can be fixed to obtain the experimental radiative and pion decay widths of the charmed vector mesons, we set $g^{2}=0.55$. The pseudoscalar meson fields appear in the combinations,

[TABLE]

where $\xi=\sqrt{U}$, $U=e^{i\phi/\sqrt{2}f}$, and

[TABLE]

For the one-loop calculation of the heavy meson masses, terms with up to two insertions of the light quark matrix $M=\mathrm{diag}(m_{u},m_{d},m_{s})$, are needed because counterterms at one-loop order are proportional to two powers of the light quark masses. Thus, the relevant lagrangian needed for the one-loop calculation is given by Jenkins (1994),

[TABLE]

where $M_{\xi}=\frac{1}{2}(\xi M\xi+\xi^{\dagger}M\xi^{\dagger})$ and $M$ is the (diagonal) light quark mass matrix. Here it is assumed $m=m_{u}=m_{d}$. If considering isosping breaking, extra terms should be added, see Eq. (2.14) of Jenkins (1994). In general, operators which respect HQSS have an expansion which starts at ${\cal O}(1)$, while the ones violating HQSS, start at ${\cal O}(\frac{1}{m_{Q}})$. The coefficients $a$ and $\sigma$ are dimensionless. These are functions of $\frac{1}{m_{Q}}$ starting at ${\cal O}(1)$ and therefore terms with these coefficients respect HQSS. These terms with $a$ and $\sigma$ give rise to the $SU(3)_{V}$ violating mass splittings among the $P^{(*)}$ mesons and to a singlet contribution to the masses, respectively. Terms accompanying the coefficients $b,c$ and $d$, correspond to chiral lagrangian terms with two insertions of the light quark mass matrix, preserving also heavy quark spin symmetry. The hiperfine splittings of the heavy mesons depend on the coefficients $(\Delta,\Delta^{(\sigma)},\Delta^{(a)})$, which start at ${\cal O}(\frac{1}{m_{Q}})$, and therefore, terms accompanying these coefficients violate HQSS. If one considers only terms with no insertions of the light quark matrix, one has $M_{P^{*}_{a}}-M_{P_{a}}=\Delta$, which is the hyperfine mass splitting at tree level. Terms with $\Delta^{(\sigma)}$ and $\Delta^{(a)}$ in Eq. (12) correspond to one insertion of the light quark mass matrix. While the former yields a light quark flavor dependent term, the latter leads to a flavor singlet contribution. Both of them contribute to counterterm coefficients in the $M_{P^{*}}-M_{P}$ splitting.

The one-loop masses of heavy mesons can be written in terms of two linear combinations, the spin average term, $\frac{1}{4}(M_{P_{a}}+3M_{P^{*}_{a}})$ and the hyperfine splitting, $M_{P^{*}}-M_{P}$, which respects and violate heavy quark spin symmetry, respectively. These, which have been calculated to one-loop in Ref. Jenkins (1994), are

[TABLE]

where the coefficients $\alpha_{l}$, $\beta_{l}^{(X)}$, $(\gamma_{l}^{(X)}-\lambda_{l}^{(X)}\alpha_{l})$, $c_{l}$, $(\gamma_{l}^{(X)}-\lambda_{l}^{(X)}\Delta)$ and $\delta c_{l}$ are given in the appendix A. The coefficients $\alpha_{l}$, $(\gamma_{l}^{(X)}-\lambda_{l}^{(X)}\alpha_{l})$, $c_{l}$, and $\delta c_{l}$ are proportional to powers of the light quark masses while $\beta_{l}^{(X)}$ and $(\gamma_{l}^{(X)}-\lambda_{l}^{(X)}\Delta)$ accompany terms proportional to $M_{X}$, the mass of the $X$ pseudoscalar meson, $\pi$, $K$ or $\eta$. Thus, $m_{H}$ and $\Delta$ can be interpreted as the spin-average mass and hyperfine splitting of the heavy mesons in the SU(3) chiral limit. The scale $\mu$ in Eqs. (13) and (14) is fixed to $770~{}$MeV.

III Fitting procedure with the LASSO method

We collect the available lattice data on low-lying charmed mesons included here in Tables 9-15, of the Appendix C. These include,

• •

Data from Ref. Kalinowski and Wagner (2015), which uses Wilson twisted mass lattice QCD with 2+1+1 dynamical quark flavors for pion masses in the range of $225-445$ MeV, $L\sim 2-3$ fm and three lattice spacings $a=0.0619,0.0815$ and $0.0885$ fm Carrasco et al. (2014) [ETMC] (Tables 9 and 10).

• •

The PACS ensembles with $N_{f}=2+1$ flavor Clover-Wilson configurations Mohler and Woloshyn (2011) for pion masses in the range $m_{\pi}\sim 150-410$ MeV, $L\simeq 2.9$ fm and a lattice spacing $a=0.0907$ MeV Aoki et al. (2009) [PACS-CS] (Table 11).

• •

The ensembles from Hadron Spectrum Collaboration for $N_{f}=2+1$ flavours of dynamical quarks with an anisotropic (clover) action Cheung et al. (2016), with $a_{s}\sim 0.12$ fm, $L\sim 1.9-3.8$ fm and $m_{\pi}\simeq 240,390$ MeV. See also Table 2 of Cheung et al. (2021) [HSC] (Table 12).

• •

The data of Table 1 of Prelovsek et al. (2021) for a $N_{f}=2+1$ simulation with Wilson dynamical fermions provided by the Coordinated Lattice Simulations (CLS) with a pion mass of $m_{\pi}=280$ MeV, $L\simeq 2.1$ fm and a lattice spacing of $a=0.08636$ fm Bruno et al. (2017) [CLS] (Table 13).

• •

Data collected in Lang et al. (2014) for $N_{F}=2$ dynamical light quarks with improved Wilson fermions, a lattice spacing $a=0.1239$, $L\simeq 2$ fm and $m_{\pi}\simeq 266$ MeV [Prelovsek et al.] (Table 14).

• •

The data of the simulation with $N_{f}=2$ improved clover Wilson sea quarks collected in Table 1 of Bali et al. (2017) for $m_{\pi}=150,290$ MeV, $L\sim 1.7-4.5$ fm and $a=0.07$ fm Bali et al. (2015) [RQCD] (Table 15).

• •

The latest studies of charmed meson masses from MILC with highly improved staggered-quark ensembles with four flavors of dynamical quarks Bazavov et al. (2018a, b), which have been done for different quark masses, $m_{\pi}\sim 130-330$, volumes $L\sim 2.5-6$ fm and lattice spacings ranging from $0.03-0.15$ fm (Table I of Bazavov et al. (2018b)) allowing for a continuous extrapolation Brambilla et al. (2018) [MILC]111We use as input the function in the continuum limit from MILC Brambilla et al. (2018) for the $D$ and $D_{s}$ masses..

Some of these data were collected in Tables VIII-XIV of Guo et al. (2018b). Indeed, the data of PACS-CS and ETMC used here correspond to Tables VIII, X and XI of Guo et al. (2018b), where a previous lattice data analysis using the chiral lagrangians developed in Lutz and Soyeur (2008); Kolomeitsev and Lutz (2004) was performed. The data of Tables XII and XIII of Guo et al. (2018b) based on MILC ensembles of the HPQCD and LHPC collaborations are not included here, because instead we include the latest MILC studies of Bazavov et al. (2018a, b) which superseed the older data in finests lattice spacing and bigger volumes that yielded a continuous extrapolation Brambilla et al. (2018).

In this work we will analyze the above data with one-loop order HH$\chi$PT given by Eqs. (13) and (14), looking for consistency between the different sets of data and trying to find a unique solution of the charmed meson masses that provides both, the heavy and light quark mass dependence of those. At present, this can only be achieved considering data from different collaborations. Regarding this, several aspects need to be discussed. First, we observe that, how the charm quark is taken to the physical point and which value is taken the observable to do so, can explain largely the deviations observed by the different collaborations, concluding that these actually come mostly from different charm quark masses at the so-called “physical point”. To compare the charm quark mass settings of the LQCD data used here, in Table 1 we compare the value of $m_{D_{s}}$ in physical units for the lightest pion mass and finest lattice spacing of every collaboration.

For the value of the lattice spacing given in Cheung et al. (2021) determinded through the $\Omega$ baryon mass, the mass of the $D_{s}$ for HSC is physical for its light pion mass, the one for RQCD and CLS is around $10$ MeV above the physical one, $20-30$ MeV above for ETMC, while being much lower for the PACS-CS, and [Prelovsek et al.] data sets in $\sim 160$ and $\sim 300$ MeV, respectively. Second, there might be discretization effects of different size. In general, we consider that the systematic effects due to discretization are larger for the charm quark mass and that these can be absorbed in the parameters $m_{H}$ and $\Delta$, as follows,

[TABLE]

where the term dependent on the lattice spacing $a$ will be different for every collaboration depending on the action used, and $m_{H}^{\prime}$, $\Delta^{\prime}$, are just the limit in the continuum of the tree level parameters. Here, we focus on fitting the $m_{H}$ and $\Delta$ parameters, introducing a new couple of parameters for every different lattice spacing, collaboration, and charm quark mass determination. There are two collaborations which perform simulations at different lattice spacings, MILC and ETMC. For MILC we perform a refit of the continuous extrapolation of Brambilla et al. (2018) that already takes into account the possible systematic effects of the staggered fermions, while for ETMC we analyze the data first and then study the discretization effects a posteriori making ${\cal O}(a^{k})\to a^{2}$ in Eqs. (15) and (16).

In the next section, first we perform individual fits for every collaboration. Then, we make an attempt to perform a global fit of data combining the different data sets. When doing so, we realize that the inclusion of the PACS-CS data set increases tremendously the $\chi^{2}$. For this reason we do not include these data in the global fit. In the case of MILC, we have a continuous function which is possible to fit in combination with the rest of data (we take $10$ data points for the $m_{D}$ mass of this function spread equally in the $m_{\pi}$ range of the lattice data, Table I of Bazavov et al. (2018b)). Since the discretization effects are embedded in the different $m_{H}$ and $\Delta$ parameters, there is no reason of concern to mix these various data sets, and indeed, as we will see, we obtain a reasonably good value of the $\chi^{2}_{\mathrm{d.o.f}}$. Thus, for the global fit, we have, $7+2n$ parameters, with $n=12$, equals to $31$ parameters, $40$ data points for $m_{D_{(s)}}$ and $30$ data for $m_{D_{(s)}^{*}}$222In the case of RQCD (Table 15) there are no data for $m_{D_{s}^{(*)}}$., and the two equations Eqs. (13) and (14). Apart from the various $(m_{H},\Delta)$ parameters introduced, the other parameters to be fitted in Eqs. (13) and (14) are, ($\frac{\sigma m_{\pi}}{B_{0}}$, $\frac{am_{\pi}}{B_{0}}$, $\frac{bm^{3}_{\pi}}{B_{0}^{2}}$, $\frac{cm^{3}_{\pi}}{B_{0}^{2}}$, $\frac{dm^{3}_{\pi}}{B_{0}^{2}}$, $\frac{\Delta^{(\sigma)}m_{\pi}}{B_{0}}$, $\frac{\Delta^{(a)}m_{\pi}}{B_{0}}$, where $m_{\pi}$ is the physical pion mass, and it is introduced to redefine the parameters in terms of adimensional quantities.

Given that the size of the data sample is not very large compared to the number of parameters, one could easily have problems with overfitting the data. In order to avoid that, we follow the procedure of Ref. Landay et al. (2017), where the Least Absolute Shrinkage and Selection Operator (LASSO) method in combination with cross validation James and Tibshirani (2009); Tibshirani (1996) and information criteria were applied for the analysis of photoproduction data with a set of models with different number of parameters. With this method, only those models which contain an optimal number of parameters able to describe the data with enough accuracy are selected. To implement the LASSO, a penalty term is added to the $\chi^{2}$

[TABLE]

In order to select the best “$\lambda$” parameter, we proceed as follows: in every iteration we divide randomly the original data set into three subsets. We keep two subsets as the training data set and one subset as the validation data set. These sets will provide the $\chi^{2}$ of the training and validation sets, $\chi^{2}_{T}$, $\chi^{2}_{V}$, respectively. As explained in Refs. James and Tibshirani (2009); Tibshirani (1996); Landay et al. (2017), while $\chi^{2}_{T}$ is a monotonously increasing function of $\lambda$, $\chi^{2}_{V}$ has a sweet spot in which the data are not underfitted (large $\lambda$) and not overfitted (small $\lambda$). The fluctuations given by the random selection of the validation set provide the error of this minimum, and also the one of the $\lambda$ value, which is chosen as the $\lambda_{\mathrm{opt}}>\lambda_{\mathrm{min}}$, compatible within errors with the minimum of $\chi^{2}_{V}(\lambda_{\mathrm{min}})$.

The result for the $\lambda_{\mathrm{opt}}$ from cross validation is compared with the one obtained from three information criteria, which are, AIC, AICc, and BIC, defined as,

[TABLE]

In the above equations, $k$ is the number of parameters that changes with $\lambda$, $n$ is the length of the data, and $L$ is the likelihood. The optimal value of $\lambda$ is given by the respective minimum of the criteria.

IV Results

In this section, we present, first, the results obtained from indiviudal fits for data sets larger enough to be fitted by means Eqs. (13), (14) (see also Eqs.  (29), (30), (31)). These data sets are: ETMC, PACS and MILC. Second, we show the result of a global fit including also the rest of the data sets, HQS, CLS, [Prelovsek et al.] Lang et al. (2014), with few data points for heavy meson masses.

All the fits are carried out with quantities in units of the lattice spacing ($am_{\pi}$, $am_{D}$…). As commented in the previous section, we assume that discretization effects are much larger for the charm quark mass and absorb the overall dependence of the charmed meson masses with the lattice spacing in the $m_{H}$ and $\Delta$ parameters (we introduce one pair of ($m_{H}$, $\Delta$) for every lattice spacing). For the data sets with more than one lattice spacing, ETMC data, we also analyze this discretization effect after the fit. We will compare this result at $a=0$ with the one from MILC (extrapolated data to the continuum) and the result of the global fit.

IV.1 ETMC

The results after performing the LASSO are collected in Fig. 1.

As can be seen in Fig. 1, in the first figure, from top to bottom, the $\chi^{2}_{\mathrm{dof}}$ is a monotonously increasing function with $\lambda$. Fig. 1 (second) shows the relevance of the parameters according to the lattice data. In this figure it is clear that the data do not constrain well the parameter $b$, which according to the analysis with LASSO, it is a superfluous degree of freedom. Next, $d$ and $\Delta_{\sigma}$ are also less relevant for the data. The $b$ and $d$ parameters correspond to chiral lagrangian terms with two insertions of the light quark mass matrix preserving HQSS, while $\Delta_{\sigma}$ violates HQSS and provides higher order corrections to the hyperfine splitting at tree level. Notice that parameters of these both types, like $c$ and $\Delta_{a}$, are, however, required by the data, i. e., not all the parameters of the same order are equally relevant.

One should have a criteria in order to decide whether one or the three parameters can be removed. We follow the criteria used in Landay et al. (2017), explained in the previous section, cross validation and information criteria. These are given in the third and fourth panels of Fig. 1, respectively. Information criteria are evaluated with all data, while for cross validation, we divide the data sets in three sets for every lattice spacing, $a$, and select $2$ sets for the training set and $1$ for the validation (test) set. As one can see, the information criteria have three minima corresponding to leaving out, one parameter, $b$, which is also the global minima, two parameters, $b$ and $d$, and three parameters, $b$, $d$ and $\Delta_{\sigma}$, respectively. The global minima coincides with the one from cross validation, in Fig. 1 (second). Thus, we select this parametrization dropping out parameter $b$ as the most predictive one.

The parameters obtained before and after the LASSO, that we call Fits I and II, respectively, are given in the two first rows of Tables 2 and 3.

We see that the result between the two different fits is very similar and that, in fact, the inclusion of the $b$ parameters, do not affect much the values of the rest of the parameters in this case.

For these data sets, one can also study the discretization effects. Assuming,

[TABLE]

with the value of $m_{H}$ and $\Delta$ for different $a$ obtained in the fit and its errors plotted in Fig. 2.

We fit the two different data from fixing the charm quark mass with the $\eta_{c}$ or $J/\psi$ mass with Eq. (18). We observe that, in these data sets, the difference between fixing the charm quark mass with the $\eta_{c}$ or $J/\psi$ mass is basically due to discretization effects, since this difference is close to zero for $a\sim 0$, see Fig. 2 (top). For the extrapolation to the continuum of $m_{H}$, we take the average of the two lines at $a\to 0$. For the hyperfine splitting, we obtain the line in Fig. 2 (bottom). Since this parameter starts its expansion at ${\cal O}(1/m_{Q})$, there is no difference between setting the charm quark mass with $\eta_{c}$ or $J/\psi$, and the error due to discretization effects of this parameter is much smaller, being much more stable when varying the lattice spacing. For $a=0$, we give the values of $m_{H}$ and $\Delta$ in the last row of Table 3.

We observe that the parameter $m_{H}$ is indeed very close to the physical value of $(m_{D}+3m_{D^{*}})/4=1973$ MeV, however, $\Delta$ differs $\sim 20\%$ from the physical hyperfine splitting, $m_{D^{*}}-m_{D}=141~{}$MeV. This motivates us to perform also fits with the expressions of the charmed meson masses at tree level, keeping terms with one insertion of the light quark mass matrix in Eq. (12). The results are given in the Appendix B. While for the one-loop expressions we obtain $\Delta=114$ MeV for $a=0$, for the tree level expressions, Eqs. (33) and (34), the continuous extrapolation of this parameter that we obtain is $\Delta=129$ MeV. The values of $m_{H}$ are quite similar for these different fits, $m_{H}\simeq 1970$. As expected, the one-loop expressions reduce the value of the hyperfine splitting in the chiral limit. We show the results of Fit II in Fig. 3. Eqs. (13) and (14) provide a very good description of the data. The $\chi^{2}_{\mathrm{dof}}$ obtained is 0.44. To evaluate the errors, here and in the following, we use automatic differentiation Ramos (2019) with the Julia ADerrors package.

Finally, we depict in Fig. 4 the extrapolation to the continuum limit of the spin average mass and the hyperfine splitting (red band) in comparison with the ETMC data set. As can be inferred, discretization effects appear to be of $\sim 2-4\%$ in the spin average mass for these data sets. However, for the hyperfine splitting, which suffers from large errors of $\sim 20-30\%$, most of the LQCD data fall inside the error band. In the continuum limit, the LQCD data agree well with the experimental value, shown with a star.

IV.2 MILC

These lattice data are based on the simulations of Bazavov et al. (2018a, b) for different pion masses and lattice spacings, which made it possible to extrapolate to the continuum the masses of the charmed mesons using EFT’s, with the theoretical framework being based on a merger of one-loop HQET and HMrAS$\chi$PT, where a generic lattice spacing dependence, as well as higher order terms, where incorporated Brambilla et al. (2018). This function has been provided to us and here we conduct a refit of it with Eqs. (13) and (14). In this case, the $\chi^{2}_{\mathrm{test}}$ is a monotonously increasing function and does not show a minimum except for $\lambda=0$. However, the information criteria, shown in Fig. 5 (top), exhibit a minimum around $\lambda\simeq 2$. Actually, for this value of $\lambda$, four parameters have been dropped out, $b$, $d$, $\Delta_{\sigma}$ and $\Delta_{a}$.

Note that the first three parameters are also the less relevant parameters for the ETMC data. The values of the parameters, before and after performing the LASSO method, Fits I and II, are given in Tables 2 and 4. In this case the data and continuous extrapolation only contain the values for the $m_{D_{(s)}}$ masses. Since the masses of the vector counterparts are absent we have to fix the value of $\Delta$. We have tried several values in the range of the quantities obtained in the other fits, and here, fixing $\Delta$ to $110$ MeV seems to provide the best result. Since four parameters are dropped after the LASSO, the values of the parameters in Fits I and II are quite different. The result of the parameters after the LASSO method has much smaller errors, telling that the method is working well. Concretely, in Fit II, the value and error obtained for $m_{H}$ is quite reasonable and also compatible with the one from the ETMC fit. The $\chi^{2}_{\mathrm{dof}}$ obtained in this case is small, 0.04, because of the refit of the function. In Fig. 6 we depict the charmed meson masses in comparison with the continuous extrapolation of the ETMC data that we did previously and the result of the physical extrapolation for the global fit presented in Sec. IV.4. We observe that the ETMC data show less dependence of the charmed meson masses with the pion mass than MILC data and the result from the global fit, but the error bands overlap each other and with the experimental charmed meson masses, shown with stars. As could be expected, $m_{D_{s}}$ does not depend much with the pion mass, since the strange quark mass is close to physical.

We can also compare our result for the charmed meson mass in the SU(3) chiral limit with the one obtained in Bazavov et al. (2018a); Brambilla et al. (2018). In Bazavov et al. (2018a) the value of $m_{D}$ in the SU(3) chiral limit obtained is $1842.7\pm 5.4$ MeV, while our fit results in $m_{D}^{SU(3)}=1877.0\pm 39.7$ MeV. These two results are consistent.

IV.3 PACS-CS

In this case the information criteria did not show a clear minimum, except for $\lambda=0$. However, cross validation gives a clear minimum in $\chi^{2}_{\mathrm{test}}$ around $\lambda\simeq 2$, as shown in Fig. 7 (top), which corresponds to dropping out the $b$ parameter, see Fig. 7 (bottom).

Nevertheless, the $\chi^{2}_{\mathrm{dof}}$ obtained is very high (much larger than $1$). This is clearly not consistent with the previous data sets, where very good values of the $\chi^{2}_{\mathrm{dof}}$ were obtained. We conclude that this lattice data set for the charmed meson masses is in tension with the rest of the lattice data sets. We notice that the high value of the $\chi^{2}$ for this collaboration is caused by one light pion mass. If we eliminate this point and we refit the data we obtain a $\chi^{2}_{\mathrm{dof}}$ of $4.4$. We show the plots of the charmed meson masses as a function of the pion mass in lattice units in Fig. 8,

were the conflicting data point is shown in red color. The values of the parameters obtained are given in Tables 2 and 5. Still, the parameters obtained in Table 2 are not so different from the other data sets. However, $m_{H}$ in Table 5 is significantly lower, suggesting that the charm quark mass in these simulations is much lower than the experimental, as also suggested in Table 1.

IV.4 Global fit

This fit includes the data of ETMC, HSC, CLS, [Prelovsek et al.] and RQCD given in Tables 9, 10, 12, 14 and 15. The data of PACS-CS are not included because of the tensions observed with the other data sets. In this case the information criteria, shown in Fig. 9 (top), have two minima, for $\lambda\simeq 0.5$ (global minima) and $\lambda\simeq 2.5$. This corresponds again to leaving out one parameter, $b$, or three, $b$, $d$ and $\Delta_{\sigma}$, see Fig. 9 (bottom). The results of the parameters before and after the LASSO are given in Tables 2 and 6.

In Figs. 10 and 11 we plot the results for the charmed meson masses in comparison with the data. In Fig. 10 we plot the spin averaged masses and the hyperfine splittings of the charmed mesons as a function of the pion mass in physical units. We leave out ETMC and [Prelovsek et al.] data from this plot for a better view. In general, the error bands for ETMC data overlap with those from HSC, and the band for [Prelovsek et al.] data results in much smaller charmed masses. We show the results for these two data sets separately in Fig. 11 in units of the lattice spacing. As can be seen all data sets included in the global fit are very well described. The $\chi^{2}_{\mathrm{dof}}$ obtained is $0.85$. In Fig. 10 we also show the result for the global fit fixing the $m_{H}$ and $\Delta$ parameters according to the experimental data. As can be seen, the physical trajectories for the spin averaged charmed meson masses are closer to the ones of RQCD data, which includes an almost physical pion mass. Still, the hyperfine splittings of the charmed mesons in the lattice data are smaller than the physical ones in $\sim 8-15\%$.

We also performed an analysis with the tree level expressions, Eqs. (33) and  (34), and the results are given in Appendix B. The $\chi^{2}_{\mathrm{dof}}$ obtained is a bit larger than with the one-loop expressions, $1.11$, but it also admits a good fit. In general, the values of $m_{H}$ and $\Delta$ are larger than the ones obtained with the one-loop expressions, as expected. The values obtained for $\Delta$ are more reasonable for the one-loop expressions than for the tree level expressions, as can be seen comparing Table 6 and Table 8, where for ETMC, larger values of $\Delta$ than the physical hyperfine splitting, are obtained using the tree level expressions.

It is interesting to see how the $m_{H}/M_{\mathrm{av}}^{c\bar{c}}$ and $\Delta$ parameters, with $M_{\mathrm{av}}^{c\bar{c}}=(m_{\eta_{c}}+3m_{J/\psi})/4$, $m_{H}$ and $\Delta$ being the spin average mass and hyperfine splitting in the chiral limit, depend on the lattice spacing. This is shown in Fig. 12.

From Fig. 12, one obtains that $m_{H}\simeq 0.63M_{\mathrm{av}}^{c\bar{c}}$, where we expect this relation to be valid for charm quark masses beyond the physical one, but not very far from the range of the charmed quark masses of the LQCD simulations studied here333At the physical point, this relation deviates $\sim 10\%$ for the bottom counterparts assuming that the physical masses behave as leading order.. Since $\Delta$ is a heavy quark spin violating parameter, one does not expect it to depend heavily on $m_{H}$. In Fig. 12, taking into account the dependence of this parameter with the lattice spacing for the ETMC data, one sees that when $a\to 0$ this would lead to a range, $\Delta\sim 90-130$ MeV. In the future, more precise LQCD data can help to determine more precisely this parameter.

V Conclusions

We have investigated the light and heavy quark mass dependence of the $D^{(*)}_{(s)}$ mesons by means of an analysis of the available lattice data from different LQCD Collaborations withing one-loop HH$\chi$PT. Model selection tools were implemented to determine the low energy constants with a higher precision.

We have found that, concerning the higher order parameters, $b,c,d,\Delta_{\sigma}$ and $\Delta_{a}$, while $c$ and $\Delta_{a}$ are required by the precision of the LQCD data, $b,d$ and $\Delta_{\sigma}$ are less relevant, and, in general, LQCD data do not restrict the $b$ parameter. We have found more precise values of the one-loop parameters, once we have applied these model selection techniques.

We also observe that, the largest difference in the various simulations comes mainly from the setting of the charm quark mass to the physical point, and that this difference becomes also larger due to discretization effects. Improvement on the methods to set the physical charm quark mass in the LQCD simulations are required at present, according to the present analysis. However, once these differences are taken into account, most of the LQCD data (all except for PACS) can be very well described together by one-loop HH$\chi$PT, considering different values of the tree level parameters, $m_{H}$ and $\Delta$. The data can be described well also with the tree level expressions taking into account one insertion of the light quark matrix instead of two, however, too high values of $\Delta$ are obtained in some cases.

The description with one loop involving pseudoscalar Goldstone bosons describes better the LQCD data and provides more reasonable values of the tree level parameters after the LASSO method is implemented. Once the charm quark mass is fixed within the theory as required by the experimental data, the charmed meson masses obtained from the global fit of LQCD data are consistent with the experimental data.

VI Acknowledgments

R. M. and F. Gil thank J. Komijani for providing the MILC continuous extrapolation to the physical point of the $D_{(s)}$ meson. We also thank A. Ramos and F. Romero for useful discussions. We acknowledge support from the CIDEGENT program with Ref. CIDEGENT/2019/015 and from the spanish national grants PID2019-106080GB-C21 and PID2020-112777GB-I00. This project has received funding from the European Union’s Horizon 2020 programme No. 824093 for the STRONG-2020 project.

Appendix A Coefficients in Eqs. (13) and  (14)

The coefficients $\alpha_{l}$, $\beta_{l}^{(X)}$, $(\gamma_{l}^{(X)}-\lambda_{l}^{(X)}\alpha_{l})$, $c_{l}$, $(\gamma_{l}^{(X)}-\lambda_{l}^{(X)}\Delta)$ and $\delta c_{l}$, are given by,

[TABLE]

[TABLE]

[TABLE]

[TABLE]

where we have assumed the isospin limit $m=m_{u}=m_{d}$ and $\delta_{1,2l}=\delta_{1,l}+\delta_{2,l}$. For the general formulas taking into account isospin breaking see Ref. Jenkins (1994). The light quark masses can be related to the leading order pseudoscalar masses, $M_{X}^{0}$, like

[TABLE]

For simplicity, we identify the physical masses of the light pseudoscalar meson masses with the leading order pseudoscalar masses, $M_{\pi}=M^{0}_{\pi}$ and $M_{K}=M^{0}_{K}$, $M_{\eta}^{2}=M_{\eta}^{02}=(4M_{K}^{02}-M_{\pi}^{02})/3$.

Appendix B Results for the tree level fits of the charmed meson masses

In this section we present the results of analyzing the charmed meson masses from LQCD with the tree level expressions including one insertion of the light quark mass matrix. The relations in these case are,

[TABLE]

[TABLE]

The parameters obtained in the various analyses using these relations are given in Tables 8 and 7. The $\chi^{2}_{dof}$ obtained is also good, $0.42$ for the ETMC fit and $1.11$ for the Global fit. For MILC we obtain a value close to zero because of the refitting of the continuously extrapolated function.

Appendix C Data included in the fits of IV