The mass spectra and wave functions for the doubly heavy baryons with $J^P=1^+$ heavy diquark core

The mass spectra and wave functions for the doubly heavy baryons are computed under the picture that the two heavy quarks inside a doubly heavy baryon, such as two $c$-quarks in $\Xi_{cc}$, combine into a heavy `diquark core' in color anti-triplet firstly, then the diquark core turns into a color-less doubly heavy baryon via combining the light $q$-quark inside the baryon. Namely both of the combinations, the two heavy quarks inside the baryon into a diquark core in color anti-triplet and the heavy diquark core with the light quark into the baryon, are depicted by relativistic Bethe-Salpeter equations (BSEs) with an accordingly QCD inspired kernel respectively, although in the paper only the heavy diquark cores with the quantum numbers $J^P=1^+$ are considered. Since the `second combination' is of the heavy diquark core and the light quark, so the structure effect of the diquark core to the relevant kernel of the BSE is specially considered in terms of the diquark-core wave functions. The mass spectra and wave functions for the `low-laying' doubly heavy baryons in the flavors $(ccq)$, $(bcq)$ and $(bbq)$ and in the quantum numbers $J^P=\frac{1}{2}^+$, $J^P=\frac{3}{2}^+$, achieved by solving the equations under the so-called instantaneous approximation, are presented properly and some comparisons with the others' results under different approaches in the literature are made.


I. Introduction
The LHCb collaboration reported their observation on the doubly charmed baryon Ξ ++ cc (ccu) in the Λ + c K − π + π + decay channel recently, where the mass and lifetime of the baryon are determined as 3.621 GeV [1] and 0.256 ps [2], and thus systematically theoretical researches on the kindred baryons (doubly heavy baryons) become imperative. In addition, LHCb also reported observations of five narrow Ω 0 c excited baryons in 2017 [3] and Ξ − b in 2018 [4]. In most theoretical models, the mass of the doubly heavy baryon Ξ +(+) cc is predicted in the range 3.5 ∼ 3.7 GeV [5][6][7][8][9][10][11][12][13][14][15][16][17][18]. The mass splitting between Ξ ++ cc and Ξ + cc is predicted to be several MeV due to the mass difference of the light quarks u and d. The predicted mass by lattice QCD is about 3.6 GeV [19][20][21], which is quite close to the LHCb's observation. The lifetimes of Ξ + cc and Ξ ++ cc are predicted to be quite long as 50 ∼ 250 and 200 ∼ 700 fs [22][23][24][25], respectively. The LHCb's measurement on Ξ ++ cc is very close to the lower boundary of the theoretical predictions. It is expected that more information on the doubly heavy baryon(s) will be avaiable in the near future.
Compared with mesons, baryons as three-quark bound states are much more complicated. However, for a doubly heavy baryon, the involved two heavy quarks move relatively slow, and thus are bound with each other strongly to form a color anti-triplet diquark core 1 . Thus, it is reasonable to reduce the three-body problem into two two-body ones. According to the Pauli principle, the baryon wave functions must be totally anti-symmetric under the interchange of any two quarks. For the Email addresses: liruo@nwpu.edu.cn (Qiang Li), zhangzx@itp.ac.cn (Chao-Hsi Chang), sqin@cqu.edu.cn (Si-Xue Qin), gl_wang@hit.edu.cn (Guo-Li Wang) 1 According to QCD and color confinement, the gluon exchange interaction between two quarks is attractive only when they are in a color anti-triplet. Moreover only in a color anti-triplet the diquark can combine with a third quark to form a colorless baryon. Thus throughout this work the diquark cores are meant to be in a color anti-triplet, i.e. the color part of the diquark-core wave function is always anti-symmetric. diquark core with even orbital angle momentum L D , the spin-flavor part of the wave function must be symmetric. Hence, (cc)-and (bb)-diquarks at ground states (S-wave) must have J P = 1 + , while (bc)-diquark, carrying different flavors, may have J P = 1 + or 0 + . So far, in this work, we only focus on the 1 2 + or 3 2 + baryons with 1 + doubly heavy diquark cores (cc, bc or bb). Whereas the 1 2 + baryons with 0 + (bc)-diquark, e.g., Ξ bc and Ω bc , will be presented elsewhere.
Based on the analysis above, we reduce the three-body problem for the doubly heavy baryons into two two-body problems. First, we deal with diquark cores in the framework of Bethe-Salpeter equation (BSE) with a QCD-inspired kernel. In the so-called instantaneous approximation, we solve the corresponding three-dimensional (Bethe-)Salpeter equation, and obtain the diquark spectra and wave functions. Then, the second step is to establish the BSE for the diquark core and the remaining light quark, where the structured diquark effects in the BSE kernel are described by the corresponding form factors. For solving the baryon BSE, the instantaneous approximation is also adopted. In this scheme, both steps only deal with two-body (Bethe-)Salpeter equations. In fact, Refs. [26][27][28][29][30][31][32][33][34][35][36][37][38][39] have adopted the similar strategy with different approximations or methods. It is known that the BSE framework has been successfully used to study the problems of two-body meson systems, e.g., mass spectra [40,41], hadronic transitions, electro-weak decays, and etc [42][43][44][45][46][47][48]. The universal agreement between the theoretical predictions and experimental observations is achieved, so here we would like to push the BSE approach further to the present baryon study. Namely, based on BSE to develop a precise and systematic approach to describe the doubly heavy baryons is our main motivation for this work.
This paper is organized as follows. In section II, we describe the 1 + diquark BSE in the instantaneous approximation, solve the mass spectra and wave functions, and then compute the diquark form factors with the Mandelstam formulation. In section III, we derive the three-dimensional diquarkquark (Bethe-)Salpeter equation, and solve the mass spectra and wave functions for baryons with J P = 1 2 + and 3 2 + . Especially, the diquark-core form factors play a quite important role in determining the BS kernel for the baryons. In section IV, we discuss the obtained results and then compare them with those from other sources. Finally, we provide a brief summary of this work.
II. Bethe-Salpeter equation for a heavy diquark core and the relevant form factors

