Tri-meson bound state $BBB^{\ast}$ via delocalized $\pi$~Bond

During the last decades, numerous exotic states which cannot be explained by the conventional quark model have been observed in experiment. Some of them can be understood as two-body hadronic molecules, such as the famous $X(3872)$, analogous to deuteron in nuclear physics. Along the same line, the existence of the triton leaves an open question whether there is a bound state formed by three hadrons. Since, for a given potential, a system with large reduced masses is more easier to form a bound state, we study the $BBB^{\ast}$ system with the one-pion exchange potential as an exploratory step by solving the three-body Schr\"odinger Equation. We predict that a tri-meson molecular state for the $BBB^{\ast}$ system is probably existent as long as the molecular states of its two-body subsystem $BB^*$ exist.

The same three bottomed meson system has been studied in Ref. [55] by calculating the scattering amplitudes between the Z b (10610) or the Z b (10650) and the bottomed meson. The universal bound states of three bottomed mesons from the Efimov effect has been ruled out. As in their calculation, only the contact interaction is included which might be the reason why they do not find a bound state. After including the long-range OPEP, the case might be different. Thus we solve the threebody Schrödinger equation to discuss the BBB * system by considering the OPEP. Without an assumption about its two-body subsystem, i.e. the molecular nature of the Z b (10610) or Z b (10650), we focus on the three-body bound state as a function of the binding energy of its subsystem. Hopefully, the present extensive investigations will be useful to deepen our the understanding of a system made off three heavy particles.
This paper is organized as follows. After the introduction, the formalism and the inputs for the BBB * system are presented in Sec. II. The dynamics of the two-body subsystem and the corresponding BO potential are given in Sec. III and Sec. IV, respectively. By constructing a proper interpolating wave functions for the BBB * in Sec. V, we solve the three-body Schrödinger equation in Sec. VI. Numerical results and discussions are given in the following section. The summary is presented in the last section. Some technicalities are relegated to the appendix.

II. FORMALISM AND THE INPUTS
The BOA has been successfully used in few-body system with several heavy and light particles [61,62]. For a three-body system with one light and two heavy mesons, such as the DD * K [60] system, the three-body Schrödinger equation is divided into two sub equations, one is the motion of light meson with two static sources and the other one is the equation for the two heavy mesons with the BO potential [60], which reflects the influence of the light meson on the dynamics of the two heavy mesons. If the interaction between the light meson and the heavy meson is attractive, it would make the two heavy mesons come closer, thus facilitating the formation of a bound state of the whole system. However, for a three heavy meson system, although the application of BOA is not straightforward, one can employ the underlying permutation symmetry, which means the corresponding dynamics is invariant under the interchange of any two constituents, for the system and continue to use the BOA for the calculation.
The OPEP indicates that there is only one virtual pion exchanged by any two constituents as shown in Fig. 1. One can use a, b and c to label the three mesons in the original channel, i.e. B * a B b B c . It changes into B a B * b B c via one-pion exchange (OPE) between a and b, and the channel B a B * b B c changes into B a B b B * c through the OPE between b and c. When the virtual pion arises between a and c, it returns back into the original channel B * a B b B c . Within this scenario, the virtual pion is not localized between any two constituents but rather shared by the whole system. It is very similar with the benzene molecule which has a pair of electrons shared by the six carbon atoms, which is called delocalized π bond in molecular physics. Since the three constituents have the same probability to be the B and B * mesons, one can write the system as B ( * ) a B ( * ) b B ( * ) c . Furthermore, the order of the a, b and c labels of the three mesons is artificial, as the system is invariant under the interchange of the a, b and c. This interchange symmetry will help to simplify our calculations. The point is that one can count the influence of each heavy meson on the dynamics of the other two mesons one by one. In other words, we can divide the system into three two-body subsystems a b, b c and a c. In each subsystem, one should add the BO potential from the remaining one. The existence of a negative common eigenvalue for the three subsystems may partly answer whether there is three-body bound state for the three heavy system. For simplicity, we call this method as Born − Oppenheimer potential method (BO potential method).
Before performing the calculation, we define the isospin wave functions of the BBB * systems as |I 2 , I 3 , I 3z with I 2 the isospin of the sub-BB * system. I 3 and I 3z represent the total isospin of the three-body system and its z direction, respectively. One thus Since BB * can couple with B * B * via OPE, the coupled channel effect is not negligible. We only consider the next close BB * B * channel in our calculation. If we distinguish the specific locations of the constituents as a, b and c, there are six channels in total, The Lagrangians with SU(2) chiral symmetry (we only consider OPE) and C-parity conservation read where the heavy flavor meson fields P and P * represent P = (B − ,B 0 ), P * = (B * − ,B * 0 ), respectively. Its corresponding heavy anti-meson fields P and P * represent P = (B + , B 0 ), P * = (B * + , B * 0 ). φ is the pion matrix We use the pion decay constant f π = 132 MeV [63]. The pionic coupling constant g = 0.57 is extracted from the width of D * + by assuming heavy quark flavor symmetry [64]. All the parameters and input datas are listed in Table I. Here, we neglect isospin breaking effect and use the masses of their charged particles. Under SU(2) chiral symmetry, the OPE interaction is of order O(p 0 ) for the three-body system. In this paper, we only take into account the OPEP to the O(p 0 ) order. Thus, there are four kinds of effective potentials. We use V 1 to denote the effective potential for the interaction BB * → B * B. V 2 and V 2 denote the process BB * → B * B * and its reverse, respectively. V 3 represents diagonal process B * B * → B * B * . Since the interactions are physical, the effective potentials should be unitary which gives V 2 = (V 2 ) † . r i j is used to denote the relative displacement between the i-th and j-th particles. Thus, the effective potentials of the three-body system in the channel space Its graphical illustration is shown in Fig. 2.
The solid and bold solid lines represent the B and B * meson fields, respectively. Dotted lines represent pion fields.

