Possible molecular states in $B^{(*)}B^{(*)}$ scatterings

We present that, if unitarizing the BB scattering amplitudes in the constituent interchange model, one can find two bound state poles for (Itot, Stot) = (0, 1) BB ∗ and BB system, which corresponds to two I(J ) = 0(1) doubly bottomed molecular states. Furthermore, it is noticed that the virtual states in (1, 0) BB, (1, 1) BB, (1, 0) BB, and (1, 2) BB systems could produce enhancements of the module squares of the scattering T -matrix just above the related thresholds, which might correspond to I(J ) = 1(0), 1(1), and 1(2) doubly bottomed molecular states, respectively. The calculation may be helpful for searching for the doubly bottomed molecular state in future experiments.


I. INTRODUCTION
Since the X(3872) was observed by Belle in 2003 [1], searches for exotic multiquark states beyond the conventional meson classifications have attracted intense attentions both from experimental and theoretical sides. In these years, many unconventional hidden-charm or hidden-beauty states have been observed, as reviewed in refs. [2][3][4]. Some of these states run out of the predictions of the quark potential model and can not be described by the naive quark model. As a result, new configurations such as hadron molecules, tetraquark states, hybrid states, are utilized to understand these exotic states. Many of these states are regarded as belonging to the valence quark configuration QQqq (where Q denotes c or b quarks and q denotes light quarks), which implies that QCD can form hadron states in an unconventional way [5].
After the observations of so many hidden-charmed or hidden-bottomed states, there arises a natural question whether there exist similar doubly-charmed or doubly-bottomed molecular states with valence quark configuration QqQq. Recently, LHCb reported the doubly charmed baryon state Ξ ++ cc [6], which inspired Karliner and Rosner to predict a doubly bottomed tetraquark state to be about 215 MeV below the B −B0 * threshold soon later [7], while Eichten and Quig predicted it to be 121 MeV below the B −B0 * threshold [8]. In fact, theoretical explorations of such doubly-heavy meson states have been carried on by several groups for a long time, as was reviewed in ref. [9]. In 1986, it was found that the bbūd state could be bounded in a non-relativistic potential model [10]. In 1999, Barnes et. al. found that the I = 0 BB * channel could be attractive and form a J P = 1 + bound state by solving the schrodinger equation of the BB, BB * , and B * B * systems [11]. After that, different approaches to investigate the possibility of forming doubly bottomed or doubly charmed meson states are studied through tetraquark models [12][13][14][15][16][17][18][19][20][21], meson exchange models [22][23][24][25], which usually calculate the binding energies by solving the Schrödinger Equation with the potentials obtained in such models. There are also some other calculations based on the QCD sum rules [26,27] or using the Lattice simulations [28][29][30][31][32][33]. These calculations basically paid their attentions to whether the bound states could be formed in the particular channels and the results in different calculations are not the same though most of them claim that there is a bound state in the I = 0 BB * channel with very high possibility. However, there could be cases where the near threshold structure in the lineshape could be caused by a virtial state or by a combined effect of the virtual state and the threshold, when the virtual state is very close to the threshold. In fact, there are models claiming that the famous X(3872) could be a virtual state just below the threshold[]. However, the method of solving the Schrödinger equation or using lattice calculations could only produce the bound state and would not give the virtual states. To study the possibility of the near threshold virtual states, one needs to have a nonperturbative scattering amplitude and annalytically continue it to the complex plane and study the pole structure of the amplitude.
In this paper, we just try to study the binding problem of BB, BB * , B * B * systems by investigating the existence of poles in the unitary meson-meson scattering matrix. As is well known, a bound state will appear as a pole of the partial-wave S matrix below the threshold on the physical Riemann sheet, while a virtual state will appear as a pole below the threshold on the second Riemann sheet. The unitary amplitude is obtained using the K-matrix method by unitarization of a Born approximation of the amplitude. The valence quark interchange model by Barnes and Swanson is adopted to provide the Born approximation of the meson-meson scattering amplitude [11,34]. Then, the unitarization of partial-wave amplitudes are derived and the poles of unitarized scattering amplitudes could be extracted by analytically continuing the energy to the complex plane. All the S-wave partial wave amplitudes of BB, BB * , B * B * scatterings are analyzed and related bound-state or virtual-state poles are searched for.
The paper is organized as follows: The Barnes-Swanson model is introduced and discussed in Section II. The K-matrix unitarization method for the partial wave Barnes-Swanson amplitude is derived in Section III. Numerical results and discussions are devoted in Section IV.

