$\bar{B}^{(\ast)} \bar{B}^{(\ast)}$ interactions in chiral effective field theory

In this work, the intermeson interactions of double-beauty $\bar{B}\bar{B}$, $\bar{B}\bar{B}^\ast$, and $\bar{B}^\ast\bar{B}^\ast$ systems have been studied with heavy meson chiral effective field theory. The effective potentials are calculated with Weinberg's scheme up to one-loop level. At the leading order, four body contact interactions and one pion exchange contributions are considered. In addition to two pion exchange diagrams, we include the one-loop chiral corrections to contact terms and one pion exchange diagrams at the next-to-leading order. The behaviours of effective potentials both in momentum space and coordinate space are investigated and discussed extensively. We notice the contact terms play important roles in determining the characteristics of the total potentials. Only the potentials in $I(J^P)=0(1^+)$ $\bar{B}\bar{B}^\ast$ and $\bar{B}^\ast\bar{B}^\ast$ systems are attractive, and the corresponding binding energies in these two channels are solved to be $\Delta E_{\bar{B}\bar{B}^\ast}\simeq -12.6^{+9.2}_{-12.9}$ MeV and $\Delta E_{\bar{B}^\ast\bar{B}^\ast}\simeq -23.8^{+16.3}_{-21.5}$ MeV, respectively. The masses of $0(1^+)$ $\bar{B}\bar{B}^\ast$ and $\bar{B}^\ast\bar{B}^\ast$ states lie above the threshold of their electromagnetic decay modes $\bar{B}\bar{B}\gamma$ and $\bar{B}\bar{B}\gamma\gamma$, and thus they can be reconstructable via electromagnetic interactions. Our calculation not only provides some useful information to explore exotic doubly-bottomed molecular states for future experiments, but also is helpful for the extrapolations of Lattice QCD simulations.


I. INTRODUCTION
Hunting for exotic multiquark states beyond the conventional meson and baryon configurations is a long-standing problem of QCD [1][2][3]. After the pioneering exotic state X(3872) was discovered in 2003 by the Belle Collaboration [4], a new perspective of hadronic physics was opened, and plenty of exotic hadrons called XYZ states were observed at experiments [5,6]. Lots of experimental and theoretical efforts have been paid to understand the nature of these states, but they still seem very elusive.
Most of these states cannot be assigned into the framework of conventional quark model. For example, if we treat X(3872) as the χ c (2P) charmonium, the mass of it would be about 78 MeV lower than the prediction of relativized quark model [7]. Furthermore, the isospin violation process X(3872) → J/ψρ with the substantial decay width makes the situation more complicated and mysterious. Thus various pictures and explanations were put forward to probe the inner structure of X(3872), such as molecular model [8][9][10][11][12][13], tetraquark state [14,15], traditional axial vector charmonium [16,17], Lattice QCD simulation [18,19], and other theoretical schemes [21] (for a review, see Ref. [5]). Among them, the shallowly bound molecular picture, i.e., the deuteronlike configuration, is the most popular one, because the mass of X(3872) is very close to the threshold of D 0D * 0 . Like study the interaction of proton and neutron is a sine qua non for understanding the characteristics of deuteron and nuclear force, if we want to comprehensively understand the nature * Electronic address: bo-wang@pku.edu.cn † Electronic address: liuzhanwei@lzu.edu.cn ‡ Electronic address: xiangliu@lzu.edu.cn of X(3872) and other resonances in XYZ family, to research the strong interactions between heavy mesons would be an ineluctable key issue.
Unlike a plethora of XYZ states have been observed in charmonium energy region [5,6,20], only two Z b states were reported in bottomonium spectrum [22]. There is also no unanimous conclusion on what the Z b (10610) and Z b (10650) really are [5], but one promising interpretation about the two charged bottomonium-like states is molecular conjecture which are composed of BB * and B * B * , respectively [23][24][25][26][27]. Stimulated by the continuously experimental observations of novelly hidden heavy flavor exotic states over the past decade, the interactions between heavy meson-heavy antimeson systems have been intensively investigated. For example, before Z b (10610) was observed in 2011, Liu et al had predicted a loosely bound S -wave BB * molecular state in the framework of one boson exchange (OBE) model [12], and later, the BB * and B * B * systems were again studied with OBE model [24-26, 28, 29]. In Ref. [30], the pion-mediated interaction associated with coupled-channel effect is exploited to study the BB * and B * B * systems. And some other related works such as QCD sum rule [31], Bethe-Salpeter approach [32], coupled-channel model [33], and so on.
As mentioned above, up to now, the observed exotic states are mainly composed of QQqq or QQqqq (where Q denotes b or c quark, and q stands for u or d quark) [5]. Whereas, nowadays, the exotic clusters with double charm or bottom still conceal themselves from the field of our view. Fortunately, LHCb Collaboration reported the observation of doubly charmed baryon Ξ ++ cc very recently [34], which triggered many discussions on whether the stable QQqq tetraquark states can exist in nature. For instance, in Ref. [35], Karliner et al predicted a doubly bottomed tetraquark bbūd with J P = 1 + at 10389±12 MeV, which lies far below the B −B * and electromagnetic decay B −B0 γ thresholds. Eichten et al also  [39], where I tot and S tot designates the total isospin and total spin of two B mesons, respectively. Due to the constraints of symmetry, the quantum numbers ofBB andB * B * systems must satisfy the the selection rule L + S tot + I tot + 2i = Even number, where L is the orbital angular momentum between two B mesons, and i denotes the isospin of one B meson, which is 1/2.