III. THE BREAK-UP STATE AND TWO-BODY SUBSYSTEM
For the three bottomed meson system, suppose one of the constituent is infinitely away from the remaining two mesons. The system can be divided into a two-body subsystem plus a free meson. A bound state solution of the two-body subsystem indicates a break-up state for the three-body system, i.e. a two-body bound state plus a free meson. In the OPE model, as there is no direct interaction between two B mesons, one could expect a break-up state with the subsystem BB * with quantum number J P = 1 + and a free meson B. We can detach the subsystem B a B * b first, and explore its binding solution. The Hamiltonian of the subsystem in the channel space where the T * = −(1/2µ * )∇ 2 ab and T * * = −(1/2µ * * )∇ 2 ab are the relative kinetic energy for the BB * and B * B * in their center-of-mass frame, respectively, with Here − − → L ab is the angular momentum operator between meson a and b. We also have the mass gap δM = M B * − M B . The effective potentials V 1 , V 2 , V 2 and V 3 depend on the isospin of the specific channels, thus we rewrite the BB * wave functions with fixed isospin For the specific channels Based on this, we can derive the scattering amplitude at the tree level where the T -matrix is the interaction part of the S -matrix and the M is defined as the invariant matrix element. In the second equation we have applied the first order of Born series expansion on the Lippmann-Schwinger equation with V f i being the effective potential. The relation between the scattering amplitude M f i and the potential V f i is where p f (i) and m f (i) denote the four-momentum and the mass of the final (initial) state. In the calculation, p 1 (E 1 , p) and p 2 (E 2 , − p) denote the four-momenta of the initial state particles in the center-of-mass system shown in Fig. 3, while p 1 (E 1 , p ) and p 2 (E 2 , − p ) denote the four-momenta of the final state particles, respectively.
is the transferred four-momentum. For convenience, we always use q = p 1 − p 1 and k = ( p 1 + p 1 )/2 instead of p and p in the calculations. The effective potential in coordinate space can be derived by Fourier transformation To take into account in a rough way the substructure of each vertex, a monopole form factor with m π the pion mass andΛ is used to suppress the contribution from UV energies. Here, ∆M = M * B − M B and ∆M = (M * B − M B )/2. As the parameter Λ is related to non-perturbative QCD, it cannot be well determined. Here we only explore its effect on the binding energy of the BB * with the quantum number J P = 1 + system. To solve the time-independent Schrödinger equation in coordinate space, the potential V( q, k) in momentum space can be transformed in to that in coordinate space as shown in the Appendix.  The isosinglet and isotriplet BB * potentials V BB * →BB * ( r) in coordinate space are shown in Fig. 4 (a) and (b), respectively, with Λ = 1440 MeV. The isosinglet potential V BB * →BB * ( r) is repulsive which does not indicate a bound solution. Nevertheless, the potential V BB * →B * B * ( r) for the isosinglet is attractive as shown in Fig. 5. So there is still the possibility of a binding solution. On the contrary, the isotriplet potential V BB * →BB * ( r) is attractive, while its potential V BB * →B * B * ( r) is repulsive. These potentials in coordinate space can be expressed as Table II. The c in Table II represents the C-parity of the corresponding channel. Here, c denotes the C-parity of the two-body system.
Since the tensor operatorŜ T leads to S-D wave mixing, the contributions from D-wave should be taken into account. Thus the wave function Ψ( r) has two parts with ψ S ( r) and ψ D ( r) the S -wave and D-wave functions, respectively. In the matrix method, we use Laguerre polynomials as a set of orthogonal basis with the normalization condition Thus the total wave function can be expanded as where φ S and φ D are the angular part of the spin and orbital wave function for the S -wave ( 3 S 1 ) and D-wave ( 3 D 1 ) states, respectively. a ( ) i and b ( ) i are the corresponding expansion coefficients of S-wave and D-wave, respectively. After solving the coupled-channel Schrödinger equation with the S-D wave mixing, we obtain the binding energy E b and its corresponding wave function Ψ(Λ, r ab ) for a given parameter Λ. Thus the wave function has the form Here, the wave function Ψ(Λ, r ab ) is normalized. If we choose the value of the parameter Λ = 1440 MeV for instance, one finds a loosely bound state for the isospin triplet system with a binding energy of 5.08 MeV, when the quantum number is J P = 1 + . There is also a loosely bound state for the isospin singlet system, when the quantum number is J P = 1 + . If the value of the parameter is chosen at Λ = 1107.7 MeV, the isospin singlet and triplet systems have the same binding energy of 5.08 MeV. The dependence of the binding energy on the parameter Λ will be given in Tables III-IV and discussed in Sec. VII.

