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

TL;DR

This paper investigates the interactions of double-beauty meson systems using chiral effective field theory, calculating potentials up to one-loop order, and predicts possible bound states relevant for future experiments.

Contribution

It provides a detailed chiral effective field theory analysis of double-beauty meson interactions, including one-loop corrections and predictions of bound states.

Findings

Attractive potentials in specific channels lead to bound states.

Predicted binding energies are approximately -12.6 MeV and -23.8 MeV.

Results aid future experimental searches and lattice QCD extrapolations.

Abstract

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.

Tables

Table 1. Table 1: The physically allowed B ¯ ​ B ¯ ¯ 𝐵 ¯ 𝐵 \bar{B}\bar{B} , B ¯ ​ B ¯ ∗ ¯ 𝐵 superscript ¯ 𝐵 ∗ \bar{B}\bar{B}^{\ast} , and B ¯ ∗ ​ B ¯ ∗ superscript ¯ 𝐵 ∗ superscript ¯ 𝐵 ∗ \bar{B}^{\ast}\bar{B}^{\ast} states Barnes:1999hs , where I tot subscript 𝐼 tot I_{\text{tot}} and S tot subscript 𝑆 tot S_{\text{tot}} designates the total isospin and total spin of two B 𝐵 B mesons, respectively. Due to the constraints of symmetry, the quantum numbers of B ¯ ​ B ¯ ¯ 𝐵 ¯ 𝐵 \bar{B}\bar{B} and B ¯ ∗ ​ B ¯ ∗ superscript ¯ 𝐵 ∗ superscript ¯ 𝐵 ∗ \bar{B}^{\ast}\bar{B}^{\ast} systems must satisfy the the selection rule L + S tot + I tot + 2 ​ i = E ​ v ​ e ​ n ​ n ​ u ​ m ​ b ​ e ​ r 𝐿 subscript 𝑆 tot subscript 𝐼 tot 2 𝑖 𝐸 𝑣 𝑒 𝑛 𝑛 𝑢 𝑚 𝑏 𝑒 𝑟 L+S_{\text{tot}}+I_{\text{tot}}+2i=Even\ number , where L 𝐿 L is the orbital angular momentum between two B 𝐵 B mesons, and i 𝑖 i denotes the isospin of one B 𝐵 B meson, which is 1 / 2 1 2 1/2 .

System	Total isospin	Total spin
	$I_{\text{tot}}$	$S_{\text{tot}}=0$	$S_{\text{tot}}=1$	$S_{\text{tot}}=2$
$\overline{B}\overline{B}$	1	even $L$
	0	odd $L$
$\overline{B}{\overline{B}}^{\ast }$	1		all $L$
	0		all $L$
${\overline{B}}^{\ast }{\overline{B}}^{\ast }$	1	even $L$	odd $L$	even $L$
	0	odd $L$	even $L$	odd $L$

Table 2. Table 2: The values of the products of polarization vectors under S 𝑆 S -wave interactions Liu:2008tn . For example, in our calculations, the term ( ϵ 1 ⋅ ϵ 2 ) ​ ( ϵ 3 ∗ ⋅ ϵ 4 ∗ ) ⋅ subscript italic-ϵ 1 subscript italic-ϵ 2 ⋅ superscript subscript italic-ϵ 3 ∗ superscript subscript italic-ϵ 4 ∗ (\epsilon_{1}\cdot\epsilon_{2})(\epsilon_{3}^{\ast}\cdot\epsilon_{4}^{\ast}) should be replaced by 3 3 3 , 0 0 , and 0 0 if the total spin S tot subscript 𝑆 tot S_{\text{tot}} equals to 0 0 , 1 1 1 , and 2 2 2 , respectively.

Terms	$S_{\text{tot}}=0$	$S_{\text{tot}}=1$	$S_{\text{tot}}=2$
$\left({ϵ}_1\cdot {ϵ}_2\right)\left({ϵ}_3^{\ast }\cdot {ϵ}_4^{\ast }\right)$	3	0	0
$\left({ϵ}_1\cdot {ϵ}_3^{\ast }\right)\left({ϵ}_2\cdot {ϵ}_4^{\ast }\right)$	1	1	1
$\left({ϵ}_1\cdot {ϵ}_4^{\ast }\right)\left({ϵ}_2\cdot {ϵ}_3^{\ast }\right)$	1	-1	1

Table 3. Table 3: The binding energies of 0 ​ ( 1 + ) 0 superscript 1 0(1^{+}) B ¯ ​ B ¯ ∗ ¯ 𝐵 superscript ¯ 𝐵 ∗ \bar{B}\bar{B}^{\ast} and B ¯ ∗ ​ B ¯ ∗ superscript ¯ 𝐵 ∗ superscript ¯ 𝐵 ∗ \bar{B}^{\ast}\bar{B}^{\ast} states obtained with different strategies in units of MeV. ``No 𝒪 ​ ( ϵ 2 ) 𝒪 superscript italic-ϵ 2 \mathcal{O}(\epsilon^{2}) LECs" represents the results obtained without considering the contributions from the finite part of 𝒪 ​ ( ϵ 2 ) 𝒪 superscript italic-ϵ 2 \mathcal{O}(\epsilon^{2}) contact terms in Eqs. ( 6 )-( 8 ).

Binding energy	No $𝒪\left({ϵ}^2\right)$ LECs	Strategy A	Strategy B
$\Delta E_{\overline{B}{\overline{B}}^{\ast }}$	$-{12.6}_{-12.9}^{+9.2}$	$-{10.4}_{-9.7}^{+7.2}$	$-{15.9}_{-12.7}^{+9.7}$
$\Delta E_{{\overline{B}}^{\ast }{\overline{B}}^{\ast }}$	$-{23.8}_{-21.5}^{+16.3}$	$-{20.1}_{-20.0}^{+14.5}$	$-{28.2}_{-23.6}^{+18.6}$

Equations

$${\mathcal{L}}_{4H}^{(0)}$$

$${\mathcal{L}}_{H\varphi }^{(1)}$$

$$H$$

$$\bar{H}$$

$$P=\left(B^{-},{\bar{B}}^0\right),P_{\mu }^{\ast }=\left(B^{\ast -},{\bar{B}}^{\ast 0}\right).$$

$${\Gamma }_{\mu }=\frac{i}{2}\left[{\xi }^{†},{\partial }_{\mu }\xi \right],u_{\mu }=\frac{i}{2}\left\{{\xi }^{†},{\partial }_{\mu }\xi \right\},$$

$$\phi=\sqrt{2}\left(\begin{array}[]{cc}\frac{\pi^{0}}{\sqrt{2}}&\pi^{+}\\ \pi^{-}&-\frac{\pi^{0}}{\sqrt{2}}\end{array}\right).$$

$${\mathcal{L}}_{4H}^{(2,h)}$$

$${\mathcal{L}}_{4H}^{(2,v)}$$

$${\mathcal{L}}_{4H}^{(2,q)}$$

