One-loop analysis of the interactions between doubly charmed baryons and Nambu-Goldstone bosons

The interactions between the spin-$1/2$ doubly charmed baryons and Nambu-Goldstone bosons are analyzed within a manifestly relativistic baryon chiral perturbation theory up to next-to-next-to leading order, by using the so-called extended-on-mass-shell scheme. We utilize heavy diquark-antiquark symmetry to estimate the low-energy constants in the chiral effective Lagrangians. The $S$- and $P$-wave scattering lengths are predicted. We find that those diagrams, vanishing exactly in the heavy-quark limit, do contribute slightly to the $S$-wave scattering lengths in reality. The influence of the spin-$3/2$ doubly charmed baryons, as heavy-quark spin partners of the spin-$1/2$ ones, on the scattering lengths is discussed as well. Finally, $S$-wave phase shifts for elastic scattering processes are presented in the energy region near threshold. Our results in this work will not only be very useful for performing chiral extrapolations of future lattice QCD data, but also provide us chiral inputs for the investigation of the spectroscopy of doubly heavy baryons.


Introduction
Over the past two decades, hadrons containing heavy quarks have been extensively explored, see e.g. Refs. [1][2][3][4][5][6] for recent reviews. This is not only because of the high experimental interest of completing the hadron spectroscopy, but also due to the theoretical importance in understanding the low-energy dynamics of quantum chromodynamics (QCD). Amongst them, the doubly heavy baryons, predicted by conventional quark model [7], have attracted quite a few attentions. They offer a unique platform for investigating the non-perturbative dynamics of light quarks in the environment of two heavy quarks and, consequently, can be used to test the correctness of theoretical frameworks such as the QCD-inspired quark model, the non-relativistic factorization theory [8] and so on.
Experimental attempts have been made to hunt for doubly heavy-flavored baryons since 2002 [9,10], but their existence has been under controversy till 2017. In 2017, the LHCb collaboration announced the first successful observation of such states [11]. It was reported that a doubly charmed baryon Ξ ++ cc with high significance is observed via the decay mode Ξ ++ cc → Λ + c K − π + π + . The existence of the Ξ ++ cc state was later confirmed through another decay mode Ξ ++ cc → Ξ + c π + [12]. It should be noted that searching for doubly charmed baryons in the above two decay modes was previously suggested by the theoretical work [13]. Determinations of the properties of the Ξ ++ cc baryon are conducted subsequently, e.g., measurements of its life time [14], mass [15] and production cross section [16]. On the other hand, the quark content of Ξ ++ cc is assigned to be [ccu] in the quark model. Therefore, its SU(3) partners, Ξ + cc and Ω + cc with contents [ccd] and [ccs] respectively, are expected to exist and are hoped to be observed at LHCb. Thereby, searches for the two states are carried out and are still ongoing. Yet, no evidence of observation has been found [17][18][19].
In fact, the search for Ξ + cc was initially performed by the SELEX collaboration twenty years ago [9]. However, the SELEX results of Ξ + cc reported in Refs. [9,10] are not confirmed by any other collaborations: FOCUS at Tevatron (proton-antiproton) collider [20], BaBar and Belle at electron-positron colliders [21,22], and LHCb at proton-proton collider [17,18]. Even worse, the experimental value of the Ξ + cc mass by SELEX is inconsistent with theoretical determinations obtained by relativistic quark model [23], effective potential models [24], heavy quark effective theory [25], lattice QCD [26][27][28], etc.. See also e.g. Ref. [29] for a brief review. The aforementioned facts lead to a long-standing puzzle in questioning the existence of doubly charmed baryons. The puzzle has been solved partly by the observation of Ξ ++ cc at LHCb [11,12]. And hopefully, it will be addressed completely in the foreseeable future, as the LHC record condition keeps improved so as to overcome the drawbacks caused by the short life times of Ξ + cc and Ω + cc [30]. The discovery of the Ξ ++ cc has stimulated a multitude of works on the theoretical side. Various approaches have been employed to decipher the underlying information of these doubly heavy baryons. Static properties, decays and productions of the doubly charmed baryons have been widely studied via quark model [31][32][33][34], extended chromomagnetic model [35], SU(3) symmetry method [36,37], operator product expansion [38], QCD sum rules [39][40][41][42], lattice QCD [43,44] and so on.
The doubly heavy baryons are composed of two heavy quarks (Q) and one light quark (q), where both heavy quark symmetry as m Q → ∞ and chiral symmetry as m q → 0 will manifest. Therefore, the doubly heavy baryons provide a novel platform to study heavy quark symmetry and chiral symmetry of light quarks simultaneously. In particular, the two heavy quarks in the baryons usually act as a form of compact heavy diquark. In the heavy quark limit (HQL), the heavy diquark belongs to the color3 c representation and serves as a static color source for the light quarks. The same color dynamics arises in the mesons containing a single heavy antiquark. The correspondence between a heavy diquark and a single heavy antiquark is known as heavy diquark-antiquark (HDA) symmetry. In consequence, the doubly heavy baryons can be related to the mesons with a heavy antiquark component through HDA symmetry [45]. In practice, effective field theories (EFTs) are efficient and powerful in the sense that they can easily implement all such kind of symmetries through the approach of constructing pertinent effective Lagrangians. Chiral perturbation theory (ChPT) [46][47][48] is the low-energy EFT of QCD, which is a powerful tool to explore hadron physics in the non-perturbative regime. It has achieved great triumphs in the pure mesonic sector. Meanwhile, various extensions have been developed in the past years [49][50][51]. Note that, ChPT for heavy hadrons is comprehensively reviewed in Ref. [6]. The version extended to the single-baryon sector [52] is known as baryon chiral perturbation theory (BChPT). However, when baryon fields are incorporated as explicit degrees of freedom in the theory, the notable power counting breaking (PCB) problem arises due to the non-zero mass of the baryons in the chiral limit. Many approaches have been proposed in order to remedy this issue. The most popular ones are the heavy-baryon (HB) formalism [53], the infrared regularization (IR) prescription [54] and the extended-on-mass-shell (EOMS) renormalization scheme [55]. Thereinto, the EOMS scheme not only restores the correct power counting but also respects the original analytic structure. Furthermore, it has been demonstrated that better convergence properties can be established compared to other approaches; see Refs. [50,56] for reviews.
Within the framework of BChPT, the properties such as masses and magnetic moments of the doubly charmed baryons have been investigated in Refs. [57][58][59][60][61][62][63][64]. For the scattering of pseudoscalar Nambu-Goldstone bosons (NGB) off doubly charmed baryons, Ref. [65] conducts a lead-order (LO) BChPT calculation. Unitarization of the LO chiral amplitude is carried out so as to search for possible exotic doubly charmed states. The low-lying spectrum of the double-charm baryons with negative parity was further studied in Ref. [66] by using a chiral potential of next-to-leading order (NLO). In Ref. [67], the scattering lengths of NGB and doubly charmed baryons are calculated up to next-to-next-to-leading order (NNLO) in the HB formalism. For a relativistic chiral description of the doubly charmed baryons, complete and minimal set of chiral effective Lagrangians have already been constructed up to O(p 4 ) in Ref. [68], with p collectively denoting the chiral expansion parameters. However, a systematical one-loop analysis of the interactions between the doubly charmed baryons and NGBs in manifestly relativistic BChPT at one-loop order is still lacking.
In present work, we have calculated the scattering amplitudes for the interactions of NGBs and doubly charmed baryons in the covariant BChPT up to NNLO, i.e. the leading one-loop order. The ultraviolet (UV) divergences stemming from the loops are removed by utilizing the modified minimal subtraction scheme, namely, the MS − 1 scheme. Furthermore, we have explicitly checked that the PCB terms can be absorbed exactly via a finite shifts of the low energy constants (LECs) according to the essence of EOMS scheme. In this way, we obtain the EOMS-renormalized chiral amplitude, which is renormalizationscale independent and can be used to derive physical observable of interest. Moreover, the obtained one-loop amplitudes possess proper analytic structure and correct power counting. And hence they are very suited to perform chiral extrapolation, when relevant lattice QCD data are available in future. In addition, our one-loop amplitudes can be applied to evaluate finite-volume corrections of lattice QCD results by merely substituting all the involved one-loop integrals with their finite counterparts that are uniformly formulated in Ref. [69].
For the lack of available data from experiments or lattice QCD at present, we have to estimate the unknown LECs in chiral effective Lagrangians by imposing the HDA symmetry mentioned above. Specifically, the unknown LECs in our case can be related to the ones involved in the Dφ scatterings, where D and φ stand for charmed D mesons and NGBs, respectively. Fortunately, the values Dφ LECs have been well determined by performing fits to lattice QCD data of the S-wave Dφ scatting lengths in Refs. [70][71][72][73]. In fact, two of the O(p 2 ) LECs, b 1 and b 2 , have been determined by fitting to lattice QCD data on the doubly charmed baryon masses in Ref. [63]. We find that the values of b 1 and b 2 obtained in Ref. [63] agree well with those by HDA symmetry, indicating the feasibility of the use of HDA symmetry in this work. The remaining LECs, that can not be constrained by HDA symmetry, are set to zero in line with the ansatz of naturalness of LECs.
Given that the LECs are pinned down, we predict the Sand P -wave scattering lengths for the elastic channels with definite strangeness S and isospin I. For S-wave scattering lengths, our relativistic results turn out to be qualitatively consistent with the ones obtained in the HB approach [67]. Sizable relativistic recoil corrections exist, for instance, in the channel of Ξ ccK → Ξ ccK with (S, I) = (−1, 0). In Ref. [67], the contributions from the diagrams, which should vanish exactly only in the HQL, are ignored. We have explicitly verified that those HQL-vanishing diagrams indeed contribute marginally to the S-wave scattering lengths in the real world with a finite heavy quark mass. In HQL, since the spin-3/2 and spin-1/2 doubly charmed baryons degenerate into a heavy quark spin (HQS) multiplet, they should be treated on equal footing. Therefore, we assess the influence of the spin-3/2 states on the scattering lengths by incorporating them as explicit degrees of freedom into the effective Lagrangians. As expected, they mainly affect the P -wave scattering lengths with quantum numbers J P = 3 2 + . Their contributions to the S-wave scattering lengths are negligible. As byproducts, we also discuss the so-called resonance contribution to the LECs by integrating out the spin-3/2 baryons.
For future reference, S-wave phase shifts are calculated for the channels of elastic scatterings in the energy regions close to the lowest thresholds under consideration. Future lattice QCD simulations of the low-energy interactions between doubly charmed baryons and NGBs are necessary to explore the spectrum of the doubly heavy baryons. Our S-wave phase shifts can be associated directly with the energy levels at non-zero momenta via the famous Lüscher formula [74] and its generalizations [75][76][77][78].
The layout of this manuscript is described as follows. Formal aspects for the scattering amplitude are introduced in section 2, where Lorentz decomposition, the strangeness-isospin structure and partial wave projection are briefly illustrated. Details on the calculation of the chiral amplitudes are exhibited in section 3. Chiral effective lagrangians relevant to our calculation up to NNLO are displayed in subsection 3.2. Tree-level and loop amplitudes are given in subsections 3.3 and 3.4, respectively. Section 4 complies the procedure of renormalization of the one-loop amplitudes within the EOMS scheme. In section 5.1, the values of the LECs are estimated. The S-wave and P -wave scattering lengths are calculated in subsection 5.2. Subsection 5.3 discusses the impact of the spin-3/2 doubly charmed baryons. The S-wave phase shifts for the channels of elastic scatterings are presented in subsection 5.4. Finally, section 6 comprises our summary and outlook. The β-functions concerning the UV-and EOMS-subtractions are relegated to appendices A and B, respectively. Appendix C contains the definition of loop integrals and explicit expressions of loop corrections for the masses and wave function renormalization constants. Application of HDA symmetry to the estimation of the LEC values is detailed in appendix D.
2 Formal aspects of scattering amplitude