IV. BORN-OPPENHEIMER POTENTIAL
As discussed in Sec. II, the BO potential reflects the influence of one of the mesons on the dynamics of the other two. For the BBB * (labeled as a, b and c) system, one can derive the BO potential from a for the bc system. The procedure is divided into the following three steps: • Considering that the particle b and c are static with the separation r bc , one can separate the degree of freedom of a from the three-body system.
• We assume the distance r bc is a parameter. The mesons b and c are static, and have one-pion interactions with meson a, which can be viewed as two static sources.
• We explore the dynamics for the meson a in the limit r bc → ∞, and subtract the binding energy for the break-up state which is trivial for the three-body bound state.
Within this scheme, we divide the motion of the system into two parts, one is the motion of the meson a relative to the mesons b and c. The other one is the relative motion between mesons b and c in the presence of the BO potential from a.
π ���π fixed fixed As illustrated in Figs. 6, we use r bc to denote the relative displacement between b and c. Further, r ab and r ac represent the displacement of the meson a relative to the meson b and c, respectively. One can separate the effective potentials for the meson a from the Eq. (4). The remaining part in Eq. (4) is the potential between bc. As discussed in the previous section, one can obtain the two-body binding energy The ψ( r ab ) and ψ ( r ab ) in the above equation are the eigenstate wave functions in Eq. (16).
In OPE model, as the virtual pion can only be exchanged between two of the BBB * subsystems, the wave function of a can be either 1 with pion exchanged between a and c. The final wave function for the meson a should be the superposition of these two components For simplicity, we neglect the mass difference for the BB * and B * B * in the kinetic operator, i.e. T * * ≈ T * . Then the Hamiltonian of the meson a is Accordingly, one can obtain the energy eigenvalue of the meson a E a (Λ, r bc ) = ψ( r ab , r ac )|H a |ψ( r ab , r ac ) where in the second step Eq. (19) and the symmetry between b and c are used. Since both the two-body energy eigenvalue E 2 and the wave functions ψ b , ψ c depend on the parameter Λ, E a is also a function of Λ. We take the parameter Λ = 1440 MeV as an example and plot E a for the isospin triplet of the BB * system in Figs. 7(a). As shown in the figure, the energy of the meson a has a minimum −17.65 MeV when r bc = 0, which corresponds to the limit that the mesons b and c are on top of each other and the system is reduced to the bc-a quasi-two-body system. When r bc → ∞, then E a tends to the two-body energy eigenvalue E 2 , i.e.−5.08 MeV. This corresponds to the situation that the meson b is infinitely far away from the meson c. Then the meson a can only form a two-body bound state with either b or c. It is not a three-body bound state anymore, but rather a two-body bound state plus a free meson state. In fact, this is nothing but the break-up state that we have discussed in the earlier sections. We also plot E a for the isospin triplet of the BB * system in Fig. 8(a), taking the parameter Λ = 1107.7 MeV. Similarly to the above, E a tends to the two-body energy eigenvalue −5.08 MeV. Therefore, we should subtract the limiting value E 2 when investigating the three-body bound state for the BBB * system. We define the BO potential as In other words, the BO potential between b and c is the energy eigenvalue of the meson a relative to that of the break-up state.