System
Total isospin Total spin all L B * B * 1 even L odd L even L 0 odd L even L odd L found the double-beauty state with constituent quarks bbūd is stable against strong interactions [36]. Based on the quarkdiquark symmetry and the input of the mass of Ξ ++ cc , Ref. [37] indicated the existence of a stable double-bottom tetraquark state in I = 0 channel.
Inspired by the observation of doubly charmed baryon Ξ ++ cc [34] and many theoretical works as mentioned before, in this work we use chiral effective field theory (χEFT) to study the physically allowedB ( * )B( * ) ([bq][bq]) (see Tab. I), i.e.,BB, BB * , andB * B * intermeson interactions, to see in which channel the potential is attractive and to search for the possible bound states. χEFT has been widely employed to study the interactions of nucleon-nucleon (N-N) and meson-meson (M-M) systems with flying colours. For the successful application of χEFT in N-N systems, one can see some important works in Refs. [54][55][56][57] and the references therein, or see Ref. [58] for an introduction to χEFT. Here, we give a short review about the use of χEFT in M-M systems, mainly focus on heavy meson sectors. In Ref. [10], AlFiky et al proposed an approach with respecting the heavy quark symmetry and chiral symmetry to deal with X(3872) as a molecular state of D 0D * 0 . Similarly, there are also an abundance of investigations on DD * (BB * ) systems with χEFT, concerning the molecular assumptions of X(3872) and two Z b states [59][60][61][62][63][64][65][66][67][68][69] (or see Refs. [5,6] for a review). In addition to utilizing χEFT to study heavy mesonheavy antimeson systems, applying χEFT to heavy mesonheavy meson systems (such as D ( * ) D ( * ) , D ( * )B( * ) and B ( * ) B ( * ) ) is also an interesting topic. In our two previous works, first in Ref. [70], we gave a tentative usage of χEFT inBB system, and obtained the strong interaction potentials in momentum space. Then in Ref. [71], we calculated S -wave DD * potential at one-loop level, noticed in the 0(1 + ) channel it is attractive, and obtained a bound state with the binding energy ∆E DD * = −15.6 +14.7 −19.2 MeV. We also notice that in Ref. [72], both open charm and bottom states (D ( * )B( * ) ) are studied with heavy meson chiral effective field theory (HMχEFT).
Based on Refs. [70,71], we naturally extend our study toB ( * )B( * ) systems, and for these kinds of states their typical quark configuration is [bq][bq]. To inquiry whether there exist such heavy flavor molecular states with double bottom or not, in our work we consider the S -wave (L = 0) interactions of differentB ( * )B( * ) systems. According to the selection rule listed in Tab. I, the I(J P ) numbers for the physically allowed states are: 1(0 + ) forBB; 1(1 + ) and 0(1 + ) forBB * ; 1(0 + ), 1(2 + ), and 0(1 + ) forB * B * . In order to get the interaction potentials for these different channels in coordinate space, we firstly calculate the amplitudes in momentum space with S U(2) HMχEFT up to one-loop level, where the amplitudes include the contributions from four-body contact interaction (FBCI), one-loop corrections to four-body contact interaction, one-pion exchange (OPE) contributions, one-loop corrections to one-pion exchange parts, and two-pion exchange (TPE) contributions. After subtracting the two-particle reducible (2PR) contributions in some typical Feynman diagrams, we make the Fourier transformation on the potentials in momentum space, and thus the potentials in coordinate space can be obtained. Finally, by solving the nonperturbative equations such as Schrödinger equation, Lippmann-Schwinger equation, and so on (in this paper, the Schrödinger equation is solved numerically), one can not only recover the 2PR contributions, but also get the binding energy ∆E in the attractive channels. In this way, we can predict the possible molecular states inB ( * )B( * ) systems and make a comparison with other phenomenological models [39,73]. This paper is organized as follows. After the introduction, we show the effective Lagrangians used in this work and the Weinberg's formalism in Sec. II. The calculations of effective potentials and the analyses are represented in Sec. III, the contributions of the low energy constants at the next-to-leading order are estimated in Sec. IV, the results are summarized in Sec. V, and some needful formulas are given in the Appendix.

