Hidden-charm and hidden-bottom molecular pentaquarks in chiral perturbation theory

The newly observed $P_c(4312)$, $P_c(4440)$ and $P_c(4457)$ at the LHCb experiment are very close to the $\Sigma_c\bar{D}$ and $\Sigma_c\bar{D}^\ast$ thresholds. In this work, we perform a systematic study and give a complete picture on the interactions between $\Sigma_c^{(\ast)}$ and $\bar{D}^{(\ast)}$ systems in the framework of heavy hadron chiral perturbation theory, where the short-range contact interaction, long-range one-pion-exchange contribution, and intermediate-range two-pion-exchange loop diagrams are all considered. We first investigate the three $P_c$ states without and with considering the $\Lambda_c$ contribution in the loop diagrams. It is difficult to simultaneously reproduce the three $P_c$s unless the effect of $\Lambda_c$ is included. The coupling between $\Sigma_c^{(\ast)}\bar{D}^{(\ast)}$ and $\Lambda_c\bar{D}^{(\ast)}$ channels is crucial for the formation of these $P_c$s. Our calculation supports the $P_c(4312)$, $P_c(4440)$ and $P_c(4457)$ to be the $S$-wave hidden-charm $[\Sigma_c\bar{D}]_{J=1/2}^{I=1/2}$, $[\Sigma_c\bar{D}^\ast]_{J=1/2}^{I=1/2}$ and $[\Sigma_c\bar{D}^\ast]_{J=3/2}^{I=1/2}$ molecular pentaquarks, respectively. Our calculation disfavors the spin assignment $J^P=\frac{1}{2}^-$ for $P_c(4457)$ and $J^P=\frac{3}{2}^-$ for $P_c(4440)$, because the excessively enhanced spin-spin interaction is unreasonable in the present case. We obtain the complete mass spectra of the $[\Sigma_c^{(\ast)}\bar{D}^{(\ast)}]_J$ systems with the fixed low energy constants. Our result indicates the existence of $[\Sigma_c^{\ast}\bar{D}^{\ast}]_J~(J=\frac{1}{2},\frac{3}{2},\frac{5}{2})$ hadronic molecules. The previously reported $P_c(4380)$ might be a deeper bound one. Additionally, we also study the hidden-bottom $\Sigma_b^{(\ast)}B^{(\ast)}$ systems, and predict seven bound molecular states, which could serve as a guidance for future experiments....


