Studying the $D_1D$ molecule in the Bethe-Salpeter equation approach

We study the possible bound states of the $D_1D$ system in the Bethe-Salpeter (BS) formalism in the ladder and instantaneous approximations. By solving the BS equation numerically with the kernel containing one-particle exchange diagrams and introducing three different form factors (monopole, dipole, and exponential form factors) at the vertices, we investigate whether the isoscalar and isovector $D_1D$ bound states may exist, respectively. We find that $Y(4260)$ could be accommodated as a $D_1D$ molecule, whereas the interpretation of $Z_2^+(4250)$ as a $D_1D$ molecule is disfavored. The bottom analog of $Y(4260)$ may exist but that of $Z_2^+(4250)$ does not.

Y (4260) and Z + 2 (4250) provide a great opportunity for understanding the strong interaction dynamics inside a hadron with moleculer inner structure assumptions since their masses are close to the D 1 D threshold. We will systematically study the D 1 D molecular state in the Bethe-Salpeter (BS) equation approach with three different form factors at the interaction verties, We will investigate the S-wave D 1 D systems with isospins I =0 and 1 being both considered. We will vary the binding energy where M is mass of the bound state) in a wide range and search for all the possible solutions with the cutoff parameter Λ in the form factor in a reasonable interval. Through this process, we will naturally check whether Y (4260) and Z + 2 (4250) may exist as a S-wave D 1 D molecular state. The possible B 1 B molecular state will also been studied in our work.
In the rest of the manuscript we will proceed as follows. In Sec. 2, we will establish the BS equation for the bound state of an axial-vector meson (D 1 or B 1 ) and a pseudoscalar meson (D or B). Then the numerical results for the D 1 D and B 1 B systems will be presented in Sec. 3. In Sec. 4 we will present a summary of our results.
2 The BS formalism for D 1 D system As discussed in Ref. [24], the flavor wave functions of Y (4250) and Z + 2 (4250) are and respectively.
Based on the picture that Y (4250) and Z + 2 (4250) are composed of an axial-vector meson (D 1 ) and a meson (D), its BS wave function is defined as where D 1 (x 1 ) and D(x 2 ) are the field operators of the pseudoscalar mesons D 1 and D at space coordinates x 1 and x 2 , respectively, P = M v is the total momentum of Y (4250) or Z + 2 (4250) and v is its velocity. Let m D 1 and m D be the masses of the D 1 and D mesons, respectively, p be the relative momentum of the two constituents, and define . The BS wave function in momentum space is defined as where X = λ 1 x 1 + λ 2 x 2 is the coordinate of the center of mass and x = x 1 − x 2 . The momentum of the D 1 meson is p 1 = λ 1 P + p and that of the D meson is p 2 = λ 2 P − p.
It can be shown that the BS wave function of the D 1 D system satisfies the following BS equation [39]: where S µν D 1 (p 1 ) and S D (p 2 ) are the propagators of D 1 and D, respectively, and K νλ (P, p, q) is the kernel, which is defined as the sum of all the two particle irreducible diagrams with respect to D 1 and D mesons. For convenience, in the following we use the variables p l (= p·v) and p t (= p − p l v) as the longitudinal and transverse projections of the relative momentum (p) along the bound state momentum (P ), respectively. Then the propagator of the D 1 meson can be expressed as and the propagator of the pseudoscalar D meson has the form where The momentum of p 1 in the numerator of Eq. (2.6) are represented by p l and p t in the following: (2.8) In the BS equation approach, the interaction between D 1 and D mesons can be due to the light vector-meson (ρ and ω) and the light scalar-meson (σ) exchanges. Based on the heavy quark symmetry and the chiral symmetry, the relevant effective Lagrangian used in this work is shown in the following [24]: where a and b are represent the light flavor quark, V µ is a 3 × 3 Hermitian matrix containing ρ, ω, K * , and φ: The coupling constants involved in Eq. (2.9) are related to each other as follows [24]: where f π = 132 MeV, g π = 3.73 and g A = 0.6 [40]. As in Ref. [41], we take |g σ | = |g σ | and |h σ | = |h σ | approximately when performing the numerical analysis. The parameters β 2 g V and λ 2 g V are given by 2g ρN N and 3 10m N (g ρN N +f ρN N ), respectively, where g 2 ρN N /4π = 0.84 and f ρN N /g ρN N = 6.10 [42]. As to the two parameters ζ 1 and µ 1 involved in the coupling constants, the information about them is very scarce and these two parameters have not been determined. However in the heavy quark limit, we can roughly assume that the , respectively. The parameters µ = 0.1 GeV −1 and ζ = 0.1 are taken in Ref. [43]. In our calculations, we will vary µ 1 from 0.05 to 0.5 GeV −1 and ζ 1 from 0.05 to 0.5, while searching for possible solutions of the D 1 D bound states.
In the following we will given the kernel for the BS equation in the ladder approximation. There have been some studies on the legitimacy of applying the ladder approximation in the BS equation [44][45][46]. In Ref. [44] it was shown that including only ladder graphs in the scalar-scalar system cannot lead to the correct one-body limit, and to solve these problems, at least crossed ladder graphs should be included. In addition, from the naive perspective, for a large coupling constant, the ladder approximation is not legitimate [45]. However, there is a significant difference between our work and that studied in Refs. [45], in which the mass of the exchanged particle (µ) is very small compared to the mass of the constituent particle (m) with µ/m = 0.15. The exchanged particles in our work are σ, ρ and ω. In Ref. [46], the authors studied the KK bound states in the BS equation, they found that when ρ is exchanged the ratio of the contribution from the crossed graph to that from the ladder one is less than 15% and the result is almost the same when ω is exchanged. Therefore, the ladder approximation is a good one which should not affect our qualitative conclusions.
where m V represent the masses of the exchanged light vector mesons ρ and ω , c I is the isospin coefficient: c 0 = 3, 1, 1 and c 1 = −1, 1, 1 for ρ, ω, and σ, respectively, ∆ and ∆ µν represent the propagators for the scalar and vector mesons, respectively. In order to describe the phenomena in the real world, we should include a form factor at each interacting vertex of hadrons to include the finite-size effects of these hadrons. For the meson-exchange case, the form factor is assumed to take the following form [47]: where Λ, m and k represent the cutoff parameter, the mass of the exchanged meson and the momentum of the exchanged meson, respectively. The value of Λ is near 1 GeV which is the typical chiral symmetry breaking scale. In general, for an axial-vector meson (D 1 ) and a pseudoscalar meson (D) bound state, the BS wave function χ µ P (p) has the following form: where f i (p) (i = 0, 1, 2, 3) are Lorentz-scalar functions and µ represents the polarization vector of the bound state. After considering the constraints imposed by parity and Lorentz transformations, it is easy to prove that χ µ P (p) can be simplified as where the function f (p) contains all the dynamics.
In the following derivation of the BS equation, we will apply the instantaneous approximation, in which the energy exchanged between the constituent particles of the binding system is neglected. In our calculation we choose the absolute value of the binding energy E b of the D 1 D system (which is defined as E b = M − m 1 − m 2 ) less than 60 MeV. In this case the exchange of energy between the constituent particles can be neglected.

