Charmonium resonances with $J^{PC}=1^{--}$ and $3^{--}$ from $\bar DD$ scattering on the lattice

TL;DR

This lattice QCD study investigates charmonium resonances with specific quantum numbers near the open-charm threshold, extracting their masses and widths, and compares findings with experimental data, including recent discoveries.

Contribution

First lattice QCD analysis of charmonium resonances with J^{PC}=1^{--} and 3^{--} considering strong D D̄ transitions, providing resonance parameters and exploring quark mass effects.

Findings

Identified the ψ(3770) resonance with mass 3780 MeV

Observed the ψ(2S) as a bound state at 3666 MeV

Found the X(3842) consistent with a 3^{--} resonance

Abstract

We present a lattice QCD study of charmonium resonances and bound states with $J^{PC}=1^{--}$ and $3^{--}$ near the open-charm threshold, taking into account their strong transitions to $\bar DD$. Vector charmonia are the most abundant in the experimentally established charmonium spectrum, while recently LHCb reported also the first discovery of a charmonium with likely spin three. The $\bar DD$ scattering amplitudes for partial waves $l=1$ and $l=3$ are extracted on the lattice by means of the L\"uscher formalism, using multiple volumes and inertial frames. Parameterizations of the scattering amplitudes provide masses and widths of the resonances, as well as the masses of bound states. CLS ensembles with 2+1 dynamical flavors of non-perturbatively $O(a)$ improved Wilson quarks are employed with $m_\pi\simeq 280$ MeV, a single lattice spacing of $a\simeq0.086$ fm and two lattice spatial extents of $L=24$ and $32$. Two values of the charm quark mass are considered to examine the influence of the position of the $\bar{D}D$ threshold on the hadron masses. For the lighter charm quark mass we find the vector resonance $\psi(3770)$ with mass $m=3780(7)$ MeV and coupling $g=16.0(^{+2.1}_{-0.2})$ (related to the width), both consistent with their experimental values. The vector $\psi(2S)$ appears as a bound state with $m=3666(10)$ MeV. The charmonium resonance with $J^{PC}=3^{--}$ is found at $m=3831(^{+10}_{-16})$ MeV, consistent with the $X(3842)$ recently discovered by LHCb. At our heavier charm-quark mass the $\psi(2S)$ as well as the $\psi(3770)$ are bound states and the $X(3842)$ remains a resonance. We stress that all quoted uncertainties are only statistical, while lattice spacing effects and the approach to the physical point still need to be explored. This study of conventional charmonia sets the stage for more challenging future studies of unconventional charmonium-like states.

Tables

Table 1. Table 1: List of two-meson interpolators O D ¯ ​ D subscript 𝑂 ¯ 𝐷 𝐷 O_{\bar{D}D} used in each lattice irrep Λ ( P ) ​ C superscript Λ 𝑃 𝐶 \Lambda^{(P)C} . For brevity we use the simplified notation D ¯ ​ ( p 1 2 ) ​ D ​ ( p 2 2 ) ¯ 𝐷 superscript subscript 𝑝 1 2 𝐷 superscript subscript 𝑝 2 2 \bar{D}(p_{1}^{2})D(p_{2}^{2}) , which refers to two-meson interpolators, with good charge conjugation C, built from a D ¯ ¯ 𝐷 \bar{D} and D 𝐷 D meson operator with momentum p → 1 subscript → 𝑝 1 \vec{p}_{1} and p → 2 subscript → 𝑝 2 \vec{p}_{2} , respectively in units of 2 ​ π / L 2 𝜋 𝐿 2\pi/L . The full expressions for these interpolators are much longer and are omitted for brevity.

	${\left|\overrightarrow{P}\right|}^2$	${\Lambda }^{\left(P\right)C}$	$O_{\overline{D}D}=\overline{D}\left(p_1^2\right)D\left(p_2^2\right)$
$O_h$	$0$	$T_1^{--}$	$\overline{D}\left(1\right)D\left(1\right)$
$Di{c}_4$	$1$	$A_1^{-}$	$\overline{D}\left(1\right)D\left(0\right)$
			$\overline{D}\left(2\right)D\left(1\right)$
$Di{c}_2$	$2$	$A_1^{-}$	$\overline{D}\left(2\right)D\left(0\right)$
			$\overline{D}\left(3\right)D\left(1\right)$

Table 2. Table 2: Lattice energies of D 𝐷 D meson for different momenta p 2 superscript 𝑝 2 p^{2} compared to the continuum values from Eq. ( 15 ). Above p 2 superscript 𝑝 2 p^{2} are given in units of ( 2 ​ π / L ) 2 superscript 2 𝜋 𝐿 2 (2\pi/L)^{2} . The spin-averaged charmonium mass of Eq. ( 6 ) is also given.

		$a{E}_D\left(p^2=0\right)$	$a{E}_D\left(p^2=1\right)$	$a{E}_D\left(p^2=2\right)$	$a{M}_{\text{av}}$
${\kappa }_c=0.12522,N_L=24$	lat	0.7732(9)	0.8132(10)	0.8516(11)	1.2318(20)
	cont	-	0.8163(9)	0.8573(8)	-
${\kappa }_c=0.12522,N_L=32$	lat	0.7711(6)	0.7942(5)	0.8166(6)	1.2348(13)
	cont	-	0.7957(9)	0.8296(8)	-
${\kappa }_c=0.12315,N_L=24$	lat	0.8453(13)	0.8814(10)	0.9164(65)	1.3559(20)
	cont	-	0.8849(13)	0.9228(12)	-
${\kappa }_c=0.12315,N_L=32$	lat	0.8433(6)	0.8642(5)	0.8844(6)	1.3587(14)
	cont	-	0.8659(8)	0.8878(8)	-

Table 3. Table 3: Summary of masses for resonances and bound states with J P ​ C = 1 − − superscript 𝐽 𝑃 𝐶 superscript 1 absent J^{PC}=1^{--} and 3 − − superscript 3 absent 3^{--} studied in our lattice simulation and compared to experiment Tanabashi  et al. ( 2018 ) . The dimensionless coupling g 𝑔 g parametrizes the width Γ = g 2 ​ p 3 / ( 6 ​ π ​ s ) Γ superscript 𝑔 2 superscript 𝑝 3 6 𝜋 𝑠 \Gamma=g^{2}p^{3}/(6\pi s) of the resonance decay ψ ​ ( 3770 ) → D ¯ ​ D → 𝜓 3770 ¯ 𝐷 𝐷 \psi(3770)\to\bar{D}D from Eq. ( 35 ). Lattice results are given for two charm quark masses, where the corresponding values of m D subscript 𝑚 𝐷 m_{D} are given at the top of the table. The splittings m − M av 𝑚 subscript 𝑀 av m-M_{\textrm{av}}~{} and m − 2 ​ m D 𝑚 2 subscript 𝑚 𝐷 m-2m_{D} are obtained taking solely lattice or solely experimental values of masses. The final mass m 𝑚 m resulting from our simulation is quoted according to Eq. ( 34 ). The last column contains the lattice results from the two-flavor simulation of Ref.  Lang  et al. ( 2015 ) based on ”fit (i)” described in Section V.6 . For the masses marked with (*) we quote the kinetic masses ( M 2 subscript 𝑀 2 M_{2} in Eq. (3.1) of Ref. Lang  et al. ( 2015 ) ) which are the appropriate masses for the Fermilab approach El-Khadra  et al. ( 1997 ); Oktay and Kronfeld ( 2008 ) pursued in Ref. Lang  et al. ( 2015 ) . These differ from the rest masses by discretization effects.

	$J^{PC}$	lat (present work)	lat (present work)	exp	lat Lang et al. (2015)
		${\kappa }_c=0.12315$	${\kappa }_c=0.12522$	${\overline{D}}^0{D}^0/D^{+}{D}^{-}$
$m_D$ [MeV]		$1927\left(2\right)$	$1762\left(2\right)$	${\overline{m}}_D\simeq 1867$ MeV	$1763\left(22\right){\left(18\right)}^{\ast }$
$m_{D_s}$ [MeV]		$1981\left(1\right)$	$1818\left(1\right)$	$1968.34\left(7\right)$
$M_{\text{av}}$ [MeV]		$3103\left(3\right)$	$2820\left(3\right)$	$3068.6\left(2\right)$	$3119\left(9\right){\left(33\right)}^{\ast }$
$m_{\pi }$ [MeV]		280	280	${\overline{m}}_{\pi }\simeq 137$MeV	266
$𝝍\left(\mathrm{𝟑𝟕𝟕𝟎}\right)$	$1^{--}$	bound st.	resonance	resonance Tanabashi et al. (2018)	resonance
g		$18.9({}_{-0.7}{}^{+0.8})$	$16.0({}_{-0.2}{}^{+2.1})$	$18.7\left(\mathrm{𝟗}\right)$	13.2(1.2)
$m-M_{\text{av}}$[MeV]		707(7)	711(7)	704.25(35)	715(7)
$m-2m_D$[MeV]		-43(8)	9(7)	38.52(35)
$𝒎$[MeV]		3776(7)	3780(7)	3773.13(35)121212We consider PDG fit estimate as the experimental value for the mass of the $\psi \left(3770\right)$ resonance throughout this article.	3784(7)
$𝝍\left(\mathrm{𝟐}𝑺\right)$	$1^{--}$	bound st.	bound st.	bound st. Tanabashi et al. (2018)	bound st.
$m-M_{\text{av}}$[MeV]		596(9)	597(10)	617.347(25)	605(6)
$m-2m_D$[MeV]		-154(10)	-105(11)	-48.383(25)
$𝒎$[MeV]		3665(9)	3666(10)	3686.097(25)	3674(6)
$𝑿\left(\mathrm{𝟑𝟖𝟒𝟐}\right)$	$3^{--}$	resonance	resonance	resonance Aaij et al. (2019)
$m-M_{\text{av}}$[MeV]		754(${}_{-7}{}^{}{}_{}{}^{+4}$)	762(${}_{-16}{}^{}{}_{}{}^{+10}$)	773.9(2)
$m-2m_D$[MeV]		4(${}_{-3}{}^{}{}_{}{}^{+9}$)	59(${}_{-16}{}^{}{}_{}{}^{+11}$)	108.2(2)
$𝒎$[MeV]		3822$({}_{-\mathrm{𝟕}}{}^{+\mathrm{𝟒}})$	3831$({}_{-\mathrm{𝟏𝟔}}{}^{+\mathrm{𝟏𝟎}})$	3842.7(2)

