Study on the possible molecular state composed of $D^*_s\bar D_{s1} $ within the Bethe-Salpeter framework

Recently a vector charmonium-like state $Y(4626)$ was observed in the portal of $D^+_sD_{s1}(2536)^-$. It intrigues an active discussion on the structure of the resonance because it has obvious significance for gaining a better understanding on its hadronic structure with suitable inner constituents. It indeed concerns the general theoretical framework about possible structures of exotic states. Since the mass of $Y(4626)$ is slightly above the production threshold of $D^+_s\bar D_{s1}(2536)^-$ whereas below that of $D^*_s\bar D_{s1}(2536)$ with the same quark contents as that of $D^+_s\bar D_{s1}(2536)^-$, it is natural to conjecture $Y(4626)$ to be a molecular state of $D^{*}_s\bar D_{s1}(2536)$, as suggested in literature. Confirming or negating this allegation would shed light on the goal we concern. We calculate the mass spectrum of a system composed of a vector meson and an axial vector i.e. $D^*_s\bar D_{s1}(2536)$ within the framework of the Bethe-Salpeter equations. Our numerical results show that the dimensionless parameter $\lambda$ in the form factor which is phenomenologically introduced to every vertex, is far beyond the reasonable range for inducing an even very small binding energy $\Delta E$. It implies that the $D^*_s\bar D_{s1}(2536)$ system cannot exist in the nature as a hadronic molecule in this model, so that we may not think the resonance $Y(4626)$ to be a bound state of $D^*_s\bar D_{s1}(2536)$, but something else, for example a tetraquark and etc.