Lorentz structure of the amplitude
The scattering process of ψ 1 (p)φ 1 (q) → ψ 2 (p )φ 2 (q ), with momenta indicated in parentheses, is described by the Lorentz-invariant amplitude, which can be decomposed as Here, ψ 1,2 ∈ {Ξ ++ cc , Ξ + cc , Ω + cc } represent the incoming and outgoing doubly charmed baryons respectively, while φ 1,2 ∈ {π ± , π 0 , K ± , K 0 ,K 0 , η} are the incoming and outgoing Goldstone bosons in order. The symbols σ and σ denote the spins of the corresponding baryons. The Mandelstam variables are defined by which satisfy the constraint s + t + u = m 2 The Lorentz decomposition of the scattering amplitude is not unique, an alternative form is given by .
Similar to the situation in pion-nucleon scattering, the decomposition in terms of D and B is more suited to perform chiral expansion [79], while the other is more practically convenient for the extraction of the structure functions A and B.

Amplitudes for given strangeness and isospin
The doubly charmed baryons and the pseudoscalar NGBs fill the representations of the SU(3) flavor group of the light (u, d, s) quarks. Namely, they belong to the SU(3) triplets and octets, respectively. In consequence, there exist a multitude of physical scattering amplitudes corresponding to the various charge states showing up in the multiplets. These scattering amplitudes can be classified, thanks to the conservation of strangeness (S) and isospin (I). As a result, the general meson-baryon scattering processes can categorized into 7 independent channels of interactions: 4 single-channel processes and 3 coupled-channel ones. In practice, the amplitudes with definite (S, I) quantum numbers can be expressed in terms of the physical amplitudes. Explicit relations are listed in the following.