V. THE CONFIGURATIONS OF THE THREE-BODY SYSTEMS
In the OPE model, there is only one pion exchanged between any two constituents in the BBB * system. The constituents will change themselves from vector mesons into pseudo-scalar mesons or vice verse when they exchange one pion. Each constituent has the same probability to be a vector meson or a pseudo-scalar meson. Thus, the symbol * is shared among them. Since only one virtual pion occurs in the BBB * molecule, the virtual pion also be shared by the three mesons. We can thus write the BBB * as B ( * ) a B ( * ) b B ( * ) c . The BO potential can describe the contribution for the one meson on the dynamics of the two remaining mesons as we have discussed in the last section. Assuming that the meson b and c are much heavier than the meson a, then we can use the Born-Oppenheimer approximation to separate the degree of freedom of a from the three-body system. In other words, it is a kind of an adiabatic approximation that we divide the degrees of freedom of the three-body system into a light one and a heavy one. The motion of the light degree of freedom is the motion of meson a relative to the three-body centre of mass. The motion of the heavy degree of freedom is the relative motion between meson b and c. When exploring the dynamics for the meson a, we can assume the meson b and c are static with the distance r bc . Then the three-body system can be simplified as a two-body system consisting of mesons b and c but with an additional BO potential generated by the meson a. Overall, only the meson a can be separated from the system due to the fact that this meson is much lighter than the other ones. A separation in this way can be a good approximation for this system. With the same procedure that we derived Eq. (20), we obtain the wave functions ψ( r ab , r ac ) for the meson a. The remaining degree of freedom is the relative motion between meson b and c that can be described by a wave function assumed as Φ( r be ), to be determined from the Schrödinger equation. Then the total wave function of the system then has the form Ψ T = Φ( r be )ψ( r ab , r ac ).
Nevertheless, the true system B ( * ) a B ( * ) b B ( * ) c is that the three mesons have little mass difference. Every meson can be considered to be a lighter one and separated from the three-body system. Thus, the system has the three basic simplification schemes. That is we can divide the system B ( * ) a B ( * ) b B ( * ) c into three kinds of two-body subsystems, i.e., B ( * ) a B ( * ) b with the BO potential created by the meson c, B ( * ) b B ( * ) c with the BO potential created by the meson a and B ( * ) a B ( * ) c with the BO potential created by the meson b as shown in Figs. 9. These three simplification schemes can be regarded as three kinds of basic configurations. The eigenstates of the three-body system should be combinations of them. As the most simplest combination, one might expect the threebody eigenstate should be the superposition of the three kinds of basic configurations. We use the ψ / a , ψ / b , ψ / c to denote these three configurations. The configuration wave function ψ / a represent the configuration that we omit the meson B ( * ) a and add the corresponding BO potential instead. Similarly, ψ / b , ψ / c denote the configurations with the BO potentials provided by the mesons B ( * ) b and B ( * ) c , respectively. Taking the configuration function ψ / a as an example, we separate the motion of the B ( * ) a relative to the other mesons B ( * ) b and B ( * ) c where their relative displacement r bc is regarded as a parameter as shown in Fig. 9(a). The wave function of the B ( * ) a has been discussed in the last section and can be written as ψ( r ab , r ac ). The remaining degree of freedom is the relative motion between B ( * ) b and B ( * ) c , which can be taken as Φ( r bc ). Thus we have the configuration function ψ / a = Φ( r bc )ψ( r ab , r ac ). The other two wave functions ψ / b and ψ / c can be obtained analogously, i.e. ψ / b = Φ( r ac )ψ( r ab , r bc ), ψ / c = Φ( r ab )ψ( r bc , r ac ) which correspond to the Fig. 9(b) and Fig. 9(c), respectively. If we regard the three configuration functions as a set of basis states, then the basis constitute a configuration space {ψ / a , ψ / b , ψ / c }. The three-body eigenstate expressed as a superposition of the three kinds of basic configurations can be described as a state vector in this configuration space. Thus, as an interpolating wave function, the three-body wave functions can be written as Ψ T = αΦ( r bc )ψ( r ab , r ac ) + βΦ( r ac )ψ( r ab , r bc ) + γΦ( r ab )ψ( r bc , r ac ) where the Φ( r bc ), Φ( r ac ) and Φ( r ab ) are undetermined functions that need to be solved. The α, β and γ are the expansion coefficients. According to Eq. (20), we rewrite the three basic configuration functions in the channel space which can be expanded as a set of Laguerre polynomials ψ / a = i φ i ( r bc )ψ( r ab , r ac ), Here the subscript i is the order of Laguerre polynomials. We define the i th order of the configuration functions as ψ i / a = φ i ( r bc )ψ( r ab , r ac ), ψ i / b = φ i ( r ac )ψ( r ab , r bc ) and ψ i / c = φ i ( r ab )ψ( r bc , r ac ). Further, C is a normalisation constant. We expect the three-body bound state that we seek for can be expressed as a state vector in the configuration space {ψ / a , ψ / b , ψ / c }. However, the configuration functions in Eq. (23) where the x i j is a parameter matrix which will be determined later. The N i are normalization coefficients. The parameter matrix x i j in the three configuration functions are the same due to the interchange symmetry for the B ( * ) a B ( * ) b B ( * ) c system. Since the i th order configuration functionψ i / a should be orthogonal with the any order of the other configuration functionψ j / a , one can gets the orthogonalization condition which gives This equation will determine the parameter matrix x i j . Considering the normalization of the i th order configuration functionψ i one can obtain the normalization equation for the N i as After solving the equations for x i j and N i , we obtain an orthonormalized configuration basis. This basis constitutes a orthonormalized configuration space. Then the eigenvector for the three-body system B ( * ) a B ( * ) b B ( * ) c can be written as a vector in the configuration space {ψ / a ,ψ / b ,ψ / c }. Therefore, we have where theα i ,β i andγ i are the i th order expansion coefficients.

