Masses of doubly heavy baryons in the Bethe-Salpeter equation approach

A doubly heavy baryon can be regarded as composed of a heavy diquark and a light quark. In this picture, we study the masses of the doubly heavy diquarkes in the Bethe-Salpeter (BS) formalism first, which are then used as one of the inputs in studying the masses of the doubly heavy baryons in the quark-diquark model. We establish the BS equations for both the heavy diquarks and the heavy baryons with and without taking the heavy quark limit, respectively. These equations are solved numerically with the kernel containing the scalar confinement and one-gluon-exchange terms. The mass of the doubly charmed baryon $\Xi_{cc}^{(\ast)}$ is obtained in both approaches, $3.60\sim3.65\,\rm GeV$ ($\Xi_{cc}^{(\ast)}$) under the heavy quark limit, $3.53\sim3.56\,\rm GeV$ for $\Xi_{cc}$ and $3.61\sim3.63\,\rm GeV$ for $\Xi_{cc}^\ast$ without taking the heavy quark limit. The masses of $\Xi_{bc}^\prime$, $\Xi_{bc}^{(\ast)}$, $\Xi_{bb}^{(\ast)}$, $\Omega_{cc}^{(\ast)}$, $\Omega_{bc}^\prime$, $\Omega_{bc}^{(\ast)}$ and $\Omega_{bb}^{(\ast)}$ are also predicted in the same way. We find that the corrections to the results in the heavy quark limit are about $-0.02\,\rm GeV\sim-0.11\,\rm GeV$ for the masses of the doubly heavy baryons.