Numerical results
In this subsection, we will solve the BS equations numerically for the D 1 D systems with I = 0 and I = 1 based on the formulas presented in Sec. 2 and study whether the S-wave D 1 D molecular states exist. We first need to reduce the BS equation (2.17) to a onedimensional from. By choosing the appropriate contour and performing the integration over p l on both sides through applying the residue theorem, we can reduce the BS equation (2.17) to a three-dimensional from. The corresponding BS wave function is in fact rotationally invariant, i.e.f (p t ) (wheref (p t ) = dp l f (p)) depends only on the norm of the three momentum, |p t |. Therefore, after completing the azimuthal integration, we can obtain the one-dimensional BS equation. The BS wave functions for D 1 D systems with I = 0 and I = 1 were solved numerically in our previous work by discretizing the integration region (0,∞) into n pieces [52] (n is chosen to be sufficiently large and we use n-point Gauss quadrature rule to evaluate the integrals). Then the BS wave function can be written as an n-dimention vector. The coupled integral BS equation becomes a matrix equation where A(|p t | n , |q t | n ) corresponding to the coefcients in Eq. (2.17). Generally, |p t | varies from 0 to +∞ and f (|p t |) will decrease to zero when |p t | → ∞. To apply the Gaussian quadrature rule, we need to convert the Gaussian integration nodes into the physical values for |q t |, which can be done using the following equation: where is a parameter introduced to avoid divergence in numerical calculations, w and y are parameters used in controlling the slopes of wave functions and finding the proper solutions for these functions, and t varies from -1 to 1. One can then obtain the numerical results of f (|p t |) by requiring the eigenvalue of the eigenvalue equation to be 1. It can be seen from Eq. (2.17) that there is only one free parameter in our model, the cutoff Λ, which can not be uniquely determined and has various forms phenomenologically. It contains the information about the nonpoint interaction due to the structures of hadrons. The value of Λ is near 1 GeV which is the typical scale of nonperturbative QCD interaction. In this work, we shall treat the cutoff Λ in the form factors as a parameter varying in a much wider range 0. 8-4.8 GeV to see if the BS equation has solutions. We also vary the parameters ζ 1 and µ 1 in the reasonable range, in order to check if the results are sensitive to the effective coupling constants.