II.1. Interaction kernel
Since the QCD-inspired interaction kernel for doubly heavy diquark cores has the same root as that of doubly heavy mesons, we briefly review the latter first. In this work, it is assumed that the instantaneous approximation (IA) works well since the considered systems always involve the heavy quark. In the IA, the kernel has no dependence on the time-component of exchanged momentum, and can be expressed as a 'revised one-gluon exchange' form as where V G ( q ) is the usual one-gluon exchange potential and V S ( q ) is a phenomenological screened potential, they behave specifically in the Coulomb gauge as [49]: where 4 3 is the color factor; a 1(2) is introduced to avoid the divergence in the region of small transfer momenta; the screened potential V S ( q ) is phenomenologically introduced and characterized by the string constant λ and the factor a 2 . The potential herein originates from the famous Cornell potential [50,51], which expresses the one-gluon exchange Coulomb-type potential at short distance and a linear growth confinement one at long distance. In order to incorporate the screening effects [52,53], the confinement potential is modified as the aforementioned form [54][55][56], where V 0 is a free constant fixed by fitting data. Note that a screened potential containing time-like vector type has been widely used and discussed in many works [57][58][59][60][61], and could achieve well description on the meson systems. Here we assumed V S also arises from the gluon exchange and hence carries the same color structure with V G . The strong coupling constant α s has the form, where Λ QCD is the scale of the strong interaction, N f is the active flavor number, and a = e is a constant regulator. For later convenience, we split V M ( q ) into two parts, where V M1 is a constant and V M2 contains all dependence on q, i.e., To construct the BSE kernel for the diquark, we need to transform an anti-fermion line into a fermion one by the charge conjugation. Considering the fact that the quark-antiquark pair in a meson is a color singlet, while the quark-quark pair inside a baryon is a color anti-triplet, and thus the corresponding color factors are 4 3 and − 2 3 , respectively, hence the kernel for the diquark core is simply assumed to be, This 'half rule' used here is widely adopted in previous works involved the baryons and quark-quark bound problems [8,9,22]. It is exact in the one-gluon exchange limit and has been found to work well in the measurement of the Ξ ++ cc mass. Its successful extension beyond weak coupling implies that the heavy quark potential factorizes into a color dependent and a space-dependent part, with the latter being the same for quarkquark and quarkantiquark pairs. The relative factor of 1 2 then results from the color algebra, just as in the weak-coupling limit [62]. Also note that, for heavy quarks, compared with the temporal component (α = 0) of the vector-vector γ α ⊗ γ α interaction kernel, the spatial components (α = 1, 2, 3) are suppressed by a factor v ∼ | p| M and will be ignored in numerical calculations. This is also consistent with the analysis in Ref. [63].   First, we briefly review the Bethe-Salpeter equation for mesons, as shown in Fig. 1(a), which in the momentum space reads