I. INTRODUCTION
In the last few decades, heavy baryons which contain one heavy quark have been studied in many references, both theoretically and experimentally [1][2][3][4][5][6][7]. Whether the theoretical models for such baryons also work for baryons containing two or more heavy quarks needs to be investigated [6,[8][9][10][11]. Ξ ++ cc is the first doubly heavy baryon observed in experiments [12]. Although the SELEX Collaboration did report the observation of Ξ + cc at a mass of 3519 ± 2 MeV years ago [13,14], other collaborations like FOCUS [15], BaBar [16] and Belle [17], on the other hand, failed to observe such a state with the properties SELEX observed.
Recently, LHCb confirmed the existence of the doubly heavy baryon Ξ ++ cc , and the mass was measured as 3621.40 ± 0.72(stat) ± 0.27(syst) ± 0.14(Λ + c ) MeV [12]. It is expected that more doubly heavy baryons will be observed in the further, hence theoretical studies of doubly heavy baryons are urgent. In this paper, we will study the masses of the doubly heavy baryons systematically. Since baryons consist of three quarks they are much more complicated than mesons. Because bb, bc and cc 1 have good spin and isospin quantum numbers, we will regard them as diquarks in the baryon, where it can be treated as a single subatomic particle with which the third quark interacts via strong interaction. Using this quark-diquark model can greatly simplify the study of these baryons which are two-body bound states [18].
The bound state equation we employ in our work is the Bethe-Salpeter (BS) equation, which has been extensively used to investigate the properties of heavy mesons and baryons in recent years [7,[19][20][21][22][23][24][25][26][27][28][29]. The kernel for the BS equation is motivated by the potential model, which includes the scalar confinement and one-gluon-exchange terms in the form of the covariant instantaneous approximation [7,25,27,29,30], and this approximation is appropriate since the energy exchange between constituents inside the heavy baryons is expected to be small when we use the constituent quark masses in the BS equation. Thus we employ this approximation in this work for the purpose of simplifying our calculations. It should be noted that solving the BS equation in the Minkowski space has been done before by other researchers [31,32], and we will also try to accomplish that by modifying our model in our future work. We will establish BS equations for doubly heavy diquarks first, then the BS equations for doubly heavy baryons. The masses of heavy diquarks will be solved out first and then used as an input to obtain the masses of doubly heavy baryons.
The masses of heavy baryons are studied with different approaches [6,10,21,[33][34][35][36][37], and most of them are done under the heavy quark limit, which works well for the doubly bottom diquark and doubly bottom baryons since m b is much larger than the QCD scale, Λ QCD . However, one may expect lager 1/m c corrections for heavy diquarks and baryons including the charm quark. Therefore, we will also establish the BS equations for doubly heavy diquarks and doubly heavy baryons without taking the heavy quark limit. We will compare the results in the two approaches. This paper is organized as follows, In Sec. II, the BS equations for the heavy diquarks will be presented under the heavy quark limit. In Sec. III, we will establish the BS equations for the ground-state doubly heavy baryons under the heavy quark limit. Then, we will present our study for the heavy scalar and axial-vector diquarks without using the heavy quark limit, the formalism of which will be presented in Sec. IV, After that, the formalism for the doubly heavy baryons without taking the heavy quark limit will be presented in Sec. V.
Then in Sec. VI, we will show our results derived from the former sections. Finally, the summary will be given in Sec. VII.

A. The BS equations for doubly heavy diquarks
In the following, we will present the BS equations for the scalar diquark composed of c and b quarks, and axial-vector diquarks composed of cc, bb and bc. First, the BS wave function for the scalar diquark composed of two fermions can be written as follows: where ψ 1 (x 1 ) and ψ 2 (x 2 ) are the quark fields, α and β are spinor indices, i, j, k are color indices, P D is the momentum of the diquark, X = λ 1 x 1 + λ 2 x 2 (λ 1 = m 1 m 1 +m 2 , λ 2 = m 2 m 1 +m 2 , m 1 (2) is the mass of the first (second) heavy quark) is the center-of-mass coordinate of the diquark, x = x 1 − x 2 is the relative coordinate and p represents the relative momentum between the two heavy quarks. Furthermore, the momenta of two heavy quarks can be written with P D and p as p 1 = λ 1 P D + p, p 2 = λ 2 P D − p, respectively. Define p l = v · p and p t = p − p l v, where v is the velocity of the diquark. The BS equation for the scalar diquark has the following matrix form [20,21]: (2) where S(p 2 ) and S(p 1 ) are the propagators for two heavy quarks, and χ P D (p) = Cχ P D (p) is the conjugate form of χ P D (p), K (1g) and K (cf ) denote the one-gluon-exchange and confinement terms, respectively, which have the following forms in the diquark picture: where α s and κ are coupling parameters in the meson case, and µ is introduced to avoid the infrared divergence in later numerical calculations, the limit µ → 0 will be taken in the end of our calculations. As for the forms of propagators, they can be constructed as follows in the leading order of 1/m Q (Q = b, c) expansion: where M D is the mass of the scalar diquark and ω 1(2) = m 2 1(2) − p 2 t . It can be seen from Eqs. (2) and (6) Then, after considering the relation in Eq. (7) and other restrictions from Lorentz convariance and parity transformation on the form of χ T P D , we obtain the following form for the BS wave function of the scalar diquark: where g 1 is a scalar function of p 2 t and p l . Defining V (1g) = −i K (1g) , V (cf ) = −i K (cf ) , g 1 = dp l 2π g 1 , and substituting Eqs. (5)-(8) into Eq. (2), we obtain the following equation for the BS wave function of the scalar diquark under the heavy quark limit: On the other hand, for the form of the BS wave function of the axial-vector diquark, we have where ξ (r) µ is the polarization vector of the axial-vector diquark, g 2 is a scalar function of p 2 t and p l . Following the same technique in the scalar diquark case, we have the following equation for the BS wave function of the axial-vector diquark under the heavy quark limit: which shows that the BS wave function for the axial-vector diquark satisfies the same equation as that for the scalar diquark when taking the heavy quark limit.

III. DOUBLY HEAVY BARYONS UNDER THE HEAVY QUARK LIMIT
In this section, we will construct the formalism for the doubly heavy baryons under the heavy quark limit.
A. Baryons composed of a heavy scalar diquark and a light quark The BS wave functions of baryons composed of a heavy scalar diquark and a light quark, like Ω bc and Ξ bc , can be defined as where ψ l (x 1 ) and φ D (x 2 ) are the field operators of the light quark and the scalar diquark, , m l and m D are the masses of the light quark and the heavy diquark, respectively.) is the center-of-mass coordinate of the baryon and x is the relative coordinate, P is the total momentum of the baryon, while k is the relative momentum between the heavy scalar diquark and the light quark. Then, k 1 = η 1 P +k, k 2 = η 2 P −k are the momenta of the light quark and the heavy scalar diquark, respectively. For the BS equation of the doubly heavy baryon, we have where G(P, k, q) is the interaction kernel, S l and S D are the propagators for the light quark and scalar heavy diquark, respectively, which can be expressed as the following forms under the heavy quark limit: where we define the variables and M is the mass of the doubly heavy baryon. For the kernel G(P, k, q), we will take the form that has been used in Ref. [21]: where the vertex Γ µ has the form: describing the interaction between diquarks and the gluon (Q 0 is a parameter which freezes F S (Q 2 ) as Q 2 → 0), and V 1 and V 2 are the scalar confinement and one-gluon-exchange terms that have the following forms in the covariant instantaneous approximation: where α sef f and κ are the coupling parameters in the baryon case, and µ is introduced to avoid infrared divergence in numerical calculations. It should be noted that the confinement parameter in the baryon case is different from that in the diquark or the meson cases, so we replace κ with κ. Since κ is around 0.2 GeV 2 in the meson case [27,30], and considering Λ QCD is the only parameter related to confinement, we expect κ around 0.04 GeV 3 due to the relation κ ∼ Λ QCD κ [7].
One the other hand, after writing down all the possible Dirac structures and considering the Lorentz and P-parity transformations, we have the following form for χ P (k): where u(v, s) is the Dirac spinor of the doubly heavy baryon with the velocity v and the helicity s, which satisfies the relation v /u(v, s) = u(v, s), and h 1(2) is a scalar function with respect to k 2 t and k l . Substituting Eqs. (14), (15), (16) and (19) into Eq. (13), we obtain the coupled integrating equations as follows, where we have defined h 1(2) = dq l 2π h 1 (2) , and V 1(2) = V 1(2) | k l =q l . The following equations based on lorentz invariance have been used while obtaining the above two coupled integrals: B. Baryons composed of a heavy axial-vector diquark and a light quark Now we consider the doubly heavy baryons composed of a heavy axial-vector diquark and a light quark. First, we give the definition of the BS wave function for such a baryon: where A µ D denotes the axial-vector field and the BS equation for the baryon in momentum space has the following form: where S µρ D (k 2 ) is the propagator of the heavy axial-vector diquark, which can be written as the following form under the heavy quark limit (m D → ∞): where O(1/m D ) terms in the propagator are omitted due to the heavy quark limit. On the other hand, motivated by the potential model [30], the following form for the kernel G ρν will be used [21] iG ρν (P, k, q) = g ρν I ⊗ where is the coupling vertex of heavy diquarks and a gluon with F A (Q 2 ) being the form factor describing the effect of the diquark structure, we will take F A (Q 2 ) = F S (Q 2 ) for simplicity [38,39]. Meanwhile, we have the following form for χ µ P (k) under the heavy quark limit: where b 1(2) is a scalar function of k 2 t , q 2 t and k t · q t , B µ (v) is the spinor of the doubly heavy baryon, which has the following explicit forms: for spin-1 2 and spin-3 2 baryons, respectively, u(v) denotes the Dirac spinor, and u µ (v) is the Rarita-Schwinger spinor. Furthermore, it can be shown that B µ m (m = 1, 2) satisfies the conditions [19]: Substituting Eqs. (14), (27), (28) and (29) into Eq. (26), and integrating q l in the integral equations, the terms like (2π) 4 q µ t q ν t can be transformed to the desired forms with the help of Eqs. (22)- (24). After all of these steps and some calculations, we arrive at: In this section, we will present the BS equation for a scalar diquark without taking the heavy quark limit m Q → ∞. To obtain the form of the BS wave function for the scalar diquark with the masses of heavy quarks being finite, we first write down the general decomposition of the following Dirac fields [20]: Then from the definition of the BS wave function of the scalar diquark, and imposing the P-parity transformation on the BS wave function corresponding to each term on the righthand-side of Eq. (34), we have where d i (i = 1, 2, 3, 4) are Lorentz scalar functions of p 2 t and p l . Now, we will check if these four functions are independent of each other. To achieve this goal, we introduce the projecting operators [40]: It can be shown that these projecting operators satisfy the conditions: With these operators, we can express the full propagators of quarks as follows: Besides, we have the following equations for the projecting operators: and Now, by multiplying a certain combination of (39) and (40) on both sides of Eq. (2), one obtains the following constrains: which may help reduce the number of independent scalar functions in Eq. (35). We find . However, when the masses of the two quarks in the diquark are different, we find no relation among these four scalar wave functions. Substituting Eqs. (35), (37), (38), K (cf ) = iV (cf ) and K (1g) = iV (1g) into Eq. (2), and completing the integration dq l 2π on both sides of Eq.
(2), we obtain four coupled integral equations: where we have defined d i (p t ) = . It can be shown that Eqs. (43) and (44) are consistent with the equations in Ref. [20] when the masses of the two quarks are equal.

B. The BS equation for doubly axial-vector diquarks
In this subsection, we will present the formalism for the axial-vector diquark. The general form for the BS equation of the axial-vector diquark is similar to that of the scalar diquark: where r stands for the index of the polarization vector of the axial-vector diquark. The BS equation for the axial-vector diquark in momentum space has the form: where Γ ρ indicates the Dirac structure of χ (r)T P D , and ε (r) ρ is the r − th polarization state of the axial-vector diquark, which satisfies the relation ε (r) ρ P ρ D = 0. In general, Dirac field operators for the axial-vector diquark can be decomposed into several terms like the scalar diquark case, and the BS wave function corresponding to each term in Eq. (34) can be further expanded as a linear combination of p t , P D and g µρ . After writing down all the terms contained in Γ ρ , we use the constrains ε where H i (i = 1, · · · , 8) are scalar functions. When the masses of the two heavy quarks are equal, we follow the same routine used in the scalar diquark case to reduce the number of the scalar functions. Consequently, we have the constrains for these scalar functions H i (i = 1, · · · , 8): where the masses of the two quarks in the diqaurk are equal (m 1 = m 2 = m). Therefore, in this case, in the eight scalar functions there are only four independent wave functions H i .
With the help of the constrains in Eq. (50), we are able to obtain four integral equations for scalar functions in the case of the axial-vector diquark: where H i (p t ) = dq l 2π H i (p). On the other hand, for the doubly heavy axial-vector diquark composed of two quarks with different masses, e.g., bc with J P = 1 + , we also have similar relations for the scalar functions in the BS wave function for the axial-vector diquark bc without taking the heavy quark limit: In this case, we can follow the same routine used before, substituting Eqs. (55)-(58) into Eq. (49), which is then substituted into the BS equation for the axial-vector diquark. We will not present the explicit forms of the integral equations for the scalar functions here because the calculations are long and demanding, we will present the numerical results for the axial-vector diquark bc in Sec. VI instead.
It is noteworthy that there are four independent scalar functions in both scalar and axialvector diquark cases, which is consistent with the physical picture since each of these two quarks has two degrees of freedom. In contrast, there is only one independent scalar function in the scalar and axial-vector diquark cases when the heavy quark limit is applied. In this subsection, we continue to use the variables defined in Sec. III. First, we can start from writing down the general form of the BS equation for baryons comprised of a heavy scalar diquark and a light quark, which has the same form as Eq. (13): except that for the propagator of the scalar diquark we will not take heavy quark limit. In this case, we have the following form for the propagator of the heavy scalar diqaurk: For the form of the BS wave function χ P (k), after considering the constrain on the Dirac spinor of the baryon u(v, s) (v /u(v, s) = u(v, s)) we write down all the possible combinations of Dirac structures and momentum (P, k t ), where Y i (i = 1 · · · 5) are scalar wave functions of k 2 t , k l and k t ·q t . Furthermore, each term in the expansion of χ P (k) transforms exactly in the way that χ P (k) transforms under P-parity and Lorentz transformations, which can help us simplify the form of χ P (k). Finally, we arrive at following the same routine we have used in the last few sections. A 1 and A 2 are scalar functions of k 2 t , k l and k t · q t . Define A 1(2) = dq l 2π A 1 (2) . Then, after we substitute Eqs. (14), (16), (60) and (62) into Eq. (59), we obtain two coupled equations for the scalar functions for the baryons composed of a doubly scalar diqaurk and a light quark as follows: It can be seen from these two equations that the corrections are generated from both the V 1 and V 2 terms compared with the situation under the heavy quark limit, which will be the reason why the masses calculated from Eqs. (63) and (64) are different from those obtained through Eqs. (20) and (21). One can prove that Eqs. (20) and (21) can be obtained by taking the heavy quark limit (m Q → ∞).

B. Baryons composed of a heavy axial-vector diquark and a light quark
The BS equation formalism for the baryons comprised of a heavy axial-vector diquark and a light quark without taking the heavy quark limit will be established in this subsection.
In this situation, the degenerate states like {Ξ cc , Ξ * cc } in the heavy quark limit will spilt into two states with different masses: spin-1 2 and spin-3 2 . We will construct different forms of BS wave functions for spin-1 2 and spin-3 2 states of these doubly heavy baryons, respectively. The most difficult thing in deriving the BS equation for this type of baryons, like Ξ ++ cc , is to find the correct form of the wave function for them without the help from the heavy quark symmetry. In the heavy quark limit, the internal dynamics is determined by the light degrees of freedom due to the SU (2) s × SU (2) f symmetry, then we can easily decompose the wave function into the following form [19]: where ψ and A µ represent the light quark and the axial-vector diquark fields, respectively, u(v) denotes the Dirac spinor of the light quark, ξ ν stands for the polarization vector of A µ , and u(v)ξ ν can be decomposed into spin-1 2 and spin-3 2 baryons B mν (m = 1, 2). The tensor ζ µν includes all the dynamics and contains two independent scalar functions in Eq. (29).
When the masses of the heavy quarks are finite the form of Eq. (65) still holds due to the LSZ reduction formula [41]. However, ζ µν , which includes internal dynamics, becomes more complicated in this scenario and will be different for spin-1 2 and spin-3 2 doubly heavy baryons. After taking into account all the possible Dirac structures and the combinations of Dirac gamma matrices and momenta (k t , P ), and imposing the constrains from P-parity and Lorentz transformations we find that the BS wave function for all the spin-1 2 doubly heavy baryons, χ µ P (k) 1/2 , can be written as the following form: where the scalar functions E i (i = 1, · · · , 6) may have different dimensions. As before, we define E(k t ) = dk l 2π E(k), i = 1, · · · , 6.
As for the propagator of the axial-vector diquark, we need its full form, It can be seen that in comparison with the propagator under the heavy quark limit, Eq. (27) [7,19,29,43], where Γ µρν is the vertex of a gluon with two axial-vector diquarks, which has the following form: where g s is the coupling constant for strong interaction in the vertex, and F A (Q 2 ) is the form factor describing the inner structure of the axial-vector diquark.
Using the same techniques as in the last few sections, we obtain coupled integral equations for E i (i = 1, · · · , 6) for spin-1 In the same way, considering all the possible Dirac structures, the combinations of Dirac gamma matrices and momenta (k t , P ), and the restrictions from P-parity and Lorentz transformations, we obtain the explicit form of BS wave function for spin-3 2 doubly heavy baryons, which has the following form: where C i (i = 1, · · · , 8) are scalar functions of k 2 t , q l and k t · q t , which correspond to eight different Dirac structures. Using the same techniques as before, we obtain the coupled

VI. NUMERICAL RESULT
In this section, we will solve all the coupled integral equations emerged in previous sections numerically. The integration region in each integral will be discretized into n pieces, with n being sufficiently large. In this way, the integral equation will be converted into an n × nmatrix equation, and the scalar wave functions of each equation will now be regarded as an n-dimensional vector. It can be shown that the integral equations, for instance, Eqs. (A1)-(A6), can be illustrated as: where E i (i = 1, · · · , 6) is an n-dimensional vector, and e ij (i, j = 1, · · · , 6) is an n × n matrix function of k t and q t . To make it more clear, we will take Eq. (9) as an example here.
In Eq. (9), the integral equation can be rewritten in the following way once we integrated the angular parts, where |p t | and |q t | are the absolute values for the relative momentum between quarks and diquarks in the initial and final states, respectively. Since we choose to integrate |q t | in Eq. (74) with the Gaussian quadrature rule, we need to convert the Gaussian integration nodes into physical values for |q t |, which can be done using the following equation where t is the ordinary Gaussian integration nodes ranging from −1 to 1, is a parameter introduced to avoid divergence in numerical calculations, w and y are parameters used in controlling the slope of wave functions and finding the proper solutions for these functions.
One can then obtain the numerical results of the g 1 (p t ) by requiring the eigenvalue of the eigenvalue equation to be 1. Similar methods can be applied to evaluate other integral equations in our work.

A. Doubly heavy diqaurks
In this subsection, the numerical results based on the theoretical framework of doubly heavy diquarks will be presented. From the potential model and some former works on the BS analysis, κ is usually taken to be around 0.20 GeV 2 for the diquark system [7,20,21,27,42]. For example, in Ref. [21], the authors used κ = 0.20 GeV 2 to obtain the masses of heavy diquarks. In the present work, we will let κ vary in a reasonable range, 0.18 − 0.20 GeV 2 , to see to what extent it will affect the masses of the doubly heavy diquarks cc, bb and cb. In our work, we will employ the values of coupling parameters (α s and κ ) from Ref. [43] and the heavy quark masses (m b = 5.02 GeV, m c = 1.58 GeV) from Ref. [27] for the calculation of the masses of the doubly heavy diquarks. These parameters, e.g., κ and α s were obtained through the investigations of heavy mesons and heavy quarkonia and led to results consistent with the experimental data [43,44].
In Table I, we present the results for the masses of diquarks under the heavy quark limit, and their counterparts with finite masses for heavy quarks. It can be seen from the table that the values under the heavy quark limit are generally bigger than their counterparts,  to the masses of diquarks show a mild increase with the increase of κ . Furthermore, the data obtained in the heavy quark limit is less independent on the parameter κ . In contrast, the masses of diquarks calculated without taking the heavy quark limit remain almost unchanged with the variation of the parameter κ . On top of that, we find that the corrections to the masses of heavy diquarks are all negative, which suggests that the contribution from 1/m Q and higher order terms cause a reduction in the masses of the diquarks.
In Figs. 1, 2    are quite similar, but the scalar wave functions of the diquarks decrease more rapidly as the diquark mass increases. This is because the radius of the diquark decreases with the increase of the diquark mass, leading to a larger average transverse momentum within the diquark [23]. Meanwhile, we plot the scalar wave functions for the scalar diquark bc in Fig. 4, where the scalar functions decrease more rapidly compared with the axial-vector bc plotted   in Fig 2. In Fig. 5, we plot the BS scalar wave functions of cc, bc and bb diquarks under the heavy quark limit, where the shapes of the wave functions are steeper when the masses of the diquarks are smaller due to the same reason as that in the case of finite heavy quark masses [23]. It should be noted that the axial-vector and scalar bc diquarks have the same BS scalar functions as shown in Eqs. (9) and (11). Obviously, all the wave functions decrease to zero when |p t | is larger than about 3.0 GeV, which is resulted from the confinement interaction.

B. Doubly heavy baryons
Now, we will use the obtained results for masses of the doubly heavy diquarks in last subsection to calculate the masses of the doubly heavy baryons in two scenarios: m Q → ∞ and m Q being finite. First, we will find the values of the coupling parameters κ and α sef f , and this is achieved by fitting the experimental data of baryons in these two scenarios (with and without taking the heavy quark limit), respectively. We take m s = 0.458 GeV, m u,d = 0.350 GeV in fitting the experimental data [27].
In the baryon case, the interaction kernel is different from that in the meson case. In general, the confinement term is still due to scalar interaction, thus the form of V 1 is still the same as that in the meson case, only the confinement parameter κ needs to be replaced with κ, which is the confinement parameter in the baryon case. We take κ to be in the range 0.02 GeV 3 ∼ 0.08 GeV 3 due to the relation between κ and κ (κ O(Λ QCD )κ ) [7]. Meanwhile, since the diquark is not a point-like particle, we introduce the form factor to describe the vertex between two diquarks and the gluon, where Q 2 0 is a parameter which freezes F (Q 2 ) when Q 2 is small. By analyzing the electromagnetic form factors of the proton [23,45], Q 2 0 = 3.2 GeV 2 can result in a good agreement with the experimental data [38].
There exists a constraint between the parameters κ and α sef f when we solve the integral equations numerically. In this way, we can obtain the corresponding value of α sef f for each κ in the range of 0.02 ∼ 0.08 GeV 3 . By fitting the experimental data for the masses of Λ Q , Ξ Q , Σ Q [46]. In this way, we can obtain the proper α sef f values, and we assume these parameters obtained from baryons composed of a light diquark and a heavy quark can also be applied to doubly heavy baryons since the strong interaction is flavor independent [7,19,21]. Following the same techniques in obtaining the values in Table II, we have the values of α sef f for heavy baryons under the heavy quark limit in Table III.  With the numerical values for the parameter Q 2 0 in the form factor and the masses of diquarks cc, bb and bc obtained before, there is only one parameter, κ, left free in our model for doubly heavy baryons. By solving the integral equations for each doubly heavy baryon numerically, we obtain the results for the masses of doubly heavy baryons under the heavy quark limit and the results for spin-1/2 and spin-3/2 baryons respectively, without taking the heavy quark limit. The results are illustrated in Table IV when κ varies between 0.02 and 0.08 GeV 3 . In Table IV, our results are different from the previous results in Ref. [21] due to different form factors applied. Furthermore, the coupling parameters κ and α sef f are obtained by fitting the known experimental data, and the masses of diquarks are also obtained within the BS formalism. The mass of Ξ cc is obtained as 3.60 ∼ 3.65 GeV under the heavy quark limit, which is in a good agreement with the experimental data 3621.40 ± 0.72 ± 0.27 ± 0.14 MeV [12]. Meanwhile, the mass for Ξ cc obtained without taking the heavy quark limit is 3.53 ∼ 3.56 GeV, which is about 0.10 GeV smaller than the result obtained under the heavy quark limit and also consistent with the experimental data. As for Ξ * cc , which is degenerate with Ξ cc under the heavy quark limit, we obtain 3.61 ∼ 3.63 GeV for its mass without taking the heavy quark limit. Considering the big mass difference between the data obtained by SELEX (3519 ± 2 MeV for Ξ + c ) and LHCb (3621.40 ± 0.72 ± 0.27 ± 0.14 MeV for Ξ ++ cc ) (where the masses of those two states are supposed to be very close due to very small isospin split between Ξ + c and Ξ ++ cc ), and based on our theoretical results listed in Table IV, we therefore suggest that those two states obtained in experiments may be spin-1/2 and spin-3/2 states of Ξ cc , respectively. Of course, this will need further confirmation from the measurement of spin-parity of Ξ cc .
In general, we can see from Table IV that corrections to the results in the heavy quark limit gradually decrease as the masses of doubly heavy baryons increase. For example, the mass of Ω cc is calculated to be about 3.73 ± 0.02 GeV under the heavy quark limit and 3.64 ± 0.01 GeV without taking the heavy quark limit, which gives the correction of about 0.10 GeV. However, in the case of Ω bb or Ξ bb , the corrections are basically 0.03 ∼ 0.06 GeV.
It is consistent with the expectation of the heavy quark effective theory and shows the importance of corrections to those in the heavy quark limit when dealing with doubly heavy baryons containing c-quark. On the other hand, it can be seen from the results that in both scenarios the results have little dependency on the coupling parameter κ, and the results obtained under the heavy quark limit are generally larger than those in the case of m Q being finite, which is similar to the results we obtained for the doubly heavy diquarks in the previous subsection. bb under the heavy quark limit with κ fixed at 0.02 GeV 3 . It can be seen form these figures that wave functions of Ξ ( * ) cc decrease more rapidly than those of Ξ * bb , which is quite similar to the situation for heavy diquarks. Also, it can be seen form the figures that when the heavy quarks have finite masses (Figs. 8 and 9) the curves of Ξ cc are clearly steeper than those of Ξ bb due to the mass difference between b and c quarks. In Figs. 8 and 9, we can see some of the curves, for example, E 6 (|p t |) (about 10 5 times bigger than E 3 (|p t |) at zero relative momentum) of Ξ cc is significantly bigger in the vicinity of zero relative momentum and decreases much more rapidly than the rest of the scalar functions, and six scalar functions ( C 3 -C 8 ) in Fig. 9 are extremely large when |p t | approaches zero. The exact values of these scalar functions are given below the figures. In Figs. 10 and 11, we also plot the scalar functions for Ξ bb and Ξ * bb without taking the heavy quark limit, which are similar to but broader than those in Figs. 8 and 9 due to smaller radii of Ξ bb and Ξ * bb .   In our calculations where the heavy quark has finite mass, Ξ * cc lies about 70 ∼ 80 MeV above Ξ cc , which is very close to the splitting between Σ c and Σ * c (64.44 MeV [46]). This is FIG. 8. The BS wave functions of Ξ cc (spin-1/2) without taking the heavy quark limit when κ = 0.02 GeV 3 ( E 6 (|p t |) is 6305 when |p t | is 0).
FIG. 10. The BS wave functions of Ξ bb (spin-1/2) without taking the heavy quark limit when consistent with the Gell-Mann−Okubo (GMO) mass relation [47,48] The similar relation is the following with u, d quarks being replaced by s quark: From Table IV, we have the splitting around 60 ∼ 70 MeV on the left-hand side, which is also relatively close to the one between Ω * c and Ω c (70.7 MeV). From the splitting information extracted here, we are able to analyze some decay possibilities. For example, since the splitting between two doubly charmed baryons Ξ * cc and Ξ cc (and Ω * cc and Ω cc ) is not big enough for the π meson emission of Ξ * cc (Ω * cc ), Ξ * cc (Ω * cc ) can only transfer to Ξ cc (Ω cc ) via γ emission.
Similarly, we have the following relation for doubly bottomed baryons: The splitting on the right-hand side of Eq. (78) is 20.2 MeV [46], in comparison to our result, about 10 ∼ 20 MeV, on the left-hand side, which is very close to the one on the FIG. 11. The BS wave functions of Ξ * bb (spin-3/2) without taking the heavy quark limit when κ = 0.02 GeV 3 ( C 3 (|p t |), C 4 (|p t |), C 5 (|p t |), C 6 (|p t |), C 7 (|p t |) and C 8 (|p t |) are 0.3, 3274, 4687, 21, 2 and 17 respectively when |p t | is 0). right. Moreover, it is noted that the splitting between doubly bottomed baryons is much smaller than that between doubly charmed baryons, which is due to the fact that b quark is much heavier than c quark, and the hyperfine splitting is reciprocal to the mass of the heavy quark.

VII. SUMMARY AND DISCUSSION
It is expected that more doubly heavy baryons will be observed in the future experiments, so theoretical studies for doubly heavy baryons are urgent at present. In this paper, we studied the masses of doubly heavy baryons in the BS formalism. We derived the BS equations for different diquark and baryon systems in the covariant instantaneous approximation, with and without taking the heavy quark limit, respectively. Then we discretized the integral equations and solved the eigenvalue equations numerically with the kernel containing the color confinement and one-gluon-exchange terms. We first calculated the masses of the doubly heavy diquarks cc, bc and bb, and gave the corrections to the masses in the heavy quark limit. After that, by fitting the experimental data for some known heavy baryons, we obtained the proper values for the parameter α sef f in the two scenarios (with and out taking the heavy quark limit) corresponding to κ ranging from 0.02 GeV 3 to 0.08 GeV 3 , where the masses of the light diquarks (uu, ud, ss) were calculated from BS equations for them and were taken as the inputs during the fitting.
We then used the masses of the doubly heavy diquarks and the coupling parameters as the inputs to calculate the masses of doubly heavy baryons. We obtained the masses of Ξ ( * ) cc , which are 3.60 ∼ 3.65 GeV in the heavy quark limit, 3.53 ∼ 3.56 GeV for Ξ cc and 3.61 ∼ 3.63 GeV for Ξ * cc when c quark has finite mass. These results are consistent with the data from SELEX (3519 ± 2 MeV for Ξ + cc ) and LHCb (3621.40 ± 0.72 ± 0.27 ± 0.14 MeV for Ξ ++ cc ) [12][13][14]. Furthermore, with the increase in the masses of the doubly heavy baryons, we found the corrections to the results in the heavy quark limit become smaller. Besides, we found the results obtained with m Q being finite depend hardly on parameters and less sensitive to the variation of α sef f in our model for the masses of both diquarks and baryons, and because of this, our predictions for these doubly heavy baryons could help locate more doubly heavy baryons in the experiments with more accuracy. On top of that, in the situation where the heavy quarks have finite masses, we obtained the hyperfine splittings between spin-1/2 and spin-3/2 doubly heavy baryons, which satisfy the GMO mass relations quite well.