II. THE MODEL
In the constituent interchange model developed by Barnes and Swanson [34,35], the meson-meson scattering amplitude is calculated by the (anti)quark-(anti)quark interactions by assuming the one-gluon-exchange (OGE) color Coulomb interaction, spin-spin interaction, and linear scalar confinement interaction. In the coordinate space, the effective interaction Hamiltonian is where i, j represent valence quark or anti-quark in different initial hadrons. The color generator F a = λ a /2 for quarks and F a = −λ aT /2 for anti-quarks. Actually, to make the scheme consistent, the quark-quark interactions have also been used in determining the meson spectroscopy and their wave functions, so the model has a small free parameter space. It is more convenient to construct the scattering amplitude in the momentum space. First, the meson state is defined by the mock state to represent its wave function as 12 and ω 12 are the spin wavefunction, flavor wave function and the color wave function, respectively. p 1 (p 2 ) and m 1 (m 2 ) are the momentum and mass of the quark (anti-quark). P = p 1 + p 2 is the momentum of the center of mass, and p = m2 p1−m1 p2 m1+m2 is the relative momentum. ψ nLML is the wave function for the meson, n being the radial quantum number. The normalization is P ′ , λ ′ | P , λ = δ λλ ′ δ( P ′ − P ), where λ represents the quantum numbers such as nJM and particle species. Then, at the lowest order the scattering amplitude for the process AB → CD is just the matrix element of interaction H I between the initial state |A, B and the final C, D| which is formed just by rearrangements of the quarks(anti-quarks) in the initial states. Thus, the S matrix can be written down as where i, f denotes the initial and final states, and the matrix element is expressed as where the C denotes the momentum of particle C, similarly for others. Since the Hamiltonian is the interaction between the constituent (anti-)quarks, the matrix element would be the integration of the product of the spatial wave functions of the (anti-)quarks in the mesons and the constituent (anti-)quark scattering amplitude. To the Born order of the (anti-)quark scattering amplitude, there are four kinds of qq − qq scattering diagrams which are labeled according to which pairs of the constituents are interacted. These are q-q interaction diagrams, i.e. "capture 1 " (C 1 ) and "capture 2 " (C 2 ), and q-q (q-q) interaction diagrams, i.e. "transfer 1 " (T 1 ), and "transfer 2 " (T 2 ), as shown in Fig. 1.
To evaluate the contributions of these diagrams, it is more convenient to redefine the momentum variables as in Fig. 2 Then, the T -matrix element M f i is contributed by the sum of all four kinds of diagrams, and the contribution of every diagram could be written down as the product of signature, flavor, color, spin, and space factors, represented as h = I sign I f lavor I color I spin I space . The signature factor is (−1) for all diagrams because of interchanging three quark lines as shown in Fig.1. The color factor is (−4/9) for two capture diagrams (C 1 and C 2 ) and (4/9) for two transfer diagrams (T 1 and T 2 ). The spin factors for C 1 , C 2 , T 1 , and T 2 diagrams in different interaction terms are listed in TABLE I. Compilation of the spin factors for C1, C2, T1, and T2 diagrams in the spin-spin hyperfine, color Coulomb, and linear potential terms [34]. Table. I. The space factor for each diagram is an overlap integral of the meson wave functions times the underlying quark transition amplitude T pot f i . The overlap integrals of wave functions could be written down explicitly as The contributions to the quark-quark amplitude T pot f i by spin-spin, color Coulomb, and Linear confinement interactions read as We will deal with the scattering amplitudes of BB, BB * , and B * B * systems. For simplicity, the u and d quarks are assumed to have the same mass. To study the physically allowed system with a certain total isospin quantum number listed in Table.  The latter two terms of the second line in Eq.(6) are related to the "symmetrized diagrams" with the quark-line exchange rather than the anti-quark line exchange [34], which has the effects of interchanging the final C and D mesons. Including both kinds of diagrams also keeps the amplitudes for BB scattering to be Bose-symmetric. Similarly, for I = 1 BB scattering, the amplitude could be also written down as The B * B * scatterings are similar to the BB scatterings, so there is no need to write the relations down explicitly. For the BB * → BB * scatterings, when I=0, the scattering amplitude could be expressed as When I=1, the scattering amplitude reads as Since the incoming two mesons are not identical, there is no Bose symmetry in this scattering amplitude. One have to notice that the latter term B 0 B − * |T |B 0 B − * , which is contributed by the diagrams with quark line interchanged, is not just the amplitude with the two final mesons interchanged in the first term.