II.2. Bethe-Salpeter equation for the heavy diquark cores
where Γ M denotes the BS vertex; S(s 1 ) and S(−s 2 ) are the propagators of the quark and antiquark, respectively; the corresponding BS wave function is defined as ψ M (p, s) ≡ S(s 1 )Γ M (p, s)S(−s 2 ). The quark internal momenta s and u are defined respectively as where λ i ≡ µ i µ 1 +µ 2 (i = 1, 2) and µ i is the constituent quark mass. The BSE normalization condition can be generally expressed as As shown in Fig. 1(b), by performing the charge conjugation transformation, the Bethe-Salpeter equation for diquark cores reads where Γ c ≡ Γ D C −1 , and C ≡ iγ 2 γ 0 denotes the charge conjugation operator; Γ D is the diquark vertex; and K c = 1 2 K M . Similarly, the diquark BS wave function can be defined as ψ c (p, s) ≡ S(s 1 )Γ c (p, s)S(−s 2 ). Note that Eq. (7) and Eq. (5) share exactly the same form with the only difference that the diquark bound strength is halved due to the color factor. Therefore, Eq. (7) can be easily solved to obtain mass spectra and wave functions of the diquark cores.
II.3. The Salpeter equation and wave function for the 1 + diquark core As pointed out in the Introduction, according to the Pauli principle, the heavy diquark core (bc) in color anti-triplet at the ground state (S-wave) may have J P = 1 + or 0 + , while (cc) and (bb) can only be in J P = 1 + . In this work, we restrict ourselves to 1 + diquark cores (cc), (bc) and (bb). Following the standard procedures in Ref. [64], we define the 3-dimensional Salpeter wave function ϕ c (p, s ⊥ ) ≡ i dsp 2π ψ c (p, s), where s p = s ·p, s ⊥ = s − s pp , andp = p µ . Then we can obtain the 3-dimensional Salpeter equation as where µ is the bound state mass; e i = µ 2 i − s 2 i⊥ (i = 1, 2) denotes the quark kinematic energy; H(s i⊥ ) ≡ 1 e i (s α i⊥ γ α + µ i )γ 0 are the usual Dirac Hamilton divided by e i ; W (s ⊥ ) ≡ γ 0 Γ c (p, s ⊥ )γ 0 is the potential energy part; and the 3-dimensional vertex is expressed as Also we obtain a constraint condition for the Salpeter wave function Accordingly, the normalization condition now becomes Since a diquark consists of two quarks, its parity is just opposite to its meson partner. Then, the wave function of the diquark with J P = 1 + can be decomposed as where s = −s 2 ⊥ andŝ ⊥ = s ⊥ s ; e is the polarization vector fulfilling the Lorentz condition e α p α = 0. Note that there are eight radial wave functions f i (i = 1, 2, · · · , 8), but only four of them are independent due to the constraint condition Eq. (10), i.e., where A ± ≡ s(e 1 ±e 2 ) µ 1 e 2 +µ 2 e 1 . Note that the meson wave functions with J P C = 1 −− share the same form of Eq. (12). Inserting the decomposition into Eq. (8) and taking the Dirac trace, we can obtain four coupled eigenvalue equations, which are explicitly shown in appendix A.1. The normalization can then be simply expressed as, Solving the coupled equations numerically, the mass spectra and wave functions of the diquark cores are obtained. The parameters a = e = 2.7183, λ = 0.21 GeV 2 , Λ QCD = 0.27 GeV, a 1 = a 2 = 0.06 GeV, and the constituent quark masses used in this work are just the same as those adopted in the previous studies and determined by fitting to the heavy and doubly heavy mesons spectra [41,[43][44][45][46][47][48]. The free parameter V 0 in the diquark system is fixed by fitting the mass of its meson partner to the experimental data, e.g., V 0 of the 1 + (cc)-diquark is determined by the J/ψ. In this work, we obtain V 0 = −0.221 GeV for 1 + (cc), V 0 = −0.147 GeV for 1 + (bc), and V 0 = −0.026 GeV for 1 + (bb). The obtained mass spectra Tab. I: Mass spectra of J P = 1 + color anti-triplet diquark cores (cc), (bc) and (bb) in units of GeV.
GeV q  where ψ V and ψ T are related to the Salpeter radial wave functions by The obtained numerical values of ψ V and ψ T are listed in Tab. II.
Tab. II: Wave functions at origin of 1 + diquark cores (cc), (bc) and (bb) for the ground states (in GeV 2 ).  In the BSEs of mesons or diquarks, the constituent quarks (antiquarks) are the point-like particles. But in the BSEs of baryons in terms of quark and diquark, the structures of the non-pointlike diquarks have to be considered. The coupling between diquarks and gluons is modified by form factors. Thus, in this subsection, we will discuss the relevant form factor within the Mandelstam formulation.
The Feynman diagrams of a doubly heavy diquark core coupling to a gluon are shown in Fig. 3. Since the coupling vertex Σ αβµ , as a 'vector current matrix element', is conserved, it can be decomposed by three independent form factors σ i (i = 1, 3, 5) as where σ i only depend on the transfer momentum t 2 ≡ (p − p ) 2 . Note that, in the instantaneous approximation, p 0 = p 0 and hence t 2 = (p ⊥ − p ⊥ ) 2 is always space-like. Considering the fact that the contributions of σ 3 and σ 5 to the BSEs of baryons is small compared with σ 1 , we will keep the dominant σ 1 only and omit the subscript for simplicity.
Corresponding to the Feynman diagrams in Fig. 3, the form factor consists of two terms where the factor 1 2 is due to the normalization condition and guarantees the form factor σ(t 2 ) = 1 at zero transfer momentum (t 2 = 0). Explicitly, for example, the amplitude Σ αβµ 1 corresponding to Fig. 3(a) reads where the contour integration over s 0 is performed and only the dominant contribution is kept. And the amplitude of Fig. 3(b) can be easily obtained by the relation Σ αβµ . Inserting the 1 + Salpeter wave function Eq. (12) into Eq. (16), we can obtain σ for 1 + (cc)diquark in the ground and first excited states, which are shown in Fig. 4. For convenience, we can parameterize the obtained numerical form factor σ as For example, we obtain A = 0.162, κ 1 = 0.109, κ 2 = 0.312 for the ground state.   III. Doubly heavy baryon as the bound state of the diquark core and a light quark In this section, we first construct the Bethe-Salpeter equation of a 1 + diquark core and a light quark, and then derive the 3-dimensional (Bethe-)Salpeter equation in the instantaneous approximation. At last, we compute mass spectra and wave functions for baryons with the quantum numbers J P = 1 2 + and 3 2 + . III.1. The Bethe-Salpeter equation of a baryon with the 1 + heavy diquark core The BS equation for a baryon with the 1 + heavy diquark core is schematically depicted by Fig. 5, which can be expressed using the matrix notation as where Γ α (P, q, r) is the BS vertex of the baryon; P , the baryon momentum with P 2 = M 2 and M denoting the baryon mass; r, the baryon spin state; (−i)K(p 1 , k 1 ; p 2 , k 2 ) represents the effective diquark-quark interaction kernel; and D βγ ( is the free diquark propagator (axialvector particle) with mass m 1 , and here we do not consider the self-energy correction of the doubly heavy diquark; the internal momenta q and k are defined as with α i = m i m 1 +m 2 and m 2 representing the constituent mass of the third quark. From now on, the symbols P and r in the BS vertex Γ α (P, q, r) would be omitted unless it is necessary.
In terms of the diquark form factor, the baryon effective interaction kernel can be simply expressed as In principle, the form factors must be calculated for off-shell diquarks. In this work, we simply generalize the on-shell form factors obtained in the previous section to the off-shell ones, namely, Σ αβµ = σ(t 2 )g αβ h µ and h = k 1 + p 1 . Notice that the diquark-quark potential is now modified by the form factor, which reflects the structure effect of the nonpoint-like diquark. The BS wave function B α (P, q) of the baryon can be defined as with the constraint condition P α B α (P, q) = 0. Accordingly, the normalization condition can be expressed as, where the operator I αβ (P, k, q) has the following form,

III.2. BSE of the doubly heavy baryons in instantaneous approximation
In the instantaneous approximation, the factor h 0 = 2(α 1 M +q P ) depends on q P and M explicitly, which is different from the case of mesons or diquarks. As we will see, the factor h 0 plays an important role in the derivation of the three-dimensional Salpeter equation of baryons. First, we introduce the instantaneous baryon kernel as and then have K αβ (p 1 , k 1 ) = h 0 g αβ K(q ⊥ − k ⊥ ). Then we define the baryon Salpeter wave function as with the constraint condition P α ϕ α = 0. Now the baryon BS wave function can be expressed as where the three-dimensional vertex Γ a (q ⊥ ) is expressed by the Salpeter wave function as Due to the constraint condition P α ϕ α (P ) = 0, the components of D αβ (p 1 ) parallel to P vanish. Hence we can rewrite D αβ (p 1 ) as Note that in Ref. [27], the term i is simply neglected. We follow Salpeter's method in Ref. [64]. Performing the contour integral over q P on both sides of Eq. (23) (see appendix A.2 for details), we eventually obtain the Schrödinger-type Salpeter equation for baryons with 1 + diquark cores as where the baryon mass M behaves as the eigenvalue; the first term expresses the contribution to the mass from the kinetic energy, and the second term that from the potential energy. Correspondingly, the normalization condition can be expressed as (see appendix A.3 for details) with ω q ≡ α 2 ω 1 − α 1 ω 2 and d αβ ≡ −g αβ − p α 1⊥ p α In this subsection, we construct the baryon wave functions directly from the good quantum number J P . Then the D-wave components are automatically included and the possible S-D mixing can also be obtained directly by solving the corresponding BSE.
For ground states with L = 0, where L denotes the diquark-quark orbital angular momentum, the 1 + diquark and the 1 2 + quark can form a baryon doublet ( 1 2 , 3 2 ) + . For the wave function of J P = 1 2 + baryons, analyzing its total angular momentum, Lorentz and Dirac structures, we can decompose it as ϕ α (P, q ⊥ , r) = g 1 + g 2 / q ⊥ q ξ 1α γ 5 u(P, r) + g 3 + g 4 / q ⊥ q q ⊥α γ 5 u(P, r) where we have defined A α to factor out the spinor; ξ 1α = (γ α + Pα M ), and q = −q 2 ⊥ ; the radial wave functions g i (i = 1, 2, 3, 4) explicitly depend on q; and u(P, r) is the Dirac spinor with momentum P and spin state r. The conjugation is defined as usualφ α (P, q ⊥ , r) = ϕ † α γ 0 . Note that in Ref. [34] only g 1 and g 2 are included. In fact, the full Salpeter wave function should have four independent terms, the 'small components' g 3 and g 4 correspond to part of the P -and D-wave components respectively, and play certain roles for the exact solutions of the J P = 1 2 + baryons, especially, they will become important for the excited states. To see the partial-wave components more clear, we rewritten the Salpeter wave function of 1 2 + in terms of the spherical harmonics as, where the coefficients are , and Γ αi ≡ (g 2 γ i ξ 1α + g 3 g αi ) with i = 1, 2, 3; g α± = ∓ 1 √ 2 (g α1 ± ig α2 ) with g αβ denoting the Minkowski metric tensor; Y m l is the usual spherical harmonics. From the above decomposition Eq. (28), we can conclude that the Salpeter wave function for 1 2 + baryon contains the S-, P -, and D-wave components, namely,: the g 1 part corresponds to the S-wave; g 2 and g 3 parts corresponds to the P -wave; and g 4 part makes contribution to both the Sand D-wave components. The numerical details will be presented in the next section.
With the specific wave function in Eq. (27), the Salpeter Eq. (25) can be expressed as Projecting the radial wave functions out, we arrive at four coupled eigenvalue equations. The details are presented in appendix A.4. Solving the equations numerically, mass spectra and radial wave functions can be obtained. Accordingly, the normalization condition Eq. (26) can be expressed in terms of radial wave functions as ; and the spinor summation r u(r)ū(r) = / P + M is used.
Similarly, for J P = 3 2 + baryons, the Salpeter wave function can be constructed on the basis of the Rarita-Schwinger spinor u α (P, r) as where we have defined tensor A αβ for later convenience; ξ 3α (P ) = (γ α − Pα M ); and t i (i = 1, · · · , 6) are radial wave functions depending on q. According to the Rarita-Schwinger equation, we have P α ϕ α = 0. Note that the full Salpeter wave function has six independent terms, compared with Ref. [34], where only the first two of them, i.e., t 1 and t 2 , are considered. The decomposition in terms of spherical harmonics is similar with the case in J P = 1 2 + and we would not repeat the details. It is evident that t 3 , t 4 , t 5 , and t 6 correspond to the P , 2 D, 4 D, and F partial waves respectively. Lacking of them means that only the largest S-wave and part of the P -wave components are considered. Notice that the D(F )-wave would always be mixed with the S(P )-wave components, namely, the t 4(5) would always contribute to the S-wave besides the D-wave components, and t 6 would also contribute to the P -wave besides the F -wave components. Further analysis of their importance will be presented later.
Now we obtain the following mass eigenvalue equation which is equivalent to six coupled eigenvalue equations similar with that for the 1 2 + case, and we will not copy details again. The normalization condition can be written as , and the following completeness relation of the Rarita-Schwinger spinor [65] has been used:

IV. Numerical results and discussions
For solving the Salpeter equation of baryons, we need to specify the parameter V 0 in the potential for the diquark core and the light quark. Unlike doubly heavy mesons, experimental data of doubly heavy baryons is very limited. So one cannot fix V 0 by fitting experimental data. Alternatively, here we compute V 0 by taking the spin-weighted average of the corresponding mesons' V 0 (see details in appendix A.5). The fixed values of V 0 are listed in Tab. III. Note that all the parameters in the model are now fixed by the meson sector.
In this work, specifically, we calculate Ξ

+26
The mass spectra for J P = 1 2 + doubly heavy baryons with flavors (ccq), (bcq) and (bbq) are presented in Tab. IV, including the diquark cores at the excited 2S and 1D states. Notice the baryon masses of 1 2 S 1/2 (1D) are comparable with that of the 1 4 D 1/2 (1S) state, namely, the different D-wave orbital angular momentum contribute similarly to the total baryon mass. The wave functions of Ξ ++ cc , as an example, are presented in Fig. 6. From the curves in the figures one may see the correspondence of the radial quantum number n B and the node number of the wave functions clearly; g 1 and g 4 correspond to the 2 S and 4 D partial wave components respectively; g 2 and g 3 represent the slight P -wave components with different radial number. Comparing the weights of the partial wave components, one can realize that the ground state n = 1 is dominated by the S-wave and has a negligible D-wave, whereas, the excited states n = 2, 3, 4 . . . have sizable D-wave components. It means that for excited states, to ignore the components g 3,4 would destroy the results. Therefore, as mentioned in the previous section, the full structures of the wave functions should be well considered.
The mass spectra for J P = 3 2 + doubly heavy baryons are presented in Tab. V, and the radial wave functions of Ξ * cc are shown in Fig. 7 as an example. For J P = 3 2 + baryons, besides the 4 S 3/2 partial waves, they also contain two different D-waves, i.e., 4 D 3/2 and 2 D 3/2 , which belongs to the quartet ( 7 2 , 5 2 , 3 2 , 1 2 ) + and the doublet ( 5 2 , 3 2 ) + , respectively. The baryon masses with the excited 2S and 1D diquarks are listed in the last two rows to make an comparison with the mass spectra with ground diquarks. Note that the masses of 1 4 S 3/2 (1D) states are close to those of 1 2 D 3/2 (1S) states, but much lower than those of the 1 4 D 3/2 (1S) states. The similar analysis on the wave functions shows that almost every state contains all the possible S-, P -, D-, and F -wave components. And the analysis on the weight of the partial waves suggests that the n = 1, 2 states are mainly characterized by 1 4 S and 2 4 S components, respectively, with D-wave mixed slightly. However, for the n = 3, 4 states, the 1 2 D and 1 4 D components become dominant, respectively. Again, these features emphasize the  Fig. 6: BS radial wave functions of the Ξ ++ cc with the energy level n = 1, · · · , 4; g 1 and g 4 correspond to the 2 S and 4 D components respectively, and g 2(3) corresponds to the P -wave; n B is one more than the number of the node; and almost every state contains the S-, P -, and D-wave components.  importance to involve the full structures of the wave functions.
To see the sensitivity of the mass spectra on the model parameters, we calculate the theoretical uncertainties by varying the parameters λ, and a 1 (2) in the quark-diquark bound potential by ±5% simultaneously, and then searching the parameter space to find the maximum deviation. The obtained theoretical errors are also listed in above Tab. IV and Tab. V, where the relative errors induced from λ and a 1 (2) are about 1% roughly. The dependence on the constituent quark mass is almost linear and not included in the error analysis.
Finally, we give a comparison between the predicted spectra of ground states and those in literatures, which is collected in Tab. VI. We notice that the mass spectra agree well with the lattice QCD calculations in Refs. [21,66,67], especially for Ξ  Fig. 7: BS radial wave functions of the Ξ * cc with the energy level n = 1, · · · , 4; t 1 , t 4 and t 5 correspond to the 4 S, 2 D and 4 D components respectively, t 2(3) and t 6 correspond to the P and F -wave respectively; and also almost every baryon state contains all the 4 S, 2 D and 4 D components.
of errors in predictions of Tab. VI with that in Ref. [21] is about 18 MeV.
Tab. VI: Comparisons of the predictions for the ground state masses (in GeV) of the doubly heavy baryons, where the results of the third [21] and forth [66,67] columns are from the lattice QCD calculations. Note that here Ξ bc , Ω bc denote the 1 2 + (bcq)-baryons with 1 + (bc)-diquark core only. In summary, within the diquark core picture, we have built a theoretical framework to deal with doubly heavy baryons. Based on the diquark picture, we reduce the three-body problem of the doubly heavy baryons into two two-body bound problems and then resolved by the corresponding BSE in the instantaneous approximation. We obtain the mass spectra and wave functions of J P = 1 2 + and 3 2 + doubly heavy baryons with 1 + diquark cores. Roughly speaking, the predicted mass spectra for the doubly heavy baryons agree with those in the literature, especially those by lattice QCD computation. Analyzing the wave functions, we found that 1 2 + and 3 2 + baryons, especially their excited states, involve unneglectable P -and D-wave components, and even F -wave for the latter, which means that the non-relativistic wave functions only consisting of S or D partial waves cannot work well to describe these baryons. In the near future, the obtained wave functions can be used to calculate lifetimes, productions, and decays of doubly heavy baryons, which can also be tested by (1 + cos 2 θ)(µ 1 e 2 + µ 2 e 1 ) , where cos θ = s· u su . These four equations can be solved numerically to obtain the diquark mass spectra and corresponding Salpeter wave functions.
A.2. Derivations of the baryon Salpeter equation with 1 + heavy diquark core To reach Eq. (25), firstly, we split out q P from the propagators, S(p 2 ) and D αβ (p 1 ) then can be expressed as, where the projector operators are defined as Λ ± (p 2⊥ ) = 1 2 1 ±Ĥ(p 2⊥ ) γ 0 ; ζ ± 1,2 are defined as, Performing the integral over q P on both sides of Eq. (23), we obtain the three-dimensional baryon Salpeter equation as By using the projector operators Λ ± , we can get the positive and negative energy Salpeter wave functions as, which are the coupled baryon Salpeter equations with 1 + diquark core, and then can be rewritten as the Schrödinger-type Eq. (25).

A.3. The normalization condition of the baryon Salpeter wave function
To obtain the normalization of ϕ α (q ⊥ ), we need the inverse of the propagators D αβ (p 1 ), which is given by, , and fulfills D −1 αγ (p 1 )D γβ (p 1 ) = d αγ ϑ γβ = δ β α . Note that d αβ (p 1⊥ ) does not explicitly depend on P 0 . The BS vertex can also be expressed by the inverse of the propagators as Γ α (q) = S −1 (p 2 )D −1 (p 1 )d αβ B β (q).