VI. THREE-BODY SCHRÖDINGER EQUATION
As we have discussed in previous sections, if the three-body binding energy is below the break-up threshold, the three-body system will disintegrate into a two-body system and a free meson. Since we only focus on the three-body bound state, we could make an energy shift and remove the energy eigenvalue E 2 for the break-up state and define a reduced Hamiltonian for the three-body system as The total Hamiltonian for the three-body system in the configuration space {ψ / a ,ψ / b ,ψ / c } can be written as

The explicit form of H is
There are two independent matrices φ i ac ψ bc |V cb 2 |φ j bc ψ ab + φ i ac ψ ab |V cb 1 |φ j bc ψ ac to be determined, where we have used the abbreviations φ i ab , φ i bc , φ i ac , ψ ab , ψ bc , ψ ac for φ( r ab ) i , φ( r bc ) i , φ( r ac ) i , ψ( r ab ), ψ( r bc ), ψ( r ac ) respectively. The expression for the H / c/ a can be easily obtained by the replacement c → b, b → c on the expression for the H / b/ a . Similarly, the expression for the H / c/ b is obtained by the replacement c → b, b → a, a → c on the expression for the H / b/ a . In fact, interchange invariance for the B ( * ) a B ( * ) b B ( * ) c system can simplify the calculation, i.e. H / c/ a = H / c/ b = H / b/ a and H / a/ a = H / b/ b = H / c/ c . Based on the above discussion, the three-body Schrödinger equation can finally be written as where the energy eigenvalue E 3 is the reduced three-body energy eigenvalue. The total energy eigenvalue relative to the BBB * mass threshold is E T = E 3 + E 2 . Solving the three-body Schrödinger equation may partly answer whether the three-body system has a loosely bound state or not.