A. Effective Lagrangians
In the framework of HMχEFT, the scattering amplitudes can be expanded order by order with a small parameter = q/Λ χ , where q is either the momentum of Goldstone bosons or the residual momentum of heavy mesons, and Λ χ represents either the chiral breaking scale or the mass of a heavy meson. Except forBB, the scattering amplitudes of bothBB * andB * B * channels at leading order O( 0 ) have the contributions from FBCI and OPE, which are described by the leading effective Lagrangians in Eq. (1) [10,61,70] and Eq. (2) [74][75][76], respectively.
Tr Hγ µH Tr Hγ µH +D b Tr Hγ µ γ 5H Tr Hγ µ γ 5H +E a Tr Hγ µ τ aH Tr Hγ µ τ aH where D a , D b , E a , and E b are four independent low energy constants (LECs), which are determined later. τ a represents the pauli matrix.
where v = (1, 0) is the four-velocity of heavy mesons, and H denotes the degenerated B and B * doublet in heavy quark limit, which can be expressed as: The last term in Eq. (2) accounts for the mass shift of B and B * , which will not vanish in chiral limit, and ∆ is the mass difference of (B, B * ) doublet. The tensor σ µν is defined as i[γ µ , γ ν ]/2. Besides, the chiral connection Γ µ and axial vector current u µ are illustrated as follows, where ξ = exp(iφ/2 f ). f is the pion decay constant, and the triplet pion field φ is defined as As indicated in Ref. [70], other forms of FBCI at the leading order with different Lorentz structures are not independent, such as Tr Hγ µH Hγ µH , Tr HH Tr HH , and Tr Hσ µνH Tr Hσ µνH , which can be expressed as the linear combinations of the terms involved in Eq. (1). Some terms like Tr Hγ 5H Tr Hγ 5H and Tr Hγ 5 τ aH Tr Hγ 5 τ aH vanish in heavy quark limit.