Introduction
The charmonium physics is one of the most charming and interesting sectors in quantum chromodynamics (QCD). On the one hand, the charmonium spectra deepen our understanding on the nonperturbative QCD and serve as a good platform to develop multifarious potential models. On the other hand, the discoveries of the exotic XY Z states challenge the conventional hadron spectra [1], since these states cannot be easily reconciled with the predictions of the conventional quark models. Furthermore, the heavy quark symmetry in the charm sector is not good enough, thus the heavy quark symmetry breaking effect would manifest itself and lead to some novel phenomena sometimes.
In 2015, two pentaquark candidates P c (4380) and P c (4450) were observed by the LHCb Collaboration in the J/ψp invariant mass spectrum via the weak decay process of Λ 0 b → J/ψpK − [2]. The discovery of these two exotica triggered many discussions on the their internal structures (for some related reviews, see refs. [1,[3][4][5][6][7]), among which, the molecular interpretation is the most favored one. In ref. [8], these two states are interpreted as the deeply bound Σ cD * and Σ * cD * molecular states in the framework of one-pion-exchange model. Whereas in ref. [9], they are regarded as the Σ * cD and Σ * heavy hadron reduction formalism to integrate out the large mass scale [40][41][42]. For the loop diagrams generated by the two-pion-exchange interactions, we will encounter another trouble, which also destroys the power counting rule. Considering the one-loop Feynman diagrams illustrated in figure 1, the scattering amplitude at the leading order of the nonrelativistic expansion is badly divergent because of the pinch singularity [32,39]. Although the problem of divergence can be solved by including the kinetic energies of Σ ( * ) c andD * at the leading order (see some more detailed discussions in refs. [32,37,39]), the amplitude would be finally enhanced by a large factor M/|p| (M could be the mass of Σ ( * ) c orD * ), which will destroy the power counting as well. This strong enhancement is the manifestation of the nonperturbative nature of the nuclear force, which is responsible for the existence of the bound pentaquark states. In other words, a nonperturbative treatment is required.
In the two seminal works [38,39], Weinberg pointed out that we shall focus on the effective potential, i.e, the contributions from two-particle-irreducible (2PI) graphs. The two-particle-reducible (2PR) part, that originates from the on-shell intermediate Σ ( * ) c andD * , should be subtracted. On the other hand, the 2PR part can be automatically recovered when the one-pion-exchange potential is inserted into the nonperturbative iterative equation, such as the Schrödinger equation or Lippmann-Schwinger equation. Therefore, the 2PI parts in the diagrams of figure 1 that contribute to the effective potentials can still be calculated perturbatively. We just need to solve a nonperturbative iterative equation with the obtained effective potential eventually.
For the Σ ( * ) cD ( * ) systems, the mass splittings in the spin doublets (Σ c , Σ * c ) and (D,D * ) do not vanish in the chiral limit, which only vanish in the strict heavy quark limit. Therefore, except for the two particular diagrams in figure 1, the intermediate states in the loops can also be their spin partners. In this case, the loop integral is well defined, and we do not need to make the 2PR subtraction, unless the inelastic one-pion-exchange couple channel is included.
In this work, we try to reproduce the newly observed P c (4312), P c (4440) and P c (4457) after simultaneously considering the leading order contact interaction and one-pion-exchange contribution, as well as the next-to-leading order two-pion-exchange diagrams. The mass splittings are kept in the loop diagrams. If these P c states are shallow bound hadronic molecules, they would be very sensitive to the subtle changes of the effective potentials. Furthermore, the nonanalytic structures, such as the terms with the logarithmic and square root functions, would emerge from the loop diagrams, which may enhance the two-pion-exchange potential to some extent. In particular, the mass difference δ betweenD * andD is larger than the pion mass m π . The heavy quark spin symmetry breaking effect has been noticed for the charmed sectors in some works [11,43]. Besides, one shall not neglect the role of Λ c , since the Λ c π couples strongly with the Σ ( * ) c . Therefore, we also include the contribution of Λ c in the loop diagrams. We will see the dramatic influences of Λ c on the Σ and with the Λ c , and an investigation on interchanging the spins of P c (4440) and P c (4457). In section 5, we study the hidden-bottom systems and predict their mass spectra. In section 6, we give a detailed examination of the heavy quark symmetry breaking effect in the hidden-charm and -bottom systems. In section 7, we conclude this work with a short summary. In the appendices A, B, C and D, we display the definitions and expressions of the loop integrals, the detailed elucidation on how to remove the 2PR contributions with the mass splittings being kept, the derivation of the spin-spin terms in the potentials, and a tentative parameterization of the effective potentials from the quark model, respectively.
We then adopt the heavy baryon reduction formalism [46] to get rid of the large baryon masses in eq. (2.3), where the heavy baryon field is decomposed into the light and heavy components by the projection operators (1 ± / v)/2, where ψ i denotes the relativistic heavy baryon field ψ 1 , ψ 3 and ψ 3 * , M i is their masses, and v µ = (1, 0) represents the four-velocity of a slowly moving heavy baryon. B i and H i are the corresponding light and heavy components, respectively. H i disappears at the leading order expansion.
Consequently, the eq. (2.3) can then be reexpressed with the nonrelativistic form as where S µ = i 2 γ 5 σ µν v ν denotes the spin operator for the spin-1 2 particle. We adopt the mass MeV [47].
Recall that the (ψ 3 , ψ 3 * ) form the spin doublet in the heavy quark limit. Thus eq. (2.7) can be rewritten as a compact form by introducing the super-field [48,49], where the super-fields ψ µ andψ µ are defined as [42,50] Expanding eq. (2.8) and comparing them with the terms in eq. (2.7), one can get the relations among the different coupling constants, The values of g 2 and g 4 can be calculated with the partial decay widths of Σ c → Λ c π and Σ * c → Λ c π [47], respectively. The other axial couplings g 1 , g 3 and g 5 can be obtained by their relations with g 2 in the framework of the quark model [51][52][53], which yields The leading order chiral Lagrangians for the interactions between the anticharmed mesons and light pseudoscalars read [54,55] where X stands for the trace of X in spinor space. The covariant derivative D µ = ∂ µ + Γ µ , δ b = mD * − mD is the mass splitting betweenD * andD. g = −0.59 represents the axial coupling constant, and its value is extracted from the partial decay width of D * + → D 0 π + [47], while the sign is determined by the quark model. We use theH to denote the super-field for the anticharmed mesons, which reads

Contact interactions
We then construct the leading order Lagrangians that account for the interactions between Σ ( * ) c andD ( * ) at the short range. We also use the super-field representations for Σ ( * ) c andD ( * ) to reduce the numbers of the LECs, which read [11] where the D a , D b , E a and E b are four independent LECs. The contact terms contain the residual contributions from the heavy degrees of freedom, which are integrated out and invisible at the low energy scale. Their values can be delicately determined from the experimental data [11] or roughly estimated with the theoretical models [35,37]. D a and D b contribute to the central potential and spin-spin interaction, respectively. E a and E b are related with the isospin-isospin interaction and contribute to the central and spin-spin interaction in spin space, respectively . At the next-to-leading order, we need the O(ε 2 ) LECs to absorb the divergences of the loop diagrams. These O(ε 2 ) contact Lagrangians shall be proportional to the m 2 π , q 2 , δ 2 a and δ 2 b . As demonstrated in ref. [36], there exist a large number of contact terms at O(ε 2 ). It is very difficult to fix all these LECs at present. Therefore, in our work, we try to combine some contributions of the O(ε 2 ) LECs with the leading ones by fitting the experimental data (At least the ones that proportional to m 2 π , δ 2 a and δ 2 b can be absorbed by renormalizing the O(ε 0 ) LECs. The ones correlated with q 2 can be largely compensated by the cutoff).
3 Analytical expressions for the effective potentials of the Σ The effective potential in momentum space can be obtained from the scattering amplitude in the following way [26], where the M 1,2 and M 3,4 are the masses of the initial and final particles. The scattering amplitude M(q) is calculated by expanding the Lagrangians in eqs. (2.7), (2.12) and (2.14). Recall that at the leading order of the nonrelativistic expansions, there are the relations [55] where ψ(p) andH(p) are the relativistic fields. χ(v) is the two-component spinor. TheH(v) is the field in eqs. (2.12) and (2.14). We then make the Fourier transformation on V(q) to get the potential V(r) in the coordinate space, where the Gauss regulator F(q) = exp −q 2n /Λ 2n is introduced to regularize V(q) nonperturbatively [56,57]. As in refs. [35][36][37]58], we set n = 2. The cutoff value Λ is commonly chosen to be smaller than ρ meson mass in the N -N system [58]. We adopt a moderate value Λ = 0.5 GeV as in ref. [11].