Equations

$$\det [{\tilde{K}}_l^{-1}(E_{cm}){\delta }_{l^{\prime }l}-B_{l^{\prime }l}^{\vec{P};\Lambda }(E_{cm})]=0.$$

$$(\rho t_l{)}^{-1}=2i(S_l-1{)}^{-1}=p^{-2l-1}{\tilde{K}}_l^{-1}-i,$$

$${\tilde{K}}_l^{-1}(E_{cm})=p^{2l+1}\cot {\delta }_l(p).$$

$$p^2=\frac{s}{4}-m_D^2,$$

$$p\cot {\delta }_1(p)-\frac{2{\mathcal{Z}}_{00}(p^2)}{\sqrt{\pi }L}=0,$$

$$M_{\text{av}}=\frac{1}{4}(m_{{\eta }_c}+3m_{J/\psi }).$$

$$M_{\text{av}}-2m_D$$

$$\tau (t^{\prime },{\mu }^{\prime },i^{\prime };t,\mu ,i)=⟨v_{i^{\prime }}^{{\mu }^{\prime }}(t^{\prime })|(D^{-1}|v_i^{\mu }(t)⟩),$$

$$k_n(t,i,\mu )=\frac{1}{\sqrt{2}}(\pm 1\pm i),$$

$$|w_n(t)⟩=\sum _{j,\nu }{k}_n(t,j,\nu )|v_j^{\nu }(t)⟩.$$

$$\tau (t^{\prime },{\mu }^{\prime },i^{\prime };t,\mu ,i)=\frac{1}{N_s}\sum _n{k}_n^{†}(t,i,\mu )⟨v_{i^{\prime }}^{{\mu }^{\prime }}(t^{\prime })|(D^{-1}|w_n(t)⟩),$$

$$\frac{1}{N_s}\sum _{n=1}^{N_s}{k}_n^{†}(t,i,\mu )k_n(t,j,\nu )\simeq {\delta }_{ij}{\delta }_{\mu \nu }.$$

$$C_{ij}(t)=⟨O_i(t)O_j^{†}(0)⟩=\sum _n{Z}_i^n{Z}_j^{n\ast }{e}^{-t{E}_n^{\text{lat}}},$$

$${\Lambda }^{PC}=T_1^{--}\text{mbox}intherestframe|\vec{P}{|}^2=0,$$

$${\Lambda }^C=A_1^{-}\text{mbox}inthemovingframe|\vec{P}{|}^2=1\text{mbox}and$$

$${\Lambda }^C=A_1^{-}\text{mbox}inthemovingframe|\vec{P}{|}^2=2,$$

$$O^{\bar{D}D}({\vec{p}}_1+{\vec{p}}_2)$$

$$-\bar{D}(R{\vec{p}}_2)D(R{\vec{p}}_1)\},$$

$$C_{ij}(t)v_n^i(t)$$

$${\lambda }_n(t)$$

$$E_D^{\text{cont}}(\vec{p})=\sqrt{m_D^2+|\vec{p}{|}^2}.$$

$$(\Delta E{)}_s=(E^{\text{lat}}{)}_s-(E_{D({\vec{p}}_1)}^{\text{lat}}{)}_s-(E_{D({\vec{p}}_2)}^{\text{lat}}{)}_s,$$

$$(E^{\text{calc}}{)}_s=(\Delta E^{\text{lat}}{)}_s+{\left(E_{D({\vec{p}}_1)}^{\text{cont}}\right)}_s+{\left(E_{D({\vec{p}}_2)}^{\text{cont}}\right)}_s$$

$$\frac{p^3\cot ({\delta }_1)}{\sqrt{s}}={\left(\frac{G_1^2}{m_1^2-s}+\frac{G_2^2}{m_2^2-s}\right)}^{-1},$$

$$\frac{p^3\cot ({\delta }_1)}{\sqrt{s}}=A+Bs+C{s}^2,$$

$$\frac{p^7\cot ({\delta }_3)}{\sqrt{s}}=\frac{m_3^2-s}{g_3^2}.$$