B. Weinberg's formalism
Before formally carrying out our calculations, we give a brief introduction to Weinberg's power counting scheme [77,78] since it is the core of the theoretical calculations in this work.
FIG. 1: Two types of 2PR Feynman diagrams of N-N scattering, which will also appear in this work.
We take the two-nucleon scattering process as an example. Considering the one-loop Feynman diagrams displayed in Fig.  1, in the non-relativistic limit the scattering amplitude of Fig.  1(a) is proportional to Inspecting the above loop integral we can find that the integration over 0 is ill-defined because it has poles above and below the real axis at 0 = ±iε (which is always called as pinch singularity [56,78]). This problem can be cured by including the kinetic energy of the nucleon in the leading terms but not treating it as a perturbation. Then the pole positions of the two nucleons' propagators will be shifted to where E = p 2 /2M N , and M N is the mass of the nucleon. As is pointed out by Weinberg [77,78], after the 0 integration being performed with residue theorem the contribution from the pole of nucleon is enhanced by a large factor M N /| p|, and this strong enhancement explicitly breaks the naive power counting with which the 0 integral should have been of O(1/| p|).
With the formalism given in Ref. [78], we concentrate on the effective potential (i.e., the irreducible graphs) but not directly calculate the scattering amplitudes. Namely, removing the 2PR part contributions that originate from intermediate on-shell nucleon states ("infrared enhancement"), the irreducible diagrams that make up the effective potential can then be calculated perturbatively. The 2PR part will be automatically recovered when the effective potential (or kernel) is inserted into the Schrödinger equation (or Lippmann-Schwinger equation).
We will encounter the same problem as discussed above when studying the interactions of two B mesons with HMχEFT. We follow the ideas in studying N-N system [78], i.e., we calculate the effective potentials of two B mesons first. One can find more details of calculations in Appendix A.
III. EFFECTIVE POTENTIALS OFB ( * )B( * ) SYSTEMS A.BB system TheBB system has been studied in Ref. [70] with HMχEFT, but in Ref. [70] the four LECs were unevaluated and only the TPE contributions in momentum space were shown, and thus the behavior of the total potential is still ambiguous. Therefore, in this part, we revisit this system and give a more intuitional form of the potential in coordinate space to see whether a bound state can exist inBB system.
At O( 2 ), the amplitudes receive both the one-loop diagrams for FBCI corrections and TPE diagrams which are illustrated in Fig. 2 and Fig. 3, respectively. The amplitudes correspond to the diagrams in Fig. 2 are listed in Eqs. (12)-(16), where the coefficients C x are flavor and LECs dependent, their concrete values are: . (17) Various scalar functions J y x are defined in Appendix A 1, m π is the mass of pion, D is the dimension where the loop integral is performed and approaches 4 at last. E represents the residual energy ofB andB * mesons, and is defined as E = EB( * ) − MB( * ). {X} r denotes the finite part of X [70], which is equivalent to make use of the modified minimal subtraction (MS) scheme. The amplitudes of the diagrams F 1.1 , Here we need to stress that in rigorous heavy quark limit, i.e., when the mass difference ∆ = 0, the amplitudes induced by the diagrams h 1.1 , h 1.2 , and B 1.1 in Fig. 2 and Fig. 3 would approach to infinity. Thus in order to get the correct potential, we employ Weinberg's formalism to subtract the 2PR contributions that stem from the double poles of the two intermediate heavy vector mesons. The contributions from the B * meson poles are enhanced in the 1/M B expansion and repeat the results of the iterated OPE diagrams, and thus we only need to consider the contribution from pion poles. The related details are shown in Appendix A 2 and A 3.
The effective potential V can be derived from the 2PI am-plitude Y 2PI by the relation where the factor −1/4 comes from OBE model [28,73], which originates from Breit approximation. The parameters used in the calculations are: ∆ = 0.045 GeV, m π = 0.139 GeV, f = 0.092 GeV, g = 0.52 [79][80][81], and the renormalization scale λ = 4π f . For simplicity, as in the OBE model [12], the transferred energy q 0 and residual energy E of heavy mesons are both set to be zero.
By calculating Eqs. (12)- (16), we obtain the result of one-loop corrections to FBCI, We can see the contribution of D a vanishes, while D b and E b emerge at O( 2 ). In addition, the result in Eq. (25) is about one order of magnitude smaller than that of Eq. (11), which manifests the convergence of chiral correction is very good.
Because the FBCI and one-loop corrections to FBCI describe the short distance interactions, and they are independent of the transferred momentum | q|. Thus the | q| dependent part can only come from TPE. The result of TPE in momentum space is displayed in Fig. 4(a). We can see that the potential is negative when | q| is varied from 0 to 300 MeV, which indicates the TPE part supplies the attractive potential, and this potential would be deeper if | q| becomes larger.
Because the potential in coordinate space can give more intuitional information, so we make the Fourier transformation on the potential in momentum space by the following way, We should note that, because the potential V(q) is expanded as the power series of the transferred momentum q in χEFT, so we need to regularize V(q) to suppress the contributions from high momentum to avoid the ultraviolet divergence of the integral, which is the manifestation of the fact that χEFT is valid only for q Λ χ . Different approaches in Refs. [82][83][84][85][86][87][88][89][90][91][92][93][94][95][96] were developed to regularize/renormalize V(q) non-perturbatively.
In this work, as in Refs. [71,[97][98][99], we adopt a simple Gauss regulator F (q) = exp(−q 2n /Λ 2n ), and we set n = 2 as in Ref. [99]. The value of the cutoff parameter Λ is commonly below the mass of ρ meson in the N-N case [100], so we use a moderate value Λ = 0.7 GeV as adopted in Ref. [71] to give predictions.
If we want to get the numerical results of the effective potential in coordinate space, we have to know the concrete values of the four LECs (see Eq. (1)). The values of the LECs have been determined in Ref. [71] by exploiting the resonance saturation model [101][102][103], which read (all in units of GeV −2 ), With the preparations above, the effective potential of 1(0 + ) BB system in coordinate space is given in Fig. 4(b). From the figure, we can read that although the TPE potential is attractive, the contribution from FBCI is dominantly repulsive, and thus the total potential is undoubtedly repulsive. This result indicates that it is impossible to find a bound state solution in 1(0 + )BB system, which is in coincidence with the calculations of Lattice QCD [50].

B.BB * system
For the S -waveBB * system, there are two different isospin states, i.e., 1(1 + ) and 0(1 + )BB * states. At the leading order, both FBCI and OPE diagrams contribute to the scattering amplitudes, which are illustrated in Fig. 5, at the next-to-leading order, the amplitude is composed of one-loop corrections to FBCI and OPE diagrams, and TPE diagrams, which are shown in Fig. 6, Fig. 7, and Fig. 8, respectively.
According to the effective Lagrangians given in Eq. (1) and Eq. (2), we can easily get the amplitudes of the elastic scattering processB(p 1 )B * (p 2 ) →B(p 3 )B * (p 4 ) from the diagrams H 2.1 and X 2.1 in Fig. 5, where the subscript I denotes the isospin of the channel, q is the transferred momentum carried by pion, 2 and * 4 designate the polarization vectors of initialB * and finalB * , respectively. In this section we only give the amplitudes of I = 1 channel, and the amplitudes of I = 0 channel are shown in Appendix B.
For the amplitudes of the one-loop corrections to FBCI of BB * system in Fig. 6, similar toBB case, we also write them out in the following,  where the values of C x are and the amplitudes of diagrams g 2.5 , g 2.6 , h 2.3 , and h 2.4 can be obtained by the following relations, The amplitudes of one-loop corrections to OPE diagram ( Fig. 7) are illustrated as follows, where Σ(m π ) = m 2 π /(24π 2 f 2 ) ln(m 2 π /λ 2 ) [58], the amplitudes of diagrams c 2.2 , c 2.5 , and c 2.6 can be obtained by the relations Here we list the amplitudes of the diagrams F 2.1 ∼B 2.3 in Fig. 8. The amplitudes of the diagrams R 2.1 ∼R 2.3 can be related to the ones of B 2.1 ∼B 2.3 , respectively. For conciseness, the terms that involve q 0 have been omitted since q 0 is set to be zero in our calculations. One can find the unabridged forms in Ref. [71].
. Notations same as in Fig. 2.
When q 0 = 0, the amplitude structures of diagrams In order to get the numerical results, some procedures are the same as we calculate theBB system, i.e., Eq. (24) should be used, and the 2PR part in diagrams h 2.1 , h 2.3 , and B 2.2 must be subtracted with the aid of Weinberg's formalism. However, some additional steps still needed, such as the pion decay constant f its bare value 0.086 GeV [71] should be adopted. Since we only consider the S -wave interaction, the terms ( 2 · * 4 ) and (q · 2 )(q · * 4 ) appeared in above equations should be replaced by [28,71,73] We first list the results of one-loop corrections to FBCI in I = 1 and I = 0 channel, respectively.
By comparing with Eq. (28) and Eq. (B1), the results again show the convergence of chiral correction is very good. The effective potentials of 1(1 + ) and 0(1 + )BB * systems in momentum space and coordinate space are shown in Fig. 9 and Fig. 10, respectively. From Fig. 9(a) and Fig. 10(a), we can see that the O( 0 ) OPE potential is attractive for both I = 1 and I = 0 channels, and the O( 0 ) OPE potential in I = 0 channel is more attractive than that of I = 1 channel. The O( 2 ) potential generated by one-loop corrections to OPE diagram is very small compared with the O( 0 ) OPE potential, which demonstrates that the convergence of chiral corrections to OPE diagram is also good. The behavior of TPE potential is totally different between these two channels. In I = 1 channel, TPE potential is attractive and it intends to be more attractive when | q| becomes larger. In I = 0 channel, TPE potential is repulsive, and would be more repulsive when | q| becomes larger.
Because the OPE and TPE potentials are all attractive in I = 1 channel, the total potential of I = 1 channel is attractive without the doubt. Although OPE provides an attractive potential in I = 0 channel, the TPE potential is repulsive and the absolute value of TPE potential is larger than OPE potential, which gives rise to a repulsive total potential in I = 0 channel.
Although the sum of OPE and TPE potentials of I = 1 channel in momentum space is attractive, we still can not conclude that there may exist a bound state in I = 1 channel, because the contribution from the contact term is independent on | q|, and not included in Fig. 9(a) and Fig. 10(a). The results in Fig. 9(b) and Fig. 10(b) are very interesting, because the inclusion of the contributions from FBCI totally reversed the sign of the results in momentum space. Although both the OPE and TPE potentials are attractive in I = 1 channel, the FBCI potential is largely repulsive, and thus the total potential of I = 1 channel in coordinate space is repulsive. For the I = 0 channel, we notice that much of potential generated from FBCI is canceled by the one of TPE approximately, and thus the residual part of FBCI potential plus the contribution from OPE play very important role, since both of them are attractive and thus the total potential is also attractive. Moreover, we find a bound state in the 0(1 + )BB * system by solving the Schrödinger equation numerically. We obtain the binding energy is ∆EBB * −12.6 +9.2 −12.9 MeV, and the corresponding root-mean-square radius is 1.0 +0.5 −0.2 fm. Analogously, the OBE model [73], quark model potential [39], and Lattice QCD [52] all demonstrated that 0(1 + )BB * system can form the molecular state. Our results are qualitatively consistent.
There are some other remarks on Fig. 10(b). One can notice that there exists a subtle cancellation between TPE and FBCI potentials in Fig. 10(b), and this phenomenon makes OPE alone be enough to bindB andB * . Although TPE is model independent, the LECs involved in the FBCI potential are determined with resonance saturation model, thus the cancellation in Fig. 10(b) is model dependent more or less. However, this result should be reasonable.BB * interactions could be related to those of BB * if we ignore the annihilation effect in BB * channel at very short distance [104]. The authors in Refs. [61,62,105] investigated DD * and BB * systems with χEFT extensively. They noticed that O( 0 ) contact-range interactions are crucial, and the inclusion of OPE only gives rise to minor modifications of the numerical results. As we have seen in Fig. 10(b), this is indeed the case, i.e., the attractive potential provided by the FBCI is stronger than that of OPE. The consistence tells that the values of the LECs we adopted in this work are reasonable, and the predictions given with moderate cutoff value Λ = 0.7 GeV are reliable. In the following, we talk about the Λ dependence of the total potentials.
We also need to stress that in Eq. (26), we introduce a Gauss regulator F (q) to regularize the potential V(q), and the function F (q) contains a cutoff parameter Λ. In theory, the binding energy should be independent of the regularization schemes adopted in Eq. (26), the scale dependence in Eq. (26) can be compensated by that of LECs. However, in practical calculations, it is difficult to remove the influences of Λ on observable. Therefore, we investigate the dependence of the total potentials for the 0(1 + )BB * system on cutoff parameter Λ, and the results are illustrated in Fig. 11. We can read that the binding in the short range would be deeper when the value of Λ is increased, but the results are not very sensitive to Λ. MeV, respectively. The result indicates that our predictions obtained at Λ = 700 MeV are stable, i.e., there exists a bound state in 0(1 + )BB * system when the cutoff is near by m ρ . Bound 0(1 + )BB * means this state is stable against strong interactions, it can not decay into its components due to the constraints of phase space. However, we notice that the mass of this state lies above the threshold of electromagnetic decay modeBBγ, thus it can be constructed via electromagnetic interactions at experiments.
C.B * B * system LikeBB system, the quantum numbers ofB * B * system must be constrainted by the selection rules given in Tab. I. With S wave, the physically allowedB * B * states are 1(0 + ), 1(2 + ) and 0(1 + ), respectively. The diagrams at the O( 0 ) are displayed in Fig. 12, and the diagrams at O( 2 ) are illustrated in Figs. 13-15. For the scattering processB * (p 1 )B * (p 2 ) →B * (p 3 )B * (p 4 ), we first list the amplitudes of the diagrams in Fig. 12 (we only give the amplitudes of I = 1 channel in this section, and the amplitudes of I = 0 channel can be found in Appendix C), where we define and The amplitudes of the diagrams in Fig. 13 are given as where the forms of the coefficients C x are The amplitudes of diagrams g 3.4 , h 3.4 , h 3.5 , and h 3.6 can be obtained by the relations The amplitudes of the diagrams in Fig. 14 are shown as follows, The amplitudes of the TPE diagrams in Fig. 15 read, . 15: Two-pion exchange diagrams ofB * B * system at O( 2 ). Notations same as in Fig. 2.
For succinctness, we also define Analogous to theBB * case (see Eq. (58)), the unlisted amplitudes of the diagrams R 3.1 ∼R 3.3 can be achieved by the following relations, To obtain the numerical results, some terms like ( i · j )( k · l ), ( i · j )(q · k )(q · l ) and (q · i )(q · j )(q · k )(q · l ) (where i , · · · , l denotes either the polarization vector of initial state or final state) appeared in Eqs. (61)-(88) must be reexpressed as a constant or a function of q. Resembling the OBE model [13], under S -wave interactions, the values of the term ( i · j )( k · l ) have been given in Tab. II. For the terms containing q, one can make the following substitutions, The 2PR contributions in diagrams h 3.1 , h 3.4 and B 3.1 must be removed, and other parameters are the same as we calculating theBB * system. The resulting potentials of I = 1 and I = 0B * B * system are depicted in Figs. 16-18, respectively. For the 1(0 + )B * B * system, from Fig. 16 we can see that the OPE potential is repulsive and the TPE potential is slightly attractive, but the contribution from FBCI is repulsive and largely dominant, which generates a fully repulsive potential for the 1(0 + )B * B * system. Thus no bound state can be found in the 1(0 + )B * B * channel in our calculations.
case. The line-shape of OPE potential is totally reversed, that is, the OPE potential is shallowly attractive for 1(2 + )B * B * system. Like the 1(0 + ) case, the potential resulted from FBCI is also repulsive and dominant. Thus the total potential of 1(2 + )B * B * system is also repulsive. So we can conclude that no bound state ofB * B * can exist in I = 1 channel.
At last, we analyze the potential of 0(1 + )B * B * system shown in Fig. 18. In rigorous heavy quark limit,B andB * mesons would be totally degenerate, so the potential ofBB * andB * B * systems should be similar if they carry the same quantum number. By comparing Fig. 18 and Fig. 10, we can see that the potentials of 0(1 + )BB * and 0(1 + )B * B * indeed have the same behaviors. The OPE and TPE potentials are attractive and repulsive for the 0(1 + )B * B * system, respectively, and the FBCI supplies a dominantly attractive potential, so the total potential is also attractive. By solving the schrödinger equation, we find a bound state in 0(1 + )B * B * channel with the binding energy ∆EB * B * −23.8 +16.3 −21.5 MeV, and the corresponding root-mean-square radius is 0.81 +0. 33 −0.13 fm.
We get the mass of 0(1 + )B * B * system also being above the threshold ofBBγγ. Therefore, it is detectable through electromagnetic interactions. In addition, we could resort to its weak decay modes as well, such as J/ψB + K 0 , with J/ψ and B + being fully reconstructable from J/ψ → + − and B + →D 0 π + .  Fig. 9.  Fig. 9.

D. The results in strict heavy quark limit
It is also very interesting to study the behaviors of effective potentials in strict heavy quark limit. As is discussed below Eq. (23), if ∆ = 0, the diagrams h 2.1 ∼h 2.4 and B 2.1 ∼B 2.3 in Fig. 6 and Fig. 8 would contribute to the 2PR parts of thē BB * scattering amplitudes. Similarly, forB * B * interactions, the diagrams h 3.1 ∼h 3.6 and B 3.1 ∼B 3.3 in Fig. 13 and Fig. 15 also contain both 2PR and 2PI contributions. In order to get the potentials, the contributions from the 2PR parts must be eliminated by using Weinberg's formalism (see Appendix A 3 for more details).
Under heavy quark limit, we find the main features of the potentials in different channels remain unchanged, so for simplicity, we only give the results of two representative states, i.e., 0(1 + )BB * andB * B * . The corresponding results are plotted in Fig. 19. Comparing the potentials with and without ∆ = 0 for 0(1 + )BB * andB * B * systems, we can find the corresponding potentials are only marginally shifted, this phenomenon indicates the heavy quark symmetry holds well for B mesons. Additionally, the effect caused by ∆ only happens at O( 2 ), and our calculations shown before have given the convergence of chiral corrections is very good. Furthermore, the mass difference ofB andB * lies far below the mass of pion m π , thus the influence of ∆ is largely suppressed. However, for D and D * systems, the mass shift ∆ > m π , one can expect more explicit effects of ∆ in DD * and D * D * systems.
B andB * would degenerate in heavy quark limit, so we anticipate the appearances of the potentials inBB * andB * B * systems should be identical if they have the same quantum number. This is vividly reflected from Fig. 19(a) and 19(b). We can see that the model independent parts, namely, OPE and TPE potentials are exactly the same for the two systems. At last, we list the binding energies and root-mean-square radii of 0(1 + )BB * andB * B * systems in the heavy quark limit, respectively. The result is in agreement with the patterns shown by cc and bb spectra, i.e., the binding would be deeper when the mass of the components is increased.

IV. ESTIMATION OF THE O( 2 ) LECS CONTRIBUTIONS
In Sec. III, we discard the contributions of O( 2 ) LECs since no enough data are available now to fit all of them. In this section, we pick these LECs up to see the influences on our numerical results. For this purpose, two different strategies are put to use. Although the estimations are rough and we are not trying to provide the very accurate LECs, they still can tell us some useful information about the reliability of the obtained results and how much the errors from the uncertainty of LECs at O( 2 ) are.

A. Strategy A
We first adopt the nonanalytic dominance approximation to give an estimation of LECs contributions [106,107], This approximation is based on the fact that the scattering matrix can be decomposed into analytic and nonanalytic part in chiral perturbation theory. The nonanalytic contributions can only be obtained from loop diagrams. However, the analytic terms can originate from both tree and loop graphs, and they are the polynomials of the expanding parameter . The analytic terms from loop diagrams can be absorbed by the LECs at the same order. Therefore, we can use the nonanalytic contributions to discuss the convergence of potentials since the O( 2 ) LECs cannot be fixed now.
In the nonanalytic dominance approximation, we assume there is a large cancellation between the analytic terms of loop diagrams and the finite parts of O( 2 ) LECs in Eqs. (6)- (8). With this approximation, we redo the calculation and give the binding energies of 0(1 + )BB * andB * B * with Λ = 0.7 GeV in the third column of Tab. III. The result given with the nonanalytic dominance approximation shows chiral expansion works well. However, one also needs to note that, in this strategy, we perhaps underestimate the effects of O( 2 ) LECs.

B. Strategy B
The O( 2 ) scattering amplitudes of 0(1 + )BB * andB * B * generated from the Lagrangians in Eq. (6) can be written as By dimensional analysis, we can naively get where Λ χ 1 GeV. Expanding the Lagrangians in Eq. (7) and Eq. (8), we notice the contributions from these two equations will be proportional to E 2 and k 2 (where E is the residual energy, and k is the residual momentum), respectively. In our calculations, the incoming and outgoing heavy mesons are on-shell, so we can assume the contributions of O( 2 ) tree diagrams mainly come from the ones governed by Eq. (6).
With the relations in Eq. (94) and the values given in Eq. (27), we can roughly estimate the numerical values of D h a , · · · , E h b to examine the influence of the LECs at O( 2 ). The corresponding results are listed in the last column of Tab. III.
In Tab. III, the difference between the results of "Strategy A" (or "Strategy B") and "No O( 2 ) LECs" can be regraded as an evaluation of the errors resulted from the uncertainty of LECs at O( 2 ). Although the two strategies used above are crude, the results in Tab. III show that the contributions of O( 2 ) LECs should be small. Neglecting the O( 2 ) tree diagrams is safe within the allowable range of errors at the early stage of the study. These LECs at O( 2 ) would be determined in a complete analysis in future when more experimental data are available.

V. SUMMARY
In summary, we have systematically investigated the intermeson interactions of double-beautyBB,BB * , andB * B * systems with HMχEFT. In addition to the O( 0 ) FBCI and OPE diagrams, we also include the O( 2 ) TPE diagrams and oneloop corrections to FBCI and OPE diagrams. The effective potentials are calculated with Weinberg's formalism [77,78], i.e., we do not calculate the scattering matrix directly since the 2PR contributions will spoil the correct power counting, and instead, we only take into account the 2PI parts of Feynman diagrams to derive the effective potentials. Moreover, with the aid of simple Gauss cutoff, we make the Fourier transformation to obtain the effective potentials with more visualized form in coordinate space, and then by iterating the potentials into Schrödinger equation, we not only can look for the bound state solutions but also can recover the 2PR contributions that we subtract before.
We only consider the S -waveB ( * )B( * ) systems in this work, and thus based on the selection rules in Tab. I, the physically allowed states are: 1(0 + )BB; 1(1 + ) and 0(1 + )BB * ; 1(0 + ), 1(2 + ), and 0(1 + )B * B * . The potentials of these six channels with different I(J P ) quantum numbers are studied and discussed in detail. We find the convergences of chiral corrections are all very good in these six channels, and the FBCI, which depicts the short-range interaction, plays the crucial role in determining the behaviors of the total potentials. For 1(0 + )BB andB * B * systems, the total potentials are both repulsive, and thus we conclude that no bound states can exist in these two channels. The same situation also holds for the 1(1 + )BB * and 1(2 + )B * B * systems.
While for the 0(1 + )BB * andB * B * systems, our results are very interesting because the total potentials of these two channels are both attractive. By solving the Schrödinger equation, we find two bound states in these two channels with the corresponding binding energy ∆EBB * −12.6 +9. Our result can qualitatively confirm the conclusions drawn by the OBE model [73], quark model potential [39], and Lattice QCD [52].
It will be an intriguing topic to search for the exotic doublebeauty molecular states in 0(1 + ) channels experimentally. These two molecular candidates cannot directly decompose into their components due to the constraints of phase space. However, theB * meson can decay intoBγ via electromagnetic interaction. We find the masses of 0(1 + )BB * andB * B * states both lie above the thresholds for decaying intoBBγ and BBγγ, and thus they are reconstructable in their electromagnetic decay modes.
Our results and above typical decay modes provide important information to future search in experiments. Additionally, the analytical chiral structures of the effective potentials between B mesons are very useful for the extrapolations of the B meson pair interactions in Lattice QCD calculations.

Calculations of J functions
For the loop integral J c x and J F x defined in Eq. (A1) and Eq. (A4), one can easily obtain their results with dimensional regularization [108]. Such as where Here, γ E is the Euler-Mascheroni constant 0.5772157. For J F 0 , just one more step, i.e., Feynman parameterization, is needed, where ∆ = x(x−1)q 2 +m 2 −iε. Performing a shift of integration variables l → l − xq so that there remain no terms linear in l in the denominator, and repeating the same step as in Eq. (A7), we can get Next, we outline the deductions of J a 0 (Eq. (A2)) which serves as a starting point for more complicated loop integrals, such as J T x . The calculations are slightly different with the integrals containing only pion propagators. We should use the following Feynman trick to combine the denominators [58], Setting A = v · l + ω + iε, B = l 2 − m 2 + iε, and n = 1, we get The terms linear in l and y in above equation can be eliminated via shifting the integration variables: (1) performing a shift l → l − yv, and making use of v 2 = 1, (2) shifting the integration variable y → y + ω. Finally, we obtain We should notice that the range of y integration now changes to [−ω, ∞), and we split this range into two parts, [−ω, 0) and [0, ∞). The first part can be easily obtained with dimensional regularization, which reads where∆ = y 2 + m 2 − ω 2 . For the second part, we first use β function to integrate out the y variables, then use dimensional regularization to integrate out l, and thus the result reads 1 8π Finally, we can get the full form of J a 0 , If one wants to obtain the results of J a 11 , J a 21 , J a 22 · · · , the shifts on integration variables l and y made above should be kept in mind, because the integration variables that appeared in the numerator also have to be shifted accordingly.
The calculations of J T x are very similar to J a x as illustrated above, the concrete procedures are: (1) Using Feynman's trick to combine the denominators of two pions propagators, then shifting l → l − xq, and making use of v · q = 0.
(2) Using Eq. (A11) to combine the denominators of heavy meson and the combined pion propagators.
(3) Then just repeating the subsequent procedures what we have done for calculating J a x . At last, we list the expressions of the used J functions in our calculations, and these functions are calculated numerically in this work.
where ∆ = y 2 + A, A = x(x − 1)q 2 + m 2 − ω 2 . The 2PR part will be recovered when the effective potential is inserted into Schrödinger equation, and thus this part has to be removed when we calculate the potentials. In the following part, we illustrate how to subtract the 2PR contributions by utilizing the Weinberg's formalism (one can find a more convenient approach in Ref. [109]), and take the J B 21 as an example.
where M 2 = x(x − 1)q 2 + m 2 , and we have used v · q = 0. The kinetic term ( p 2 − l 2 )/2M i (M i is the mass of heavy meson) is included in the denominators of heavy mesons propagators for avoiding the pinch singularity when ω = δ = 0 (cf. discussions in Ref. [78]). In Eq. (A31), the pion poles contribute, and there is also a contribution from the double heavy meson poles. The contribution from the latter one is enhanced in the 1/M B expansion and gives the result of the iterated O( 0 ) OPE diagram, which is just the 2PR part that we need to eliminate. In other words, only the contribution from the pion poles accounts for the effective potential. The pion poles are located at Picking out l α l β term in the second line of Eq. (A31) and substituting it with g αβ (l 2 0 − l 2 )/D, then we can get the scalar function J B 21 defined in Eq. (A6). Closing the l 0 contour integral in the lower half-plane with the pole of interest located at l 0 = l 2 + M 2 − iε, the 2PI part of J B 21 can be easily obtained by using residue theorem and dimensional regularization in D − 1 dimensions, the corresponding result is written as In the following, we write out the used 2PI parts of J h x and J B x functions in this calculations.
J h

2PI
= m 2 L + m 2 64π 2 1 + 2 ln The I = 0 amplitudes that come from Fig. 5 are: (B2) The I = 0 amplitudes that come from Fig. 6 are: Y g 2.1 I=0 = C g 2.1 g 2 f 2 ( 2 · * 4 ) J g Y g 2.7 I=0 = C g 2.7 g 2 f 2 ( 2 · * 4 ) J g 22 r (m π , E + ∆, E + ∆), (B7) where Note that the coefficients C x appeared here should not be confused with those for the I = 1 channels. The amplitudes of diagrams g 2.5 , g 2.6 , h 2.3 , and h 2.4 can be obtained by the relations The I = 0 amplitudes that come from Fig. 7 are: The amplitudes of diagrams c 2.2 , c 2.5 , and c 2.6 can be obtained by the relations The I = 0 amplitudes that come from Fig. 8 are: When q 0 = 0, the amplitudes of the diagrams R 2.1 ∼R 2.3 can be obtained by the relations The I = 0 amplitudes that come from Fig. 15 are: (m π , E, E + ∆, q), (C25) (C26)