III. PARTIAL WAVE DECOMPOSITION AND UNITARIZATION
In the previous section, we only calculated the lowest order scattering amplitude, which contains different partial wave component and does not satisfy the nonperturbative unitarity. Thus it needs partial wave decomposition and unitarization. Here we briefly describe the partial-wave decomposition and unitarity relation in our convention.
In general, consider the scattering process of 12 → 1 ′ 2 ′ , where all particle are massive. We use k(k ′ ) to denote the momentum of partical 1 (1 ′ ) in the initial(final) states in the center of mass system and E 1 , E 2 (E ′ 1 , E ′ 2 ) to denote the energies for the two initial (final) particles with the total energy E = E 1 + E 2 . The spins and the third components are denoted as s 1 , σ 1 for particle 1 and so on. By using the convention of [36], the amplitude could be expanded in partial waves as If we only considered the spin-spin, color Coulomb, and linear confinement interactions here, the total orbital angular momentum and the total spin are conserved separately. The partial wave amplitude with total angular momentum j will be whereθ is the angle between k and k ′ . That means, in this calculation, the partial-wave elastic unitarity condition is as simple as that of the scalar particles Im[M j ls ] = πM j * ls M j ls , where we have omitted the repeated subscript ls for brevity and Im[ ] means the imaginary part of the related function. If we redefine the amplitude M j ls = 1 π ρt l by extracting the kinematic factor ρ = 2| k| E , one will obtain a familiar form similar to the elastic unitarity condition of partial waves for scalar particles as or in a more concise form In the constituent interchange model, only the Born term is calculated, so there is no elastic cut in the scattering amplitude and it does not obey the unitarity relation. One need to restore it by adopting a suitable unitarization scheme. Here we use the K-matrix unitarization method, which could be regarded as summing over all the bubble chains. Then, the unitarized partial-wave S matrix element could be represented as and the scattering T -matrix element is The pole of S-matrix element below the threshold on the real axis of the first Riemann sheet, satisfying 1 − iρt l = 0, represents a bound state. When the unitarity relation is satisfied, the scatteing S-matrix on the second Riemann sheet is the inverse of that on the first Riemann sheet, S 2nd sheet = 1/S 1st sheet , that means, the zero point of the first Riemann sheet corresponds to the pole of the second sheet. Thus, the zero point satisfying 1 + iρt l = 0 below the threshold represents a virtual state.

IV. NUMERICAL CALCULATIONS AND DISCUSSIONS
The parameters in the calculation is provided by the Godfrey-Isgur (GI) model [37] because its interactions are similar to those in the Barnes-Swanson model and it presented a generally successful prediction to the meson mass spectrum. In the GI model, the wave functions of mesons are expanded in a series of a very large number of harmonics oscillator (HO) wave functions, which make it difficult to decompose the amplitude in the angular momentum in an analytical form, so we approximate the meson wave function by a HO wave function carrying the same radial quantum number and orbital angular momentum as the meson, with its effective radius obtained by the rms radius r rms of the related meson state in the GI model.
The running coupling function is parameterized as α s (q 2 ) = 0.25e −q 2 + 0.15e − q 2 10 + 0.20e − q 2 1000 to saturate the result of perturbation theory calculation in the large q 2 region and avoid the divergence in the low q 2 region, and the quark masses used here are m u = m d = 0.22GeV, m b = 4.977GeV. The strength coefficient of the confinement linear potential is b = 0.18. We used the HO wave functions, with the oscillator parameters of the bottomed mesons as β B = 0.579GeV, β B * = 0.542GeV (β is defined in the HO wave function by ψ( r) ∼ (polynomial) × e −β 2 r 2 /2 ).
In this calculation, only the partial S-wave scatterings of B ( * ) B ( * ) are investigated. The scattering systems are labeled by their total isospins and total spins as (I tot , S tot ). Thus, the B ( * ) B ( * ) systems discussed here, which could have non-vanishing partial S-waves, are (1, 0) BB, (1, 1) BB * , (0, 1) BB * , (1, 0) B * B * , (0, 1) B * B * , (1, 2) B * B * , as listed in Table. II. Using the standard parameters listed above, there are a bound state found in (0, 1) BB * system and (0, 1) B * B * system respectively, and one virtual pole is found in each of the other systems, whose pole positions are listed in Table.III. All the bound states and virtual states are just near the thresholds of the related channels.
Usually, if there is a virtual state close enough to the threshold, it will produce an enhancement for the absolute square of the scattering amplitude |T | 2 just above the related thresholds, which could be discerned in experiments. Whether the threshold enhancement could be recogonized as a state also depends on the interplay between the wave functions and the threshold. As shown in Fig.3, although there is a near threshold virtual state in each of the (1, 0) BB, (1, 1) BB * , (1, 0) B * B * , and (1, 2) B * B * systems, the near threshold peaks of |T | 2 in (1, 1) BB * and (1, 2) B * B * systems are more obvious than the other two even though the virtuals states in the latter are closer to the related thresholds than the former. III. The pole positions of the physically allowed systems for B ( * ) B ( * ) . The subscript "v" denotes a virtual state and "b" a bound state.