$${\tilde{\chi }}_{\pm }={\chi }_{\pm }-\frac{1}{2}\text{Tr}[{\chi }_{\pm }],{\chi }_{\pm }={\xi }^{†}\chi {\xi }^{†}\pm \xi \chi \xi ,\chi =m_{\pi }^2.$$

$$i\int \frac{d^D\ell }{(2\pi {)}^D}\frac{1}{(-{\ell }_0+i\epsilon )({\ell }_0+i\epsilon )({\ell }^2-m_{\pi }^2+i\epsilon )}.$$

$${\ell }_0=\pm \left(\mathcal{E}-{\vec{\ell }}^2/2M_N\right)\pm i\epsilon ,$$

$${\mathcal{Y}}_{I=1}^{(0)\text{FBCI}}=8(D_a+E_a).$$

$${\mathcal{Y}}^{g_{1.1}}$$

$${\mathcal{Y}}^{g_{1.2}}$$

$${\mathcal{Y}}^{h_{1.1}}$$

$${\mathcal{Y}}^{h_{1.2}}$$

$${\mathcal{Y}}^{z_{1.1}}$$

$${\mathcal{C}}^{g_{1.1}}$$

$${\mathcal{C}}^{h_{1.1}}$$

$$\{X{\}}_r=\underset{D\to 4}{\lim }\left(X-L\frac{\partial }{\partial L}X\right)+\frac{1}{16{\pi }^2}\underset{D\to 4}{\lim }\left(\frac{\partial }{\partial D}\frac{\partial }{\partial L}X\right),$$

$${\mathcal{Y}}^{F_{1.1}}$$

$${\mathcal{Y}}^{T_{1.1}}$$

$${\mathcal{Y}}^{T_{1.2}}$$

$${\mathcal{Y}}^{B_{1.1}}$$

$${\mathcal{Y}}^{R_{1.1}}$$

$$\mathcal{V}=-\frac{1}{4}{\mathcal{Y}}_{\text{2PI}},$$

$${\mathcal{Y}}_{I=1}^{(2)\text{FBCI}}=4(-0.18D_b+0.1E_a-0.08E_b).$$

$$\mathcal{V}(r)=\int \frac{d^{3}q}{(2\pi {)}^3}{e}^{-i\mathbf{q}\cdot \mathbf{r}}\mathcal{V}(\mathbf{q})\mathcal{F}(\mathbf{q}).$$

$$D_a$$

Full text

I introduction