VII. NUMERICAL RESULTS
As discussed in the above sections, there is only one free parameter Λ in the monopole form factor introduced in Sec. III, which reflects in a rough way the internal structure of the interacted hadrons. In other words, the size of hadron is proportional to 1/Λ, which is still unknown from the fundamental theory. For the deuteron case, the parameter Λ is within the range 0.8 ∼ 1.5 GeV. One would expect that the size of heavier bottom system is smaller than the size of deuteron, leading to a larger Λ. Thus, we vary the parameter Λ from 0.9 GeV to 1.6 GeV to study whether the BBB * system is bound or not.
In order to show the properties of the two-body interactions for the BB * , we first present the numerical results for the break-up state in Tables III-IV. We plot the effective potentials for the BB * in Fig. 4 and Fig. 5, where the regularization parameter is fixed at 1440 MeV. In these figures, (a) and (b) correspond to the isospin I = 0 and I = 1 cases, respectively. After carefully solving the coupled-channel Schrödinger equation with the treatment of the S-D wave mixing, we find loosely bound states for both cases.
For the isospin triplet case, i.e. I 2 = 1, the dependence of the binding energy of the two-body BB * system on the regularization parameter Λ is shown in Table III. The energy threshold of the break-up state for the BBB * is just the two-body energy eigenvalue of the BB * plus the mass of the three static free meson. We use E 2 denote the energy eigenvalue of the BB * . When the parameter Λ varies from 1380 MeV to 1560 MeV, there is a bound state solution with the binding energy 2.11 ∼ 15.59 MeV and the root-mean-square radius 1.51 ∼ 0.67 fm. The S -wave component takes over 99.13% ∼ 99.36% comparing to the value 0.87% ∼ 0.64% for D-wave. The BB * and B * B * channels have probabilities 91.91% ∼ 76.50% and 8.09% ∼ 23.50%, respectively. The proportion of the B * B * channel and D-wave component are relatively small. However, as the value of the regularization parameter Λ increases, the B and B * interacting with pion is more like a point particle, the proportion of the B * B * channel increases greatly, while the D-wave component decreases. We plot the radial wave functions of the S -wave and D-wave Fig. 10(a) for the system BB * , obviously, the bound state we have found is ground state. III: Bound state solutions of the BB * system with the isospin I 2 = 1. Λ is the parameter in the form factor. E 2 is the energy eigenvalue. The binding energy is −E 2 . r rms is the root-mean-square radius. α and β are the probabilities for the components BB * and B * B * , respectively. For the isospin singlet case, i.e. I 2 = 0, the dependence of the binding energy of the two-body BB * system on the regularization parameter Λ is shown in Table IV. We also use E 2 to denote the energy eigenvalue of the BB * . When the parameter varies from 1040 MeV to 1220 MeV, there is a bound state solution with binding energy 1.88 ∼ 14.78 MeV and the root-mean-square radius 2.03 ∼ 0.95 fm. The S -wave component is 86.09% ∼ 77.21% compared to the value 13.91% ∼ 22.79% for D-wave. The BB * and B * B * channels have probabilities 92.05% ∼ 75.09% and 7.95% ∼ 24.91%, respectively. The proportion of the B * B * channel and D-wave component are relatively small which is similar with the case of isospin triplet. As the value of the regularization parameter Λ increases, the proportion of the B * B * channel increases greatly. Different from the case of isospin triplet the S-wave component decrease and D-wave increase as Λ increases. As the parameter Λ increases, all of the effective potentials become stronger. The S-wave potential increases faster than the D-wave potential for the isospin triplet case, while it is reverse for the isospin singlet case. In order to check whether the bound state we have found is the ground state, we also plot the radial wave functions of the S -wave and D-wave Fig. 10(b) for the system BB * . TABLE IV: Bound state solutios of the BB * with the isospin I 2 = 0. Λ is the parameter in the form factor. E 2 is the energy eigenvalue. The binding energy is −E 2 . r rms is the root-mean-square radius. α and β are the probabilities for the component BB * and B * B * , respectively. In Sec. II, we have listed the isospin wave functions of the BBB * which are expressed as |I 2 , I 3 , I 3z . After solving the threebody Schrödinger equation via the method of Sec. V, we find that all of these isospin eigenstates have bound state solutions. As long as the two-body system BB * has a loosely bound state, the three-body system BBB * is most likely to have a loosely bound state, too. We have collected the dependence of the three-body bound state solutions on the two-body binding energy in Tables V-VI. The bound state solutions for the state |1, 3 2 , ± 1 2 (± 3 2 ) are shown in Table V. The three-body binding energy relative to their break-up states is 5.67 MeV, when the parameter Λ is chosen at 1440 MeV and the two-body binding energy of their subsystems BB * is 5.08 MeV. To search for the dependence on the binding energy of the two-body system E 2 , we change the parameter Λ. It turns out that if the value of the E 2 varies from −0.18 MeV to −17.97 MeV, then the reduced three-body energy eigenvalue E 3 decreases from −0.19 MeV to −18.99 MeV and the total three-body energy eigenvalue E T decreases from −0.38 MeV to −36.95 MeV. The structure of the three-body bound state is a regular triangle with the root-mean-square length of one side decreasing from 3.98 fm to 0.65 fm. In order to illustrate the strength of the BO potential, we also collect its minimum V BO (0) in the table within the range of -3.43∼-37.24 MeV as the E 2 increases. As E 2 increases, the effective attraction between B and B * becomes stronger, the BO potential becomed deeper, so then the three-body system becomes tighter and has a larger binding energy. From the results in the table, we can also see that the dominant wave between any two B ( * ) in the BBB * is S wave and the dominant channel is the BBB * instead of the BB * B * channel. For comparison, we plot the wave functions for any two B ( * ) in the BBB * system and that for the two-body BB * system in Fig. 10(a) with Λ = 1440 MeV. The shapes of these exhibit little difference. From another perspective, one more B meson has little effect on the size of the system but greatly increases the binding energy.
For the state |0, 1 2 , ± 1 2 , we also find a loosely bound solution, which is shown in Table VI. The three-body binding energy relative to their break-up states is 7.18 MeV, when the parameter Λ is chosen at 1107.7 MeV and the two-body binding energy of their subsystems BB * is 5.08 MeV. In order to show the dependence on the binding energy of the two-body system E 2 , we also change the parameter Λ. We find that if the value of the E 2 varies from −0.19 MeV to −17.10 MeV, then the reduced three-body energy eigenvalue E 3 decreases from −0.32 MeV to −20.96 MeV and the total three-body energy eigenvalue E T decreases from −0.51 MeV to −38.06 MeV. The structure of the three-body bound state is a regular triangle with the root-mean-square length of one side decreasing from 3.89 fm to 0.93 fm. As an illustration for the strength of the BO potential, we also list its minimum V BO (0) in the table within the range of −2.15∼ −30.15 MeV as E 2 increases. Similar to the I 3 = 3 2 case, the dominant wave between any two B ( * ) in the BBB * is S wave and the dominant channel is the BBB * one. In order to show that one more B meson has little effect on the size of the system, we also plot the wave functions for any two B ( * ) in the BBB * system and that for the two-body BB * system in Fig. 10(b) with Λ = 1107.7 MeV. Here we chose the parameter Λ = 1107.70 MeV for a better comparison with the |1, 3 2 , ± 1 2 (± 3 2 ) case, since both cases have the same two-body binding energy 5.08 MeV. Bound state solutions of the BBB * with the isospin I 3 = 3/2. E 2 is the energy eigenvalue of its subsystem BB * with the isospin I 2 = 1. E 3 is the reduced three-body energy eigenvalue relative to the break-up state of the BBB * system. E T is the total three-body energy eigenvalue relative to the BBB * threshold. V BO (0) is the minimum of the BO potential. r rms represents the root-mean-square radius of any two B in the BBB * system. The S -wave and D-wave represent the probabilities for S -wave and D-wave components in any two B in the BBB * . The α and β denote the probabilities for the BBB * and BB * B * components, respectively.  is the energy eigenvalue of its subsystem BB * with the isospin I = 0. E 3 is the reduced three-body energy eigenvalue relative to the break-up state of the BBB * system. E T is the total three-body energy eigenvalue relative to the BBB * threshold. V BO (0) is the minimum of the BO potential. r rms represents the root-mean-square radius of any two B in the BBB * system. The S -wave and D-wave represent the probabilities for S -wave and D-wave components in any two B in the BBB * . The α and β denote the probabilities for the BBB * and BB * B * components, respectively. The state |1, 1 2 , ± 1 2 also has a loosely bound solution. However, in our calculation the states |1, 1 2 , ± 1 2 and |1, 3 2 , ± 1 2 (± 3 2 ) are degenerate. This is due to the fact that in the OPE model only two-body interactions are considered, and the two-body interaction is only depends on the total isospin of the two interacting mesons. The states |1, 1 2 , ± 1 2 and |1, 3 2 , ± 1 2 (± 3 2 ) have the same two-body interaction but may have different three-body interaction. If we further consider the calculation to the next-tonext leading order, this degeneracy may disappear. The calculation that contains three-body interactions via pion exchange is quite complicated, which is left for the further work.
The numerical results show that the three-body binding energy |E 3 | increases as the two-body binding energy |E 2 | increases. One may wonder whether there is a critical value of the |E 2 |, below which the three-body system has no bound state solution. After lots of calculations, it turns out that all of the isospin eigenstates have no such critical value. That is to say, no matter what a small value for the two-body binding energy |E 2 |, as long as the two-body system BB * has a loosely bound state, the three-body system BBB * probably has a loosely bound state. To show this conclusion explicitly, we plot the dependence of the three-body binding energy on the variety of two-body binding energy in Fig. 11. The isospin eigenstate |1, 3 2 , ± 1 2 (± 3 2 ) and |1, 1 2 , ± 1 2 cases are shown in Fig. 11(a), where E BBB * I=3/2 and E BB * I=1 denote the reduced three-body binding energy and two-body binding energy, respectively. When the two-body binding energy E BB * I=1 approaches 0 MeV, the reduced three-body binding energy E BBB * I=3/2 approaches a small value of about 0.06 MeV. Similarly, we also plot the dependence curve for the |0, 1 2 , ± 1 2 case in Fig. 11(b), where E BBB * I=1/2 and E BB * I=0 denote the reduced three-body binding energy and two-body binding energy in this case, respectively. As shown in the figure, the reduced three-body binding energy E BBB * I=1/2 also has a small value of about 0.12 MeV, when the two-body binding energy E BB * I=0 approaches zero.
FIG. 10: Plot of various wave functions. The blue lines represent the wave functions for any two constituents in the BBB * . The red lines denote the wave functions for its subsystem BB * . (a) corresponds to the isospin states |1, 3 2 , ± 1 2 (± 3 2 ) and |1, 1 2 , ± 1 2 cases. (b) corresponds to the isospin state |0, 1 2 , ± 1 2 case. Here we chose the parameter Λ = 1440 MeV in t (a) and Λ = 1107.7 MeV in (b) for a better comparison of all the cases, since they have the same two-body binding energy of 5.08 MeV.  11: Dependence of the reduced three-body binding energy on the two-body binding energy of its subsystem BB * . The red point is the critical point which indicates the lower limit of the required binding energy of the isotriplet BB * to form a three-body bound state. (a) corresponds to the isospin states |1, 3 2 , ± 1 2 (± 3 2 ) and |1, 1 2 , ± 1 2 cases, while (b) corresponds to the isospin state |0, 1 2 , ± 1 2 case.