Σ cD system
Since theDDπ vertex is forbidden by the parity conservation law, the leading order effective potential for the Σ cD system only arises from the contact terms [diagram (X 1.1 ) in the figure 2]. One can readily get where I 1 and I 2 represent the isospin operators of the Σ Figure 2. The leading order Feynman diagrams that account for the O(ε 0 ) effective potentials of the Σ cD (X 1.1 ), Σ cD * (X 2.1 , H 2.1 ), Σ * cD (X 3.1 ) and Σ * cD * (X 4.1 , H 4,1 ) systems. We use the thin line to denote thē D meson, and other notations are the same as those in figure 1. Figure 3. The two-pion-exchange diagrams of the Σ cD system at O( 2 ). These diagrams are classified as the football diagram (F 1.1 ), triangle diagrams (T 1.i ), box diagrams (B 1.i ) and crossed box diagrams (R 1.i ). The internal heavy baryon lines in diagrams (T 1.3 ), (B 1.1 ) and (R 1.1 ) can also be the Λ c . The notations are the same as those in figure 2. Figure 4. The next-to-leading order Feynman diagrams that contribute to the vertex corrections and wave function renormalizations. Each graph denotes the one type of diagrams with the same topological structure.
The analytical expressions of the two-pion-exchange diagrams in figure 3 read When the contribution of Λ c is included, it will appear in the graphs (T (3.14) In above equations, the loop functions J y x are defined in appendix A. d is the dimension where the loop integral is performed and approaches four at last. E represents the residual energies of the Σ ( * ) c andD ( * ) , which is defined as ). E is set to zero in our calculations.
The leading order potential for the Σ cD * system stems from the contact interaction and onepion-exchange diagrams [graphs (X 2.1 ) and (H 2.1 ) in figure 2], which reads where σ is the Pauli matrix. The spin operator S 1 of Σ c satisfies S 1 = 1 2 σ. The operator T = iε * ×ε (ε and ε * are the space components of polarization vectors of the initial and finalD * meson) is correlated with the spin operator S 2 of theD * meson by the relation S 2 = −T . Thus the σ · T term represents the spin-spin interaction (see appendix C). Since only the S-wave interaction is considered, one can use the following replacement rules in the potentials, The two-pion-exchange diagrams for the Σ cD * system are shown in figure 5. The potentials from these graphs read Considering the contribution of Λ c : (3.31) From the above equations we see that, in the S-wave interactions, only the central terms and spinspin interactions survive in the effective potentials.

Σ * cD system
Like the Σ cD system, the leading order potential for the Σ * cD system only stems from the contact terms [diagram (X 3.1 ) in figure 2]. The expression reads We see that the contact potential of the Σ * cD system equals to the one of the Σ cD system in eq. (3.4), because the O(ε 0 ) contact Lagrangian is constructed in the heavy quark limit. The heavy quark breaking effect will be manifested in the loop diagrams when the mass splittings are considered in the propagators of the heavy matter fields. 3 The two-pion-exchange diagrams are illustrated in figure 6. The analytical results for these diagrams are given as Including the contribution of Λ c : The leading order diagrams for Σ * cD * system are the graphs (X 4.1 ) and (H 4.1 ) in figure 2.
The potentials from these two graphs read where the operator σ rs is related to the spin operator S 1 of the Σ * c with S 1 = 3 2 σ rs (see the detailed derivations in appendix C), so the σ rs · T term represents the spin-spin interaction as well. We see the O(ε 0 ) potentials for Σ * cD * resemble the ones for Σ cD * in eqs. (3.15) and (3.16).
The two-pion-exchange diagrams are displayed in figure 7. The potentials originate from these graphs read Unlike the two-pion-exchange potentials of the Σ cD * system, there exists a very simple relation represents that only the coefficient of J B 21 in the square brackets should be replaced with the one of J R 21 , while the other terms remain unchanged. For example, the Including the contribution of Λ c : s are equal to the ones in eq. (3.55) for i = 3 and i = 4, respectively. In the above equations, we notice that a new spin-spin structure (σ rs · T ) 2 arises in the box and crossed box diagrams, which is the characteristic interaction structure for the high spin particle systems. Such a structure cannot appear in the two-body potentials with spin-1 2 particle, such as the Σ cD * system.
Due to the constraints of the commutation and anticommutation relations of the Pauli matrix, the spin operator of a spin-1 2 particle appears at most once. On the other hand, the (σ rs · T ) 2 terms do not emerge at the tree level, where the heavy quark symmetry is satisfied. In other words, this structure is also the reflection of the heavy quark symmetry breaking effect at the one-loop level, which indeed disappears if we set the mass splittings in the loops to be zeros (this structure will persist for the diagrams with Λ c as the intermediate state, since the mass splittings δ c and δ d do not vanish even in the rigorous heavy quark limit).