Hunting for exotic multiquark states beyond the conventional meson and baryon configurations is a long-standing problem of QCD GellMann:1964nj ; Zweig:19xx ; Jaffe:1976ig . After the pioneering exotic state $X(3872)$ was discovered in 2003 by the Belle Collaboration Choi:2003ue , a new perspective of hadronic physics was opened, and plenty of exotic hadrons called $XYZ$ states were observed at experiments Chen:2016qju ; Guo:2017jvc . 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 $\chi_{c}^{\prime}(2P)$ charmonium, the mass of it would be about 78 MeV lower than the prediction of relativized quark model Godfrey:1985xj . Furthermore, the isospin violation process $X(3872)\to J/\psi\rho$ 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 Swanson:2003tb ; Wong:2003xk ; AlFiky:2005jd ; Fleming:2007rp ; Liu:2008fh ; Liu:2008tn , tetraquark state Maiani:2004vq ; Chen:2010ze , traditional axial vector charmonium Barnes:2003vb ; Kong:2006ni , Lattice QCD simulation Chiu:2006hd ; Padmanath:2015era , and other theoretical schemes Kang:2016jxw (for a review, see Ref. Chen:2016qju ). Among them, the shallowly bound molecular picture, i.e., the deuteron-like configuration, is the most popular one, because the mass of $X(3872)$ is very close to the threshold of $D^{0}\bar{D}^{\ast 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 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 Chen:2016qju ; Guo:2017jvc ; Esposito:2016noz , only two $Z_{b}$ states were reported in bottomonium spectrum Belle:2011aa . There is also no unanimous conclusion on what the $Z_{b}(10610)$ and $Z_{b}(10650)$ really are Chen:2016qju , but one promising interpretation about the two charged bottomonium-like states is molecular conjecture which are composed of $B\bar{B}^{\ast}$ and $B^{\ast}\bar{B}^{\ast}$, respectively Bondar:2011ev ; Sun:2011uh ; Dias:2014pva ; Ohkoda:2011vj ; Kang:2016ezb . 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 $B\bar{B}^{\ast}$ molecular state in the framework of one boson exchange (OBE) model Liu:2008fh , and later, the $B\bar{B}^{\ast}$ and $B^{\ast}\bar{B}^{\ast}$ systems were again studied with OBE model Sun:2011uh ; Dias:2014pva ; Ohkoda:2011vj ; Sun:2012zzd ; Zhao:2015mga . In Ref. Liu:2016kqx , the pion-mediated interaction associated with coupled-channel effect is exploited to study the $B\bar{B}^{\ast}$ and $B^{\ast}\bar{B}^{\ast}$ systems. And some other related works such as QCD sum rule Chen:2015ata , Bethe-Salpeter approach He:2014nya , coupled-channel model Coito:2016ads , and so on.

As mentioned above, up to now, the observed exotic states are mainly composed of $Q\bar{Q}q\bar{q}$ or $Q\bar{Q}qqq$ (where $Q$ denotes $b$ or $c$ quark, and $q$ stands for $u$ or $d$ quark) Chen:2016qju . 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 $\Xi_{cc}^{++}$ very recently Aaij:2017ueg , which triggered many discussions on whether the stable $QQ\bar{q}\bar{q}$ tetraquark states can exist in nature. For instance, in Ref. Karliner:2017qjm , Karliner et al predicted a doubly bottomed tetraquark $bb\bar{u}\bar{d}$ with $J^{P}=1^{+}$ at $10389\pm 12$ MeV, which lies far below the $B^{-}\bar{B}^{\ast}$ and electromagnetic decay $B^{-}\bar{B}^{0}\gamma$ thresholds. Eichten et al also found the double-beauty state with constituent quarks $bb\bar{u}\bar{d}$ is stable against strong interactions Eichten:2017ffp . Based on the quark-diquark symmetry and the input of the mass of $\Xi_{cc}^{++}$, Ref. Mehen:2017nrh indicated the existence of a stable double-bottom tetraquark state in $I=0$ channel.

Actually, the exploration on the existence of stable $QQ\bar{q}\bar{q}$ tetraquark states is always an intriguing topic. Many meaningful works, including quark model potential, QCD sum rule, color-magnetic interaction, and Lattice QCD, have been utilized to study $bb\bar{u}\bar{d}$ systems. Date back to 1985, in the framework of non-relativistic potential model, Ref. Zouzou:1986qh found the state with quark contents $bb\bar{u}\bar{d}$ to be bound. Later, Barnes et al Barnes:1999hs investigated $\bar{B}\bar{B}$, $\bar{B}\bar{B}^{\ast}$, and $\bar{B}^{\ast}\bar{B}^{\ast}$ intermeson interactions in quark model potential, and after solving two-meson Schrödinger equations they found $I=0$ $BB^{\ast}$ ($[b\bar{q}][b\bar{q}]$) channel is attractive (for other related works, see Refs. Vijande:2009kj ; Ebert:2007rn ; Zhang:2007mu ; Du:2012wp ; Luo:2017eub ). In Refs. Michael:1999nq ; Detmold:2007wk ; Bicudo:2012qt ; Brown:2012tm ; Bicudo:2015vta ; Bicudo:2015kna ; Francis:2016hui ; Bicudo:2016ooe , the potential between two $B$ mesons as a function of two static $b$ quark distance is calculated in quenched Lattice QCD (for a review, see Ref. Beane:2008dv ), and especially in Refs. Bicudo:2015kna ; Francis:2016hui ; Bicudo:2016ooe the binding energy region of $bb\bar{u}\bar{d}$ system is evaluated to be $20-200$ MeV, which is compatible with phenomenological models' predictions Karliner:2017qjm ; Eichten:2017ffp .

Inspired by the observation of doubly charmed baryon $\Xi_{cc}^{++}$ Aaij:2017ueg and many theoretical works as mentioned before, in this work we use chiral effective field theory ($\chi$EFT) to study the physically allowed $\bar{B}^{(\ast)}\bar{B}^{(\ast)}([b\bar{q}][b\bar{q}])$ (see Tab. 1), i.e., $\bar{B}\bar{B}$, $\bar{B}\bar{B}^{\ast}$, and $\bar{B}^{\ast}\bar{B}^{\ast}$ intermeson interactions, to see in which channel the potential is attractive and to search for the possible bound states.

$\chi$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 $\chi$EFT in $N$-$N$ systems, one can see some important works in Refs. Bernard:1995dp ; Epelbaum:2008ga ; Machleidt:2011zz ; Meissner:2015wva and the references therein, or see Ref. Scherer:2002tk for an introduction to $\chi$EFT. Here, we give a short review about the use of $\chi$EFT in $M$-$M$ systems, mainly focus on heavy meson sectors. In Ref. AlFiky:2005jd , 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^{0}\bar{D}^{\ast 0}$. Similarly, there are also an abundance of investigations on $D\bar{D}^{\ast}(B\bar{B}^{\ast})$ systems with $\chi$EFT, concerning the molecular assumptions of $X(3872)$ and two $Z_{b}$ states Fleming:2008yn ; Baru:2011rs ; Valderrama:2012jv ; Nieves:2012tt ; Meng:2014ota ; Baru:2015tfa ; Wang:2013kva ; Baru:2015nea ; Jansen:2015lha ; Braaten:2015tga ; Baru:2016iwj (or see Refs. Chen:2016qju ; Guo:2017jvc for a review). In addition to utilizing $\chi$EFT to study heavy meson-heavy antimeson systems, applying $\chi$EFT to heavy meson-heavy meson systems (such as $D^{(\ast)}D^{(\ast)}$, $D^{(\ast)}\bar{B}^{(\ast)}$ and $B^{(\ast)}B^{(\ast)}$) is also an interesting topic. In our two previous works, first in Ref. Liu:2012vd , we gave a tentative usage of $\chi$EFT in $\bar{B}\bar{B}$ system, and obtained the strong interaction potentials in momentum space. Then in Ref. Xu:2017tsr , we calculated $S$-wave $DD^{\ast}$ potential at one-loop level, noticed in the $0(1^{+})$ channel it is attractive, and obtained a bound state with the binding energy $\Delta E_{DD^{\ast}}=-15.6^{+14.7}_{-19.2}$ MeV. We also notice that in Ref. Abreu:2015jma , both open charm and bottom states ($D^{(\ast)}\bar{B}^{(\ast)}$) are studied with heavy meson chiral effective field theory (HM$\chi$EFT).

Based on Refs. Liu:2012vd ; Xu:2017tsr , we naturally extend our study to $\bar{B}^{(\ast)}\bar{B}^{(\ast)}$ systems, and for these kinds of states their typical quark configuration is $[b\bar{q}][b\bar{q}]$. 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 different $\bar{B}^{(\ast)}\bar{B}^{(\ast)}$ systems. According to the selection rule listed in Tab. 1, the $I(J^{P})$ numbers for the physically allowed states are: $1(0^{+})$ for $\bar{B}\bar{B}$; $1(1^{+})$ and $0(1^{+})$ for $\bar{B}\bar{B}^{\ast}$; $1(0^{+})$, $1(2^{+})$, and $0(1^{+})$ for $\bar{B}^{\ast}\bar{B}^{\ast}$. In order to get the interaction potentials for these different channels in coordinate space, we firstly calculate the amplitudes in momentum space with $SU(2)$ HM$\chi$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 $\Delta E$ in the attractive channels. In this way, we can predict the possible molecular states in $\bar{B}^{(\ast)}\bar{B}^{(\ast)}$ systems and make a comparison with other phenomenological models Barnes:1999hs ; Li:2012ss .

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.

II effective lagrangians and Weinberg's formalism

II.1 Effective Lagrangians

In the framework of HM$\chi$EFT, the scattering amplitudes can be expanded order by order with a small parameter $\epsilon=q/\Lambda_{\chi}$, where $q$ is either the momentum of Goldstone bosons or the residual momentum of heavy mesons, and $\Lambda_{\chi}$ represents either the chiral breaking scale or the mass of a heavy meson. Except for $\bar{B}\bar{B}$, the scattering amplitudes of both $\bar{B}\bar{B}^{\ast}$ and $\bar{B}^{\ast}\bar{B}^{\ast}$ channels at leading order $\mathcal{O}(\epsilon^{0})$ have the contributions from FBCI and OPE, which are described by the leading effective Lagrangians in Eq. (1) AlFiky:2005jd ; Valderrama:2012jv ; Liu:2012vd and Eq. (2) Burdman:1992gh ; Wise:1992hn ; Yan:1992gz , respectively.

[TABLE]

where $D_{a}$, $D_{b}$, $E_{a}$, and $E_{b}$ are four independent low energy constants (LECs), which are determined later. $\tau^{a}$ represents the pauli matrix.

[TABLE]

where $v=(1,\vec{0})$ is the four-velocity of heavy mesons, and $H$ denotes the degenerated $B$ and $B^{\ast}$ doublet in heavy quark limit, which can be expressed as:

[TABLE]

[TABLE]

The last term in Eq. (2) accounts for the mass shift of $B$ and $B^{\ast}$, which will not vanish in chiral limit, and $\Delta$ is the mass difference of $(B,B^{\ast})$ doublet. The tensor $\sigma^{\mu\nu}$ is defined as $i[\gamma_{\mu},\gamma_{\nu}]/2$. Besides, the chiral connection $\Gamma_{\mu}$ and axial vector current $u_{\mu}$ are illustrated as follows,

[TABLE]

where $\xi=\exp(i\phi/2f)$. $f$ is the pion decay constant, and the triplet pion field $\phi$ is defined as

[TABLE]

As indicated in Ref. Liu:2012vd , other forms of FBCI at the leading order with different Lorentz structures are not independent, such as $\textnormal{Tr}\left[H\gamma_{\mu}\bar{H}H\gamma_{\mu}\bar{H}\right]$, $\textnormal{Tr}\left[H\bar{H}\right]\textnormal{Tr}\left[H\bar{H}\right]$, and $\textnormal{Tr}\left[H\sigma_{\mu\nu}\bar{H}\right]\textnormal{Tr}\left[H\sigma^{\mu\nu}\bar{H}\right]$, which can be expressed as the linear combinations of the terms involved in Eq. (1). Some terms like $\textnormal{Tr}\left[H\gamma_{5}\bar{H}\right]\textnormal{Tr}\left[H\gamma_{5}\bar{H}\right]$ and $\textnormal{Tr}\left[H\gamma_{5}\tau^{a}\bar{H}\right]\textnormal{Tr}\left[H\gamma_{5}\tau_{a}\bar{H}\right]$ vanish in heavy quark limit.

At the next-to-leading order, i.e., $\mathcal{O}(\epsilon^{2})$, the scattering amplitudes can be decomposed into four parts, one-loop corrections to FBCI and OPE, TPE, and the tree diagrams governed by $\mathcal{O}(\epsilon^{2})$ Lagrangians. In our calculations, we introduce the $\mathcal{O}(\epsilon^{2})$ FBCI Lagrangians to renormalize the $\mathcal{O}(\epsilon^{2})$ loop diagrams, which read Liu:2012vd ; Xu:2017tsr ,

[TABLE]

[TABLE]

[TABLE]

where

[TABLE]

The large amounts of LECs appeared in Eqs. (6)-(8) contain both the finite and infinite parts, and the infinite parts can be used to cancel the divergences of the loop diagrams at $\mathcal{O}(\epsilon^{2})$. Nevertheless, the finite parts can not be fitted due to the lack of enough data right now, and thus we neglect the contributions from $\mathcal{O}(\epsilon^{2})$ tree diagrams that governed by the Lagrangians in Eqs. (6)-(8) in our computations for the moment. However, an estimation of the errors caused by the $\mathcal{O}(\epsilon^{2})$ tree graphes is presented in Sec. IV.

II.2 Weinberg's formalism

Before formally carrying out our calculations, we give a brief introduction to Weinberg's power counting scheme Weinberg:1990rz ; Weinberg:1991um since it is the core of the theoretical calculations 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

[TABLE]

Inspecting the above loop integral we can find that the integration over $\ell_{0}$ is ill-defined because it has poles above and below the real axis at $\ell_{0}=\pm i\varepsilon$ (which is always called as pinch singularity Weinberg:1991um ; Machleidt:2011zz ). 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

[TABLE]

where $\mathcal{E}=\vec{p}^{2}/2M_{N}$, and $M_{N}$ is the mass of the nucleon. As is pointed out by Weinberg Weinberg:1990rz ; Weinberg:1991um , after the $\ell_{0}$ integration being performed with residue theorem the contribution from the pole of nucleon is enhanced by a large factor $M_{N}/|\vec{p}|$, and this strong enhancement explicitly breaks the naive power counting with which the $\ell_{0}$ integral should have been of $\mathcal{O}(1/|\vec{p}|)$.

With the formalism given in Ref. Weinberg:1991um , 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$\chi$EFT. We follow the ideas in studying $N$-$N$ system Weinberg:1991um , 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 of $\bar{B}^{(\ast)}\bar{B}^{(\ast)}$ systems

III.1 $\bar{B}\bar{B}$ system

The $\bar{B}\bar{B}$ system has been studied in Ref. Liu:2012vd with HM$\chi$EFT, but in Ref. Liu:2012vd 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 in $\bar{B}\bar{B}$ system.

From Tab. 1, one can see that only $1(0^{+})$ $\bar{B}\bar{B}$ state survives under $S$-wave interaction, which corresponds to the $B^{-}(p_{1})B^{-}(p_{2})\to B^{-}(p_{3})B^{-}(p_{4})$ elastic scattering process. At the leading order $\mathcal{O}(\epsilon^{0})$, OPE is excluded and the amplitude only receive contribution from the FBCI, by expanding the Lagrangian in Eq. (1), one can easily get

[TABLE]

At $\mathcal{O}(\epsilon^{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),

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

where the coefficients $\mathcal{C}^{x}$ are flavor and LECs dependent, their concrete values are:

[TABLE]

Various scalar functions $\mathcal{J}_{x}^{y}$ are defined in Appendix A.1, $m_{\pi}$ is the mass of pion, $D$ is the dimension where the loop integral is performed and approaches 4 at last. $\mathcal{E}$ represents the residual energy of $\bar{B}$ and $\bar{B}^{\ast}$ mesons, and is defined as $\mathcal{E}=E_{\bar{B}^{(\ast)}}-M_{\bar{B}^{(\ast)}}$. $\{X\}_{r}$ denotes the finite part of $X$ Liu:2012vd ,

[TABLE]

which is equivalent to make use of the modified minimal subtraction ($\mathrm{\overline{MS}}$) scheme.

The amplitudes of the diagrams $F_{1.1}$, $T_{1.1}$, $T_{1.2}$, $B_{1.1}$, and $R_{1.1}$ in Fig. 3 are illustrated in Eqs. (19)-(23), respectively.

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

Here we need to stress that in rigorous heavy quark limit, i.e., when the mass difference $\Delta=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^{\ast}$ 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 $\mathcal{V}$ can be derived from the 2PI amplitude $\mathcal{Y}_{\text{2PI}}$ by the relation

[TABLE]

where the factor $-1/4$ comes from OBE model Sun:2012zzd ; Li:2012ss , which originates from Breit approximation.

The parameters used in the calculations are: $\Delta=0.045$ GeV, $m_{\pi}=0.139$ GeV, $f=0.092$ GeV, $g=0.52$ Patrignani:2016xqp ; Ohki:2008py ; Detmold:2012ge , and the renormalization scale $\lambda=4\pi f$. For simplicity, as in the OBE model Liu:2008fh , the transferred energy $q_{0}$ and residual energy $\mathcal{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,

[TABLE]

We can see the contribution of $D_{a}$ vanishes, while $D_{b}$ and $E_{b}$ emerge at $\mathcal{O}(\epsilon^{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 $|\vec{q}|$. Thus the $|\vec{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 $|\vec{q}|$ is varied from [math] to $300$ MeV, which indicates the TPE part supplies the attractive potential, and this potential would be deeper if $|\vec{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,

[TABLE]

We should note that, because the potential $\mathcal{V}(\bf q)$ is expanded as the power series of the transferred momentum $q$ in $\chi$EFT, so we need to regularize $\mathcal{V}(\bf q)$ to suppress the contributions from high momentum to avoid the ultraviolet divergence of the integral, which is the manifestation of the fact that $\chi$EFT is valid only for $q\ll\Lambda_{\chi}$. Different approaches in Refs. Fleming:1999ee ; Phillips:1997xu ; Frederico:1999ps ; Birse:2005um ; Beane:2001bc ; PavonValderrama:2005gu ; Nogga:2005hy ; PavonValderrama:2005wv ; PavonValderrama:2007nu ; Entem:2007jg ; Yang:2007hb ; Beane:2008bt ; Valderrama:2007ja ; Cordon:2009pj ; Gamermann:2009uq were developed to regularize/renormalize $\mathcal{V}(\bf q)$ non-perturbatively. In this work, as in Refs. Xu:2017tsr ; Ordonez:1995rz ; Epelbaum:1999dj ; Ren:2016jna , we adopt a simple Gauss regulator $\mathcal{F}(\mathbb{q})=\exp(-\mathbb{q}^{2n}/\Lambda^{2n})$, and we set $n=2$ as in Ref. Ren:2016jna . The value of the cutoff parameter $\Lambda$ is commonly below the mass of $\rho$ meson in the $N$-$N$ case Epelbaum:2014efa , so we use a moderate value $\Lambda=0.7$ GeV as adopted in Ref. Xu:2017tsr 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. Xu:2017tsr by exploiting the resonance saturation model Epelbaum:2001fm ; Ecker:1988te ; Du:2016tgp , which read (all in units of GeV*-2*),

[TABLE]

With the preparations above, the effective potential of $1(0^{+})$ $\bar{B}\bar{B}$ 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^{+})$ $\bar{B}\bar{B}$ system, which is in coincidence with the calculations of Lattice QCD Bicudo:2015kna .

III.2 $\bar{B}\bar{B}^{\ast}$ system

For the $S$-wave $\bar{B}\bar{B}^{\ast}$ system, there are two different isospin states, i.e., $1(1^{+})$ and $0(1^{+})$ $\bar{B}\bar{B}^{\ast}$ 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 process $\bar{B}(p_{1})\bar{B}^{\ast}(p_{2})\to\bar{B}(p_{3})\bar{B}^{\ast}(p_{4})$ from the diagrams $H_{2.1}$ and $X_{2.1}$ in Fig. 5,

[TABLE]

where the subscript $I$ denotes the isospin of the channel, $q$ is the transferred momentum carried by pion, $\epsilon_{2}$ and $\epsilon_{4}^{\ast}$ designate the polarization vectors of initial $\bar{B}^{\ast}$ and final $\bar{B}^{\ast}$, 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 $\bar{B}\bar{B}^{\ast}$ system in Fig. 6, similar to $\bar{B}\bar{B}$ case, we also write them out in the following,

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

where the values of $\mathcal{C}^{x}$ are

[TABLE]

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,

[TABLE]

The amplitudes of one-loop corrections to OPE diagram (Fig. 7) are illustrated as follows,

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

where $\Sigma(m_{\pi})=\left[m_{\pi}^{2}/(24\pi^{2}f^{2})\right]\ln(m_{\pi}^{2}/\lambda^{2})$ Scherer:2002tk , the amplitudes of diagrams $c_{2.2}$, $c_{2.5}$, and $c_{2.6}$ can be obtained by the relations

[TABLE]

Here we list the amplitudes of the diagrams F_{2.1}$$\sim$$B_{2.3} in Fig. 8. The amplitudes of the diagrams R_{2.1}$$\sim$$R_{2.3} can be related to the ones of B_{2.1}$$\sim$$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. Xu:2017tsr .

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

When $q_{0}=0$, the amplitude structures of diagrams R_{2.1}$$\sim$$R_{2.3} resemble the ones of B_{2.1}$$\sim$$B_{2.3}, correspondingly. Hence we can get the amplitudes of diagrams R_{2.1}$$\sim$$R_{2.3} by multiplying the proper coefficients to the amplitudes of diagrams B_{2.1}$$\sim$$B_{2.3}, and substituting the loop integral $\mathcal{J}_{x}^{B}$ with $\mathcal{J}_{x}^{R}$ accordingly. One can readily verify the following correlations,

[TABLE]

In order to get the numerical results, some procedures are the same as we calculate the $\bar{B}\bar{B}$ 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 Xu:2017tsr should be adopted. Since we only consider the $S$-wave interaction, the terms $(\epsilon_{2}\cdot\epsilon_{4}^{\ast})$ and $(q\cdot\epsilon_{2})(q\cdot\epsilon_{4}^{\ast})$ appeared in above equations should be replaced by Sun:2012zzd ; Li:2012ss ; Xu:2017tsr

[TABLE]

We first list the results of one-loop corrections to FBCI in $I=1$ and $I=0$ channel, respectively.

[TABLE]

By comparing with Eq. (28) and Eq. (138), the results again show the convergence of chiral correction is very good.

The effective potentials of $1(1^{+})$ and $0(1^{+})$ $\bar{B}\bar{B}^{\ast}$ 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 $\mathcal{O}(\epsilon^{0})$ OPE potential is attractive for both $I=1$ and $I=0$ channels, and the $\mathcal{O}(\epsilon^{0})$ OPE potential in $I=0$ channel is more attractive than that of $I=1$ channel. The $\mathcal{O}(\epsilon^{2})$ potential generated by one-loop corrections to OPE diagram is very small compared with the $\mathcal{O}(\epsilon^{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 $|\vec{q}|$ becomes larger. In $I=0$ channel, TPE potential is repulsive, and would be more repulsive when $|\vec{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 $|\vec{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^{+})$ $\bar{B}\bar{B}^{\ast}$ system by solving the Schrödinger equation numerically. We obtain the binding energy is $\Delta E_{\bar{B}\bar{B}^{\ast}}\simeq-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 Li:2012ss , quark model potential Barnes:1999hs , and Lattice QCD Bicudo:2016ooe all demonstrated that $0(1^{+})$ $\bar{B}\bar{B}^{\ast}$ 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 bind $\bar{B}$ and $\bar{B}^{\ast}$. 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. $\bar{B}\bar{B}^{\ast}$ interactions could be related to those of $B\bar{B}^{\ast}$ if we ignore the annihilation effect in $B\bar{B}^{\ast}$ channel at very short distance Hippchen:1991rr . The authors in Refs. Valderrama:2012jv ; Nieves:2012tt ; Guo:2013sya investigated $D\bar{D}^{\ast}$ and $B\bar{B}^{\ast}$ systems with $\chi$EFT extensively. They noticed that $\mathcal{O}(\epsilon^{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 $\Lambda=0.7$ GeV are reliable. In the following, we talk about the $\Lambda$ dependence of the total potentials.

We also need to stress that in Eq. (26), we introduce a Gauss regulator $\mathcal{F}(\mathbf{q})$ to regularize the potential $\mathcal{V}(\mathbf{q})$, and the function $\mathcal{F}(\mathbf{q})$ contains a cutoff parameter $\Lambda$. 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 $\Lambda$ on observable. Therefore, we investigate the dependence of the total potentials for the $0(1^{+})$ $\bar{B}\bar{B}^{\ast}$ system on cutoff parameter $\Lambda$, 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 $\Lambda$ is increased, but the results are not very sensitive to $\Lambda$. By solving the Schrödinger function, we find when $\Lambda=600$ MeV, 700 MeV, and 770 MeV ($m_{\rho}$), the corresponding binding energies are $\Delta E_{\bar{B}\bar{B}^{\ast}}\simeq-8.1$ MeV, $-12.6$ MeV, and $-15.6$ MeV, respectively. The result indicates that our predictions obtained at $\Lambda=700$ MeV are stable, i.e., there exists a bound state in $0(1^{+})$ $\bar{B}\bar{B}^{\ast}$ system when the cutoff is near by $m_{\rho}$.

Bound $0(1^{+})$ $\bar{B}\bar{B}^{\ast}$ 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 mode $\bar{B}\bar{B}\gamma$, thus it can be constructed via electromagnetic interactions at experiments.

III.3 $\bar{B}^{\ast}\bar{B}^{\ast}$ system

Like $\bar{B}\bar{B}$ system, the quantum numbers of $\bar{B}^{\ast}\bar{B}^{\ast}$ system must be constrainted by the selection rules given in Tab. 1. With $S$ wave, the physically allowed $\bar{B}^{\ast}\bar{B}^{\ast}$ states are $1(0^{+})$, $1(2^{+})$ and $0(1^{+})$, respectively. The diagrams at the $\mathcal{O}(\epsilon^{0})$ are displayed in Fig. 12, and the diagrams at $\mathcal{O}(\epsilon^{2})$ are illustrated in Figs. 13-15.

For the scattering process $\bar{B}^{\ast}(p_{1})\bar{B}^{\ast}(p_{2})\to\bar{B}^{\ast}(p_{3})\bar{B}^{\ast}(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),

[TABLE]

where we define

[TABLE]

and

[TABLE]

The amplitudes of the diagrams in Fig. 13 are given as

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

where the forms of the coefficients $\mathcal{C}^{x}$ are

[TABLE]

The amplitudes of diagrams $g_{3.4}$, $h_{3.4}$, $h_{3.5}$, and $h_{3.6}$ can be obtained by the relations

[TABLE]

The amplitudes of the diagrams in Fig. 14 are shown as follows,

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

The amplitudes of the TPE diagrams in Fig. 15 read,

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

For succinctness, we also define

[TABLE]

Analogous to the $\bar{B}\bar{B}^{\ast}$ case (see Eq. (III.2)), the unlisted amplitudes of the diagrams R_{3.1}$$\sim$$R_{3.3} can be achieved by the following relations,

[TABLE]

To obtain the numerical results, some terms like $(\epsilon_{i}\cdot\epsilon_{j})(\epsilon_{k}\cdot\epsilon_{l})$, $(\epsilon_{i}\cdot\epsilon_{j})(q\cdot\epsilon_{k})(q\cdot\epsilon_{l})$ and $(q\cdot\epsilon_{i})(q\cdot\epsilon_{j})(q\cdot\epsilon_{k})(q\cdot\epsilon_{l})$ (where $\epsilon_{i},\cdots,\epsilon_{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 $\vec{q}$. Resembling the OBE model Liu:2008tn , under $S$-wave interactions, the values of the term $(\epsilon_{i}\cdot\epsilon_{j})(\epsilon_{k}\cdot\epsilon_{l})$ have been given in Tab. 2. For the terms containing $q$, one can make the following substitutions,

[TABLE]

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 the $\bar{B}\bar{B}^{\ast}$ system.

The resulting potentials of $I=1$ and $I=0$ $\bar{B}^{\ast}\bar{B}^{\ast}$ system are depicted in Figs. 16-18, respectively.

For the $1(0^{+})$ $\bar{B}^{\ast}\bar{B}^{\ast}$ 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^{+})$ $\bar{B}^{\ast}\bar{B}^{\ast}$ system. Thus no bound state can be found in the $1(0^{+})$ $\bar{B}^{\ast}\bar{B}^{\ast}$ channel in our calculations.

Wh show the effective potential of $1(2^{+})$ $\bar{B}^{\ast}\bar{B}^{\ast}$ system in Fig. 17. The behavior of TPE potential is similar to the $1(0^{+})$ case. The line-shape of OPE potential is totally reversed, that is, the OPE potential is shallowly attractive for $1(2^{+})$ $\bar{B}^{\ast}\bar{B}^{\ast}$ system. Like the $1(0^{+})$ case, the potential resulted from FBCI is also repulsive and dominant. Thus the total potential of $1(2^{+})$ $\bar{B}^{\ast}\bar{B}^{\ast}$ system is also repulsive. So we can conclude that no bound state of $\bar{B}^{\ast}\bar{B}^{\ast}$ can exist in $I=1$ channel.

At last, we analyze the potential of $0(1^{+})$ $\bar{B}^{\ast}\bar{B}^{\ast}$ system shown in Fig. 18. In rigorous heavy quark limit, $\bar{B}$ and $\bar{B}^{\ast}$ mesons would be totally degenerate, so the potential of $\bar{B}\bar{B}^{\ast}$ and $\bar{B}^{\ast}\bar{B}^{\ast}$ 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^{+})$ $\bar{B}\bar{B}^{\ast}$ and $0(1^{+})$ $\bar{B}^{\ast}\bar{B}^{\ast}$ indeed have the same behaviors. The OPE and TPE potentials are attractive and repulsive for the $0(1^{+})$ $\bar{B}^{\ast}\bar{B}^{\ast}$ 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^{+})$ $\bar{B}^{\ast}\bar{B}^{\ast}$ channel with the binding energy $\Delta E_{\bar{B}^{\ast}\bar{B}^{\ast}}\simeq-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^{+})$ $\bar{B}^{\ast}\bar{B}^{\ast}$ system also being above the threshold of $\bar{B}\bar{B}\gamma\gamma$. Therefore, it is detectable through electromagnetic interactions. In addition, we could resort to its weak decay modes as well, such as $J/\psi B^{+}K^{0}$, with $J/\psi$ and $B^{+}$ being fully reconstructable from $J/\psi\to\ell^{+}\ell^{-}$ and $B^{+}\to\bar{D}^{0}\pi^{+}$.

III.4 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 $\Delta=0$, the diagrams h_{2.1}$$\sim$$h_{2.4} and B_{2.1}$$\sim$$B_{2.3} in Fig. 6 and Fig. 8 would contribute to the 2PR parts of the $\bar{B}\bar{B}^{\ast}$ scattering amplitudes. Similarly, for $\bar{B}^{\ast}\bar{B}^{\ast}$ interactions, the diagrams h_{3.1}$$\sim$$h_{3.6} and B_{3.1}$$\sim$$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^{+})$ $\bar{B}\bar{B}^{\ast}$ and $\bar{B}^{\ast}\bar{B}^{\ast}$. The corresponding results are plotted in Fig. 19. Comparing the potentials with and without $\Delta=0$ for $0(1^{+})$ $\bar{B}\bar{B}^{\ast}$ and $\bar{B}^{\ast}\bar{B}^{\ast}$ 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 $\Delta$ only happens at $\mathcal{O}(\epsilon^{2})$, and our calculations shown before have given the convergence of chiral corrections is very good. Furthermore, the mass difference of $\bar{B}$ and $\bar{B}^{\ast}$ lies far below the mass of pion $m_{\pi}$, thus the influence of $\Delta$ is largely suppressed. However, for $D$ and $D^{\ast}$ systems, the mass shift $\Delta>m_{\pi}$, one can expect more explicit effects of $\Delta$ in $DD^{\ast}$ and $D^{\ast}D^{\ast}$ systems.

$\bar{B}$ and $\bar{B}^{\ast}$ would degenerate in heavy quark limit, so we anticipate the appearances of the potentials in $\bar{B}\bar{B}^{\ast}$ and $\bar{B}^{\ast}\bar{B}^{\ast}$ 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. For the model dependent parts, i.e., contact interactions, at $\mathcal{O}(\epsilon^{0})$ the contributions for $\bar{B}\bar{B}^{\ast}$ and $\bar{B}^{\ast}\bar{B}^{\ast}$ are equal, and the only difference comes from the loop corrections. Because at $\mathcal{O}(\epsilon^{2})$, different contact structures with various combinations of LECs are involved into the loop diagrams.

At last, we list the binding energies and root-mean-square radii of $0(1^{+})$ $\bar{B}\bar{B}^{\ast}$ and $\bar{B}^{\ast}\bar{B}^{\ast}$ systems in the heavy quark limit, respectively.

[TABLE]

The result is in agreement with the patterns shown by $c\bar{c}$ and $b\bar{b}$ spectra, i.e., the binding would be deeper when the mass of the components is increased.

IV Estimation of the $\mathcal{O}(\epsilon^{2})$ LECs Contributions

In Sec. III, we discard the contributions of $\mathcal{O}(\epsilon^{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 $\mathcal{O}(\epsilon^{2})$ are.

IV.1 Strategy A

We first adopt the nonanalytic dominance approximation to give an estimation of LECs contributions Liu:2012uw ; Bijnens:1995yn , 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 $\epsilon$. 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 $\mathcal{O}(\epsilon^{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 $\mathcal{O}(\epsilon^{2})$ LECs in Eqs. (6)-(8). With this approximation, we redo the calculation and give the binding energies of $0(1^{+})$ $\bar{B}\bar{B}^{\ast}$ and $\bar{B}^{\ast}\bar{B}^{\ast}$ with $\Lambda=0.7$ GeV in the third column of Tab. 3. 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 $\mathcal{O}(\epsilon^{2})$ LECs.

IV.2 Strategy B

The $\mathcal{O}(\epsilon^{2})$ scattering amplitudes of $0(1^{+})$ $\bar{B}\bar{B}^{\ast}$ and $\bar{B}^{\ast}\bar{B}^{\ast}$ generated from the Lagrangians in Eq. (6) can be written as

[TABLE]

By dimensional analysis, we can naively get

[TABLE]

where $\Lambda_{\chi}\simeq 1$ GeV.

Expanding the Lagrangians in Eq. (7) and Eq. (8), we notice the contributions from these two equations will be proportional to $\mathcal{E}^{2}$ and $k^{2}$ (where $\mathcal{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 $\mathcal{O}(\epsilon^{2})$ tree diagrams mainly come from the ones governed by Eq. (6).

With the relations in Eq. (94) and the values given in Eq. (III.1), we can roughly estimate the numerical values of $D_{a}^{h}$, $\cdots$, $E_{b}^{h}$ to examine the influence of the LECs at $\mathcal{O}(\epsilon^{2})$. The corresponding results are listed in the last column of Tab. 3.

In Tab. 3, the difference between the results of Strategy A" (or Strategy B") and ``No $\mathcal{O}(\epsilon^{2})$ LECs" can be regraded as an evaluation of the errors resulted from the uncertainty of LECs at $\mathcal{O}(\epsilon^{2})$. Although the two strategies used above are crude, the results in Tab. 3 show that the contributions of $\mathcal{O}(\epsilon^{2})$ LECs should be small. Neglecting the $\mathcal{O}(\epsilon^{2})$ tree diagrams is safe within the allowable range of errors at the early stage of the study. These LECs at $\mathcal{O}(\epsilon^{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-beauty $\bar{B}\bar{B}$, $\bar{B}\bar{B}^{\ast}$, and $\bar{B}^{\ast}\bar{B}^{\ast}$ systems with HM$\chi$EFT. In addition to the $\mathcal{O}(\epsilon^{0})$ FBCI and OPE diagrams, we also include the $\mathcal{O}(\epsilon^{2})$ TPE diagrams and one-loop corrections to FBCI and OPE diagrams. The effective potentials are calculated with Weinberg's formalism Weinberg:1990rz ; Weinberg:1991um , 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$-wave $\bar{B}^{(*)}\bar{B}^{(*)}$ systems in this work, and thus based on the selection rules in Tab. 1, the physically allowed states are: $1(0^{+})$ $\bar{B}\bar{B}$; $1(1^{+})$ and $0(1^{+})$ $\bar{B}\bar{B}^{\ast}$; $1(0^{+})$, $1(2^{+})$, and $0(1^{+})$ $\bar{B}^{\ast}\bar{B}^{\ast}$. 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^{+})$ $\bar{B}\bar{B}$ and $\bar{B}^{\ast}\bar{B}^{\ast}$ 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^{+})$ $\bar{B}\bar{B}^{\ast}$ and $1(2^{+})$ $\bar{B}^{\ast}\bar{B}^{\ast}$ systems.

While for the $0(1^{+})$ $\bar{B}\bar{B}^{\ast}$ and $\bar{B}^{\ast}\bar{B}^{\ast}$ 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 $\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. We predict their masses to be

[TABLE]

Our result can qualitatively confirm the conclusions drawn by the OBE model Li:2012ss , quark model potential Barnes:1999hs , and Lattice QCD Bicudo:2016ooe .

It will be an intriguing topic to search for the exotic double-beauty 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, the $\bar{B}^{\ast}$ meson can decay into $\bar{B}\gamma$ via electromagnetic interaction. We find the masses of $0(1^{+})$ $\bar{B}\bar{B}^{\ast}$ and $\bar{B}^{\ast}\bar{B}^{\ast}$ states both lie above the thresholds for decaying into $\bar{B}\bar{B}\gamma$ and $\bar{B}\bar{B}\gamma\gamma$, 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.

Acknowledgments

B. W is very grateful to Prof. Shi-Lin Zhu for very helpful discussions. B. W also thanks Hao Xu and Rui Chen for some useful discussions. This project is supported by the National Natural Science Foundation of China under Grants No. 11175073 and No. 11705072. This work is also supported by the Fundamental Research Funds for the Central Universities. Xiang Liu is also supported in part by the National Program for Support of Top-notch Young Professionals.

Appendix A Loop integral

A.1 Definitions of $\mathcal{J}$ functions

The various $\mathcal{J}$ functions that appeared in the amplitudes of the previous sections are defined in the following, and they can be obtained by calculating the loop integrals in $D$ dimensions.

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

where the representation $X\vee Y\vee Z\vee\cdots$ denotes the symmetrized tensor structure of $X^{\alpha}Y^{\beta}Z^{\gamma}\cdots+\cdots$, and can be written as,

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

A.2 Calculations of $\mathcal{J}$ functions

For the loop integral $\mathcal{J}_{x}^{c}$ and $\mathcal{J}_{x}^{F}$ defined in Eq. (95) and Eq. (98), one can easily obtain their results with dimensional regularization Peskin:1995ev . Such as

[TABLE]

where

[TABLE]

Here, $\gamma_{E}$ is the Euler-Mascheroni constant $0.5772157$.

For $\mathcal{J}_{0}^{F}$, just one more step, i.e., Feynman parameterization, is needed,

[TABLE]

where $\Delta=x(x-1)q^{2}+m^{2}-i\varepsilon$. Performing a shift of integration variables $l\to l-xq$ so that there remain no terms linear in $l$ in the denominator, and repeating the same step as in Eq. (101), we can get

[TABLE]

Next, we outline the deductions of $\mathcal{J}_{0}^{a}$ (Eq. (96)) which serves as a starting point for more complicated loop integrals, such as $\mathcal{J}_{x}^{T}$. The calculations are slightly different with the integrals containing only pion propagators. We should use the following Feynman trick to combine the denominators Scherer:2002tk ,

[TABLE]

Setting $A=v\cdot l+\omega+i\varepsilon$, $B=l^{2}-m^{2}+i\varepsilon$, and $n=1$, we get

[TABLE]

The terms linear in $l$ and $y$ in above equation can be eliminated via shifting the integration variables: (1) performing a shift $l\to l-yv$, and making use of $v^{2}=1$, (2) shifting the integration variable $y\to y+\omega$. Finally, we obtain

[TABLE]

We should notice that the range of $y$ integration now changes to $[-\omega,\infty)$, and we split this range into two parts, $[-\omega,0)$ and $[0,\infty)$. The first part can be easily obtained with dimensional regularization, which reads

[TABLE]

where $\tilde{\Delta}=y^{2}+m^{2}-\omega^{2}$. For the second part, we first use $\beta$ function to integrate out the $y$ variables, then use dimensional regularization to integrate out $l$, and thus the result reads

[TABLE]

Finally, we can get the full form of $\mathcal{J}_{0}^{a}$,

[TABLE]

If one wants to obtain the results of $\mathcal{J}_{11}^{a}$, $\mathcal{J}_{21}^{a}$, $\mathcal{J}_{22}^{a}\cdots$, 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 $\mathcal{J}_{x}^{T}$ are very similar to $\mathcal{J}_{x}^{a}$ as illustrated above, the concrete procedures are:

(1) Using Feynman's trick to combine the denominators of two pions propagators, then shifting $l\to l-xq$, and making use of $v\cdot q=0$.

(2) Using Eq. (105) 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 $\mathcal{J}_{x}^{a}$.

At last, we list the expressions of the used $\mathcal{J}$ functions in our calculations, and these functions are calculated numerically in this work.

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

where $\Delta=y^{2}+A$, $A=x(x-1)q^{2}+m^{2}-\omega^{2}$.

[TABLE]

A.3 Removing the 2PR contributions from $\mathcal{J}_{x}^{h}$ and $\mathcal{J}_{x}^{B}$ functions

We can not directly calculate the full expressions of $\mathcal{J}_{x}^{h}$ and $\mathcal{J}_{x}^{B}$ functions as defined in Eq. (A.1) and Eq. (A.1), respectively, because they contain both 2PR and 2PI contributions. 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. Zhu:2004vw ), and take the $\mathcal{J}_{21}^{B}$ as an example.

[TABLE]

where $\mathcal{M}^{2}=x(x-1)q^{2}+m^{2}$, and we have used $v\cdot q=0$. The kinetic term $(\vec{p}^{2}-\vec{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 $\omega=\delta=0$ (cf. discussions in Ref. Weinberg:1991um ). In Eq. (A.3), 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 $\mathcal{O}(\epsilon^{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

[TABLE]

Picking out $l^{\alpha}l^{\beta}$ term in the second line of Eq. (A.3) and substituting it with $g^{\alpha\beta}(l_{0}^{2}-\vec{l}^{2})/D$, then we can get the scalar function $\mathcal{J}_{21}^{B}$ defined in Eq. (A.1). Closing the $l_{0}$ contour integral in the lower half-plane with the pole of interest located at $l_{0}=\sqrt{\vec{l}^{2}+\mathcal{M}^{2}}-i\varepsilon$, the 2PI part of $\mathcal{J}_{21}^{B}$ can be easily obtained by using residue theorem and dimensional regularization in $D-1$ dimensions, the corresponding result is written as

[TABLE]

In the following, we write out the used 2PI parts of $\mathcal{J}_{x}^{h}$ and $\mathcal{J}_{x}^{B}$ functions in this calculations.

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

Appendix B The amplitudes of $0(1^{+})$ $\bar{B}\bar{B}^{\ast}$ state

The $I=0$ amplitudes that come from Fig. 5 are:

[TABLE]

The $I=0$ amplitudes that come from Fig. 6 are:

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

where

[TABLE]

Note that the coefficients $\mathcal{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

[TABLE]

The $I=0$ amplitudes that come from Fig. 7 are:

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

The amplitudes of diagrams $c_{2.2}$, $c_{2.5}$, and $c_{2.6}$ can be obtained by the relations

[TABLE]

The $I=0$ amplitudes that come from Fig. 8 are:

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

When $q_{0}=0$, the amplitudes of the diagrams R_{2.1}$$\sim$$R_{2.3} can be obtained by the relations

[TABLE]

Appendix C The amplitudes of $0(1^{+})$ $\bar{B}^{\ast}\bar{B}^{\ast}$ state

The $I=0$ amplitudes that come from Fig. 12 are:

[TABLE]

The $I=0$ amplitudes that come from Fig. 13 are:

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

where

[TABLE]

and the amplitudes of diagrams $g_{3.4}$, $h_{3.4}$, $h_{3.5}$, and $h_{3.6}$ can be obtained by the relations

[TABLE]

The $I=0$ amplitudes that come from Fig. 14 are:

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

The $I=0$ amplitudes that come from Fig. 15 are:

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]