VIII. SUMMARY AND DISCUSSION
In the present paper, we have performed an extensive study on the possibility of the BBB * system to form tri-meson molecules. Based on the Born-Oppenheimer potential method as well as the OPE scheme, we derived the three-body Schrödinger equation for the system BBB * . Since the regularization parameter Λ is hard to be determined, we choose the parameter in the range of 0.9∼1.6 GeV and show various bound state solutions of the BBB * system. After careful treatments of the S-D wave mixing and the coupled-channel BB * B * , we found all of the isospin eigenstates expressed by the |I 2 , I 3 , I 3z have bound state solutions. For instance, in the states |1, 3 2 , ± 1 2 (± 3 2 ) and |1, 1 2 , ± 1 2 , the three-body binding energy relative to their break-up states is 5.67 MeV, when the parameter Λ is chosen at 1440 MeV and the two-body binding energy of their subsystems BB * is 5.08 MeV. In the state |0, 1 2 , ± 1 2 , the three-body binding energy relative to their break-up states is 7.18 MeV, when the parameter Λ is chosen at 1107.7 MeV and the two-body binding energy of their subsystems BB * is also 5.08 MeV. After careful calculations, we find no critical value for the two-body binding energy, which indicates the lower limit of the required binding energy of their subsystem BB * to form a three-body bound state. That is to say, no matter however small the two-body binding energy is, as long as the two-body subsystem BB * has a loosely bound state, the three-body system BBB * is most likely to have a loosely bound state.
The BO potential method we have used in this paper is an adiabatic approximation that we divide the degrees of freedom of the motion for the three-body system into a light one and a heavy one. Then we simplify the three-body system into a two-body system only with heavy degrees of freedom but with an additional BO potential generated by the relative light meson. Since the system B ( * ) a B ( * ) b B ( * ) c has little mass difference on its constituents, the motion of every constituent can be regarded as a light degree of freedom. Therefore, the eigenstates of the three-body system should be the combinations of all of the possible cases. As the most simplest combination, one might expect the three-body eigenstate should be a superposition all of the possible cases. In other words, the three-body bound state solutions we have listed in the last section are approximate solutions. Maybe the strict solutions will be a more complicated combinations. To answer this question requires further research.
Our calculations are based on the OPE scheme, which leading order in the chiral power counting (negelcting contact interactions). Since only one virtual pion occurs in the BBB * molecule, the virtual pion is also shared by the three mesons. Therefore, the three constituents in the BBB * system share one virtual pion which corresponds to a delocalized pion bond. It is attractive and strong enough to make them form a three-body molecular state.
As a short summary, with the delicate efforts of the long-range one-pion exchange, the S-D wave mixing and coupled-channel effects, we have investigated the existence of the loosely bound tri-meson molecules BBB * and find that it is very easy to form a tri-meson molecular state as long as its two-body subsystem BB * has a molecular state. Hopefully, the present extensive investigations will be useful to the understanding of the few-body hadronic systems and the future well-developed experiments on hadron collisions will provide us with a platform to seek out the tri-meson molecules.