Partial wave projection
The partial wave expansion of the scattering amplitude takes the form [80,81] T (S,I) where σ = (σ 1 , σ 2 , σ 3 ) is a vector with components being Pauli matrices, and χ (χ † ) is the spinor of the incoming (outgoing) baryon. Furthermore,q = q/|q| andq = q /|q | are the directions of the incoming and outgoing Goldstone mesons, respectively. P (z) is the Legendre polynomial and P (z) is its derivative. z = cos θ with θ the scattering angle. Here is the orbital momentum, and the notation ± indicates the total angular momentum is J = ± 1 2 . In Eq. (2.21), the amplitude T (S,I) ± is the so-called dimensionless partial-wave amplitude, which possesses definite quantum numbers of strangeness S, isospin I and total angular momentum J. It is popular to redefine a new partial wave amplitude f which is more preferable in investigations of the analytic properties [82]. The inversion of Eq. (2.21) yields f (S,I) where A (S,I) and B (S,I) are given by In the above equations, E ψ 1 (E ψ 2 ) is the energy of the incoming (outgoing) doubly charmed baryon, E φ 1 (E φ 2 ) denotes the energy of the incoming (outgoing) meson, and q (q ) stands for the three momentum of the incoming (outgoing) meson. In the center-of-mass (CM) frame, one has For an elastic scattering process ψφ → ψφ, the relevant partial wave amplitude f (S,I) ± (s) in Eq. (2.23) can be simplified to [83,84] f (S,I) is the CM energy of doubly charmed baryon and |q| = λ 1/2 (s, m 2 ψ , m 2 φ )/(2 √ s) refers to the modulus of the CM momentum.

Scattering length
Scattering length is one of the most important quantities characterizing the properties of strong interaction. In what follows, we derive the formulae for the calculation of the scattering lengths in the elastic scattering channels.
In the vicinity of threshold, the amplitudes f (2.28) The coefficients on the right hand side of the above equation are referred to as threshold parameters. Specifically, the coefficients of the first three terms, a On the other hand, for small t, A (S,I) (s, t) and B (S,I) (s, t) can be expressed as , a (S,I) , a (S,I) . (2.33) Given that the expressions of the A and B functions are known, the above formulae can be readily applied to obtain the scattering lengths analytically.

Power counting
The amplitudes for the processes of on-shell scatterings are multivariate functions of masses and Mandelstam variables. Since the baryon masses are non-zero in the chiral limit, the chiral expansion of the amplitudes in the vicinity of threshold can be organized in powers of the following quantities, with Λ χ denoting the chiral symmetry breaking scale. Accordingly, the power counting rules for those parameters are set as where p is a collective symbol representing the small parameters. An important feature in BChPT is that each Feynman diagram under consideration is characterized by a chiral dimension D. Namely, the importance of the diagram is regarded to be the order of (p/Λ χ ) D . In our case of one-baryon sector, the chiral dimension D for a given diagram can be determined by the naive power counting rule, with L the number of loops, V n the number of the n th -order vertices, I φ the number of internal pion lines, and I ψ the number of internal doubly charmed baryon lines. However, there exist pieces in the loop amplitudes, originating from the diagrams with internal baryon lines, which violates the above power counting rule (3.3). These pieces are known as PCB terms, as mentioned in the Introduction. The emergence of PCB terms is due to the fact that the masses of doubly charmed baryons do not vanish in the chiral limit, as pointed out by Ref. [52]. We will address this issue by using the EOMS scheme in section 4.