Numerical results without and with the Λ c
The newly observed three P c states, P c (4312), P c (4440) and P c (4457) have been studied with the same framework in our previous paper [11], in which we did not include the contribution of the Λ c . There are three scenarios in ref. [11]. In scenario I, the LECs are estimated from the N -N data, but the result is not good, because we cannot reproduce the P c (4457). In scenario II, the LECs are determined by fitting the data of the three P c s, yet the result is still unsatisfactory. In scenario III, the P c s are simultaneously reproduced in a relatively small parameter region when the couple channel effect is included. In this part, we revisit these states without and with the Λ c contribution, and give a comparison with the result in scenario II of ref. [11].

The three P c states without the Λ c
Up to now, the four LECs in eq. (2.14) are still unknown. But we do not have to determine each of them since the forms of the O(ε 0 ) contact potentials are homogeneous for definite isospin states. There are only two independent LECs in nature if the isospin-isospin terms are absorbed into the relevant LECs with the following redefinitions, Thus the O(ε 0 ) contact potentials of the Σ ( * ) cD ( * ) systems can be rewritten as 1 The masses and widths of the newly observed three P c states and the previously reported P c (4380) are displayed in table 1. The closest thresholds, binding energies as the Σ  Table 1. The experimental and theoretical information of the P c (4312), P c (4440), P c (4457) [10], and P c (4380) [2]. The corresponding binding energies are obtained with the thresholds of Σ  With the above preparations, as in ref. [11], we vary the D 1 and D 2 in the ranges [−100, 150] GeV −2 and [−100, 100] GeV −2 respectively to search for the possible region where the three P c states can coexist. We mainly focus on the I = 1 2 states, because these P c states are observed in the mass spectra of J/ψp. To make a comparison, we present both the results without and with the Λ c in figures 8(a) and 8(b), respectively 2 . We assume that the P c (4312), P c (4440) and molecular states, respectively. We use three colored bands to denote the region of parameters with binding energy [−30, 0] MeV for each system, respectively. Considering the hadronic molecules are loosely bound states, we set −30 MeV as the lower limit of the bindings. The intersection point of two black solid lines designates the coordinate value (D 2 , D 1 ) where the corresponding two P c s can coexist. Ideally, the three The line-shape of the effective potentials for the three P c s in this case have been given in ref. [11], where a set of parameters D 1 = 42 GeV −2 and D 2 = −12.5 GeV −2 in the overlap region are adopted. Here, we use these two values to calculate the binding energies of the [Σ

Role of the Λ c
As mentioned above, we cannot give a good description for the P c s if we only consider the spin partners of Σ ( * ) c andD ( * ) in the two-pion-exchange diagrams. In this part, we are going to include the contributions of Λ c in the loops. Since both the Σ c and Σ * c can decay into Λ c π, the strong couplings between Σ ( * ) c and Λ c π should not be neglected.   The result with the Λ c being included is illustrated in figure 8(b), from which we find that there exists a very large overlap among the three colored bands. The small triangle surrounded by three straight lines just locates in the overlap. Besides, the intersection points between two of the three solid lines are very close to each other, and they almost meet at a point if we consider the experimental errors. In other words, the three P c can be synchronously reproduced in this case. The result in figure 8(b) is in good agreement with the experimental data.

∆E
We system might correspond to the previously reported P c (4380) [2]. Therefore, we urge the experimentalists to reanalyze the data to see whether P c (4380) is the most deeply bound one in the [Σ respectively. In the following, we analyze their behaviors in detail. Σ cD ( * ) systems: The results in figures 9(a), 9(b) and 9(c) all demonstrate that the contact term supplies the very strong attractive potential. From eq. (4.2) we know that O(ε 0 ) contact term for the Σ cD system only contains the central potential, while the spin-spin contact term appears for the Σ cD * system. Thus their difference is mainly caused by the spin-spin interaction. Meanwhile, the small difference between their O(ε 0 ) contact potentials indicates that the spin-spin interaction is rather weak and only serves as the hyperfine splittings.
There is no one-pion-exchange potential for the Σ cD due to the vanishingDDπ vertex.  of the two-pion-exchange potentials, which are repulsive at the short range, but become weakly attractive at the intermediate range. This is the typical feature of the nuclear force [32].  table 2, the result of Σ * cD is about eight times larger than that of the Σ cD . These two systems have the same O(ε 0 ) contact potentials [e.g., see eq. (4.2)]. The one-pion-exchange contribution vanishes for both systems. Thus the difference can only arise from the two-pion-exchange potentials, as shown in figure 10(a). One can notice the behaviors of the two-pion-exchange potential for the Σ * cD is attractive at the short-range and weakly repulsive at the intermediate-range, which is in contrast to that of the Σ cD [e.g., see figure 9(a)]. Therefore, if one only considers the O(ε 0 ) contribution, it is unlikely to obtain the significant difference between the Σ cD and Σ * cD systems. So we eagerly hope the future analysis at LHCb can help us confirm this observation. system as an example. The two-pion-exchange potential is attractive if we do not consider the Λ c , while it becomes repulsive when the Λ c is involved. This can well explain why the binding of the [Σ cD ] 1 2 state is much deeper without the Λ c (see table 2). The magnitude of the change from the minimum to the maximum in these two cases is about 120 MeV, which is even larger than the minimum of the total potential [see figure 9(a)]. The enhancement is mainly generated by the accidental degeneration of the Σ cD and Λ cD * systems, since the contribution of the box diagram (B 1.1 ) is MeV is tiny. Another reason that may cause the enhancement is the contributing diagrams with the Λ c are only (T 1.3 ) and (B 1.1 ). Unlike the Σ cD * system, the accidental cancelations among several diagrams cannot happen. In other words, the Λ c indeed plays a crucial role in the formation of the P c (4312).
For the [Σ cD * ] 1 2 system, since the whole contribution of the two-pion-exchange potential is much weaker than the O(ε 0 ) contact term [see figure 9(b)], the influence of Λ c on this state is not so apparent as in Σ cD . However, it is still very important to the existence of P c (4457) and the possible [Σ * cD * ] J bound states (e.g., see the data in table 2).
In figure 11, we also show the dependence of the two-pion-exchange potentials on the mass splitting δ c . One can see that they are very sensitive to the δ c . The loop integrals generally contain two structures. One is the analytic term, which is the polynomials of the m 2 π , q 2 , δ 2 , etc.. Another one is the nonanalytic term, which comprises the typical multivalued functions, such as log X and √ X (X is the polynomials of the m 2 π , q 2 , δ 2 .). The physical value of the δ c is about 168 MeV, which is larger than the pion mass m π . We then decrease its value to 100 MeV and 65 MeV. One can anticipate the dependence on δ is regular if the terms that make up the potential are only polynomials, but the variation trend in figure 11 is irregular. This phenomenon indicates the nonanalytic terms can distort the O(ε 2 ) potentials, which are vital to the formations of the P c states. The contributions of the nonanalytic terms incorporate the complicated light quark dynamics, which are almost impossible to estimate from quark models.
After the above discussions, one may wonder whether it is possible to reproduce the three

An episode: interchanging the spins of P c (4440) and P c (4457)
The J P quantum numbers of the P c (4312), P c (4440) and P c (4457) are not determined yet [10]. The theoretically favored J P for P c (4440) and P c (4457) in this paper and some previous works [11][12][13][14][15]  The result of interchanging the spin assignment of P c (4440) and P c (4457) is given in figure 12(b). Marvelously, the result is comparable with the one in figure 8(b), i.e., it seems this assignment can well describe the experimental data, likewise. However, one shall note that the values of the LECs (D 1 , D 2 ) in the center of the small triangle are (58, −31) GeV −2 , while these values in figure 8(b) are (52, −4) GeV −2 . The shift of D 1 in these two cases is small, but the D 2 in the first case is about eight times larger than that of the latter one. One has to largely enhance the contribution of the O(ε 0 ) spin-spin interaction to reverse the canonical order 3 of the spins of P c (4440) and P c (4457). 13. We notice the correspondence and consistence with the N -N system. The leading order contact Lagrangian for the N -N system reads [32],

The binding energies of the [Σ
where C S and C T are two independent LECs. One would see that they respectively correspond to the D 1 and D 2 in our work if we write out the O(ε 0 ) contact potential of the N -N system, The values of C S and C T have been precisely determined by fitting the np scattering phase shift at the next-to-next-to-next-to-leading order of chiral perturbation theory [62]. For the np system, which gives 4 If absorbing the minus sign of eq. (4.3) into the C S and C T , one would see the redefined C S and C T share the same sign with the D 1 and D 2 , correspondingly. Meanwhile, the ratio of the absolute values for C S and C T gives R ST = |C S |/|C T | 18, which is compatible with the R 12 for D 1 and D 2 . However, this ratio for the case of interchanging the spins of P c (4440) and P c (4457) is about 1.9, which is one order of magnitude smaller than the R ST , because of the spin-spin term in the contact potential is immoderately enhanced. On the one hand, from the point of potential model, the spin-spin term is suppressed by the factor 1/(m Σ ( * ) c m D ( * ) ) (e.g., see appendix D). On the other hand, one can build a mandatory connection between the contact terms of chiral effective field theory and the one-boson-exchange model with the help of resonance saturation model [63,64]. As the heavy fields, ρ, ω, f 0 , a 0 , etc., which are equally treated in one-boson-exchange model, are integrated out in chiral effective field theory, and their contributions are packaged into the LECs. The (ω, f 0 ) and (ρ, a 0 ) mesons account for the isospin-isospin unrelated D a and related E a , respectively [e.g., see eq. (2.14)] [35]. Meanwhile, the ω and ρ mesons couple to the matter fields via the P -wave interaction due to the parity conservation. Each vertex contains one momentum. In other words, the ω and ρ mesons are responsible for the momentum-dependent spin-spin interaction, which cannot be matched with the O(ε 0 ) D b and E b . Therefore, the momentum-independent contributions for D b and E b can only come from the axial-vector mesons, such as (h 1 , f 1 ) and (b 1 , a 1 ). The masses of these states reside around 1.2 GeV, which are much heavier than those of ω and ρ, and suppress the value of D 2 .

Hidden-bottom molecular pentaquarks
The above study for the hidden-charm pentaquarks can be extended to the hidden-bottom case, once the coupling constants and mass splittings are replaced by the bottomed ones. The coupling constants g 2 and g 4 for the bottom baryons can be calculated with the partial decay widths of Using the average values of the decay widths of Σ + b → Λ 0 b π + and Σ * + b → Λ 0 b π + [47], we get g 2 = −0.51, g 4 = 0.91. The other couplings can then be obtained with the relations in eq. (2.11), which yield, The axial coupling g of the B mesons cannot be directly derived from the experiments due to absence of phase space for B * → Bπ, so we adopt the average value g = −0.52 from the lattice calculations [65,66]. Similarly, the mass splittings are correspondingly given by where the masses of the Σ ( * )+ b and B ( * )0 are used [47].
The small scale expansion [67] is used in eqs. (2.8) and (2.12), i.e., the mass splitting δ is treated as another small scale in the Lagrangians. This expansion works well for the systems with one heavy matter field [43]. The loop integrals in these systems are the polynomials of δ, thus the convergence of the chiral expansion is not affected as long as the δ ∼ m π or smaller than m π . But the situation becomes different for the systems with two heavy matter fields. The loop integral of the box diagram is proportional to 1/(δ x + δ y ). If δ x + δ y is of the order of the pion mass, the convergence of the expansion could still be good. For example, for the Σ ( * ) cD ( * ) systems [11], MeV, which is much smaller that the pion mass 5 . Therefore, if we still adopt the same procedure as used in the Σ systems, the amplitudes of some typical box diagrams would be largely amplified, which results in extremely strong attractive or repulsive potential. This is unphysical and mainly caused by the poles of the heavy matter fields. In some previous works [11,37], the mass splittings are discarded in the box diagrams to subtract the 2PR contributions. Here, we develop a method to remove the heavy matter field poles in the box diagrams with the mass splittings being kept (see appendix B for more details).
In order to predict the possible P b states, we also need to know the LECs D 1 and D 2 for the Σ ( * ) b B ( * ) systems. In principle, they should be fixed from experimental data or the results from lattice QCD, which are not available at present. Thus, we estimate the ranges of D 1 and The masses of the hidden-bottom molecules are all above 11 GeV. Like their P c partners, they may be observed from the Υ(1S)N and Υ(2S)N final states. We hope future experiments to hunt for these P b states. We conclude this section by borrowing one of the famous phrases from R. P. Feynman: "There is plenty of room at the 'bottom'." [68]

Heavy quark symmetry breaking effect
The QCD Lagrangian has heavy quark symmetry (HQS) when the heavy quark mass m Q → ∞. For a heavy hadron containing one single heavy quark, the strong interaction would be independent of the heavy flavors in this limit. Meanwhile, the heavy quark will decouple with the light degrees of freedom. The multiplet associated with the heavy quark spin would be degenerate in the heavy quark limit. However, the physical masses of the heavy quarks are finite, such as m c ∼ 1.5 GeV, m b ∼ 5 GeV. Therefore, the effects of the heavy quark flavor symmetry breaking and spin symmetry breaking are explicit. For example, the axial coupling g for B * Bπ is about 17% smaller that that of the D * Dπ. Thus the value 17% can be roughly regarded as the breaking size of the heavy quark flavor symmetry. In addition, the heavy quark spin symmetry (HQSS) breaking is more obvious, such as the mass splittings of (B * , B) and (D * , D) are about 45 MeV and 142 MeV, respectively.
HQS can be used to relate the coupling constants to one another, such as the axial coupling constants in the O(ε 0 ) Lagrangians. Under HQS, the heavy quarks only serve as the spectators. The interaction between Σ ( * ) c andD ( * ) is mediated by their inner light degrees of freedom, i.e., the light diquarks in Σ ( * ) c and the light quark inD ( * ) . Therefore, the S-wave effective potentials between Σ ( * ) c andD ( * ) at the quark level can be parameterized as [11] V HQS quark−level = V c + V s l 1 · l 2 , (6.1) where V c and V s denote the central term and spin-spin term, respectively. l 1 and l 2 are the spins of the light degrees of freedom of the Σ ( * ) c andD ( * ) , respectively. With the potentials at the quark level, one can build the relations between different channels at the hadron level by parameterizing the hadron level potentials as where S 1 and S 2 are the spin operators of the Σ ( * ) c andD * , respectively. One can easily verify 6 3) The leading order potentials obviously satisfy the above relations obtained from HQS. One can also testify the one-loop level analytical expressions satisfy the above relations as well when d → 4 and δ a,b → 0.
The HQS breaking effect would manifest itself in the loop diagrams if δ a,b = 0 7 . When δ a = δ b = 0, all the box diagrams would become the 2PR ones, thus we have to remove the 2PR contributions. In order to compare with the cases of δ a = δ b = 0, we also subtract the 2PR contributions from the (δ a , δ b ) = 0 cases (see appendix B). The 2PI two-pion-exchange potentials for the Σ (m π , ω, q), (6.4) where A and B are the positive numbers. The corresponding J T ij functions generally contain two structures, the odd function of ω and the even one (see appendix A). The odd part is proportional to where F (x, . . . ) and G (x, y, . . . ) denote the integrands which are the functions of (x, . . . ) and (x, y, . . . ), respectively. These two terms vanish in the heavy quark limit, i.e., when δ a = δ b = 0. The even one is proportional to Only this term contributes when δ a = δ b = 0. Therefore, one would see a different scenario when the nonzero mass splittings are considered, since the two terms in eq. (6.5) also contribute. For the diagrams (T 1.1 ) and (T 3.1 ), ω = −δ b , while for the diagrams (T 2.2 ) and (T 4.2 ), ω = δ b . Thus the HQS breaking effect is totally different for the Σ ( * ) cD and Σ ( * ) cD * systems, because the eq. (6.5) is the odd function of ω, which is very sensitive to the sign of the ω. In addition, the integrands F , G and H always have the nonanalytic structures, such as the logarithmic and square root terms. So the variations of the graphs (T 1.2 ) and (T 3.3 ) are not so dramatic as those of the (T 1.1 ) and (T 3.1 ), because δ a < m π , whereas δ b > m π . The HQS breaking effect expounded above issues from the loop diagrams, which is the quantum physics of the light degrees of freedom at the low energy, and cannot be modified by any unknown physics that happens at the high energy.

Summary and conclusion
In the April of this year, the LHCb collaboration reported the observation of the three pentaquark states P c (4312), P c (4440) and P c (4457) [10]. They were subsequently interpreted as the molecular states by many theoretical works [11][12][13][14][15] due to the proximities to the Σ cD and Σ cD * thresholds. In this paper, we have systematically investigated the interactions between the charmed baryons Σ ( * ) c and anticharmed mesonsD ( * ) in the framework of chiral perturbation theory. To this end, we have simultaneously considered the short-range contact interaction, long-range onepion-exchange contribution, intermediate-range two-pion-exchange loop diagrams, as well as the influence of the mass splittings on the effective potentials.
When we fix the total isospin as I = 1 2 , the original four independent LEC can be reduced to two. These two LECs can be fitted using the binding energies of the P c (4312), P c (4440) and P c (4457) as inputs. We first attempt to reproduce the newly observed three P c s via only considering the spin partners of the Σ cD ( * ) in the loops. But we fall into the same dilemma as in the scenario II of our previous work [11], i.e., it is nearly impossible to reproduce the three P c s synchronously in this case. Considering the strong couplings between Σ ( * ) c and Λ c π, we then include the contribution of Λ c in the loop diagrams. Three P c s are simultaneously reproduced at this point. This indicates the Λ c plays a very important role for the formation of these P c s. We also notice that only considering the Λ c cannot describe the P c s either. The subtle interplay between the channels with Λ c and the ones with Σ 3/2 molecular state, is a deeper bound hadronic molecule in our calculation. This is mainly caused by the important contribution from the two-pion-exchange diagrams, which is the essential difference with the predictions from the quark model and leading order effective field theory. These two approaches do not contain the nonanalytical terms, such as the powers of log q 2 and q 2 , which are irregular and may give the enhanced contributions sometimes. These terms cannot be predicted accurately from the aspects of the quark model.
We also study the hidden-bottom Σ 1/2 1/2 . The hidden-bottom ones might be observed from the ΥN final states. We give a complete picture on the mass spectra of the hidden-charm and hidden-bottom molecular pentaquarks, and there are overall fourteen bound states in our calculations. The discovery of P c s at the LHCb is just the beginning for the community to search for the exotic multiquark matters.
The heavy quark symmetry is always exploited to predict the mass spectra of the hiddencharm and hidden-bottom systems. Since m b is much larger than m c , so the predictions from the HQS in the bottom sector is more reliable because the correction from the next-to-leading order heavy quark expansion is very small. But the reliability of the HQS in the charm sector is still questionable. So we examine the HQS breaking effect in the loop diagrams by considering the mass splittings in the propagators of the intermediate states. As expected, the HQS in the hiddenbottom systems is much better than that in the hidden-charm cases. Besides, for some accidental reasons, the HQS as an approximation in the Σ cD and Σ * cD systems is not as good as in the others. The two-pion-exchange potentials become totally different with the mass splittings or not. One reason is the mass difference between the initialD and intermediateD * is −δ b and some triangle diagrams are very sensitive to the sign of the mass difference. Another reason is δ b > m π , so the nonanalytic structures, e.g., logarithmic and square root terms in the loop functions would be enhanced to distort the potentials. This enlightens us that the HQS breaking effect shall not be ignored if we want to give a comprehensive description of the effective potentials, especially for the interactions between the charmed hadrons.
We hope the lattice QCD simulations on the hidden-charm and hidden-bottom pentaquark systems could be carried out in the future, which can help us to get a deeper insight into the inner structures of these exotica. The analytical expressions derived in this work can also be used to perform the chiral extrapolations.
where the notation X ∨ Y ∨ Z ∨ · · · represents the symmetrized tensor structure of X α Y β Z γ · · · + · · · , which are given as g ∨ q 2 ≡ q α q β g γδ + q α q δ g βγ + q α q γ g βδ + q γ q δ g αβ + q β q δ g αγ + q β q γ g αδ , These J functions can be directly calculated with the dimensional regularization in d dimensions, or by an iterative way as shown in ref. [11]. Their detailed expressions read The L is defined as where γ E is the Euler-Mascheroni constant 0.5772157. We adopt the MS scheme to renormalize the loop integrals.

B Removing the 2PR contributions
Sometimes, we need to subtract the 2PR contributions from the box diagrams, which can be recovered by inserting the one-pion-exchange potentials into the iterative equations. For the case of ω = δ = 0, the 2PR part must be discarded due to the pinch singularity. This can be easily done by using the simple derivative relation given in eq. (A.15). In this part, we develop a new method to make such a subtraction with the help of the principal-value integral method. In this way, we can subtract the 2PR part in a diagram with nonvanishing mass splittings, which has no pinch singularity. Considering the loop integral of a box diagram with the following form, where the Lorentz structure L µν···α (l) ≡ l µ l ν · · · l α . This integral can be straightforwardly disassembled into two parts through the following way, The principal-value integral method tells that If we replace the x with the v · l + ω + i and −v · l + δ + i , the integral can be divided into two parts, the principal-vale part and the Dirac delta part. The Dirac delta part is the pole contribution of the matter fields, which corresponds to the 2PR part in the time-ordered perturbation theory.The principal-value part is just the 2PI contribution. In other words, the 2PI part of the integral I can be written as As long as we can derive the form of I 2PR , we could obtain I 2PI , since the complete form of I has been given in appendix A. The calculation of the I 2PR is simple due to the special property of the delta function. We take the calculation of the I 2PR part of J B 21 as an example. We first show the concrete form of the I 2PR , where we have used the Feynman parameterization to combine the denominators of the propagators of the light pseudoscalars, and M 2 = x(x − 1)q 2 + m 2 . Besides, we have also utilized the approximation v · q 0 in the two delta functions. Choosing L µν···α (l − xq) to be L µν (l − xq) we would be in the position to calculate the 2PR part of the J B 2i (denoted by J B 2i 2PR ). For the J B 21 2PR , we have where the function N (ω) = x(x − 1)q 2 + m 2 − ω 2 − i . Following the same procedure given above, we can get all the 2PR parts of the J B ij functions.
One can avoid the lengthy and tedious calculations by adopting another trick. The loop integrals of the box diagrams can be constructed from the ones of the triangle diagrams [e.g., see eq. (A.15)], the finite part of the loop functions J T ij that make up the J B ij actually contains two types of functions, one is the odd function of ω, and the other one is the even function of ω. Therefore, the renormalized J T ij can be written as where O T ij (ω) and E T ij (ω) represent the odd and even parts of the J T ij (ω), respectively. The other two variables m and q are omitted for simplicity. It can be proved that O T ij (ω) and E T ij (ω) account for the 2PI and 2PR parts of the J B ij , respectively. For example, we find the − 1 16π 1 0 dx N (ω) in eq. (B.7) is just the opposite of the E T 12 (ω). With the simple properties of the odd and even functions, we can readily obtain When ω and δ approach to zero, this formula evolves into the derivative relation in eq. (A.15). One can easily testify the remainder ones indeed satisfy the eq. (B.9), likewise.

C Spin transition operators
In calculating the loop diagrams of the Σ * cD * system, we encountered some intractable scalarproducts, such as (ū · ε * )(u · ε) and (ū · ε)(u · ε * ), where the u µ denotes the spinor-vector of the spin-3 2 Rarita-Schwinger field ψ µ , and ε µ represents the polarization vector of the spin-1 fieldP * µ .ū µ and ε * µ are their conjugations, respectively. We notice that these structures involving polarizations can be transformed into the spin-spin interaction terms by introducing the so-called spin transition operators for the spin-3 2 and spin-1 fields, respectively.

C.1 Vector field
In the rest frame of a vector particle, the space components of the polarization vectors with different helicity λ = 0, ±1 read, The ε(λ) can be obtained with the following relation, where S t is the spin transition operator for the spin-1 field. The matrix form of the S t is One can easily verify that D A tentative parameterization of the effective potential from the quark model Assuming a pair of c andc quarks are produced in the high energy colliding process, and they are surrounded by the largely separated light quarks u and d. At the very short c andc separation r, the c andc quarks interact with the perturbative one-gluon-exchange Coulomb potential. There is essentially no screening of the cc interaction due to the much farther separated u and d quarks. Before the hadronization occurs, the effective potential at this size scale can be written as [69,70] V ij (r i , s i , r j , s j ) = −C α s 4 1 |r i − r j | − δ 3 (r) 8π 3m i m j s i · s j + · · · + · · · , (D.1) where we only show the central term and spin-spin interaction. Other terms such as the tensor force and spin-orbit interaction are omitted for the S-wave case. The C denotes the color factor. α s is the strong coupling constant. r i , s i and m i represent the position, spin, and mass of the i-th quark, respectively. We need the cc color singlet to supply an attractive core, thus C = 16 3 . In order to avoid the c andc pair to rapidly move far away from each other with large velocity, we assume that the cc pair is produced near the threshold. When the distance between the slowly moving c andc increases, the light quarks u and d start to screen the color interaction at this point. Then the five quarks form two weakly interacting color singlet clusters Σ ( * ) c andD ( * ) . The force between them is nothing but just the residual color interaction similar to the van der Waals force between neutral molecules. At this size scale, the attractive core from c andc still works, but attenuates rapidly with the increase of the separation r. At the same time, the heavy quark spin decouples, and the spin-spin interaction is transferred to their inner light degrees of freedom. If ignoring other higher order contributions, one could roughly parameterize the potential as follows, where A and B are two independent constants with the same dimension, which can be determined by fitting the data. l 1 and l 2 denote the spins of the inner light degrees of freedom of Σ ( * ) c andD ( * ) , respectively 8 . d ∈ [1,2] fm stands for the characteristic size of a hadronic molecule, we choose the upper limit d = 2 fm. x is always chosen to be 1.5 or 2 for some phenomenological considerations. Here we use x = 2 as in ref. [71]. Obviously, the strength of the spin-spin term is suppressed by the factor 1/(m Σ ( * ) c m D ( * ) ). By fitting the binding energies of the P c (4312), P c (4440) and P c (4457), we obtain A 2.45 We see the newly observed three P c s can be simultaneously reproduced, and other four systems all have the binding solutions. The ∆E for [Σ * cD ] 3 2 in the quark model is smaller than that of the chiral perturbation theory. In addition, the bindings of the [Σ * cD * ] J systems are larger than 8 The matrix element of l1 · l2 can be found in ref. [11].