system(threshold) Total isospin
Total spin Itot In Ref. [35], the authors extracted approximated local potentials from a part of the Born amplitudes which are calculated using the BS model. Then by solving the two-meson schrodinger equation with these potentials, they found that only the I = 0 BB * channel is attractive enough to form a I(J P ) = 0(1 + ) bound state with a binding energy of −5.5 MeV, which is consistent with our result. We also found that the (0, 1) B * B * system could also form another I(J P ) = 0(1 + ) bound state. Since our unitarization approach is using the full Born amplitude and the wave functions and coupling function used in the calculation are also different from the method In Ref. [35], the differences between the results may not be surprising. In addition, our method provide more informations about the appearance of the virtual states. Another approach using the heavy meson chiral effective theory in Ref. [22] also found that the I(J P ) = 0(1 + ) BB * and B * B * systems are attractive and obtain two bound states with their binding energies to be ∆E BB * ∼ −12.6 +9.2 −12.9 MeV and ∆E B * B * ∼ −23.8 +16.3 −21.5 MeV, respectively. This may indicate a strong interaction in both of the (0, 1) BB * and the (0, 1) B * B * system, which is similar in our result.
There is a simple qualitative argument for the appearance of the virtual states in different channels in our calculation. We have stated that the virtual state would produce a near threshold enhancement in the |T | 2 to be observed in the experiments, which means that there is a large scattering length in these processes. In fact, from the theoretical point of view, it is the large scattering length calculated from BS model that causes the existence of the virtual states. In the K-matrix formalism, the virtual state pole position is the solution to 1 + iρt l = 0. If the perturbative scattering amplitude t l at the threshold, which is just proportional to the scattering length, is positive and large enough, 1 + iρt l = 0 will be saturated near the threshold, because the kinematic factor iρ is negative below the threshold and approaches to zero at the threshold. We can also see that the larger the scattering length is, the nearer the position of solution is to the threshold. Because of the large scattering length, |T | 2 will usually present a peak just above the threshold, which is observed in the experiments as a state. Similarly, if t l is negative and its absolute value is large, it will imply a near-threshold bound state. But for the deep bound state, it will depend on the detailed behavior of t l below the threshold. Thus, in general, as long as a model can produce a scattering length with a large A further remark about the two I(J P ) = 0(1 + ) states generated in BB * and B * B * systems is in order. Since only the S-wave amplitudes are considered in our calculation and different systems are decoupled, the two states in (I tot = 0, S tot = 1) BB * and B * B * actually have the same quantum numbers but are different states. If these two channels are coupled to each other in reality, they could appear in both channels and may affect each other.

V. SUMMARY
In this paper, we use the Barnes-Swanson constituent interchange model to provide the Born term of B ( * ) B ( * ) scattering amplitude and then unitarize them using the K-matrix method. By analytically continuing the unitarized scattering amplitudes to the complex energy plane, the dynamically generated bound states or virtual states could be extracted. Two bound states with I(J P ) = 0(1 + ) are found at about 10600.9 MeV in BB * scattering and at about 10648.6 MeV in B * B * scattering, respectively. We also found that several near threshold virtual states could be found in (1, 0) BB, (1, 1) BB * , (1, 0) B * B * , and (1, 2) B * B * systems, which might produce threshold enhancements and correspond to doubly bottomed molecular states.
In comparison, the other methods used in discussing the states in these systems, such as solving the schrödinger equation [35] using the potentials extracted from the amplitude, or the heavy meson chiral effective theory discussion [22], could only produce the near threshold bound state information and can not say anything about virtual states. Besides the bound states, our method also predicts the near-threshold virtual states in these systems which may have observable effects and could also inspire more interests in searching for threshold enhancements in experimental explorations. We also provide a simple qualitative explanation that it is the large scattering length calculated from the BS model that is essential in generating these virtual states which may be a more general result. From these different approaches, we expect that the existence of the near threshold composite states may be a general result in the S-wave BB, BB * and B * B * systems.