Chiral effective Lagrangian
The chiral effective Lagrangian, which is relevant to our calculation of the meson-baryon scattering amplitude up to O(p 3 ), takes the form with the superscripts '2i' and 'j' in the brackets representing the chiral dimensions. For the purely mesonic sector, we need the following terms [48] where · · · stands for the trace in the flavor space, F is the pion decay constant in the SU(3) chiral limit [51,85], and L k (k = 4, 5, 6, 8) denote the mesonic LECs. The NGB fields are collected in U , which reads Furthermore, in the isospin limit m u = m d =m, the chiral operator χ can be written as where B 0 = − 0|qq|0 /3F 2 , and 0|qq|0 is the quark condensate in the chiral limit [47].
For the baryonic sector, we proceed with the SU(3) triplet ψ in which the doubly charmed baryons are contained. The physical states Ξ ++ cc , Ξ + cc and Ω + cc form the baryon triplet ψ, which reads The full set of the chiral operators up to and including O(p 4 ), describing the interactions between Goldstone bosons and doubly charmed baryons, is constructed in Ref. [68]. In our current case, the required terms are given as follows: where Here, m is the baryon mass in the chiral limit and g denotes the bare axial-vector coupling constant. The baryonic LECs b j (j = 1, · · · , 7) and c k (k = 11, · · · , 20) are unknown parameters in units of GeV −1 and GeV −2 , respectively. H.c. stands for terms obtained by hermitian conjugation.

Tree amplitudes
According to the aforementioned power counting rule, the tree-level Feynman diagrams contributing to the meson-baryon scattering amplitude up to O(p 3 ) are shown in Figure 1. Note that the corresponding crossed diagrams are not displayed in Figure 1. The LO tree-level amplitude takes the following form with the functions F and G defined by where m is the mass of the intermediate doubly charmed baryons in the chiral limit. All the involved coefficients are given in Table 1. The Weinberg-Tomazawa (WT) term, accompanied by the coefficient C WT , stems from diagram (b) in Figure 1. The s-channel diagram (a) gives the terms proportional to C (1) S , while its crossed partner yields the pieces with coefficient C The coefficients C (2) i (i = 1, · · · , 4) are compiled in Table 2. All of these coefficients are obtained from diagram (c) of Figure 1. Table 1. Coefficients appearing in the tree amplitudes of O(p 1 ). The exchanged doubly charmed baryons are indicated in the square brackets. The scattering processes are categorized by strangeness S and isospin I. Table 2. Coefficients in the tree amplitudes of O(p 2 ).
The tree-level amplitude at O(p 3 ) can be written as The coefficients are shown in Tables 3, 4

Leading one-loop contributions
is used for brevity. Abbreviations: is used for brevity.
Abbreviations: amplitudes, we do not need to distinguish the physical masses and the bare masses, since the caused difference is of higher order beyond our accuracy. Note that the contributions of diagrams corresponding to one-loop corrections on the external legs are incorporated via wave function renormalization, which will be discussed in the next section. Let us begin with the baryonic sector. The dressed propagator of the doubly charmed baryons is defined as

Masses and wave function renormalization constants
where m denotes the bare baryon mass and Σ ψ ( / p) refers to the baryon self-energy up to leading one-loop order. The sum of one-particle-irreducible diagrams contributing to the two-point function is denoted by −iΣ ψ ( / p), which comprises contact and one-loop diagrams shown in Figure 3. Namely, one has where ζ ψ are unit vectors in the SU(3) flavor spaces, Furthermore, the chiral results for the self energies in Eq. (4.2) are given by The pole position of the dressed propagator (4.1) defines the physical mass m ψ of the baryon. That is, where δm loop Ξcc and δm loop Ωcc are given in Eq. (C.5). The wave function renormalization constant is defined as the residue of the pole term of the dressed propagator, where m ψ is the physical baryon mass as specified in Eq. (4.6). In view of Eq. (4.1), the wave function renormalization constant Z ψ is given by where a prime means performing derivative with respect to / p. Explicit expression of Z ψ can be readily obtained by substituting Eq. (4.4) into the above equation. The readers are referred to Appendix C for the final results for Z Ξcc and Z Ωcc .
Likewise, one can derive the wave function renormalization constants for the Goldstone bosons. Nevertheless, they have been extensively calculated elsewhere [48,51,86]. For completeness, we quote the results in the following