$$\frac{(ap)^{2l+1}\cot(\delta)}{\sqrt{a^{2}s}}=\left\{\begin{array}[]{ll}\left(\frac{[0.63(33)]^{2}}{[1.4966(30)]^{2}-a^{2}s}+\frac{[3.69(37)]^{2}}{[1.5457(32)]^{2}-a^{2}s}\right)^{-1}&~{}~{}~{}l=1\\ \frac{[1.568(11)]^{2}-a^{2}s}{[0.07(3)]^{2}}&~{}~{}~{}l=3\\ \end{array}\right.$$

$$\frac{(ap)^{2l+1}\cot(\delta)}{\sqrt{a^{2}s}}=\left\{\begin{array}[]{ll}\left(\frac{[0.40(50)]^{2}}{[1.63(1)]^{2}-a^{2}s}+\frac{[4.15(13)]^{2}}{[1.6745(17)]^{2}-a^{2}s}\right)^{-1}&~{}~{}~{}l=1\\ \frac{[1.6883(85)]^{2}-a^{2}s}{[0.025(17)]^{2}}&~{}~{}~{}l=3\\ \end{array}\right.\,,$$

$$\frac{(ap)^{2l+1}\cot(\delta)}{\sqrt{a^{2}s}}=\left\{\begin{array}[]{ll}{\scriptstyle-0.601({}^{+200}_{-240})+0.554({}^{+200}_{-170})(a^{2}s)-0.1267({}^{+380}_{-440})(a^{2}s)^{2}}&~{}~{}~{}l=1\\[12.0pt] \frac{[1.5639({}^{+83}_{-14})]^{2}-a^{2}s}{[0.058({}^{+94}_{-55})]^{2}}&~{}~{}~{}l=3\\ \end{array}\right.$$

$$\frac{(ap)^{2l+1}\cot(\delta)}{\sqrt{a^{2}s}}=\left\{\begin{array}[]{ll}{\scriptstyle-0.08({}^{+11}_{-13})+0.107({}^{+100}_{-80})(a^{2}s)-0.028({}^{+15}_{-17})(a^{2}s)^{2}}&~{}~{}~{}l=1\\[12.0pt] \frac{[1.6869({}^{+33}_{-32})]^{2}-a^{2}s}{[0.042({}^{+53}_{-22})]^{2}}&~{}~{}~{}l=3\\ \end{array}\right.$$

$$\Omega (\mu ,A(s))=\frac{\det (A)}{\det (({\mu }^2+A{A}^{†}{)}^{\frac{1}{2}})},$$

Full text

††preprint: MITP/19-029

I Introduction

The charmonium spectrum determined from experiment challenges our theoretical understanding of the internal dynamics of mesons containing charm quarks. While charmonia below the open-charm decay threshold are found to fit within a quark-antiquark picture, the properties of several resonances above the threshold, generally referred to as $XYZ$ states, call for exotic interpretations. Various suggestions have been made based on phenomenological approaches, effective field theories, potential models, etc. These include compact diquark-antidiquark states, mesonic molecules, hybrid mesons with gluonic excitations, or resonances generated through coupled channel scattering. The theoretical approaches are motivated by phenomenology and approximate certain regimes of the strong interaction. However, a complete description of the observed charmonium spectrum from Quantum ChromoDynamics (QCD) also requires complementary first principles non-perturbative calculations, which can be performed using lattice QCD. Ab-initio determinations of ground state charmonia, such as in Refs. Burch et al. (2010); Donald et al. (2012); Bečirević and Sanfilippo (2013); Briceño et al. (2012); Yang et al. (2015); Galloway et al. (2014); Mathur et al. (2018); DeTar et al. (2019), illustrate the power of these non-perturbative approaches in precisely determining masses and splittings. Lattice calculations of other conventional charmonium bound states and resonances aim to satisfactorily describe well understood excitations and demonstrate the efficacy of the methods employed. This gives confidence in investigations of unconventional charmonium-like states and of those for which the experimental situation is less clear.

To date, several lattice QCD calculations of the excited charmonium spectrum with isospin zero have been performed within the single hadron approach, where one employs only mesonic interpolators of $\bar{c}c$ type and assumes that this qualitatively captures the single meson spectrum in a finite volume Dudek et al. (2008); Bali et al. (2011); Liu et al. (2012); Mohler et al. (2013); Cheung et al. (2016). So far, only one lattice study has gone beyond the single-hadron approach, taking into account the strong decays of charmonium resonances above the open-charm threshold.111The strong decay of the $Z_{c}(3900)$ resonance with isospin 1 has been studied in the HALQCD approach Ikeda et al. (2016); Ikeda (2018). In that study an exploratory determination of the $\bar{D}D$ scattering amplitude in the $l=1$ and $l=0$ partial waves was performed Lang et al. (2015). The low-lying vector and scalar charmonium spectra were calculated in a single lattice volume in the center of momentum frame (CMF). The authors extracted the resonance and bound state parameters in the vector channel, while the conclusions in the scalar channel indicated different possible interpretations and left several unresolved puzzles for both theory and experiment.

Note that all previous lattice studies of charmonium were limited to the CMF Dudek et al. (2008); Bali et al. (2011); Liu et al. (2012); Mohler et al. (2013); Prelovsek and Leskovec (2013); Padmanath et al. (2015); Lang et al. (2015); Cheung et al. (2016). However, lattice calculations in the moving frames (i.e. with non-zero total momentum) can provide additional information on the relevant two-meson scattering amplitudes (for a review, see Ref. Briceño et al. (2018)). The challenge in spectroscopy using lattice techniques is that, due to the reduced symmetry on the lattice, several $J$ can contribute (in the CMF) to any given lattice irrep. In the moving frames, the situation is much worse as the finite volume spectrum in any lattice irrep gets a contribution from states with different $J^{P}$ (in their rest frame), as they are not good quantum numbers already in the infinite volume continuum and also on the lattice. As the first step, we extracted the excited charmonium spectrum within the single hadron approach on the lattice in multiple inertial frames and identified the continuum quantum numbers Padmanath et al. (2019). Such a spin identified single-hadron spectrum tells us where to expect effects from resonances in different partial waves, which in turn provides valuable information for the parameterization of scattering amplitudes. Thus this study allows one to fully exploit the data in moving frames in the present and related future studies.

The current article presents a lattice QCD investigation of charmonium resonances taking into account their most important strong decay, which goes beyond the single hadron approach. We focus on the scattering of a pair of $\bar{D}$ and $D$ mesons in the $l=1$ partial wave, which features charmonium resonances and bound states with $J^{PC}=1^{--}$, to determine the respective masses as well as the couplings with the $\bar{D}D$ scattering channel. We also consider the $\bar{D}D$ scattering in partial wave $l=3$ yielding information on the lowest charmonium resonance with $J^{PC}=3^{--}$, for which LHCb discovered a candidate named $X(3842)$ Aaij et al. (2019).

The $\bar{D}D$ scattering channel in partial wave $l=1$ is the dominant hadronic decay mode of the $\psi(3770)$ (with branching ratio $Br~{}=~{}93\pm 9\%$). It is a well-established vector resonance generally accepted to be predominantly a conventional $\bar{c}c$ state. We therefore assume elastic $\bar{D}D$ scattering, neglecting the influence of all other allowed decay modes of the $\psi(3770)$ and of higher thresholds. The determination of the relevant scattering amplitudes and resonance parameters, such as the mass and decay width of the $\psi(3770)$, will serve as a demonstration of our realization of Lüscher’s finite volume treatment. In contrast to the only previous study Lang et al. (2015), the $\bar{D}D$ scattering amplitude in $p$-wave is determined utilizing the finite volume spectrum from two lattice ensembles with different physical volumes, each in three different inertial frames. The study is performed for two charm quark masses, one below and one above the physical $m_{c}$. In this way we explore the dependence of the spectra on the position of the $\bar{D}D$ threshold.

The recent LHCb discovery of the charmonium $X(3842)$ with likely quantum numbers $J^{PC}=3^{--}$ Aaij et al. (2019) in part motivates our investigation of this narrow state with $\Gamma\simeq 3~{}$MeV. This resonance is also considered since it inevitably affects the spectrum of vector excitations in a finite hyper-cubic volume. While we are able to determine the mass of this charmonium resonance, the finite volume spectrum in our lattice setup is not dense enough to reliably estimate the very narrow width of this state.

The layout of this paper is as follows. In Section II, the relevant quantization condition relating the infinite volume scattering amplitudes to the finite volume energy spectrum is discussed. The details of the lattice setup and related technical information are briefly outlined in Section III. In Section IV, our procedure for the determination of the excited charmonium spectrum and our approach for dealing with discretization errors in hadron observables involving charm quarks is discussed. We summarize our findings in Section V and conclude in Section VI.

II Scattering in a finite box: Lüscher formalism

In a Euclidean box of finite size, the QCD spectrum is discrete. Energy levels corresponding to bound states are affected by exponentially suppressed finite (spatial) volume effects, while a power law dependence on the size of the box is expected for energy levels related to two-particle states above the strong decay threshold. Lüscher’s finite volume formalism and its extensions Lüscher (1986, 1991a, 1991b); Briceño (2014) (see the review Ref. Briceño et al. (2018) to find more references) show that the energy dependence of infinite volume scattering amplitudes determines the discrete energy spectrum in a finite-sized box with periodic boundary conditions.

We are interested in $\bar{D}D$ scattering, i.e. single-channel scattering of spin-zero particles with equal masses, in inertial frames with zero and non-zero total momenta $\vec{P}$. In this case the scattering matrix $S_{l}(p)=e^{2i\delta_{l}(p)}$ can be expressed in terms of a single phase shift $\delta_{l}(p)$ in the partial-wave $l$, where $p$ is the momentum of each hadron in the center of momentum frame. The spin of the system is zero, since the scattering particles are spin-less and the total angular momentum equals the orbital angular momentum ($J=l$).

The quantization condition connecting the infinite volume scattering amplitudes with the finite volume discrete energy $E$ of the eigenstates can be expressed in our case as

[TABLE]

Here the determinant is computed over all partial waves $l$ and $E_{cm}$ is the eigenenergy in the center of momentum frame. The real function $\tilde{K}_{l}$ parametrizes the infinite volume scattering amplitude $S_{l}$ and the transition amplitude $t_{l}$ in the partial wave $l$ as

[TABLE]

with $\rho=2p/E_{cm}$, or equivalently

[TABLE]

The kinematic variable $p$ is the momentum of each $D$ meson in its center-of-momentum frame

[TABLE]

with $s=E_{cm}^{2}$. We follow the definitions of $\tilde{K}$ in Ref. Brett et al. (2018) and the definition of $t_{l}$ used by the Hadron Spectrum Collaboration (see e.g. Ref. Dudek et al. (2014)).222However, unlike in Ref. Brett et al. (2018), we do not divide our energy levels or physical quantities by the mass of the scattering particles. The factor $p^{-2l-1}$ in Eq. (2) ensures a smooth kinematic behaviour of $\tilde{K}$ at the threshold.

The box-matrix $B^{\vec{P};\Lambda}_{l^{\prime}l}(E_{cm})$, as introduced in Ref. Morningstar et al. (2017), depends on $E_{cm}$, the lattice spatial size $L$ and the lattice irreducible representation (irrep) $\Lambda$ under investigation. Its entries are customarily written in terms of the known $\mathcal{Z}$-functions defined in Ref. Lüscher (1986). Note that $B$ also has off-diagonal entries in the partial wave indices $l^{\prime}$ and $l$, due to the reduced rotational symmetries in a hyper-cubic box. The latter leads to the finite volume spectrum for any lattice irrep containing information from scattering in multiple partial waves.

The main focus of the present work are the resonances and bound states of charmonium in the vector channel that are related to $\bar{D}D$ mesons scattering in $p$-wave. In the simplest approximation, in which the influence of higher partial waves $l>1$ are neglected, the quantization condition in Eq. (1) allows one to compute $\delta_{1}$ directly from a given measured energy level in a finite box. In this case, the determinant condition reduces to the well-known Lüscher relation Lüscher (1986)

[TABLE]

for zero total momentum.

In the case when the effects of higher partial waves are non-negligible, the properties of the scattering amplitudes $S_{l}$ are accessible only after finding and fitting a suitable parametrization of all the relevant $\tilde{K}_{l}$, so that the determinant condition in Eq. (1) is solved for each of the energy levels in the finite volume spectrum. The finite volume mixing of even and odd partial waves is not allowed if the scattering particles have equal masses Rummukainen and Gottlieb (1995); Dudek et al. (2013), as in our case. Hence the contamination of higher partial waves in the $l=1$ excitations due to broken rotational symmetry can occur only from odd partial waves $l=3,5,\dots$. The contributions from the $l>5$ partial waves are expected to be strongly suppressed in the near-threshold region explored here and it can be neglected. We will therefore investigate partial waves $l=1$ and $l=3$ (as well as their mutual influence) through the finite volume quantization condition of Eq. (1). We use the package TwoHadronsInBox to compute the determinant in Eq. (1) when fitting the parametrization of the $\tilde{K}$-matrix to our energy levels, following the “determinant residual” method Morningstar et al. (2017).

III Lattice setup

The results presented in this work have been determined from two $N_{f}=2+1$ ensembles labeled U101 and H105 generated by the Coordinated Lattice Simulation (CLS) consortium Bruno et al. (2015); Bali et al. (2016). The inverse gauge coupling is $\beta=6/g^{2}=3.4$ and the lattice spacing is $a=0.08636(98)(40)$ fm Bruno et al. (2017). The lattice volume is $T\times L^{3}=96\times 32^{3}$ for H105 and $128\times 24^{3}$ for U101. The Wilson-clover action  Sheikholeslami and Wohlert (1985) is non-perturbatively improved to remove the $O(a)$ lattice artefacts, with the clover term set equal to $c_{\textrm{sw}}=1.986246$. The strange and light hopping parameters are $\kappa_{l}=0.13697$ and $\kappa_{s}=0.13634079$, respectively, leading to pion and the kaon masses of $m_{\pi}=280(3)$ MeV and $m_{K}=467(2)$ MeV. Note that the light quark mass ($m_{l}$) is heavier than its physical value while the strange quark mass ($m_{s}$) is lighter. This is because the U101 and H105 ensembles lie on a trajectory where the physical point is approached keeping the average quark mass $2m_{l}+m_{s}$ fixed. The gauge and the fermion fields fulfill open boundary conditions in the time direction and periodic boundary conditions in the spatial directions. The total number of configurations analyzed is 255 for the U101 ensemble and 492 for H105, with one configuration taken every 20 and 16 MDUs, respectively, from two different replica in each case.

III.1 Choosing the charm quark mass

The bare charm quark mass, set by the hopping parameter $\kappa_{c}=1/(2m_{c}+8)$ for Wilson-type actions, is a parameter to be determined with care. As we aim to study the strong decay of near-threshold charmonium resonances, the particles relevant for tuning the charm quark mass are the lowest $\bar{c}c$-states, the $J/\psi$ and the $\eta_{c}$, and the decay products of their excited states, such as the $D$ and $D_{s}$ mesons. Away from the physical point, the masses of all these particles cannot be tuned simultaneously to their experimental values. However, the most critical quantity for our current investigation is the position of the strong decay threshold relative to the $J/\psi$ and the $\eta_{c}$, rather than the mass of an individual charmonium state. We aim to explore how charmonia near and above the open-charm threshold depend on its position. We also want to investigate the systematics coming from the tuning of the charm quark mass. To this end we employ two different values of $\kappa_{c}$ given by 0.12522 and 0.12315, corresponding to the $D$ meson mass being below and above its physical value, see Table 3. The main results for the position of the resonances in physical units will be presented relative to the spin average mass of the $J/\psi$ and the $\eta_{c}$,

[TABLE]

This is motivated by the fact that some of the discretization effects (which can be significant at a lattice spacing of $a\approx 0.09$ fm) are canceled when computing energy differences. The spin-average mass is equal to 3103(3) MeV for $\kappa_{c}=0.12315$ and 2820(3) MeV for $\kappa_{c}=0.12522$. The energy gap between $M_{\textrm{av}}$ and 2$m_{D}$ is

[TABLE]

Note that the relative position of the $\bar{D}D$ threshold with respect to $M_{\textrm{av}}$ is closer to the experimental value for $\kappa_{c}=0.12522$ than for $\kappa_{c}=0.12315$.

III.2 Setup of the distillation algorithm

We employ the distillation method to compute our correlation functions Peardon et al. (2009), i.e. the propagation of the quarks is computed in the distillation space of the eigenvectors of the 3D Laplacian. As a consequence, we are able to separate the construction of the interpolators from the computation of the light and charm propagators without resorting to storing full propagators. Only at a later time, the quark-lines of the Feynman diagrams are contracted and summed to obtain the entries of the correlation matrix. We employ 90 and 150 eigenvectors of the lattice Laplacian for the U101 and the H105 gauge field ensembles, respectively. We use full distillation for U101. Full distillation is also used on H105 for the charm propagators and the $\phi$-matrices (defined in Eq. (12) of Ref. Peardon et al. (2009)). The cost of full distillation for the determination of light quark propagators on the H105 ensemble would have been prohibitively expensive, as $150\times 4$ inversions of the Dirac operators on a large volume would have been required for each flavor and each time-source/sink. Therefore, we use stochastic distillation Morningstar et al. (2011) for the light quark propagators on the H105 ensemble only. The perambulator, defined as

[TABLE]

requires to compute the inverse of the Dirac matrix $D$ on a source $|v_{i}^{\mu}(t)\rangle$ which is zero everywhere except for the timeslice $t^{\prime}$ and spin index $\nu$, where it is equal to the Laplacian eigenvector $j$. The number of inversions required can be reduced by considering random sources. We introduce $N_{s}$ complex $Z_{2}$ stochastic noise vectors

[TABLE]

using full-time dilution. The stochastic source given to the multigrid solver is constructed for every time-slice $t$ as

[TABLE]

The perambulator can be approximated as

[TABLE]

which is equal to the exact expression of Eq. (7) in the limit $N_{s}\rightarrow\infty$, since

[TABLE]

Thanks to the use of full-time dilution in our stochastic distillation scheme, there is no bias for contributions to the correlation matrix involving light annihilation diagrams. However, we need to consider two independent sets of random vectors for diagrams involving light quarks propagating at different timeslices. We therefore use two sets of 20 noise vectors, and we average over the two only for light annihilation diagrams. Following this approach, we are able to reduce the number of the expensive inversions of the Dirac operator for the light quarks by a factor ten on our largest volume. At the same time, we retain the advantages of full-dilution for the purely charm contributions and employ full distillation to compute the correlation functions between $\bar{c}c$-type operators, that need to be estimated precisely for the study of charmonium excited states we are interested in.

IV Energy spectrum

IV.1 Determination of the finite volume spectrum

Many excited energy levels need to be accurately determined on several different volumes and/or in different momentum frames in order to extract the scattering matrix according to Lüscher’s method. The excited energy spectrum $E_{n}^{\textrm{{lat}}}$ is extracted from Euclidean time two-point correlation functions

[TABLE]

constructed from a basis of interpolators {$O_{i}$} with the desired quantum numbers. Here $Z_{i}^{n}=\langle 0|O_{i}|n\rangle$ refers to the operator state overlap.

Our main purpose is to investigate $\bar{D}D$-meson scattering in partial wave $l=1$. We determine the charmonium spectrum in three lattice irreducible representations that contain $l=1$ (they all also contain partial wave $l=3$):

[TABLE]

where each irreducible representation corresponds to a different inertial frame, with total momentum $\vec{P}$. The values of $|\vec{P}|^{2}$ are given in units of $(2\pi/L)^{2}$. Note that charge conjugation $C$ is a good quantum number in all frames, whereas parity $P$ is a good quantum number only in the rest frame. In order to reliably extract the relevant low energy spectrum, a large basis of interpolators {$O_{i}$} that includes single-meson as well as two-meson operators in all the irreps is required.

For the single-meson interpolators, we follow the procedure proposed in Refs. Dudek et al. (2010); Thomas et al. (2012) and construct all interpolators with up to two gauge covariant derivatives. As a first step in the analysis, the spectrum is extracted using a basis of these single-meson interpolators and, following the spin identification procedure discussed in Ref. Padmanath et al. (2019), the continuum quantum numbers $J^{P}$ of the levels in the energy region of interest are identified. In this way, we can understand the quantum numbers of the energy levels that are related to a naive $\bar{c}c$ description that enter into our phase shift analysis.

The two-meson interpolators are built using the projection

[TABLE]

where $\vec{p}_{1}$ and $\vec{p}_{2}$ are the momenta of the scattering particles such that $\vec{P}=\vec{p}_{1}+\vec{p}_{2}$, $R$ is a group element of the point symmetry group $G$ of the inertial frame with momentum $\vec{P}$, $T^{\Lambda}_{r,r}(R)$ is the representation for the group element $R$ in the irrep $\Lambda$ and $r$ is the row of the irrep $\Lambda$. As $\bar{D}D$ is the predominant decay mode of the $\psi(3770)$ resonance, we include $\bar{D}D$ interpolators for all momentum combinations within the relevant low energy region. In Table 1, we show the two-meson interpolators considered in our calculation in each of the lattice irreps being analyzed. In this study, we neglect the $p$-wave decay $J/\psi\eta$ and in addition the final states including three or more hadrons, such as $J/\psi\pi^{+}\pi^{-}$ and $J/\psi\pi^{0}\pi^{0}$.

The inclusion of two-meson interpolators in the basis brings additional Wick contractions. We compute all relevant Wick contraction diagrams for the single- and two-meson operators (see Fig. 1 of Ref. Lang et al. (2015)). We exclude diagrams with disconnected charm quark contractions, that lead to the decay of charmonium into light hadrons. The correlation functions constructed from the single and two meson operators are averaged over the spin and momentum polarization to increase the statistical precision. We also bin our data in blocks of size equal to one, two and four to correctly address the autocorrelation of our configurations. For further details, see Appendix A.

The extraction of the energy levels proceeds via a variational approach Michael (1985); Lüscher and Wolff (1990); Blossier et al. (2009) by solving a generalized eigenvalue problem (GEVP)

[TABLE]

for the correlation matrices in Eq. (12). We set $t_{0}=2$ in our analysis. The eigenvalues $\lambda_{n}(t)$ are fitted to a single or a double exponential function and the quality of the fits is compared to estimate the best fitting intervals. We remark that increasing $t_{0}$ leads to consistent results, however, the errors also increase. Corrections to these fit forms arise when using open boundary conditions in the temporal direction and the source or sink timeslice of the correlation function is close to one of the boundaries. However, all our measurements are made in the bulk of the lattice, 28 time slices away from either boundary, and no such corrections are observed.

In Figs. 1 and 2, we present the finite volume excited energy spectrum (above the $J/\psi$) in all the three frames studied for $\kappa_{c}=0.12522$ and $\kappa_{c}=0.12315$, respectively. The statistical uncertainty is smaller on the H105 ensemble, in part due to the larger number of configurations analyzed compared to U101. The color coding of the levels indicates the quantum numbers. Blue levels are related to the $\psi(2S)$ bound state and the green levels are identified to have quantum numbers $J^{PC}=3^{--}$.333This state appears twice for $\vec{P}^{2}=2$ according to Thomas et al. (2012) and Table I of Padmanath et al. (2019): one level has helicity zero and the other helicity two.. The remaining levels in red arise from charmonia with $J^{PC}=1^{--}$ and from $\bar{D}D$ scattering levels shifted due to interaction. We extract up to 4 or 5 excited states for each lattice irrep considered.

Within the energy range explored, we also observe a $2^{--}$ level in the $A_{1}^{-}(|\vec{P}|^{2}=2)$ spectrum. The helicity two components of the $2^{--}$ state appear in this irrep as can be seen from Table II of Thomas et al. (2012) and Table I of Padmanath et al. (2019). Taking the naive energy level for $\kappa_{c}=0.12522$ and using the continuum dispersion relation, we determine the mass of this state in its rest frame to be 3825(8) MeV, which is consistent with the lowest $2^{--}$ charmonium state observed by Belle Bhardwaj et al. (2013) and BESIII Ablikim et al. (2015). The coupling of this excitation with the rest of the spectrum is expected to be negligible as, in a study of vector channel scattering of two equal mass pseudoscalar mesons, $J^{P}=2^{-}$ is an unnatural quantum number and $l=2$ is an irrelevant partial wave. This level is absent when we omit the relevant $2^{--}$ interpolators from the operator basis, while the rest of the spectrum remains unaffected. Hence we exclude this level from the amplitude analysis and also do not show it in the figures.

The thresholds for the channels $D_{s}D_{s}$ and $\bar{D}D^{*}$, which are omitted in our simulations, are shown in Fig. 1 and 2. These thresholds are lower than some of the employed energy levels on the $N_{L}=24$ ensemble, however we expect that the effect of the omitted channels is very small because they appear in $p$-wave and they open at relatively high energy. Note that there are no close-by two-particle levels of the omitted channels because they appear in $p$-wave. The influence of the omitted channels on the results from the $N_{L}=32$ ensemble is expected to be negligible444We have verified that masses and couplings of the extracted charmonia from the $N_{L}=32$ ensemble alone are compatible with the ones based on combined fits of all volumes..

IV.2 Approach to dispersion relations and charm-quark discretization errors

The lattice energy levels presented in Figs. 1 and 2 are inputs to the quantization condition given in Eq. (1) (realized using the TwoHadronsInBox package of Ref. Morningstar et al. (2017)), which is based on continuum dispersion relations. However, lattice energy levels are subject to discretization effects and can deviate from continuum expectations. Such lattice artifacts are larger in magnitude for charm than for the light quarks $u,d$ and $s$, since, in dimensionless units, $1\sim am_{c}\gg am_{u/d,s}$ even on many currently available lattices. At finite lattice spacing $a$, there are non-negligible deviations of the dispersion relation $E^{\textrm{lat}}_{D}(p)$ of the $D$-meson from its continuum counterpart

[TABLE]

As is evident from Table 2, there are non-negligible deviations from $E^{\textrm{cont}}_{D}(\vec{p})$ particularly for the smaller physical volume (for which the lattice momentum is larger), and hence care is required when extracting the scattering matrix in the $\bar{D}D$ channel. In this section we describe our approach to this issue.

First, we determine the energy shift of each interacting eigenstate with total momentum $\vec{P}$ with respect to the nearest non-interacting state $\bar{D}(\vec{p}_{1})D(\vec{p}_{2})$

[TABLE]

where $\vec{p}_{1,2}=\vec{n}_{1,2}\tfrac{2\pi}{L}$ and $\vec{p}_{1}+\vec{p}_{2}=\vec{P}$. All three energies on the r.h.s. are extracted from their corresponding eigenvalues of the GEVP, separately. Here, $(E^{\textrm{lat}})_{s}$ is the energy of the interacting $\bar{D}D$ system, presented in Section IV.1, while $(E^{\textrm{lat}}_{D(\vec{p}_{1})})_{s}$ is the energy of a single $D$ meson with momentum $\vec{p}_{1}$ measured on the lattice. All are extracted on a given bootstrap or jackknife sample $s$. If the scattering matrix is equal to the identity, the phase shift is equal to zero and hence the energy shifts in Eq. (16) should be equal to zero. Given that the quantization condition is based on the continuum dispersion relation,555The Lüscher $\mathcal{Z}$-functions, that enter the box matrix $B$ in the quantization condition, have poles at energies $E=E^{\textrm{cont}}_{D(\vec{p}_{1})}+E^{\textrm{cont}}_{D(\vec{p}_{2})}$ with $\vec{p}_{1,2}=\vec{n}_{1,2}\tfrac{2\pi}{L}$ and $\vec{p}_{1}+\vec{p}_{2}=\vec{P}$, while $\tan\delta$ in the Lüscher-type relations is a sum of terms with $\mathcal{Z}$-functions in the denominator. Therefore $\delta=0$ at those energies $E$. we ensure this important constraint is satisfied by using the energies

[TABLE]

as input to Eq. (1). The energies $(E^{\textrm{cont}}_{D(\vec{p}_{1})})_{s}$ and $(E^{\textrm{cont}}_{D(\vec{p}_{2})})_{s}$ are computed from Eq. (15) using the lattice momenta $\vec{p}_{1,2}$ and the $D$-meson mass at rest (also determined from the bootstrap or jackknife sample $s$). Note that, the energies $E^{\textrm{calc}}$ are equal to $E^{\textrm{lat}}$ in the continuum limit $a\rightarrow 0$ by construction (Eqs. (16) and (17)). Furthermore, the above $E^{\textrm{calc}}$ also clearly provides a zero phase-shift when $\Delta E$ in Eq. (16) is zero, even if non-negligible heavy-quark discretization effects are present in our simulation.

V Results for charmonia with $J^{PC}=1^{--}$ and $3^{--}$

The ground state in the vector channel is the $J/\psi$, while the first excited state, the $\psi(2S)$, lies below the open charm decay threshold. Beyond these two states, there is a relatively narrow resonance, the $\psi(3770)$ with $\Gamma^{\textrm{exp}}\approx 27$ MeV, and a broader resonance, the $\psi(4040)$, both which can decay to $\bar{D}D$ in $p$-wave.

The near-threshold charmonium resonance $\psi(3770)$ is the main target in this investigation. In particular, we aim to extract the pole position of the $\psi(3770)$ and its coupling to $\bar{D}D$, which is related to the $\psi(3770)$ decay width. The study of the vector channel provides a benchmark for our simulations enabling us to assess systematic uncertainties which will also arise when considering more complex cases where several decay channels are involved (and for which the experimental situation is unclear). We assume that the $\psi(3770)$ resonance is only coupled with the $\bar{D}D$ scattering channel and hence its parameters can be completely determined from the elastic $\bar{D}D$ scattering amplitude. Experimentally, the $\psi(3770)$ decays into $\bar{D}D$ with a branching ratio of $93~{}_{+8}^{-9}~{}\%$ Tanabashi et al. (2018). In this work, we neglect the effects from other hadronic decay modes, such as $J/\psi\eta$, $J/\psi\pi^{0}\pi^{0}$ and $J/\psi\pi^{+}\pi^{-}$, which collectively have a branching ratio below 0.2%. We also neglect decays into light hadrons through charm-annihilation.

Our second aim is to study the lowest $3^{--}$ charmonium resonance. A candidate for this state referred to as the $X(3842)$666The LHCb collaboration has not identified the quantum numbers of this state but assumes $J^{PC}=3^{--}$ as the mass and width of the $X(3842)$ fit model expectations. was recently discovered by LHCb in $\bar{D}D$ decay and has width $\Gamma^{\textrm{exp}}\approx 2.8$ MeV Aaij et al. (2019).

V.1 Fits of the phase shift

We perform a combined fit of the energy levels calculated on the U101 and H105 ensembles using bootstrapping for the determination of the uncertainties (see Appendix A for details). The ground state level in each frame, corresponding to $J/\psi$, is excluded from our fits. We include the energy level corresponding to the $\psi(2S)$ to use the maximal information from our spectrum and to correctly describe the spectrum in the vicinity of the threshold. The fitting forms we have explored for the vector channel are the “double-pole”

[TABLE]

inspired by two states in the region of interest, the $\psi(2S)$ and the $\psi(3770)$777Here $\cot(\delta_{1})=0$ at $s=m_{1}^{2}$ and $s=m_{2}^{2}$. Between these two energies, $\cot(\delta_{1})$ has a pole and $T$ has a zero at $s=\frac{G_{2}^{2}m_{1}^{2}+G_{1}^{2}m_{2}^{2}}{G_{1}^{2}+G_{2}^{2}}$. The presence of a pole in $\cot(\delta_{1})$ allows to fulfill easily the consistency check in Sec. V.3 and the requirement that the parametrization should not predict unobserved energy levels.; and the “quadratic” form

[TABLE]

where $p^{2}=s/4-m_{D}^{2}$. Other parameterizations were investigated. A linear form did not describe the data well; a triple-pole form would require additional data at large momenta to enable a third pole to be resolved; a double pole form with an additional constant is described in Appendix C.

Although many partial wave amplitudes in principle contribute to the finite volume spectrum in any given irrep starting at the $\bar{D}D$ threshold, the effect of higher partial waves should be small, unless there is a resonance contributing to a higher partial wave in the energy region of interest. The spin-identified finite volume spectrum for charmonium from our previous investigation Padmanath et al. (2019) indicates a low lying $J^{PC}=3^{--}$ resonance, that contributes to the $l=3$ partial wave. In order to investigate the effects of such a low-lying $l=3$ partial wave excitation on the rest of the discrete spectrum, we parametrize it with a simple Breit-Wigner form given by

[TABLE]

Our fits include all the fourteen excited energy levels we have been able to determine on the H105 and the U101 ensembles. For both $\kappa_{c}$, the fits are more constrained by the data from the H105 ensemble. The “double pole” fit is presented for the $1^{--}$ channel in Fig. 3(a) for $\kappa_{c}=0.12522$ and in Fig. 3(b) for $\kappa_{c}=0.12315$. For this parametrization (and Eq. (20)), we have been able to fit all our data with

[TABLE]

for $\kappa_{c}=0.12522$ with a $\chi^{2}/\textrm{d.o.f.}=0.5$, while we have

[TABLE]

for $\kappa_{c}=0.12315$ with a $\chi^{2}/\textrm{d.o.f.}=0.3$.

The quadratic fit is presented for the $1^{--}$ channel in Fig. 4(a) for $\kappa_{c}=0.12522$ and in Fig. 4(b) for $\kappa_{c}=0.12315$. Compared to the double pole fit form, the minimum of the correlated $\chi^{2}$ is more unstable across different bootstrapping resamples, resulting in an asymmetric distribution of the fitted parameters with long tails. We therefore quote asymmetric uncertainties for this parametrization. The fitted quadratic functional form for $\kappa_{c}=0.12522$ is equal to

[TABLE]

with a $\chi^{2}/\textrm{d.o.f.}=1.32$, while for $\kappa_{c}=0.12315$

[TABLE]

with a $\chi^{2}/\textrm{d.o.f.}=1.15$.

The quality of our fits for the combined channel case is presented in Figs. 5, 6, 7 and 8, in terms of the $\Omega$ function,

[TABLE]

as introduced in Ref. Morningstar et al. (2017) to minimize the $\chi^{2}$ in the “determinant residual method”. Here the matrix $A(s)$ is the argument of the determinant in Eq. (1) and $\mu=8$ is a regularization parameter, see Appendix A for further details on the choice of $\mu$. Each crossing of the $\Omega$ function with the vertical axis $\Omega=0$ represents a predicted energy level from our parametrization, which can be compared with our observed energy spectrum. The number of energy levels is indeed equal to the number of solutions of $\Omega=0$ in the energy region between $E\simeq m_{\psi(2S)}$ (the blue energy level) and the highest level. In particular, there are two approximately degenerate energy levels when the fitted $\Omega$-function has a double zero (this occurs for $3^{--}$ when $P^{2}=2$). We remark that our scattering amplitude captures the energy dependence only in the energy region where we have data and is not applicable outside of this region. There is an additional crossing of $\Omega=0$ in certain irreps for $\kappa_{c}=0.12315$ in Fig. 5, however, it occurs for $E<m_{\psi(2S)}$. The amplitude is not constrained and is, therefore, not applicable there. 888This additional crossing would not occur if the slope of $p^{3}\cot\delta_{1}/\sqrt{s}$ was slightly bigger (in terms of the absolute value) in the energy region $E<m_{\psi(2S)}$ of Fig. 9(b) (see also Appendix C for details).

The parameters for the $3^{--}$ resonance are quite unstable across the bootstrap samples, resulting in large uncertainties, especially for the coupling $g_{3}$. To check the stability of the fits in the vector channel, we have also performed an alternative analysis: the levels related to $J^{PC}=3^{--}$ are identified (following the methods utilized in Ref. Padmanath et al. (2019)) and excluded from the fits. The resulting couplings and masses for the $1^{--}$ resonance are found to be compatible with the values quoted above and we can therefore conclude that the influence of the $3^{--}$ resonance on the $\psi(3770)$ is negligible in our simulation.

V.2 Bound states, virtual bound states and resonances

Experiments and lattice QCD determine the scattering amplitudes $t_{l}$ for real energies. The physical interpretation in terms of (virtual) bound states and resonances is however conventionally performed by looking at poles of the $t_{l}(s)$ in Eq. (2), continued to the complex $s$-plane

[TABLE]

The two Riemann sheets I and II, from the square root branch cut are defined to have Im$(\rho)>0$ and Im$(\rho)<0$, respectively999The square-root branch cut for $\rho$ is chosen in the conventional sense, such that it runs from the threshold $s=4m_{D}^{2}$ to $+\infty$ along the positive real axis. We neglect the effects from the left hand cut and from the freedom in choosing the real part of $i\rho$ in Eq. (26) (such as using a Chew-Mandelstam phase space), either of which results in pole structure deep below the energy region being studied.. The product

[TABLE]

is an analytic function of $s$ only, where our lattice results for $p^{2l+1}\cot\delta_{l}/\sqrt{s}$ were expressed in terms of $s$ in Section V.1.

The scattering amplitude $t_{l}$ in Eq. (26) has a pole at a given $p$ if $\cot\delta(p)=i$. A bound state appears as a pole of the $t_{l}$ matrix at $p=i\kappa_{B}$, i.e. on the real axis below the scattering threshold on the first Riemann sheet. A virtual state appears as a pole at $p=-i\kappa_{B}$, i.e. on the real axis on the second Riemann sheet. In both cases $\kappa_{B}$ denotes a real positive number. It is convenient to express these pole conditions in terms of the quantity $p^{2l+1}\cot\delta$, which is considered in the Section V.1 for $l=1,3$. Such poles appear if the conditions

[TABLE]

are fulfilled. A resonance corresponds instead to a pair of poles away from the real axis above threshold on the second sheet, influencing the physical region in the first sheet.

V.3 Consistency check for a bound state in partial wave $l$

The $S(p)=e^{2i\delta(p)}$ matrix for a bound state has a pole at $p=i\kappa_{B}$, as discussed in the previous section. The analyticity of the $S(p)$ matrix near the real energy axes implies the sign of the pre-factor in front of the pole to be

[TABLE]

where $\beta^{2}$ is a real and positive number. This functional form is derived from the general form of the solutions for $S(k)$ continued to complex $p$ in Ref. Sitenko (1971) (see Eq. (7.60)) and applied for $s$-wave in, for example, in Ref. Iritani et al. (2017). Equation (30) must be satisfied for all bound states, but it does not apply to virtual bound states. We have explicitly verified that the relation (30) with positive $\beta^{2}$ is satisfied for $s$-wave and $p$-wave bound states in 3D non-relativistic quantum mechanics for various forms of the potential (square well, Yukawa, Wood-Saxon, and several other forms of the type $V(r)=\frac{1}{r^{n}}\exp{(-Ar^{m})}$).

Let us express the condition in Eq. (30) in terms of $p^{2l+1}\cot(\delta_{l})$ in Eq. (28) which was fitted from our data. For this purpose we consider the dependence of the quantity

[TABLE]

on $p^{2}$ near the bound state pole, where it equals zero in Eq. (28). Inserting in Eq. (31) $\cot\delta=i(S+1)/(S-1)$ with $S$ from Eq. (30), and expressing $p$ with $p^{2}$ as $p=i|p|=i\sqrt{-p^{2}}$, we obtain

[TABLE]

Taking the derivative of the previous expression with respect to $p^{2}$ and evaluating it at $p^{2}=-\kappa_{B}^{2}$, we get

[TABLE]

The above condition implies that the slope of $p^{2l+1}\cot\delta$ in terms of $p^{2}$ has to be smaller101010Note that this applies for the values of these quantities, not for absolute values of these quantities. than the slope of $-(p^{2})^{l}\sqrt{-p^{2}}$ at the position of the bound state in partial wave $l$, while it does not apply for virtual bound states. Note that the position of the bound state is where these two curves cross. The condition for $s$-wave was referred to as a sanity-check in Ref. Iritani et al. (2017).

Let us now turn to the fitted phase shifts in Section V.1 and investigate which results satisfy the condition of Eq. (30) or equivalently the consistency check of Eq. (33). The slope of $p^{2l+1}\cot(\delta)$ (the red curves in Fig. 9) has to be smaller than the slope of $-(p^{2})^{l}\sqrt{-p^{2}}$ (the blue curve in Fig. 9) where these two curves cross. This is satisfied for the double-pole fits, while it is not satisfied for the quadratic fit as can be seen from Fig. 9. Note that all these slopes are negative, and in this case the absolute value of slope for $p^{2l+1}\cot(\delta)$ has to be larger than the absolute value of slope for $-(p^{2})^{l}\sqrt{-p^{2}}$.

V.4 Bound states in the vector channel, and mass and coupling of the $\psi(3770)$

The physical parameters for bound states and resonances are based on our fits of $D\bar{D}$ scattering in p-wave. The first step is to ensure that our fitted parametrizations fulfill all the physical requirements, and in particular the consistency check given by the condition in Eq. (33). The double-pole fits in Section V.1 satisfy this condition, while the quadratic fits do not, as explained in the previous subsection. Therefore, in the following, in terms of the physics conclusions we focus on the double-pole fits.

The absolute value of the resulting amplitude $t_{l=1}(s)$ (Eqs. (26) and (27)) in the complex energy plane is shown in Figs. 10 and 11 for $\kappa_{c}=0.12522$ and $0.12315$, respectively. The locations of these poles on both Riemann sheets are collected in Fig. 12. These poles are related to (virtual) bound states and resonances according to the general criteria from Section V.2. The final masses of various states will be quoted according to

[TABLE]

where the mass splitting $[m-M_{\textrm{av}}]^{\textrm{lat}}$ with respect to $M_{\textrm{av}}$ of Eq. (6) is determined from the lattice. Figure 13 and Table 3 summarize our results for various states, which are detailed below:

• •

$\mathbb{\kappa_{c}=0.12522}$ ($m_{D}=1762(2)$ MeV): For our lighter charm quark mass there is a complex conjugated pair of poles in the second sheet corresponding to the $\psi(3770)$ resonance, see Fig. 10(b). Figure 10(a) shows a single pole of the $t$-matrix on the real axis in the first Riemann sheet, corresponding to the bound state $\psi(2S)$. In addition, there is a virtual bound state on the real axis on the second sheet.

The parameters for the $\psi(3770)$ resonance are extracted from the linearized (Breit-Wigner) behavior of Eq. (21) in the resonance region

[TABLE]

The resonance width $\Gamma$ can not be directly compared to experiment since the phase space (proportional to $p^{3}$) for our unphysical $D$ meson mass is different to that in experiment. Therefore we extract the dimensionless coupling $g=16.0({}^{+2.1}_{-0.2})$ that parametrizes the $\psi(3770)\to\bar{D}D$ width. This coupling agrees within errors with the experimental value $g^{\textrm{exp}}=18.7(9)$ obtained from $\Gamma^{\textrm{exp}}$. The resonance mass equals $m_{2}$ of Eq. (18). The fitted value is given in Eq. (21) and (utilizing Eq. (34)) this yields a final mass $m_{2}=3780(7)$ MeV, which is consistent with the experimental value of $3773.13(35)$ MeV.

These parameters are verified by also considering the location of the resonance pole in the amplitude of Eq. 21, located at $m_{2}^{pole}+\tfrac{1}{2}\Gamma^{\textrm{pole}}=3780(4)-0.97(63)i$ MeV. The pole is plotted in Fig. 12. The value of $\Gamma^{pole}$ results from the formula on the right in Eq. (35) to a coupling $g=17(9)$111111The phase space $p^{2}/s$ in (35) is extracted from $m^{\textrm{pole}}\simeq 3780$ MeV and threshold position 3772 MeV in Fig. 12 (phase space in simulation gives also compatible result)., consistent with the value quoted in the previous paragraph.

The mass of the $\psi(2S)$ bound state is obtained from $p^{3}\cot(\delta)/\sqrt{s}$ of Eq. (21) by finding a pole for real energies below threshold according to Eq. (28): it corresponds to the crossing of the red and blue curves in Fig. 9(a). The resulting mass in Table 3 is about $20~{}$MeV smaller than the experimental mass.

• •

$\mathbb{\kappa_{c}=0.12315}$ ($m_{D}=1927(2)$ MeV): $\psi(3770)$ is a shallow bound state for our heavier charm quark mass, in contrast to experiment where it is a strongly decaying resonance. It corresponds to a near-threshold pole in the first sheet of Fig. 11(a), arising from the fact that the crossing of our fits with the real axis occurs on the left of the axis origin (below threshold) in Fig. 3(b). The mass of the bound state $m=3776(7)$ MeV in Table 3 is determined according to Eq. (28) and corresponds to the near-threshold crossing of the red and blue lines in Fig. 3(b). The $\psi(3770)$ cannot strongly decay into $\bar{D}D$, but we still extract the coupling $g=18.9({}^{+0.8}_{-0.7})$, which parametrizes the slope of $p^{3}\cot\delta_{1}/\sqrt{s}$ near the real axes according to Eq. (35).

The $\psi(2S)$ corresponds to the lighter of the two bound states. Its mass $m\simeq 3687~{}$MeV, determined from the bound-state condition of Eq. (28), has a large error due to the huge uncertainty of the coupling $G_{2}$ in the double-pole parametrization in Eq. (22). This mass is consistent within errors with the value $m=3665(9)~{}$MeV obtained directly from the first-excited energy-level at $\vec{P}=0$ and $L=32$. This agreement is expected as the influence of the threshold on a state which lies about 150 MeV below is very small. We quote the latter value as the final result for the $\psi(2S)$ mass at this $\kappa_{c}$.

Our results for the masses of the $\psi(3770)$ and $\psi(2S)$ states are summarized and compared to experiment in Fig. 13 and Table 3. The mass $m$ and the coupling $g$ for the $\psi(3770)$ agree with the experimental values within errors, which applies to both charm-quark masses considered. The mass of the $\psi(2S)$ is about $2\sigma$ below the experimental value, which is not unreasonable given the lack of continuum and quark-mass extrapolations to the physical values.

The above results (that are based on a combined fit to two different ensembles) are subject to a technical subtlety. Both ensembles have been generated with the same bare lattice parameters, but with different lattice sizes. Hence sub-leading exponentially suppressed finite volume effects, which are generally neglected in the Lüscher finite volume formalism, can sometimes lead to different physical situations for different box sizes. Such a delicate situation can happen when a pole singularity in the complex energy plane is expected to be very close to the threshold: a small volume may exhibit different physics from that of a larger volume. For example, in our case for the lighter $D$-meson mass, we see that our $N_{L}=24$ data have a tendency to prefer a bound state (unlike the $N_{L}=32$ data), although the data are compatible with a resonance within the large errors. In order to confirm that our combined fit does not lead to conclusions that are different from the physical situation in the larger volume (where the neglected exponential corrections are small) we performed an analysis of the $N_{L}=32$ data alone. We find the mass and the coupling from such an analysis to be $m_{2}=3782(7)$ MeV and $g=15({}^{+2}_{-1})$, which are in agreement with the results from our combined fits.

V.5 $J^{PC}=3^{--}$ resonance

The LHCb collaboration recently reported the discovery of the $X(3842)$ charmonium state in the $\bar{D}D$ invariant mass Aaij et al. (2019) with a mass $m^{\textrm{exp}}\simeq 3842~{}$MeV, a very narrow width $\Gamma^{\textrm{exp}}\simeq 2.8~{}$MeV and likely quantum numbers $J^{PC}=3^{--}$. The presence of such a $3^{--}$ resonance was explored by fitting the parametrization of the $\tilde{K}$ matrix in Eq. (1) to our energy levels including the Breit-Wigner form in Eq. (20) for $l=3$.

The resulting mass for $\kappa_{c}=0.12522$ (see Eq. (21)) is

[TABLE]

The mass obtained from our lattice simulation is consistent with the mass of the discovered $X(3842)$ state. This adds confidence to the interpretation of the $X(3842)$ as a $3^{--}$ charmonium resonance. The results for both values of the charm quark mass are compared to experiment in Table 3. Our simulation is not sensitive to the width of this very narrow resonance since there are no $\bar{D}D$ eigenstates in the very narrow energy region $m\pm\Gamma$.131313This applies for the width $\Gamma=\Gamma^{\textrm{exp}}$, while the width $\Gamma$ is even narrower on our lattice due to smaller phase space. In particular, the finite volume lattice energy by itself would lead to compatible results for the mass of the $X(3842)$. For this reason we can not compare the width or the $X(3842)\to\bar{D}D$ coupling with experiment.

V.6 Comparison to previous results

Elastic $\bar{D}D$ scattering and its relation to the $\psi(2S)$ and $\psi(3770)$ resonances have previously been investigated in a lattice study in Ref. Lang et al. (2015) for $m_{\pi}\simeq 156~{}$MeV and $266~{}$MeV. In this subsection we compare our current study to these older results, focusing on those from the $m_{\pi}\simeq 266~{}$MeV ensemble, since the pion mass is similar and the errors are smaller. The most significant improvement over the calculation in Ref. Lang et al. (2015) lies in the use of two lattice volumes at the same lattice spacing and light-quark masses, combined with calculations in moving frames. The latter is facilitated by a new approach to dealing with discretization effects described in Section IV.2 along with the spin identification techniques pioneered in Ref. Dudek et al. (2009) and employed on our ensemble(s) in Ref. Padmanath et al. (2019). The use of two volumes and several moving frames results in a significant increase in the number of energy levels obtained, which renders more information on the dependence of the $l=1,3$ scattering amplitudes on $E_{cm}$.

One important consequence is that the $J^{PC}=3^{--}$ resonance, which contributes to the same irreducible representations on a hyper-cubic lattice, can now be included in our fits. In Ref. Lang et al. (2015) it had to be assumed that after omission of an energy level related to the presence of this resonance, its effect on the remaining finite volume energies was negligible. For our current ensembles this assumption can instead be tested, and our results agree with the mass of the recently observed $X(3842)$ Aaij et al. (2019).

The $l=1$ scattering amplitude in Ref. Lang et al. (2015) was constrained only at three values of $E_{cm}$ and the results for $\psi(2S)$ and $\psi(3770)$ were summarized in Table V of Ref. Lang et al. (2015) for two types of fits. Fit (ii) in Eq. (5.6) was similar to the ”quadratic fit” used in the present work and does not satisfy the consistency check of Eq. (30) for the $\psi(2S)$ bound state (this constraint on the physical parameterizations was not checked in Ref. Lang et al. (2015)). Fit (i) was the Breit-Wigner fit of the energy region above threshold and rendered the $\psi(3770)$ resonance parameters summarized in the right column of Table 3; the mass is consistent with our present result while the coupling is somewhat smaller. The $\psi(2S)$ mass $m=3674(6)~{}$MeV obtained from the naive energy level (second level in Table III of Ref. Lang et al. (2015)) agrees within errors with our present result. The two-pole fit with four parameters could not be preformed using just three energy levels in Ref. Lang et al. (2015).

Another advantage of the present study is the use of the CLS gauge ensemble library, which in principle enables us to extend the current study to further ensembles with different lattice spacings and light quark masses, while maintaining our current setup.

VI Conclusions

We determined the elastic $\bar{D}D$ scattering amplitude in the energy region of the $\psi(2S)$ bound state and the $\psi(3770)$ resonance from lattice QCD simulations using the finite volume Lüscher method. The results we obtain are from simulations with 2+1 flavors of dynamical light quarks at unphysical masses and at a single lattice spacing. We investigate different parametrizations of the $p$-wave scattering amplitude and impose a consistency check on the extracted bound state poles to ensure the physical bound state constraint in Eq. (33). Using a “double pole” parametrization defined in Eq. (18), the bound state mass of the $\psi(2S)$ and the resonance parameters of the $\psi(3770)$ have been extracted from a simultaneous fit to the energy spectrum from three different moving frames and two lattice volumes.

The lattice study of a near-threshold charmonium state away from the physical point requires a careful choice of the charm quark mass for the strong decay to $\bar{D}D$ to be allowed. To this end, on each gauge field ensemble, we employ two charm quark mass values leading to $D$-meson masses 100 MeV below and 60 MeV above the physical value. The masses of the $\psi(2S)$ and $\psi(3770)$, measured as splittings from the spin-averaged ground-state mass, are in good agreement with experiment for these values of the heavy quark mass. The result for the coupling $g$ of the $\psi(3770)$ resonance with the $\bar{D}D$ scattering channel is compatible with the experimental value, even for the charm quark mass where the $\psi(3770)$ is a bound state. The results from both heavy quark masses employed in this work, the results from previous existing work and the corresponding experimental values are summarized in Table 3.

We also investigate the $\bar{D}D$ scattering in partial wave $l=3$ that contains $3^{--}$ resonances. One of the aims is to investigate how this partial wave contributes to the finite volume spectra. Including the lattice energy levels related to the lowest $3^{--}$ state in the amplitude analysis and parametrizing the corresponding $\bar{D}D$ scattering amplitude with a Breit-Wigner form (Eq. 20), we are able to extract the mass of this state. Nice agreement is found with the mass of the recently discovered $X(3842)$, which is interpreted as a $3^{--}$ charmonium resonance. However, we are unable to reliably determine its coupling with the $\bar{D}D$ scattering channel as the finite volume spectrum of $\bar{D}D$ scattering channels are not sufficiently dense to be sensitive to its decay width. We also find our results for the vector channel remain unchanged with the exclusion of the $l=3$ partial wave from the amplitude analysis and hence conclude the influence of the low lying $3^{--}$ resonance on the lattice estimates for the $\psi(3770)$ resonance parameters is negligible in our calculation.

The main challenge for future determinations of scattering amplitudes in the charmonium spectrum will be the presence of multiple coupled channels for the study of conventional and exotic resonances away from the $\bar{D}D$-threshold. A first preliminary analysis of the scalar channel has been presented in Ref. Collins et al. , and we plan to complete the study in a forthcoming publication. In addition it would be desirable to obtain a denser set of points by using more and/or larger volumes, to reduce the discretization effects by calculations on CLS ensembles with a smaller lattice spacing, and to investigate the approach to the physical point by utilizing ensembles with lighter pion masses. In particular, discretization errors are an important source of systematic uncertainty of our current work. Using the CLS trajectory with physical strange-quark mass Bali et al. (2016) would be particularly attractive, as the splitting between the $D$ and $D_{s}$ mesons on these ensembles is closer to the physical splitting and therefore results in a larger energy window for elastic $\bar{D}D$-scattering.

VII Acknowledgements

We thank our colleagues A. Schäfer and S. Weishäupl for help and support during the execution of the project. We thank Ben Hörz for his help and support with the package TwoHadronsInBox and C. B. Lang for his variational analysis codes. We are also thankful to Sinya Aoki, Zohreh Davoudi, Jozef Dudek, Andrew Hanlon and Akaki Rusetksy for valuable discussions. We are grateful to the Mainz Institute for Theoretical Physics (MITP) for its hospitality and its partial support during the course of this work. We use the multigrid solver of Refs. Heybrock et al. (2014, 2016); Richtmann et al. (2016); Georg et al. (2017) for the inversion of the Dirac operator. Our code implementing distillation is written within the framework of the Chroma software package Edwards and Joo (2005). The simulations were performed on the Regensburg iDataCool cluster, and the SFB/TRR 55 QPACE 2 Arts et al. (2015) and QPACE 3 machines. We thank our colleagues in CLS for the joint effort in the generation of the gauge field ensembles which form a basis for the here described computation. The Regensburg group was supported by the Deutsche Forschungs-gemeinschaft Grant No. SFB/TRR 55 and the European Union Innovative Training Network Grant No. 813942 (EuroPLEx). M. P. acknowledges support from the EU under grant no. MSCA-IF-EF-ST-744659 (XQCDBaryons). S. Prelovsek was supported by Slovenian Research Agency ARRS (research core funding No. P1-0035 and No. J1-8137).

Appendix A Error determination

A combined fit of the phase shift from the energy spectrum in various lattice irreducible representations and ensembles requires a careful consideration of all the sources of correlations and systematic uncertainties. The fitting interval used to extract the energy levels from Eq. (14) is a well-known source of systematic error. An example of our fits is presented in Fig. 14. We have performed fits for several different intervals and compared the stability of the raw energy shifts as defined in Eq. (16). We observed that the dependency of the energy shifts on the fitting interval of the scattering $D$-mesons is relatively small. There is a dependence on the starting point of the fit of the excited charmonium state, that stabilizes when the excited state contamination is under control. The raw correlators are binned to take into account the autocorrelations in our data and we have chosen a bin-size equal to two, after observing that the error is stable if larger bin-sizes are used.

Jackknife and bootstrapping are methods widely employed in lattice QCD to determine the errors on secondary quantities extracted or fitted from the raw data. Each bootstrapping resample consists of 255 and 492 randomly chosen configurations from the U101 and H105, respectively. The first step in the analysis is the determination of the covariance matrix for the correlated $\chi^{2}$ minimization of our $K$-matrix parametrization. This is estimated simply as the standard correlation between two variables as given by bootstrapping. The minimum of the correlated $\chi^{2}$, defined from the $\Omega$-function of Eq. (25) as given in Eqs. (38) and (44) of Ref. Morningstar et al. (2017), is recomputed for each bootstrapping resample. We extract from here the distribution of secondary derived quantities such as the parameters of the $K$-matrix or the positions of the poles in the complex plane. The $\Omega$-function depends on an additional regularization parameter $\mu$, and we have verified that varying $\mu$ in the range between 2 and 16 does not change the central value nor the error of the fitted parameters significantly. The crossing with the vertical axis, determining the expected spectrum, does not change significantly if $\mu$ is varied in the above interval, see Fig. 16. Note however that the value of $\mu$ determines the slope and the curvature of the $\Omega$-function at the crossing.

Given that our fits combine two different ensembles with a different number of configurations, single elimination jackknife is not an optimal choice. However, it is worth to compare bootstrapping and jackknife if we restrict the fits to a single ensemble, say U101. Looking to Fig. 15, we can see that jackknife (for which the spread of the resamples is smaller by a factor $\sqrt{N_{\textrm{conf}}}$ compared to bootstrapping) determines the error as a linear approximation of the Lüscher $\mathcal{Z}$-functions around the central value of the phase shift. It might happen however that a given energy level lies very close to a turning point of the $\mathcal{Z}$-function, as in the case of the four rightmost points of Fig. 15. In this situation, the strong non-linear nature of the phase-shift analysis could result in a systematic underestimation of the error when using jackknife. For the same reason, the bootstrap distribution of the phase shift is not symmetric and has long tails. We therefore conclude that bootstrapping provides a more reliable determination of the uncertainty compared to jackknife for phase-shift studies.

The error on the fitted parameters, as well as on the couplings and masses of the charmonium states, is determined as the standard deviation of the bootstrap samples, except for the case when the distribution is asymmetrical or has long tails. In this situation, we quote asymmetric uncertainties as determined by the interval that fits 68 % of the bootstrap samples, i.e. by an interval that excludes 16 % of the bootstrap samples from the left and right tail.

Appendix B Correlation matrices of the fitted parameters

In this appendix we include the correlation matrix $C$ between the parameters $\kappa_{i}$ of the “double pole” fits described in Sec. V.1 which fulfill the consistency checks of Sec. V.3. The correlation matrix is defined as

[TABLE]

where the expectation values appearing in the numerator and the denominator are computed directly from the bootstrapping distributions and $\bar{\kappa}_{i}=\langle\kappa_{i}\rangle$.

For the “double pole” fit form, the six parameters $\kappa_{i}$ are given in the order $m_{2}$, $m_{1}$, $m_{3}$, $G_{2}$, $G_{1}$ and $g_{3}$. The correlation matrices are

[TABLE]

for $\kappa_{c}=0.12522$, and

[TABLE]

for $\kappa_{c}=0.12315$.

Appendix C Alternative parametrization aiming to describe the full vector charmonium spectrum

The “double pole” fit form is able to reproduce the charmonium spectrum including the $\psi(2S)$, as discussed in Sec. V.1. The presence of an additional crossing of the $\Omega$-function below the $\psi(2S)$ is an interesting feature that suggests the possibility of describing the full charmonium spectrum including even the $J/\psi$. For example, by using a slightly modified “double pole” fit form with an additional linear parameter $\rho$

[TABLE]

The corresponding $\Omega$-function determined from the fits are presented in Fig. 17. While the couplings and the masses $G_{1}$, $G_{2}$, $m_{1}$ and $m_{2}$ are consistent within error with Eq. (21) and (22) after the inclusion of the $J/\psi$, the parameter $\rho$ fluctuates strongly across the bootstrapping resamples. Given that the spectrum below the $\psi(2S)$ is well separated from the two resonances we are investigating, we can safely exclude the $J/\psi$ from our main analysis. We emphasize again that our focus is the scattering amplitude in the energy region above the $\psi(2S)$ (the consistency check of Sec. V.3 would not be satisfied for the $J/\psi$ with the parameterization 38).