D 1 D system
In our calculation, we choose to work in the rest frame of the system in which P = (M ,0). We take the averaged masses of the mesons from the PDG The 3D cutoff Λ M , Λ D and Λ E graphics are presented in Figure 2. From these figures, we see that the cutoff show small variations with respect to the changes in the parameters ζ 1 and µ 1 and the cutoff Λ is more sensitive to ζ 1 than µ 1 . We can also intuitively find the cutoff Λ in the dipole form factor are larger than those in the monopole and exponential form factors while they vary in reasonable regions. Therefore, Y (4260) could be a molecular state in appropriate effective coupling constants and the cutoff.
In Figure 3,  In Ref. [26,27], the authors interpreted Y (4260) as a D 1 D bound state, but they also predicted a significantly smaller mass of about 4.22 GeV. Soon, after Y (4220) was observed by BESIII Collaboration [6]. The mass of Y (4220) is about 58 MeV below the threshold for D 1 D. Therefore, we varying the binding energy E b in the region from 0 to -60 MeV trying to find all the possible solutions.
From our calculations, the D 1 D system can not be an I = 1 bound state. Hence like Z + 2 (4250) can not be an I = 1 D 1 D molecule. The reason for that is the exchanges of ρ and ω cancel each other for the quantum number I = 1, and the small coupling constant g σ determines that including the σ contribution cannot make a substantial change to the kernel as studied in Ref. [53].   For the I = 0 D 1 D system, we find several regions for the cutoffs. The results are listed in Table 1, from which we can see the cutoffs Λ in three difference form factors are in reasonable ranges and all of them become smaller with the increase of the binding energy. We depict the 3D graphics (see Figure 4) showing the variation of the cutoffs with respect to the parameters ζ 1 and µ 1 when E b = −60 MeV. For other 3D graphics with different E b , the variation trend is the same as that of Figure 4. We can conclude that the D 1 D system could be a molecular state. Our result is consistent with the meson exchange model based on the heavy meson chiral perturbation theory [24] and the chiral quark model in which the isoscalar channel is found to be easier to bind than the isovector channel for the same components [54].

B 1 B system
The same procedure can be easily extended to study the bottom analogs of Y (4260) and Z + 2 (4250), by simply replacing the charm quark and antiquark with the bottom quark and antiquark, respectively. The masses of the bottom mesons are m B 1 = 5726.0 MeV and m B = 5279.4 MeV [5]. Due to the heavier mass of the bottom meson the kinematic term is relatively small, resulting in the bottom system to form the molecular state more easier.
We use the same set of parameters as in the D 1 D system. With these parameters the I = 0 B 1 B bound state always exists with the reasonable cutoff. The numerical results for the B 1 B system are shown in Table 2. Similar to the D 1 D system, all the cutoffs Λ M , Λ D and Λ E become smaller with the increase of the binding energy. The variations of the cutoffs with respect to the parameters ζ 1 and µ 1 when E b = −60 MeV are presented in Figure 5.

summary
In this work, we studied whether Y (4260) and Z + 2 (4250) could be a D 1 D molecular state in the Bethe-Salpeter equation approach. In our model, we applied the ladder and instantaneous approximations to obtain the kernel containing one-particle-exchange diagrams and introduced three different form factors (the monopole form factor, the dipole form factor, and the exponential form factor) since the constituent particles and the exchanged particles are not pointlike. The cutoff Λ introduced in the form factors reflects the effects of the structures of interacting particles. Since Λ is controlled by nonperturbative QCD and cannot be determined accurately, we let it vary in a reasonable range within which we try to find possible bound states of the D 1 D system. Since the two effective coupling constants µ 1 and ζ 1 have not been determined we varied them in larger ranges, i.e. [0.05, 0.5] GeV −1 and [0.05, 0.5], respectively.
Our numerical results indicated that when the parameters are within reasonable ranges the D 1 D system can form an I = 0 molecular state but cannot form an I = 1 molecular state. In other words, Y (4260) could be accommodated as a D 1 D molecule. However, the existence of Z + 2 (4250) as a molecule requires that the coupling constant g σ be enhanced by several times. Consequently, the interpretation of Z + 2 (4250) as a D 1 D molecule is disfavored. Its structure should be studied further.
The bottom analogs of Y (4260) and Z + 2 (4250) were also studied. Similar to the D 1 D system, the B 1 B system with I = 0 also can form a molecular state but cannot for I = 1. We expect forthcoming experimental measurements to test our model for the D 1 D and B 1 B systems.