The full scattering amplitude within EOMS scheme
The Lehmann-Symanzik-Zimmermann (LSZ) reduction formula indicates that the full scattering amplitude, i.e. the on-shell transition amplitude, is related to the amputated Green function T in the momentum space through where T ψ 1 φ 1 →ψ 2 φ 2 has been calculated in Section 3 and the wave function renormalization constants of the involved fields are presented in Subsection 4.1. Within our working accuracy, the full amplitude should be truncated at O(p 3 ). Therefore, the chiral expansion of the full amplitude can be written as tree + T tree + T tree + T loop + T wf , (4.15) where the numbers in the superscripts denote the chiral orders. The last term is counted as O(p 3 ) and takes the form tree , (4.16) . It is worth mentioning that this term perturbatively incorporates the effect of the multiplication of the wave function renormalization constants in Eq. (4.14), and the above procedure is usually called wave function renormalization. In another word, it actually takes into account the contribution of the diagrams with one-loop corrections on the external legs, as displayed in Figure 4. The full scattering amplitude in Eq. (4.15) has UV divergencies and PCB terms. Here we adopt the EOMS scheme [55] to address these issues. The EOMS scheme contains two steps: MS-1 subtraction and an extra finite renormalization.
In MS-1 scheme, the bare baryon mass and axial-vector coupling constant, LECs are divided into renormalized and divergent parts. The divergent ones are used to cancel the UV divergence from loop amplitudes. To be specific, these quantities are written in the following form where µ is the renormalization scale, and Then we utilize a finite renormalization to restore the correct power counting. Since the PCB terms are polynomials of chiral expansion parameters, they can be properly absorbed by the LECs in the tree amplitudes. The EOMS-renormalized parameters are defined by 1, 2, . . . , 7) . (4.18) Here, the EOMS β-functions are collected in Appendix B. Note that the O(p 3 ) LECs c r j (µ) maintain unchanged, i.e. c j = c r j . The renormalized amplitudes obtained in EOMS scheme possess of original analytic structure and respect the power counting rule.
Finally, it should be mentioned that, in practical computations, the chiral limit decay constant F is always replaced by the physical decay constants, F π , F K and F η . For the amplitude of the process ψ 1 φ 1 → ψ 2 φ 2 , this is achieved by making the following substitution It should be noted that such a substitution, when carried out for the O(p 1 ) tree amplitude, generates O(p 3 ) pieces that should be kept in our calculation. For the O(p 2 ) and O(p 3 ) amplitudes, one can merely make the replacement F 2n → (F φ 1 F φ 2 ) n , since the caused differences are of higher orders beyond the accuracy of our calculation. A merit of Eq. (4.19) is that the physical decay constants are properly chosen according to the incoming and outgoing Goldstone bosons. For instance, the F 2 in the amplitude of Ξ cc π → Ξ cc π is changed to F 2 π rather than e.g. F 2 K , while the one in the inelastic scattering amplitude of Ξ cc π → Ξ cc η is substituted by F π F η . NLO expressions of the decay constants in ChPT can be found in Ref. [48], which read with the chiral corrections being The bare parameters L i 's are related to the corresponding UV-renormalized ones by (4.24) The Γ i functions are given in Ref. [48] and are shown in Eq. (A.1) for easy reference.

Parameters
Masses and decay constants used in our work are collected in Table 6. Since the Ξ + cc and Ξ ++ cc are treated as members of an isospin doublet, their masses are degenerate in the isospin limit and we take m Ξcc = 3621.55 MeV from the PDG review [87]. Unlike the Ξ cc states, there is no experimental value for the Ω cc baryon so far, therefore we employ m Ωcc = 3738 MeV determined by lattice QCD [28].
For the decay constants of pion and kaon, we adopt their world-average values from Ref. [87]. It is known that the physical η and η meson are superpositions of the singlet η 0 and octet η 8 , that can be systematically formulated in large-N c ChPT, see e.g. Refs. [88,89]. Nevertheless, in the standard SU(3) ChPT we are using here, the singlet η 0 is absent and the octet member η 8 is approximately regarded as the physical η meson. Based on the SU(3) ChPT calculation performed in Ref. [86], one has F η 1.27F π 117 MeV. This value is more or less in agreement with the recent lattice QCD determination [90], if one assumes F η F 8 η with F 8 = 115.0 MeV. Due to the lack of relevant experimental and lattice QCD data, a big challenge we encounter is the determination of the unknown LECs in the Lagrangians (3.10), (3.11) and (3.12). Thanks to the HDA symmetry [45], the doubly charmed baryon sector can be connected with the ones in the charmed meson sector. Derivation of the relations is detailed in Appendix D.
Specifically, for the LO coupling constant, one has  Their values have been determined by fitting to lattice QCD data in Ref. [72]. Therein, four different fits were performed using the method of unitarized ChPT (UChPT). The fit denoted by UChPT-6(b ) was done by enforcing the naturalness requirements of the fitting parameters. Moreover, the validity range of UChPT is ensured by further excluding the lattice data with a pion mass larger than 600 MeV. As it was expected, the obtained Dφ LECs from the UChPT-6(b ) fit are more natural than the others in Ref. [72]. Therefore we utilize the UChPT-6(b ) outputs, collected in the last column of Table 7 for completeness, to estimate the values of b i (i = 1, · · · , 6) andc i (i = 11, 12, 20). The obtained results are shown in the 2nd column of Table 7. Importantly, those LECs determined by HDA symmetry are acceptable in the sense that they turn out be of natural size. Furthermore, the corresponding uncertainties are obtained by Monte Carlo propagating the errors of the Dφ LECs of Ref. [72]. 4 Table 7. It is worth noting thath 1 in Table 7 was fixed by the mass difference between D and D s mesons in Ref. [72]. Likewise, in our case the LECb 2 can be estimated by the mass difference of the Ξ cc and Ω cc baryons. With the help of Eq. (4.7) and Eq.  Table 6 have been used. One can see that theb 2 value in Eq. (5.4) is comparable with the one in Table 7, justifying the validity of HDA symmetry to a certain extent.

Mention that the loop corrections have been neglected and the masses in
In fact, a more reliable determination ofb 1 andb 2 has been conducted in Ref. [63]. The two LECs are pinned down by fitting to the lattice QCD data of the baryon masses, leading to 5b 1 = −0.09 ± 0.08 GeV −1 ,b 2 = −0.09 ± 0.09 GeV −1 .
(5.5) Theb 1 -b 2 parameter space allowed by 1-σ uncertainties cover the values ofb 1 andb 2 in Table 7 and Eq. (5.4). Therefore, we use Eq. (5.5) forb 1 andb 2 throughout this paper. It should be also stressed that the numbers in Eq. (5.5) are obtained at the renormalization scale µ = 1 GeV. For consistency, the same renormalization scale is chosen, during our numerical computation of one-loop corrections. Apart from the fixed parameters discussed above, there are still eight unknown LECs. Since they can not be estimated via HDA symmetry, we assume them to be zero,b 7 = 0.0 GeV −1 andc k = 0.0 GeV −2 , k ∈ {13, · · · , 19}. Such an assumption is more or less reasonable in view of the smallness for most of the LECs in Table 7. A solid determination of these parameters is expected to be done only when lattice QCD data of e.g. scattering lengths are available in the future.
Finally, the values of the LECs L r 4 and L r 5 can be found in Ref [48], L r 4 = −0.3 × 10 −3 and L r 5 = 1.4 × 10 −3 , which are obtained at the scale µ = M ρ . An recent update of the mesonic LECs is summarized in Ref. [93], and comparable results of L r 4 and L r 5 are achieved. Those values at µ = M ρ can be translated to the ones at µ = 1 GeV with the help of the following relations