Since 2003 many exotic resonances X, Y and Z bosons [9][10][11][12][13][14][15][16][17][18] have been experimentally observed, such as X(3872), X(3940), Y (3940), Z(4430), Y (4260), Z c (4020), Z c (3900), Z b (10610) and Z b (10650) (of course, not a complete list). The states have attracted attention of theorists because their structures obviously are beyond the simple qq settings for mesons. If we can firmly determine their compositions, it would definitely enrich our knowledge on hadron structures and moreover shed light on the non-perturbative QCD effects at lower energy ranges. Studies with different explanations on the inner structures [19] have been tried, such as molecular state, tetraquark or dynamical effect [20]. Anyway, all the ansatz have a certain reasonability, but a unique picture or criterion for firmly determining the inner structures is still lacking. Nowadays, the majority of phenomenological researchers conjectures what the concerned exotic states are made of, just based on the available experimental data. Then by comparing the results with new data one can check the validity degree of the proposal. If the results obviously contradict to the new measurements of better accuracy, the ansatz should be abandoned. Following this principle we explore Y (4626) by assuming it to be a molecular state of D * sD s1 (2536), and then using more reliable theoretical framework to check the scenario and see if the proposal from our intuition can be valid.
Concretely, in this work supposing Y (4626) as a D * sD s1 (2536) molecular state, we employ the Bethe-Salpeter (B-S) equation which is a relativistic equation established on the basis of quantum field theory, to study the two-body bound state [21]. Initially, the B-S equation was used to study the bound state of two fermions [22][23][24], later the method was generalized to the system of one-fermion-one-boson [25]. In Ref. [26,27] the authors employed the Bethe-Salpeter equation to study some possible molecular states, such as KK and BK system. With the same approach the bound state of Bπ, D ( * ) D ( * ) , B ( * ) B ( * ) are studied [28,29]. Recently the approach was applied to explore doubly charmed baryons [30,31] and pentaquarks [32,33]. In this work, we try to calculate the spectrum of Y (4626) composed of a vector meson and an axial vector meson.
If two constituents can form a bound state the interaction between them should be large enough to hold them into a bound state. The chiral perturbation theory tells us that two hadrons interact via exchanging a certain mediate meson(s) and the forms of the effective vertices are determined by relevant symmetries, but the coupling constants generally are obtained by fitting data. For the molecular states, since two constituents are color-singlet hadrons the exchanged particles are some light mesons with proper quantum numbers. It is noted that even though there are many possible light mesons contributing to the effective interaction between the two constituents, generally one or several of them would provide the dominant contributions. Moreover, beyond it, most of time the scenario with some other mesons exchange should also be taken into account, because even though the extra contributions are small comparing to the dominant one(s), they sometimes are not negligible, namely it would make the secondary contribution to the effective interaction. Then the effective kernel for the B-S equation is set. For the D * sD s1 (2536) system, the contribution of η [34][35][36] dominates, whereas in Ref. [34] the authors suggested σ exchange makes the secondary contribution. In our case considering the concerned quark contents of D * s and D s1 (2536), the contribution of f 0 (980) should stand as the secondary one. The effective interactions induced by exchanging η are deduced with the heavy quark symmetry [34][35][36] and we present the formulas in the appendix A. With the effective interactions we can derive the kernel and establish the corresponding B-S equation.
With all necessary parameters being beforehand chosen and input, the B-S equation is solved numerically. In some cases the equation does not possess a solution if one or several parameters are set within a reasonable range, then a conclusion must be drawn that the proposed bound state should not exist in nature. On the contraries, a solution of the B-S equation with reasonable parameters implies that the corresponding bound state is formed. In that case, simultaneously the B-S wave function is obtained which can be used to calculate the rates of strong decays, which will help experimentalist to design new experiments for further measurements.
This paper is organized as follows: after this introduction we will derive the B-S equation related to possible bound state composed of D * s andD s1 (2536) which are a vector and an axial vector meson respectively. In section III the formula for its strong decays are present. Then in section IV we will solve the B-S equation numerically. Since Y (4626) is supposed to be a molecular bound state, the input parameters must be within a reasonable range, but our results say that this mandatory condition cannot be satisfied, thus we think that such a molecular state of D * sD s1 (2536) may not exist. However, as we deliberately set the parameters to a region which is not favored by all previous phenomenological works, we can obtain the required spectrum and corresponding wavefunctions. With the wavefunction we evaluate the strong decay rate of Y (4626) and present our results by figures and tables. Section IV is devoted to a brief summary. Since the newly observed resonance Y (4626) contains hidden charms and its mass is close to the sum of the masses of D * s andD s1 where D * s −D s1 corresponds to D * + s − D − s1 or D * − s − D + s1 , a conjecture about its molecular structures composed of D * s andD s1 is favored. For a state with spin-parity being 1 − , its spatial wave function is in S wave. There are two FIG. 1: the bound states of D * sDs1 formed by exchanging η.
). We will focus on such an ansatz and try to find numerical results by solving the relevant B-S equation.
A. The Bethe-Salpeter (B-S) equation for 1 − D * sDs1 molecular state By the effective theory D * s andD s1 interact mainly via exchanging η. The Feynman diagram at the leading order is depicted in Fig. 1. To take into account the secondary contribution induced by exchanging other mediate mesons, in Ref. [34] the authors consider a contribution of exchanging σ to the effective interaction. Since there are neither u nor d constituents in D * s andD s1 , their coupling to σ would be very weak, thus the secondary contribution to the interaction may come from exchanging f 0 (980) instead. The relevant Feynman diagrams are shown in Fig. 2. The relations between relative and total momenta of the bound state are defined as where p 1 and p 2 (q 1 and q 2 ) are the momenta of the constituents, p and q are the relative momenta between the two constituents of the bound state at the both sides of the diagram, P is the total momentum of the resonance, η i = m i /(m 1 + m 2 ) and m i (i = 1, 2) is the mass of the i-th constituent meson. k is the momentum of the exchanged mediator.
A detailed analysis on the Lorentz structure [24,26,27] determines the form of the B-S wave function of the bound state comprising a vector and an axial vector mesons ( D * s and where a, b, c and d are Lorentz indices. The wave function in the momentum space can be obtained by carrying out a Fourier transformation Using the so-called ladder approximation one can get the B-S equation deduced in earlier references [21][22][23] where ∆ 1aα and ∆ 2bβ are the propagators of D * s andD s1 respectively, K αβµν (P, p, q) is the kernel determined by the effective interaction between two constituents which can be calculated from the Feynman diagrams in Fig. 1 and 2. In order to solve the B-S equation, we decompose the relative momentum p into the longitudinal component p l (= p · v) and the transverse one p µ t (= p µ − p l v µ )=(0, p T ) with respect to the momentum of the bound state P (P = Mv).
where M is the mass of the bound state Y (4626), ω i = p l 2 + m 2 i and ξ is the parameter in the gauge fixing. We set ξ to be 1 in our calculation.
By the Feynman diagrams shown in Fig. 1 and 2, the kernel K αβµν (P, p, q) is written as where m s is the mass of the exchanged meson, s D * s f 0 and g D s1 D s1 f 0 are the concerned coupling constants and ∆(k, m s ) = i/(k 2 − m 2 s ). Since the two constituents of the molecular state are not on shell, at each interaction vertex a form factor should be introduced to compensate the off-shell effect. The form factor is employed in many references [37][38][39][40], even though it has different forms. Here we set it as: where k is the three-momentum of the exchanged meson and Λ is a cutoff parameter. Indeed, the form factor is introduced phenomenologically and there lacks any reliable knowledge on the value of the cutoff parameter Λ. Λ is often parameterized to be λΛ QCD + m s with Λ QCD = 220 MeV which is adopted in some references [37][38][39][40]. As suggested, the order of magnitude of the dimensionless parameter λ should be close to 1. In our later numerical computations, we set it to be within a wider range of 0 ∼ 4 . The wave function can be written as where ǫ is the polarization vector of the bound state and f (p) is the radial wave function. The three-dimension spatial wave function is obtained after integrating over p l Substituting Eqs. (7) and (9) into Eq. (4) and multiplying ε abf g χ * g P (x 1 , x 2 )P f on both sides one can sum over the polarizations of both sides. Employing the so-called covariant instantaneous approximation q l = p l i.e. using p l to replace q l in K(P, p, q), the kernel K(P, p, q) does not depend on q 1 any longer. Then we take a typical procedure: integrating over q l on the right side of Eq. (4), multiplying dp l (2π) on the both sides of Eq. (4), and integrating over p l on the left side, to reduce the expression into a compact form. Finally we obtain (11) with While we integrate over p l on the right side of Eq. (11) there exist four poles which are By choosing an appropriate contour we only need to evaluate the residuals at Here , one can integrate out the azimuthal part and then Eq. (11) is reduced into a one-dimensional integral equation with

B. Normalization condition for the B-S wave function
In analog to the cases in Refs. [26,27] the normalization condition for the B-S wave function of a bound state should be i 6 where P 0 is the energy of the bound state which is equal to its mass M in the center of mass frame. I(P, p, q) is a product of reciprocals of two free propagators with a proper weight.
In our earlier work [29] we found that the term K abαβ (P, p, q) in brackets is negligible, so that now we ignore it as done in Ref. [41].
After performing some manipulations we obtain the normalization of the radial wave function 1 π 2

III. THE STRONG DECAYS OF THE MOLECULAR STATE Y (4626)
Now we investigate the strong decays of Y (4626) using the effective interactions.
The amplitude can be parameterized as [42] A a = g 0 The factors g 0 and g 2 are extracted from the expressions of A a . Then the partial width is expressed as The corresponding Feynman diagram is depicted in Fig. 3 (b) whereD ′ s1 denotes D s (2460) through the whole paper. The amplitudes is The amplitude can also be parameterized as where ǫ 2 is the polarizations ofD s (2460). The factors g ′ 0 and g ′ 2 can be extracted from the expressions of A b .

(c). The amplitudes is
where ǫ 1 and ǫ 2 are the polarization vectors of D s (2460) andD * s respectively. The total amplitude can be parameterized as [42] A c = g 10 ε µναβ P µ ǫ 1ν ǫ 2α ǫ * β + The factors g 10 , g 11 and g 12 are extracted from the expressions of A c . Fig. 3

(d). The amplitude is
where ǫ 1 and ǫ 2 are the polarization vectors of D * s andD s (2460) respectively. The total amplitude for the strong decay of Y (4626) → D * s (1 − ) +D s (2460)(1 + ) can also be expressed as The factors g ′ 10 , g ′ 11 and g ′ 12 are extracted from the expressions of A d .
The Feynman diagram is depicted in Fig. 3(e) whereD s2 representsD s (2572). The amplitudes is, where ǫ 2 is the polarization tensor ofD s (2572)(2 + ). The total amplitude is written as The factors g 20 can be extracted from the expressions of A e .
The amplitude is still written as The factors g ′′ 0 and g ′′ 2 are extracted from the expressions of A f .

A. the numerical results
Before we numerically solve the B-S equation all necessary parameters should be priori determined as accurate as possible. The masses m D * s , m D s0 , m D s1 , m D ′ s1 , m D s2 , m η and m ρ come from the databook [43]. The coupling constants in the effective interactions g D s1 D * s η , η and g D s1 D s2 η are taken from the relevant literatures and their values and related references are collected in the Appendix.
With these input parameters the B-S equation Eq. (12) can be solved numerically. Since it is an integral equation, an efficient way for solving it is discretizing it and then turns solving the integral equation to an algebraic equation group. Concretely, we let the variables |p T | and |q T | be discretized into n values Q 1 , Q 2 ,...Q n ( n=129 in our calculation) and the equal gap between two adjacent values as Qn−Q 1 n−1 . Here we set Q 1 =0.001 GeV and Q n =2 GeV . The n values of f (|p T |) constitute a column matrix on the left side of the equation and the n elements f (|q T |) construct another column matrix on the right side of the equation as shown below. In this case, the functions in the curl bracket of Eq. (12) multiplying |q T | 2 12M 2 (2π) 2 would be an effective operator acting on f (|q T |). It is specially noted that because discretizing the equation, even |q T | 2 12M 2 (2π) 2 turns from continuous integration variable into n discrete values which are involved in the n × n coefficient matrix. Substituting the n pre-set Q i values into those functions, the operator turns into an n × n matrix which associates the two column matrices. It is noted that Q 1 , Q 2 ,...Q n should take sequential values.
As is well known, if a homogeneous equation possesses non-trivial solutions, the necessary and sufficient condition is det|A(∆E, λ) − I| = 0 (I is the unit matrix) where A(∆E, λ) is just the aforementioned coefficient matrix. Thus solving the integral equation just turns to a sort of eigenvalue searching problem which is a familiar issue in quantum mechanics, in particular, the eigenvalue is required to be unit in this problem. Here A(∆E, λ) is a function of the binding energy ∆E = m 1 + m 2 − M and parameter λ. The following procedure is a bit tricky. Inputting a supposed ∆E, we vary λ to make det|A(∆E, λ) − I| = 0 hold. One can note that the matrix equation (A(∆E, λ) ij )(f (j)) = β(f (i)) is exactly an eigenequation.
Using the values of ∆E and λ, we seek out all possible "eigenvalues" β. Among them only β = 1 is the solution we expect. In the process of solving the equation group, the value of λ is determined, and actually it is the solution of the equation group with β = 1. Meanwhile using the obtained λ, one achieve the corresponding wavefunction f (Q 1 ), f (Q 2 )...f (Q 129 ) which just is the solution of the B-S equation.
Generally λ should be within the range around the order of unit. In Ref. [37] the authors fixed λ to be 3. In our earlier paper [40] the value of λ varied from 1 to 3. In Ref. [33] we set the value of λ within a range of 0 ∼ 4 by which as believed, a bound state of two hadrons can be formed. When the obtained λ is much beyond this range, one would conclude that the molecular bound state may not exist or at least is not a stable state. But it is really noted that the form factor is phenomenologically introduced and the parameter λ is usually fixed via fitting data, i.e. neither the form factor nor the value of λ are derived from an underlying theory, but based on our intuition (or say, a theoretical guess). Since the concerned processes are dominated by the non-perturbative QCD effects whose energy scale is about 200 MeV, we have reason to believe that the cutoff should fall within a range around a few hundreds of MeV to 1 GeV, and by this allegation one can guess that the value of λ should be close to unit. However, from other aspect, this guess does not have a solid support, further phenomenological studies and a better understanding on low energy field theory are needed to get more knowledge on the form factor and value of λ. So far, even though we believe this range for λ which sets a criterion to draw our conclusion at present, we cannot absolutely rule out the possibility that some other values of λ beyond the designated region may hold. That is why we proceed to compute the decay rates of Y (4626) based on the molecule postulate. (see below the numerical results for clarity of this point).
By our strategy, for the state Y 2 we let ∆E = 0.021 GeV which is the binding energy of the molecular state as M D * s + M D s1 (2536) − M Y (4626) . Then we try to solve the equation |A(∆E, Λ) − I| = 0 by varying λ within a reasonable range. In other words, we are trying to determine a λ whose value falls in the range of 0 to 4 as suggested in literature, to satisfy the equation.
As our result, we have searched a solution of λ within a rather large region, but unfortunately find that there is no solution which can satisfies the equation.
However, for the Y 1 state if one still keeps ∆E = 0.021 GeV but set λ = 10.18, the equation |A(∆E, λ) − I| = 0 holds while the contributions induced by exchanging both η and f 0 (980) are included. Instead, if the contribution of exchanging f 0 (980) (Fig. 2) is ignored, with the same ∆E one could get a value of λ which is very close to that without the contribution of f 0 (980). It means that the contribution from exchanging f 0 (980) is very small and can be ignored safely. Meanwhile by solving the eigenequation we obtain the wavefunction f (Q 1 ), f (Q 2 )...f (Q 129 ). The normalized wavefunction is depicted in Fig. 4 with different ∆E.
Due to existence of an error tolerance on measurements of the mass spectrum, we are allowed to vary ∆E within a reasonable range to fix the values of λ again, for the D s1D * s system, the results are presented in Tab. I. Apparently for a reasonable ∆E any λ value which is obtained by solving the discrete B-S equation is far beyond 4. Does the result imply that D s1D * s fails to form a bound state? We will further discuss its physical significance in next section.
A new resonance Y (4626) has been experimentally observed [1], and it is the fact everybody acknowledges, but what composition it has, demands a theoretical interpretation. The molecular state explanation is favored by an intuitive observation. However our theoretical study does not support the allegation that Y (4626) is the molecule of D * sD s1 . On other respect, the above conclusion is based on a requirement: λ must fall in a range of 0∼4, which is determined by phenomenological studies done by many researchers. However, λ being in 0∼4 is by no means a mandatory condition because it is not deduced form an underlying principle and lacks real foundation. Therefore even though our result does not favor the molecular structure for Y (4626), we still proceed to study the transitions Y . Using the wave function we calculate the form factors g 0 , g 2 , g ′ 0 , g ′ 2 , g 10 , g 11 , g 12 , g ′ 10 , g ′ 11 , g ′ 12 , g 20 , g ′′ 0 , g ′′ 2 defined in Eqs. (18, 21, 23, 25, 27 and 29). With these form factors we get the decay widths of (2573) and Y (4626) → D sDs2 (2536) which are denoted as Γ a , Γ b , Γ c , Γ d , Γ e and Γ f presented in Table II. Theoretical uncertainties originate from the experimental errors, namely the theoretically predicted curve expands to a band.
Of course, exchanging two η mesons can also induce a potential as the next-to-leading order (NLO) contribution, but it undergoes a loop suppression therefore, we do not consider that contribution i.e. one-boson-exchange model is employed in our whole scenario.

V. CONCLUSION AND DISCUSSION
In this work we explore the bound state composed of a vector and an axial vector within the B-S equation framework. Concretely we study the resonance Y (4626) which is assumed   found as a solution of the B-S equation, we would claim that resonance could be a molecular state. From the results shown in table I one can find that even for a small binding energy (we deliberately vary the value of the binding energy), the λ which makes the equation to hold, must be larger than 9 which is far beyond the favorable one in literature so that we tend to think the molecular state of D * sD s1 (2536) does not exist unless the coupling constants get larger than those given in Appendix.
As aforementioned discussion, the λ in the form factor at each vertex is phenomenologically introduced and does not receive a solid support from any underlying principle, therefore, we may suspect its application regime which might be the pitfall of the phenomenology. Thus we try to overcome this barrier to extend the value to a region which obviously deviates from the region favored by the previous works. As a λ value beyond 9, the solution of the B-S equation exists, and the B-S wavefunction is constructed. Just using the wavefunctions, we calculate the decay rates of Y The important and detectable issue is the decay patterns deduced above. This would compose a crucial challenge to the phenomenological scenario. If the decay patterns deduced in terms of the molecule assumption are confirmed (within an error tolerance) , it would imply that the constraint on the phenomenological application of form factor which is originating from the chiral perturbation can be extrapolated to a wider region. By contrary, if the future measurements negate the predicted decay patterns, one should claim that the assumption that Y (4626) is a molecular state of D * sD s1 (2536) fails. The resonance would be in different structure, such as tetraquark or hybrid etc.
We lay our hope on the future experimental measurements on those decay portals, which can help us to clarify the structure of Y (4626).
L DD * P = g DD * P D b (∂ µ M) ba D * µ † a + g DD * P D * µ b (∂ µ M) ba D † a + c.c., (A5) where c.c. is the complex conjugate term, a and b represent the flavors of light quarks, f 0 denotes f 0 (980), M is 3 × 3 hermitian and traceless matrix M =