Prediction of scattering lengths
Once all the involved LECs are pinned down, we are now in the position to calculate the Sand P -wave scattering lengths numerically. By definition, scattering lengths can be calculated via However, the fraction on right hand side can not be computed numerically exact at threshold for ≥ 1. That is, the fraction value becomes undefined when one has zeros in the denominator, as the modulus of the three momentum is vanishing. To avoid this issue, by making use of Eq. (2.33), we have derived analytical expressions for the Sand P -wave scattering lengths, with which numerical computations can be carried out precisely. For calculations with accuracy up to NNLO, derivatives of one-loop integrals are involved, which are handled by adopting the techniques proposed in Ref. [94].  The convergence of SU(3) ChPT has remained under debate for over several decades, see e.g. Refs. [56]. Generally speaking, the convergence in the three-flavor ChPT calculations is usually worse than the two-flavor case, due to the relatively large strange quark mass. 6 Here it would be also interesting to have a look at the convergence properties of the chiral expansion of the ψφ scattering lengths. We first discuss the processes involving pion mesons. It can be seen from Table 8 that the chiral series for the elastic Ξ cc π scatterings, with (S, I) = (0, 3/2) or (0, 1/2), converge well if one only concerns the first two orders. Namely, the O(p 1 ) contributions are significantly larger than the O(p 2 ) ones. Nevertheless, although the sums of O(p 3 ) trees and loops in the two channels are small, there are comparable to the O(p 2 ) trees. We ascribe the failure of convergence to the underestimation of the O(p 3 ) trees for poor information on the LECs. As for Ω cc π → Ω cc π, scattering length starts to 6 It should be pointed out that a one-loop calculation is not really sufficient to make a solid statement about convergence. It is well-known from both mesonic and baryonic sectors of ChPT that one needs at least two loops to make a significant statement about the chiral expansion. For example, an estimate for the convergence radius of the chiral expansion of the nucleon mass at the two-loop level was performed in Refs. [95,96], indicating a breakdown of convergence already below mπ ≈ 360 MeV. The convergence might be even worse in the SU(3) case, since the kaons and the eta are substantially heavier. contribute at O(p 2 ), since the LO term for this channel (see Eq. (3.14) and Table 1) is identical to zero exactly.
On the other hand, the scatterings of kaon and eta mesons off doubly charmed baryons have bad convergence properties, due to the emergence of large masses of kaon and eta as chiral expansion parameters. Such a non-convergent behaviour usually implies resummation of diagrams is required so that higher-order contributions can be implemented nonperturbatively. Similar situations happened for SU(3) meson-meson scatterings [86], and kaon-nucleon scatterings [97,98]. In general, the resummation procedure restores unitarity and extends the applicability range of ChPT, see Refs. [29] for a recent review. It should be emphasized that, the so-called UChPT amplitude plays a crucial role in investigating the spectrum of doubly charmed baryons. For instance, possible candidates of negative-parity doubly charmed baryons have been found on the basis of the LO [65] and NLO [66] ChPT amplitudes. It is also worth noting that the O(p 3 ) amplitudes turn out to be significant and an inclusion of their effects may improve the studies made in Refs. [65,66]. Especially for the two channels, Ξ cc K → Ξ cc K with (S, I) = (1, 1) and Ω ccK → Ω ccK with (S, I) = (−2, 1/2), since there exist large cancellations between the LO and NLO contributions, the NNLO corrections become dominating. In a word, the ChPT amplitudes up to NNLO, we obtain in this work, can be applied to systematically scrutinize the doublycharmed-baryon spectroscopy in the future, once more experimental or lattice QCD data are available.
For comparison, the HB results of Ref. [67] are shown in the last column of Table 8. For all the channels, our relativistic results are qualitatively consistent with the ones obtained in the HB formalism. Relativistic corrections are mainly responsible for the differences. Another source might be owing to the fact that the HQL-vanishing diagrams are not taken into account in the calculation of Ref. [67]. Nevertheless, their contributions are negligible, which will be illustrated below. Lastly, the spin-3/2 doubly charmed baryons are incorporated as explicit degrees of freedom in Ref. [67], which would lead to discrepancies as well. As we can see, good agreements are observed within one-sigma uncertainties for the five channels: Ω ccK with (S, I) = (−2, 1/2), Ξ cc K with (S, I) = (1, 0), Ω cc η with (S, I) = (−1, 0), Ω cc π with (S, I) = (−1, 1) and Ξ cc π with (S, I) = (0, 1/2). Such an observation implies that the net effects of relativistic corrections, contributions of HQL-vanishing diagrams and resonance-exchanging diagrams are slight for those channels in S wave.
P -wave scattering lengths with J P = 3 2 + and J P = 1 2 + are complied in Table 9 and Table 10, respectively. It is found that, at O(p 1 ), the P -wave scattering lengths are entirely saturated by the crossing partner of diagram (a). For some channels like Ξ cc π scatterings, contributions from O(p 3 ) loops to the P -wave scattering lengths turn out to be sizeable. Before ending this subsection, we intend to discuss the contributions of HQL-surviving diagrams in more detail. In the HQL, most of the diagrams in Figures 1, 2 and 4 do not contribute at threshold. This can be illustrated by performing an expansion in terms of the inverse of the baryon mass, i.e. HB projection. For the axial term in Eq. (3.10), the HB projection yieldsψ g 2 / uγ 5 ψ −→N v gS ν · u N v + α 1 m −1 + α 2 m −2 + · · · (5.8) Table 9. P -wave scattering lengths a (S,I) 1+ (J P = 3 2 + ) in units of 10 −2 fm 3 .

(S, I)
Processes  Table 10. P -wave scattering lengths a

(S, I)
Processes where the spin matrix is S µ v ≡ i 2 γ 5 σ µν v ν with the four vector v µ = (1, 0, 0, 0). In the HQL, m → ∞, all the inverse mass terms approach to zero. Only the first term that is independent of mass survives. However, it vanishes at threshold due to the feature of derivative coupling of the NGBs. Specifically, one has which corresponds to S v · q φ in the momentum space, with q φ being the momentum of Goldstone boson. At threshold, q φ = (m φ , 0), leading to S v · q φ = 0. Subsequently, all the diagrams containing axial vertices disappear in the HQL. Therefore, only the diagrams that are irrelevant to the axial coupling g survive. For easy reference, they have been labelled by boxed "HQL" in Figures 1, 2 and 4. Results of scattering lengths obtained by including only the HQL-surviving diagrams are shown in Table 11. As expected, for the S-wave scattering lengths the deviation of the HQL results in Table 11 from the ones in Table 8 is negligible. However, the results of Pwave scattering lengths change dramatically by switching off the HQL-vanishing diagrams, especially for a 1+ with J P = 3 2 + . One may probably ascribe such large variations to the absence of the HQS partners like spin-3/2 doubly charmed baryons. They are expected to contribute equally as the spin-3/2 baryons, since the HQS symmetry is exact in HQL. In the next subsection, we will evaluate the effect of the spin-3/2 doubly charmed baryons. Table 11. Results of scattering lengths obtained by taking only the HQL-surviving diagrams into consideration. The Sand P -wave scattering lengths are in units of fm and 10 −2 fm 3 , respectively.
The inclusion of spin-3/2 baryons in covariant BChPT is complicated. An appropriate power counting rule should be assigned to the new parameter, i.e. the mass difference between the spin-1/2 and spin-3/2 baryons, denoted by ∆ hereafter. Like the treatment of the ∆(1232) resonances in the traditional BChPT, here we adopt the so-called δ-counting rule that was proposed in Ref. [102]. Specifically, in the energy region near threshold the mass splitting ∆ is counted as O(p 1/2 ) in addition to Eq. (3.2). According to this power counting rule, the spin-3/2 baryon propagator is O(p −1/2 ). Loops with internal spin-3/2 baryon lines are at least of order O(p 7/2 ), and therefore their contributions are beyond the accuracy of our calculation. For the calculation up to O(p 3 ), we only need to take into account the LO and NLO tree-level contributions from explicit spin-3/2 baryons, which are O(p 3/2 ) and O(p 5/2 ), respectively. In order to avoid the emergence of too many new unknown LECs, the NLO Born-term contribution is omitted here. In fact, as pointed out in Ref. [103], the NLO Born term is redundant in the sense that its contribution can be absorbed by redefining the coupling h A in the LO Born term (see Eq. (5.12)) and the LECs in the contact terms (3.16) and (3.18).
For the purpose of our calculation the following LO Lagrangian are needed, where M and h A are the mass of the spin-3/2 baryons and the coupling constant of the ψψ φ interactions in the chiral limit, respectively. A is an arbitrary real parameter but A = −1/2. The spin-3/2 doubly charmed baryons are collected in the triplet Feynman diagrams relevant to our calculation are displayed in Figure 5. Tree-level amplitudes are derived, which can be written in the form as where the coefficients C S and C (1) U are the same as those in Eq. (3.14), which are shown in Table 1. Nevertheless, the exchanged spin-1/2 baryons indicated in Table 1 should be substituted by their spin-3/2 partners. Explicit expressions of the functions F and G are given by As argued in Ref. [101], physical quantities are independent of the parameter A. Therefore, we have set A = −1 for convenience as done in Ref. [103]. Table 12. LO contributions of the spin-3/2 ψ baryons to the scattering lengths. ∆ denotes the mass difference between the spin-3/2 and spin-1/2 baryons. The Sand P -wave scattering lengths are in units of fm and 10 −2 fm 3 , respectively. There are two unknown parameters: h A and M . We need to assign appropriate values to them, so that scattering lengths can be evaluated numerically. Thanks to HQS, the LO coupling constant h A can be related to the g in the Lagrangian (3.10) by h A = 2 √ 3g. The chiral limit mass M is replaced by the physical mass M ψ . The physical mass of spin-3/2 baryon is assumed to be M ψ = m ψ + ∆, where ∆ is the mass splitting. Since the spin-1/2 and spin-3/2 baryons differ only in the relative orientation of the quark spins, the resultant difference in their masses is attributed to the spin-spin interaction in the viewpoint of potential model. Here, following Ref. [67], we roughly take the difference to be around 100 MeV. Three values of ∆ are adopted: ∆ = 50 MeV, ∆ = 100 MeV and ∆ = 150 MeV. Contributions of spin-3/2 baryons to the Sand P -wave scattering lengths are compiled in Table 12.
It can be seen from Table 12 that for S wave the size of the contribution of spin-3/2 baryons is nearly negligible. The sand u-channel exchanges of ψ , corresponding to diagram (a) and (b) in Figure 5 respectively, mainly affect the P -wave scattering lengths. For a given strong interaction, strangeness and isospin are conserved definitely. Therefore, the s-channel exchange of spin-3/2 Ξ cc and Ω cc states contribute only to the scattering processes with (S, I) = (0, 1/2) and (−1, 0), respectively. This can also be justified by seeing the coefficients C  Table 1. Furthermore, we have checked that it dominates the a 1+ scattering lengths in the two coupled channels with (S, I) = (0, 1/2) and (−1, 0), which is in accordance with the conservation law of total angular momentum. The remaining processes get contributions entirely from the crossed diagram, i.e. diagram (b) in Figure 5.
In parallel, one may also study the resonance contribution to the LECs with Eq. (5.12). In the viewpoint of effective field theory, a low-energy effective Lagrangian contains only low-lying degrees of freedom, while resonances at hard scale have been integrated out. The information of resonance contributions is regarded to be encoded in the LECs, i.e, the coefficients of the operators in the chiral effective Lagrangian without resonances. The pertinent contributions of resonances to the LECs can be efficiently achieved by a matching procedure carried out at the level of effective Lagrangian. Such a procedure has already been used, e.g., in the analyses of the LECs in pion-nucleon scattering [49,104] and Dφ interactions [105]. In our current situation, we intend to utilize the above-mentioned technique to evaluate the influence of the ψ states on the LECs. Specifically, the ψexchange amplitudes in Eq. (5.12) are expanded in terms of ν = (s − u)/(4m), t and m φ , and then compared to the contact-term contribution in Eqs. (3.16) and (3.18). As a result, the ψ -exchange contributions to the O(p 2 ) LECs read b ψ 1,2,4, where m is the spin-1/2 baryon mass in the chiral limit. For the O(p 3 ) LECs, the ψ states contribute as For the sake of numerical estimation, we take the chiral limit mass m equal to the average of the physical masses of the Ξ cc and Ω cc , i.e., m = (m Ξcc + m Ωcc )/2 = 3679.8 MeV. Numerical results, obtained with three different mass splitting values ∆ = (50, 100, 150) MeV, are displayed in Table 13.
In Table 13, the magnitudes of the LECs become smaller as the mass splitting ∆ gets larger, which can be inferred from Eqs. (5.14), (5.15) and (5.16). Consequently, the spin-3/2 baryon contribution to the scattering lengths is expected to decrease as ∆ increases.  In Table 12, such a behavior can be clearly observed e.g. for the P -wave scattering lengths a 1− . However, anomaly happens for the scattering lengths a 1+ with (S, I) = (−1, 0) and (S, I) = (0, 1/2). This anomaly actually indicates that the spin-3/2 baryon fields can not be integrated out in the two channels. In other words, explicit inclusion of spin-3/2 doubly charmed baryons is necessary if one intends to well determine the P -wave scattering length with J P = (3/2) + and (S, I) ∈ {(−1, 0), (0, 1/2)}.

S-wave Phase shifts
With the chiral amplitudes, one can compute partial-wave phase shifts straightforwardly, which are functions of the CM energy √ s. Although it is undoable to extract phase shifts experimentally, they can be related to energy levels according to Lüscher formula [74] and its extensions [75][76][77][78], which can be computed by lattice simulations in future.
Usually, the partial-wave amplitude f (S,I) ± (s) from BChPT do not obey partial wave unitarity exactly, since they are derived perturbatively up to a certain order. The method of extracting phase shifts from perturbative amplitudes has been discussed, e.g, in Ref. [103]. Namely, one can calculate the phase shifts in the elastic scattering region by using where |q| is the modulus of the CM momentum. In the present work, we are only interested in the S-wave interactions, whose strength is expected to be stronger than that of higher partial waves. Moreover, one is allowed to ignore the effects of the spin-3/2 HQS cousins of the spin-1/2 doubly charmed baryons, since their impact on the S-wave phase shifts are negligible, as discussed in the preceding subsection. In Figure 6, the S-wave phase shifts for the processes of elastic scattering are plotted for the energy region from threshold  contributions are represented by violet dashed, orange dash-dot-dotted, green dash-dotted lines, in order. It can be seen that the convergence properties of the chiral expansion in the two elastic channels, Ξ cc π with (S, I) = (0, 1/2) and (0, 3/2) are relatively better, compared to the other channels. In fact, most of the LECs we used here are estimated by imposing HDA symmetry, while some of the LECs in the O(p 3 ) tree amplitudes are simply assumed to be zero under the requirement of naturalness. Therefore, a solid conclusion on the convergence properties can be drawn only when relevant lattice QCD data are available. For comparison, the HQL results of phase shifts are shown as well, which are obtained by setting g = 0. They are represented by the black dotted lines in the figure. It is found that the HQL-vanishing diagrams contribute slightly to the phase shifts, which is similar to the case of S-wave scattering lengths.

Summary and outlook
In this work, we have performed a NNLO calculation of the scattering amplitudes for the interactions between NGBs and doubly charmed baryons within the framework of relativistic BChPT. The EOMS scheme is employed to handle with the UV divergences and the PCB terms originating from the loops. We find that most of the unknown LECs in the chiral effective Lagrangians can be pinned down by making use of the HDA symmetry. The obtained LEC values enable us to make predictions of the Sand P -wave scattering lengths. We show that the HQL-surviving diagrams dominate the contribution to the Swave scattering lengths, while the HQL-vanishing ones contribute marginally. The influence of spin-3/2 double-charm baryons on the scattering lengths is also evaluated in detail. Their contributions to the LECs are estimated as well. For future reference, the energy-dependent S-wave phase shifts are plotted for the elastic scattering channels in the energy regions near the lowest thresholds. Our chiral results can be applied to perform chiral extrapolations of future lattice QCD data, and can also be used to investigate the spectroscopy of doubly heavy baryons systematically.

A Γ-functions and β-functions
In this appendix, the Γ-functions in Eq.
Our aim is to derive relations between the ψφ LECs in the Lagrangians (3.10), (3.11) and (3.12)