Isgur-Wise function in $B_c$ decays to charmonium with the Bethe-Salpeter method

TL;DR

This paper calculates the Isgur-Wise function for $B_c$ decays to charmonium states using the Bethe-Salpeter method, highlighting the importance of relativistic corrections for accurate form factor descriptions.

Contribution

The paper introduces a method to extract the Isgur-Wise function for $B_c$ decays to various charmonium states using the Bethe-Salpeter approach, emphasizing relativistic effects.

Findings

Relativistic corrections significantly impact the Isgur-Wise function.

The method provides numerical results for multiple charmonium final states.

Overlapping integrals of wave functions are used to derive the Isgur-Wise function.

Abstract

The Isgur-Wise function vastly reduces the weak-decay form factors of hadrons containing one heavy quark. In this paper, we extract the Isgur-Wise functions from the instantaneous Bethe-Salpeter method, and give the numerical results for the $B_c$ decays to charmonium where the final states include $1S$, $1P$, $2S$ and $2P$. The overlapping integral of the wave functions for the initial and final states is the Isgur-Wise function, as the heavy quark effective theory does. In the case of accurate calculation, describing form factors need to introduce more relativistic corrections which are the overlapping integrals with the relative momentum between the quark and antiquark to Isgur-Wise function. The relativistic corrections to Isgur-Wise function provide greater contributions especially involving the excited state, and therefore are necessary to be adopted.

Tables

Table 1. Table 1: The normalization coefficients of different processes.

	final state	${\eta }_c$	$J/\psi$	$h_c$	${\chi }_{c0}$	${\chi }_{c1}$	${\chi }_{c2}$
	$C$	$4\sqrt{M{M}^{\prime }}$	$4\sqrt{M{M}^{\prime }}$	$\frac{4M}{\sqrt{3}}$	$4M$	$4\sqrt{\frac{2}{3}}M$	$\frac{4M{M}^{\prime }}{\sqrt{3}}$

Equations

$$\frac{\mathrm{d}\Gamma (\bar{B}\to D\ell \bar{\nu })}{\mathrm{d}\omega }=\frac{G_F^2}{48{\pi }^3}|V_{cb}{|}^2(m_B+m_D{)}^2{m}_D^3({\omega }^2-1{)}^{3/2}{\xi }^2(\omega ).$$

$$T=\frac{G_F}{\sqrt{2}}{V}_{cb}{\bar{u}}_{{\nu }_{\ell }}{\gamma }^{\mu }(1-{\gamma }_5)v_{\ell }⟨(c\bar{c})(P_f)|J_{\mu }|B_c^{+}(P)⟩,$$

$$\langle(c\bar{c})|\bar{b}\gamma^{\mu}(1-\gamma^{5})c|B_{c}^{+}\rangle=\int\frac{\mathrm{d}\vec{q}}{(2\pi)^{3}}{\rm Tr}\bigg{[}\overline{\varphi}_{P_{f}}^{++}(\vec{q}\>^{\prime})\frac{\not{P}}{M}\varphi_{P}^{++}(\vec{q}\>)\gamma^{\mu}(1-\gamma^{5})\bigg{]},$$

$$\langle(c\bar{c})|\bar{b}\gamma^{\mu}(1-\gamma^{5})c|B_{c}^{+}\rangle=\int\frac{\vec{q}\>^{\prime 2}\mathrm{d}|\vec{q}\>^{\prime}|\mathrm{d}(\cos\theta)}{(2\pi)^{2}}{\rm Tr}\bigg{[}\overline{\varphi}_{P_{f}}^{++}(\vec{q}\>^{\prime})\frac{\not{P}}{M}\varphi_{P}^{++}(\vec{q}\>)\gamma^{\mu}(1-\gamma^{5})\bigg{]},$$

$$⟨P|\bar{b}{\gamma }^{\mu }(1-{\gamma }^5)c|B_c^{+}⟩=S_{+}(P+P_f{)}^{\mu }+S_{-}(P-P_f{)}^{\mu },$$

$$⟨V|\bar{b}{\gamma }^{\mu }(1-{\gamma }^5)c|B_c^{+}⟩=(t_1{P}^{\mu }+t_2{P}_f^{\mu })\frac{ϵ\cdot P}{M}+t_3(M+M_f){ϵ}^{\mu }+\frac{2t_4}{M+M_f}\mathrm{i}{\epsilon }^{\mu \nu \sigma \delta }{ϵ}_{\nu }{P}_{\sigma }{P}_{f\delta },$$

$$⟨T|\bar{b}{\gamma }^{\mu }(1-{\gamma }^5)c|B_c^{+}⟩$$

$$+\frac{2t_4}{M+M_f}\mathrm{i}{\epsilon }^{\mu \beta \sigma \delta }{ϵ}_{\alpha \beta }\frac{P^{\alpha }}{M}{P}_{\sigma }{P}_{f\delta },$$

$${ϵ}_{\mu }^{(\lambda )}(P_f)P_f^{\mu }$$

$$\sum _{\lambda =0,\pm }{ϵ}_{\mu }^{(\lambda )}(P_f){ϵ}_{\nu }^{†(\lambda )}(P_f)$$

$${ϵ}_{\alpha \beta }^{(\lambda )}(P_f)P_f^{\alpha }$$

$$\sum _{\lambda =0,\pm 1,\pm 2}{ϵ}_{\mu \nu }^{(\lambda )}(P_f){ϵ}_{\alpha \beta }^{†(\lambda )}(P_f)$$

$$\Gamma =\frac{1}{8M(2\pi {)}^3}\int \frac{|{\vec{P}}_{\ell }|}{E_{\ell }}\mathrm{d}|{\vec{P}}_{\ell }|\int \frac{|{\vec{P}}_f|}{E_f}\mathrm{d}|{\vec{P}}_f|\sum _{\lambda }|T{|}^2,$$

$${\Psi }_P(\vec{q})=(/P+M){\gamma }_{5}f(\vec{q}),$$

$${\phi }_{0^{-}}(q_{\perp })$$

$$f_3(q_{\perp })$$

$$f_4(q_{\perp })$$

$${\phi }_{0^{-}}^{++}(q_{\perp })=\left[A_1(q_{\perp })+\frac{/P}{M}{A}_2(q_{\perp })+\frac{/q_{\perp }}{M}{A}_3(q_{\perp })+\frac{/P/q_{\perp }}{M^2}{A}_4(q_{\perp })\right]{\gamma }^5,$$

$$A_1$$

$$A_2$$

$$\int \frac{\mathrm{d}\vec{q}}{(2\pi {)}^3}4f_1{f}_2{M}^2\left\{\frac{m_1+m_2}{{\omega }_1+{\omega }_2}+\frac{{\omega }_1+{\omega }_2}{m_1+m_2}+\frac{2{\vec{q}}^2(m_1{\omega }_1+m_2{\omega }_2)}{(m_2{\omega }_1+m_1{\omega }_2{)}^2}\right\}=2M.$$

$${\omega }_i\equiv \sqrt{m_i^2+\vec{q}{}^2}\approx m_i+\frac{\vec{q}{}^2}{2m_i}$$

$$\begin{array}{rl}& \int \frac{\mathrm{d}\vec{q}{}^{\prime }}{(2\pi {)}^3}{f}_1(|\vec{q}|)f_1^{\prime }(|\vec{q}{}^{\prime }|),\int \frac{\mathrm{d}\vec{q}{}^{\prime }}{(2\pi {)}^3}{f}_1(|\vec{q}|)f_2^{\prime }(|\vec{q}{}^{\prime }|),\\ & \int \frac{\mathrm{d}\vec{q}{}^{\prime }}{(2\pi {)}^3}{f}_2(|\vec{q}|)f_1^{\prime }(|\vec{q}{}^{\prime }|),\int \frac{\mathrm{d}\vec{q}{}^{\prime }}{(2\pi {)}^3}{f}_2(|\vec{q}|)f_2^{\prime }(|\vec{q}{}^{\prime }|),\end{array}$$

$${\xi }_{00}(v\cdot v^{\prime })=C\int \frac{\mathrm{d}\vec{q}{}^{\prime }}{(2\pi {)}^3}\overline{f{f}^{\prime }},$$

$$\begin{array}{rl}{\xi }_{11}& =C\int \frac{\mathrm{d}\vec{q}{}^{\prime }}{(2\pi {)}^3}\overline{f{f}^{\prime }}\frac{|\vec{q}{}^{\prime }|\cos \theta }{\sqrt{M{M}^{\prime }}},{\xi }_{20}=C\int \frac{\mathrm{d}\vec{q}{}^{\prime }}{(2\pi {)}^3}\overline{f{f}^{\prime }}\frac{\vec{q}{}^{\prime 2}}{M{M}^{\prime }},\\ {\xi }_{22}& =C\int \frac{\mathrm{d}\vec{q}{}^{\prime }}{(2\pi {)}^3}\overline{f{f}^{\prime }}\frac{\vec{q}{}^{\prime 2}{\cos }^2\theta }{M{M}^{\prime }},{\xi }_{31}=C\int \frac{\mathrm{d}\vec{q}{}^{\prime }}{(2\pi {)}^3}\overline{f{f}^{\prime }}\frac{|\vec{q}{}^{\prime }{|}^3\cos \theta }{\sqrt{(M{M}^{\prime }{)}^3}},\\ {\xi }_{33}& =C\int \frac{\mathrm{d}\vec{q}{}^{\prime }}{(2\pi {)}^3}\overline{f{f}^{\prime }}\frac{|\vec{q}{}^{\prime }{|}^3{\cos }^3\theta }{\sqrt{(M{M}^{\prime }{)}^3}},{\xi }_{40}=C\int \frac{\mathrm{d}\vec{q}{}^{\prime }}{(2\pi {)}^3}\overline{f{f}^{\prime }}\frac{\vec{q}{}^{\prime 4}}{(M{M}^{\prime }{)}^2},\\ {\xi }_{42}& =C\int \frac{\mathrm{d}\vec{q}{}^{\prime }}{(2\pi {)}^3}\overline{f{f}^{\prime }}\frac{\vec{q}{}^{\prime 4}{\cos }^2\theta }{(M{M}^{\prime }{)}^2},\end{array}$$

$$\int \frac{\mathrm{d}\vec{q}}{(2\pi {)}^3}4M{f}^2=1.$$

$$\begin{array}{rl}{S}_{+}& =-\frac{M+M^{\prime }}{2\sqrt{M{M}^{\prime }}}{\xi }_{00}+\frac{1}{4\sqrt{M{M}^{\prime }}}\left[b_1\frac{1}{m_1}+b_2\frac{1}{m_2}\right]\alpha {\xi }_{00}\\ & +\frac{1}{4P^{\prime }}\left[-b_1\frac{1}{m_1}-b_2\frac{1}{m_2}+a_1\frac{1}{m_1^{\prime }}+a_2\frac{1}{m_2^{\prime }}\right]{\xi }_{11}+\frac{(M-M^{\prime })P^{\prime 2}}{8\sqrt{M{M}^{\prime }}}\frac{1}{m_1{m}_2}{\alpha }^2{\xi }_{00}\\ & +\frac{P^{\prime }}{8}\left[(E^{\prime }+M^{\prime })\left(\frac{1}{m_1{m}_1^{\prime }}+\frac{1}{m_2{m}_2^{\prime }}\right)+(E^{\prime }-M^{\prime })\left(\frac{1}{m_1{m}_2^{\prime }}+\frac{1}{m_2{m}_1^{\prime }}\right)-(M-M^{\prime })\frac{1}{m_1{m}_2}\right]\alpha {\xi }_{11}\\ & +\frac{\sqrt{M{M}^{\prime }}}{8}\left[(M-M^{\prime })\left(\frac{1}{m_1{m}_2}-\frac{1}{m_1{m}_2^{\prime }}-\frac{1}{m_2{m}_1^{\prime }}+\frac{1}{m_1^{\prime }{m}_2^{\prime }}\right)-(M+M^{\prime })\left(\frac{1}{m_1{m}_1^{\prime }}+\frac{1}{m_2{m}_2^{\prime }}\right)\right]{\xi }_{20}\\ & +\frac{\sqrt{M{M}^{\prime }}}{8}(M-E^{\prime })\left[\frac{1}{m_1{m}_1^{\prime }}+\frac{1}{m_1{m}_2^{\prime }}+\frac{1}{m_2{m}_1^{\prime }}+\frac{1}{m_2{m}_2^{\prime }}\right]{\xi }_{22}\\ & -\frac{P^{\prime 2}}{16\sqrt{M{M}^{\prime }}}\left[b_1\frac{1}{m_1^3}+b_2\frac{1}{m_2^3}\right]{\alpha }^3{\xi }_{00}+\frac{P^{\prime }}{16}\left[b_1\frac{3}{m_1^3}+b_2\frac{3}{m_2^3}+a_1\frac{1}{m_1{m}_2{m}_2^{\prime }}+a_2\frac{1}{m_1{m}_2{m}_1^{\prime }}\right]{\alpha }^2{\xi }_{11}\\ & -\frac{\sqrt{M{M}^{\prime }}}{16}\left[b_1\left(\frac{1}{m_1^3}-\frac{1}{m_2{m}_1^{\prime }{m}_2^{\prime }}\right)+b_2\left(\frac{1}{m_2^3}-\frac{1}{m_1{m}_1^{\prime }{m}_2^{\prime }}\right)\right]\alpha {\xi }_{20}\\ & -\frac{\sqrt{M{M}^{\prime }}}{8}\left[b_1\frac{1}{m_1^3}+b_2\frac{1}{m_2^3}+a_1\frac{1}{m_1{m}_2{m}_2^{\prime }}+a_2\frac{1}{m_1{m}_2{m}_1^{\prime }}\right]\alpha {\xi }_{22}\\ & +\frac{M{M}^{\prime }}{16P^{\prime }}[b_1\left(\frac{1}{m_1^3}-\frac{1}{m_2{m}_1^{\prime }{m}_2^{\prime }}\right)+b_2\left(\frac{1}{m_2^3}-\frac{1}{m_1{m}_1^{\prime }{m}_2^{\prime }}\right)\\ & -a_1\left(\frac{1}{m_1^{{}^{\prime }3}}-\frac{1}{m_1{m}_2{m}_2^{\prime }}\right)-a_2\left(\frac{1}{m_2^{{}^{\prime }3}}-\frac{1}{m_1{m}_2{m}_1^{\prime }}\right)]{\xi }_{31}\end{array}$$

$$\begin{array}{rl}{S}_{-}& =\frac{M-M^{\prime }}{2\sqrt{M{M}^{\prime }}}{\xi }_{00}+\frac{1}{4\sqrt{M{M}^{\prime }}}\left[-c_1\frac{1}{m_1}+c_2\frac{1}{m_2}\right]\alpha {\xi }_{00}\\ & +\frac{1}{4P^{\prime }}\left[c_1\frac{1}{m_1}-c_2\frac{1}{m_2}+d_1\frac{1}{m_1^{\prime }}+d_2\frac{1}{m_2^{\prime }}\right]{\xi }_{11}-\frac{(M+M^{\prime })P^{\prime 2}}{8\sqrt{M{M}^{\prime }}}\frac{1}{m_1{m}_2}{\alpha }^2{\xi }_{00}\\ & +\frac{P^{\prime }}{8}\left[(M^{\prime }+E^{\prime })\left(\frac{1}{m_1{m}_1^{\prime }}+\frac{1}{m_2{m}_2^{\prime }}\right)+(E^{\prime }-M^{\prime })\left(\frac{1}{m_1{m}_2^{\prime }}+\frac{1}{m_2{m}_1^{\prime }}\right)+2(M+M^{\prime })\frac{1}{m_1{m}_2}\right]\alpha {\xi }_{11}\\ & +\frac{\sqrt{M{M}^{\prime }}}{8}\left[(M+M^{\prime })\left(-\frac{1}{m_1{m}_2}+\frac{1}{m_1{m}_2^{\prime }}+\frac{1}{m_2{m}_1^{\prime }}-\frac{1}{m_1^{\prime }{m}_2^{\prime }}\right)+(M-M^{\prime })\left(\frac{1}{m_1{m}_1^{\prime }}+\frac{1}{m_2{m}_2^{\prime }}\right)\right]{\xi }_{20}\\ & -\frac{\sqrt{M{M}^{\prime }}}{8}(M+E^{\prime })\left[\frac{1}{m_1{m}_1^{\prime }}+\frac{1}{m_1{m}_2^{\prime }}+\frac{1}{m_2{m}_1^{\prime }}+\frac{1}{m_2{m}_2^{\prime }}\right]{\xi }_{22}\\ & +\frac{P^{\prime 2}}{16\sqrt{M{M}^{\prime }}}\left[c_1\frac{1}{m_1^3}-c_2\frac{1}{m_2^3}\right]{\alpha }^3{\xi }_{00}+\frac{P^{\prime }}{16}\left[-c_1\frac{3}{m_1^3}+c_2\frac{3}{m_2^3}+d_1\frac{1}{m_1{m}_2{m}_2^{\prime }}+d_2\frac{1}{m_1{m}_2{m}_1^{\prime }}\right]{\alpha }^2{\xi }_{11}\\ & +\frac{\sqrt{M{M}^{\prime }}}{16}\left[c_1\left(\frac{1}{m_1^3}-\frac{1}{m_2{m}_1^{\prime }{m}_2^{\prime }}\right)-c_2\left(\frac{1}{m_2^3}-\frac{1}{m_1{m}_1^{\prime }{m}_2^{\prime }}\right)\right]\alpha {\xi }_{20}\\ & -\frac{\sqrt{M{M}^{\prime }}}{8}\left[c_1\frac{1}{m_1^3}-c_2\frac{1}{m_2^3}-d_1\frac{1}{m_1{m}_2{m}_2^{\prime }}-d_2\frac{1}{m_1{m}_2{m}_1^{\prime }}\right]\alpha {\xi }_{22}\\ & +\frac{M{M}^{\prime }}{16P^{\prime }}[-c_1\left(\frac{1}{m_1^3}-\frac{1}{m_2{m}_1^{\prime }{m}_2^{\prime }}\right)+c_2\left(\frac{1}{m_2^3}-\frac{1}{m_1{m}_1^{\prime }{m}_2^{\prime }}\right)\\ & -d_1\left(\frac{1}{m_1^{{}^{\prime }3}}-\frac{1}{m_1{m}_2{m}_2^{\prime }}\right)-d_2\left(\frac{1}{m_2^{{}^{\prime }3}}-\frac{1}{m_1{m}_2{m}_1^{\prime }}\right)]{\xi }_{31}\end{array}$$

$$\begin{array}{rl}⟨{\eta }_c|\bar{b}{\gamma }^{\mu }(1-{\gamma }^5)c|B_c^{+}⟩& =-\sqrt{M{M}_f}\left[v^{\mu }+v_f^{\mu }\right]{\xi }_{00},\\ S_{\pm }& =\mp \frac{M\pm M^{\prime }}{2\sqrt{M{M}^{\prime }}}{\xi }_{00}.\end{array}$$

$$C\int \frac{\mathrm{d}\vec{q}{}^{\prime }}{(2\pi {)}^3}\overline{f{f}^{\prime }}\frac{|\vec{q}|\cos \beta }{\sqrt{M{M}^{\prime }}},$$

$$⟨J/\psi |\bar{b}{\gamma }^{\mu }(1-{\gamma }^5)c|B_c^{+}⟩=\sqrt{M{M}_f}\left[ϵ\cdot v{v}_f^{\mu }-(v\cdot v_f+1){ϵ}^{\mu }+\mathrm{i}{\epsilon }^{\mu \nu \sigma \delta }{ϵ}_{\nu }{v}_{\sigma }{v}_{f\delta }\right]{\xi }_{00}.$$

$$0^{-}:\frac{M+/P}{2\sqrt{M}}{\gamma }^5\Psi ,1^{--}:\frac{M+/P}{2\sqrt{M}}/ϵ\Psi .$$

$$\begin{array}{rl}{0}^{++}& :\frac{/q_{\perp }}{|\vec{q}|}\frac{M+/P}{2\sqrt{M}}\Phi ,1^{++}:\mathrm{i}{\epsilon }_{\mu \nu \alpha \beta }\sqrt{\frac{3}{2}}\frac{P^{\nu }}{M}\frac{q_{\perp }^{\alpha }}{|\vec{q}|}{ϵ}^{\beta }\frac{M+/P}{2\sqrt{M}}{\gamma }^{\mu }\Phi ,\\ 2^{++}& :\sqrt{3}{ϵ}_{\mu \nu }{\gamma }^{\mu }\frac{q_{\perp }^{\nu }}{|\vec{q}|}\frac{M+/P}{2\sqrt{M}}\Phi ,1^{+-}:\sqrt{3}\frac{q_{\perp }\cdot ϵ}{|\vec{q}|}\frac{M+/P}{2\sqrt{M}}{\gamma }^5\Phi .\end{array}$$

Full text

I Introduction

Under the heavy quark effective theory (HQET), a semileptonic decay process can be related to a rotation of the heavy quark flavor or spin Neubert and Rieckert (1992); Neubert (1996). In the limit $m_{Q}\to\infty$ ($Q$ denotes the heavy quark or anti-quark), this rotation is a symmetry transformation. The form factors depend only on the Lorentz boost $\gamma=v\cdot v^{\prime}$ which connects the rest frames of the initial state and final state. The transition can be described by a dimensionless function $\xi(v\cdot v^{\prime})$. Heavy-quark symmetry reduces the weak-decay form factors of heavy hadrons to this universal function. These relations were derived by Isgur and Wise firstly Isgur et al. (1989); Isgur and Wise (1990), so called Isgur-Wise function (IWF).

HQET vastly simplifies the calculations, and plays a crucial role in extracting the values of $|V_{cb}|$ and $|V_{ub}|$. For example, the differential semileptonic decay rate for $B\to D$ in the heavy-quark limit can be model-independent described by Neubert (1996)

[TABLE]

The decay rate depends on only two quantities, $|V_{cb}|$ and $\xi(\omega)$. If the differential semileptonic decay rate is measured by experiments, one can obtain the value of $|V_{cb}|\xi(\omega)$. Exploiting the normalization of Isgur-Wise function $\xi(1)=1$, the value of $|V_{cb}|$ can be extracted. Conversely exploiting the given value of $|V_{cb}|$, the differential semileptonic decay rate can be obtained by calculating the Isgur-Wise function. But Isgur-Wise function can not be calculated by perturbation theory in principle, and can only be abtained by various phenomenological models. Obviously the latter is model-dependent. Then a great deal of efforts were directed to study the Isgur-Wise function and its applications in different frameworks. Kiselev et al. exploited the Isgur-Wise function in a covariant potential model to calculate the semileptonic $B\to D^{(*)}\ell\nu$ Decays Kiselev (1995). Olsson et al. explored the Isgur-Wise Function in the relativistic flux tube model Olsson and Veseli (1995). Huang et al. calculated the leading Isgur-Wise function parametrizing the semileptonic decays $\Lambda_{b}\to\Lambda_{c1}^{1/2,3/2}$ by using the QCD sum rules Huang et al. (2001). Lee et al. calculated the subleading Isgur-Wise function for the same processes by using the QCD sum rules Lee and Park (2003). Hassanabadi and Rahmani et al. investigated the Isgur-Wise Function for semileptonic decays of bottom and charmed hadrons in a series of models, including a Cornell potential combined a improved variational approach Rahmani and Hassanabadi (2014), the baryon hyperspherical approach Hassanabadi et al. (2014a), a combination of linear confining and Hulthén potentials Hassanabadi et al. (2014b), the three-boy Schrödinger equationHassanabadi and Rahmani (2015), a combination of Deng-Fan-type and harmonic potentials Hassanabadi et al. (2016), a quantum isotonic nonlinear oscillator potential model Rahmani and Hassanabadi (2017) and so on. Choudhury et al. had studied the renormalization scale dependence of IW function by using a wave function with linear part of the Cornell potential as perturbation Choudhury and Pathak (2014). Yazarloo et al. employed the nonrelativistic Schrödinger equation with a new potential model to calculate the semileptonic decay width and branching ratio of the $B\to D^{*}\ell\nu$ process Yazarloo and Mehraban (2016). In a potential model, the Isgur-Wise function usually corresponds to the overlapping integral of the wave functions obtained by this model.

However the lowest order result is not accurate enough due to the heavy quark approximation. The symmetry-breaking corrections are needed when the study becomes more precise, since the masses of the heavy quarks or anti-quarks are not infinite actually. In addition the radiation correction can not be ignored. The HQET provides a systematic framework to analyze these corrections. For example, Luke analyzed the $1/m_{Q}$ corrections for the more complicated case of weak decay form factors Luke (1990). Falk et al. analyzed the structure of $1/m^{2}_{Q}$ corrections for both meson and baryon weak decay form factors, and calculated the leading QCD radiative corrections Falk et al. (1992, 1990). Neubert et al. developed a method to express the universal functions, that appear in leading and next-to-leading order in the $1/m_{Q}$ expansion of current matrix elements, in terms of hadronic form factors evaluated at maximum recoil Neubert and Rieckert (1992). Leibovich et al. investigated the $1/m_{Q}$ corrections for semileptonic $\Lambda_{b}$ decays to excited $\Lambda_{c}$ at zero recoil Leibovich and Stewart (1998). Cardarelli presented an investigation of the Isgur-Wise form factor relevant for the semileptonic decay $\Lambda_{b}\to\Lambda_{c}\ell^{-}\nu_{\ell}$ performed within a light-front constituent quark model, including both radiative effects and first-order power corrections in the inverse heavy-quark mass Cardarelli and Simula (2000). Coleman et al. investigated a new effective IW function with a shift in slope, and all of $1/m_{Q}$ corrections and QCD radiative corrections are distilled Coleman et al. (2001). Other efforts of complements are too many to be list here.

The flavor-spin symmetry can be used not only in weak decays of ground hardons containing one heavy quark, but also in the processes involving orbitally and radially excited states containing one heavy quark. Isgur and Wise exploited the flavor-spin symmetry to obtain model-independent predictions in weak decays from pseudoscalar meson of a heavy quark $Q_{i}$ to P wave excited states of another heavy quark $Q_{j}$ in term of two Isgur-Wise functions Isgur and Wise (1991). Gan et al. employed the QCD sum rule approach to estimate the Isgur-Wise functions for the semileptonic decays of the ground statemeson $B_{s}$ into orbitally excited D-wave $\bar{c}s$ mesons Gan et al. (2015). Yaouanc calculated the Isgur-Wise functions that involves only heavy mesons with light cloud $j^{P}=1/2^{-}$ and their radial excitations by Bakamjian-Thomas relativistic quark model Le Yaouanc et al. (2016). Bec̆irević et al. have studied the $1/m_{Q}$ corrections of form factors for semi-leptonic transitions between $B$ mesons and radially excited $D$ mesons in the Bakamjian-Thomas relativistic quark model Bečirević et al. (2017a, b).

However, the validity of HQET is suspectable, when the systems contain two or more heavy degrees of freedom. Some work have explored the application of the heavy quark symmetry to describe the weak decays of hadrons containing two heavy quarks. Sanchis-Lozano estimated that the flavor symmetry losts, while the spin symmetry holds for the double-heavy meson Sanchis-Lozano (1995). As for the baryon containing two heavy quarks, two heavy quarks constitute a bosonic diquark whose mass could be regarded as infinite, and therefore HQET still holds true for the diquark-light quark system. But HQET is not suitable for dealing with the diquark subsystem. Ebert et al. exploited the relativistic quark model to obtain the diquark wave functions, and then the transition amplitudes of heavy diquarks $bb$ and $bc$ going respectively to $bc$ and $cc$ are expressed through the overlap integrals of corresponding diquark wave functions Ebert et al. (2004). Pathak et al. treated $B_{c}$ meson as a typical heavy-light meson like $B$ or $D$ within a QCD potential model, and the semileptonic decay rates of $B_{c}$ meson into $\eta_{c}$, $J/\psi$ are exploited Pathak and Choudhury (2013). Das et al. computed the slopes and curvatures of the form factors of semileptonic decays of heavy-light mesons included $B_{c}$ Das et al. (2016). The error of the result for $B_{c}$ is very large. Wang at al. obtained form factors for the $B_{c}$ into S-wave and P-wave charmonium by three universal Isgur-Wise functions Wang and Zhu (2018). Almost all of these studies just use the formulas in HQET directly, but the usefulness of the heavy quark symmetry to describe double heavy hadrons has been little reported.

The aim of this paper is to investigate the heavy quark symmetry in the double-heavy mesons without the use of heavy quark limit. The heavy quark symmetry, if exists, must be reflected in the results obtained by a suitable phenomenological model without the use of heavy quark limit. The instantaneous Bethe-Salpeter (BS) equation is just a very effective method to deal with double-heavy mesons. This method has a comparatively solid foundation because both the BS equation and the Mandelstam formula are established on relativistic quantum field theory. Meanwhile the instantaneous approximation is reasonable, since the quark and antiquark in double-heavy mesons are both heavy. This method can give an analytic expression, so the symmetry can be found intuitively. However, the accuracy may not be as good as Lattice QCD. We choose the semileptonic $B_{c}$ decays to charmonium , and the final mesons involve the orbitally and radially excited states. We conclude that the flavor symmetry breaks, while the spin symmetry holds for the double-heavy meson, as Sanchis-Lozano estimated Sanchis-Lozano (1995). Further, we investigate the behaviors of the Isugr-Wise function and other high-order relativistic corrections to Isgur-Wise function.

Actually, some works have been done on IWF with BS equation. Kugo et al. expanded BS equation in orders of the inverse heavy quark mass and defined the leading term in the expansion of the first form factor as IWF Kugo et al. (1994). El-Hady et al. pointed out that the IWF can be related to the overlap integral of normalized meson wave functions in the infinite momentum frame and it should be possible to calculate the form factors directly without using the heavy quark limit Abd El-Hady et al. (1995). Zoller et al. calculated the numerical IWF by multiplying quark masses with a large factor directly Zoller et al. (1995). Chang et al. obtained two universal functions in $B_{c}\to h_{c},\chi_{c}$ with the instantaneous BS method, but the wave functions they used are nonrelativistic Chang and Chen (1994); Chang et al. (2002). Nowadays the instantaneous BS method has developed to be quite covariant, and the full Salpeter equations are solved for different $J^{P(C)}$ states Cvetic et al. (2004); Chang et al. (2005); Wang (2007, 2009). So the relativistic correction which equate to symmetry-breaking correction has been taken into account. In this paper we shall not derive these corrections from HQET, but attempt to extract the IWF from the solutions of the instantaneous BS method. Note that we do not use the heavy quark limit. The results show that only the lowest order correction to IWF is not enough, and higher-order relativistic corrections need to be introduced for more accurate results.

The paper is organized as follows. In section II, we give the useful formulas for the $B_{c}$ decays to charmonium. In section III, we give the relativistic wave function for $0^{-}$ state in the instantaneous BS method. In section IV, we extract the IWF and give the analytical results. In section V, we give the numerical results and discussions. We summarize and conclude in section VI, and put the Salpeter equation and some wave functions in the appendix A.

II Form factors and semileptonic decay width

For the $B_{c}^{+}\to(c\bar{c})\ell^{+}\nu_{\ell}$ processes shown in figure 1, the transition amplitude element reads

[TABLE]

where $(c\bar{c})$ denotes charmonium; $V_{cb}$ is the Cabibbo-Kobayashi-Maskawa (CKM) matrix element; $J_{\mu}\equiv V_{\mu}-A_{\mu}$ is the charged current responsible for the decays; $P$ and $P_{f}$ are the momenta of the initial $B_{c}^{+}$ and the final charmonium, respectively.

The hadronic transition element can be written as the overlapping integral over the initial and final relativistic BS wave functions within Mandelstam formalism. We would not solve the full BS equation, but the instantaneous one, namely, the full Salpeter equation. We perform the instantaneous approximation to the transition element Chang et al. (2006) and write it as

[TABLE]

where $\varphi_{P}^{++}$ denotes the positive energy component of the instantaneous BS wave function of the initial state; $\overline{\varphi}_{P_{f}}^{++}\equiv\gamma^{0}(\varphi_{P_{f}}^{++})^{\dagger}\gamma^{0}$ is the Dirac conjugate of the positive energy component of the final state; $m^{\prime}_{1}$ and $m^{\prime}_{2}$ are the masses of quark and antiquark in the final state, respectively, and $\vec{q}\>^{\prime}=\vec{q}-\frac{m^{\prime}_{1}}{m^{\prime}_{1}+m^{\prime}_{2}}\vec{P}_{f}$ is the relative momentum between them. In this paper, we keep only the positive energy component $\varphi^{++}$ of the relativistic wave functions, because the contributions from other components are much smaller than $1\%$ in transition of $B_{c}\to(c\bar{c})$ Wang et al. (2017). This matrix element can also be written in the framework in which the momentum $\vec{q}\>^{\prime}$ is the integral argument by means of a suitable Jacobi transformation,

[TABLE]

where $\vec{q}=\vec{q}\>^{\prime}+\alpha\vec{P}_{f},\alpha=\frac{m^{\prime}_{1}}{m^{\prime}_{1}+m^{\prime}_{2}}$ and $\theta$ is the angle between $\vec{q}\>^{\prime}$ and $\vec{P}_{f}$. The Eq. (4) is more convenient, because some matrix elements we calculate in this paper involve a $P$-wave final state Chang et al. (2002).

For $B_{c}^{+}\to P\ell^{+}\nu_{\ell}$ (here $P$ denotes $\eta_{c}$ or $\chi_{c0}$), the hadronic matrix element can be written as

[TABLE]

where $S_{+}$ and $S_{-}$ are the form factors. For $B_{c}^{+}\to V\ell^{+}\nu_{\ell}$ (here $V$ denotes $J/\psi$, $h_{c}$ or $\chi_{c1}$), the hadronic matrix element can be written as

[TABLE]

where $\epsilon^{\mu}$ is the polarization vector of the final vector meson; $t_{1}$, $t_{2}$, $t_{3}$ and $t_{4}$ are the form factors. For $B_{c}^{+}\to T\ell^{+}\nu_{\ell}$ (here $T$ denotes $\chi_{c2}$), the hadronic matrix element can be written as

[TABLE]

where $\epsilon_{\alpha\beta}$ is the polarization tensor of the final tensor meson; $t_{1}$, $t_{2}$, $t_{3}$ and $t_{4}$ are the form factors.

The summation formulas for polarization of the final vector meson are

[TABLE]

The summation formulas for polarization of the final tensor meson are

[TABLE]

Where $S_{\mu\nu}=-g_{\mu\nu}+\frac{P_{f\mu}P_{f\nu}}{M_{f}^{2}}$. Finally, the semileptonic decay width can be expressed as

[TABLE]

where $\vec{P}_{\ell}$ is the three-dimensional momentum of the final lepton, and $\vec{P}_{f}$ is the three-dimensional momentum of the final meson. In this paper, we only calculate the form factors but no longer calculate the decay widths.

III Relativistic wave function

Usually, the nonrelativistic wave function for a pseudoscalar is written as Chang and Chen (1994)

[TABLE]

where $M$ and $P$ are the mass and momentum of the meson, respectively; $\vec{q}$ is the relative momentum between the quark and antiquark in the meson, and the radial wave function $f(\vec{q})$ can be obtained numerically by solving the Schrodinger equation.

But in our method, we solve the full Salpeter equation. The form of wave function is relativistic and depends on the $J^{P(C)}$ quantum number of the corresponding meson. For a pseudoscalar, the relativistic wave function can be written as the four items constructed by $P$, $q_{\perp}$ and $\gamma$-matrices Kim and Wang (2004)

[TABLE]

where $q=p_{1}-\alpha_{1}P=\alpha_{2}P-p_{2}$ is the relative momentum between quark (with momentum $p_{1}$ and mass $m_{1}$) and antiquark (momentum $p_{2}$ and mass $m_{2}$), $\alpha_{1}=\frac{m_{1}}{m_{1}+m_{2}}$, $\alpha_{2}=\frac{m_{2}}{m_{1}+m_{2}}$; $q_{\perp}=q-\frac{P\cdot q}{M^{2}}P$, in the rest frame of the meson, $q_{\perp}=(0,\vec{q})$.

All the items in the wave function Eq. (12) have the quantum number of $0^{-}$. This wave function is a general relativistic form for a pseudoscalar with the instantaneous approximation. If we set the items with $f_{3}$ and $f_{4}$ to zero, and set $f_{1}=f_{2}$, the relativistic wave function is reduced to the Schrodinger wave function Eq. (11).

Taking into account the last two equations in Eq. (39), we obtain the relations

[TABLE]

where the quark energy $\omega_{i}=\sqrt{m^{2}_{i}-q_{\perp}^{2}}=\sqrt{m^{2}_{i}+\vec{q}^{\>2}}$ ($i=1,2$). The wave function corresponding to the positive energy projection has the form

[TABLE]

where

[TABLE]

The normalization condition reads

[TABLE]

By solving the full Salpeter equation, the numerical values of wave functions $f_{1}$ and $f_{2}$ are obtained. The positive energy component Eq. (14) is brought into the Mandelstam formula Eq. (3). After the trace and integral are finished, the form factors $S_{+}$ and $S_{-}$ can be calculated numerically. In this paper, besides the wave function for $0^{-}$ state, we also need the wave functions for the states of $1^{--}$ ($J/\psi$), $1^{+-}$ ($h_{c}$), $0^{++}$ ($\chi_{c0}$), $1^{++}$ ($\chi_{c1}$) and $2^{++}$ ($\chi_{c2}$). We put the $2^{++}$ state wave function in the appendix A and the others can be referred to Geng et al. (2019).

IV Isgur-Wise function

Since $B_{c}$ and charmonium are the weak-binding states, the approximation

[TABLE]

is taken in this paper. This approximation requires the three-dimensional relative momentum $|\vec{q}|$ between quarks much less than the masses of quarks, but its contribution will be suppressed by the wave function $f_{i}(\vec{q})$ in the large $|\vec{q}|$ interval. After performing this approximation and the trace on the matrix element Eq. (4), the dependence of all the form factors on the overlapping integrals of the wave functions for the initial state and the final state becomes transparent. For instance, one type of overlapping integrals are

[TABLE]

where $f_{i}$ denotes the wave function of initial state, and $f^{\prime}_{i}$ denotes the wave function of final state. Two wave functions from the same meson is very close numerically, i.e., $f_{1}\approx f_{2}$ and $f^{\prime}_{1}\approx f^{\prime}_{2}$. So the four overlapping integrals in Eq. (18) are approximately equal, and for convenience they are replaced by their average which is denoted as

[TABLE]

where $C$ is the normalized coefficient; $v,v^{\prime}$ are the four dimensional velocities of the initial state and final state respectively, and $\overline{ff^{\prime}}=\frac{f_{1}f^{\prime}_{1}+f_{1}f^{\prime}_{2}+f_{2}f^{\prime}_{1}+f_{2}f^{\prime}_{2}}{4}$. There are other overlapping integrals with the relative momentum $\vec{q}\>^{\prime}$ being inserted. They may be the relativistic corrections to the function $\xi_{00}$. We denote them as $\xi_{qx}$, where subscript $q$ denotes the power of the relative momentum $\vec{q}\>^{\prime}$, subscript $x$ denotes the power of $\cos\theta$, and $\theta$ is the angle between $\vec{q}\>^{\prime}$ and $\vec{P}_{f}$, i.e.,

[TABLE]

and so on. When the final state is S-wave meson, we keep the first six functions and abandon the higher order $\mathcal{O}(q^{4})$; When the final state is P-wave meson whose wave function includes a $\vec{q}$, $\xi_{00}$ disappears, thus we reserve the first eight functions and abandon the higher order $\mathcal{O}(q^{5})$. The normalized coefficients based on the normalized formulas are shown in Table. 1 for each process. Taking the process $B_{c}\to\eta_{c}$ as an example, the initial and final states are both $0^{-}$ state. With the approximations $\vec{q}=0$, $f_{1}=f_{2}$ and $\omega_{i}=m_{i}$, Eq .(16) is deduced as

[TABLE]

So the normalized wave function of $0^{-}$ state is $2\sqrt{M}f$, and the normalized coefficient is $4\sqrt{MM^{\prime}}$ for the process $B_{c}\to\eta_{c}$.

The form factors of semileptonic decay $B_{c}\to\eta_{c}\ell\nu_{\ell}$ can be written as

[TABLE]

where $a_{1}=E^{\prime 2}-E^{\prime}M+E^{\prime}M^{\prime}-MM^{\prime}=M^{\prime}(E^{\prime}-M)(\omega+1)$, $a_{2}=E^{\prime 2}-E^{\prime}M-E^{\prime}M^{\prime}+MM^{\prime}=M^{\prime}(E^{\prime}-M)(\omega-1)$, $b_{1}=MM^{\prime}-E^{\prime}M^{\prime}+E^{\prime}M-M^{\prime 2}=M^{\prime}(M-M^{\prime})(1+\omega)$, $b_{2}=MM^{\prime}-E^{\prime}M^{\prime}-E^{\prime}M+M^{\prime 2}=M^{\prime}(M+M^{\prime})(1-\omega)$, $b_{1}=\vec{P}_{f}^{2}-a_{1}$, $b_{2}=a_{2}-\vec{P}_{f}^{2}$.

[TABLE]

where $c_{1}=E^{\prime}M+E^{\prime}M^{\prime}+MM^{\prime}+M^{\prime 2}$, $c_{2}=E^{\prime}M-E^{\prime}M^{\prime}-MM^{\prime}+M^{\prime 2}$, $d_{1}=E^{\prime 2}+E^{\prime}M+E^{\prime}M^{\prime}+MM^{\prime}$, $d_{2}=E^{\prime 2}+E^{\prime}M-E^{\prime}M^{\prime}-MM^{\prime}$, $c_{1}=d_{1}-\vec{P}_{f}^{2}$, $c_{2}=d_{2}-\vec{P}_{f}^{2}$.

The function $\xi_{00}$ may be directly related to the Isgur-Wise function appearing in HQET for $0^{-}\to 0^{-}$ decays. Because the form factors in this process will degenerate into those in the nonrelativistic limit if only the function $\xi_{00}$ is considered Neubert (1996),

[TABLE]

The other functions are the relativistic corrections ($1/m_{i}$ corrections) to the leading order IWF $\xi_{00}$, where $i$ denote a quark or anti-quark in the initial and final mesons. The number of $\vec{q}\>^{\prime}$ contained in the function $\xi_{qx}$ (subscript $q$) corresponds to the order of the correction. Note that there should have been another type of overlapping integrals with the relative momentum $\vec{q}$ in the initial state. For example,

[TABLE]

where $\beta$ is the angle between $\vec{q}$ and $\vec{P}_{f}$. Due to the relation $\vec{q}=\vec{q}\>^{\prime}+\alpha\vec{P}_{f},\alpha=\frac{m^{\prime}_{1}}{m^{\prime}_{1}+m^{\prime}_{2}}$, this overlapping integral Eq. (25) is decomposed into $\xi_{11}+\alpha|\vec{P}_{f}\>|\xi_{00}$. So the item involving $\alpha\xi_{00}$ should be considered as relativistic correction of the same order as $\xi_{11}$. Generally, the item involving $\alpha^{n}\xi_{qx}$ is the $q+n$ order relativistic correction ($1/m_{i}^{q+n}$ correction) which can be confirmed in Eq. (22) and (23). The process $0^{-}\to 1^{--}$ is the same as above case. The leading order result is agree with HQET Neubert (1996), i.e.,

[TABLE]

It is very natural that the leading order results in this paper are entirely consistent with HQET for $0^{-}\to 0^{-}$ or $1^{--}$ processes. Because for the leading order results, the terms involving $\not{q}$ disappear and $\omega_{i}=m_{i}$, so that the BS wave functions degenerate into the nonrelativistic case, i.e.,

[TABLE]

A pseudoscalar meson and its corresponding vector have the same radial wave function $\Psi$ in the nonrelativistic limit. But in this paper, the radial wave functions are obtained by solving BS equation, and the $\Psi$ in Eq. (27) corresponds to the normalized wave function $2\sqrt{M}f_{i}$. They are not exactly the same numerically. And then, the numerical results of IWF is not close to HQET and contains part of the relativistic correction.

For P-wave meson as the final state, the nonrelativistic wave functions are usually written as

[TABLE]

and these states have the same radial wave function $\Phi$. Similarly in this paper, the radial wave functions are obtained by solving BS equation, and the $\Phi$ in Eq. (28) corresponds to the normalized wave function. In these cases, $\xi_{00}$ disappears and $\xi_{11}$ is IWF. Due to the change of orbital angular momentum, this weak decay process can not correspond to a scattering process simply. We give the leading order results in the case that only the function $\xi_{11}$ is considered,

[TABLE]

These results are not agree with Ref. Wang and Zhu (2018). The latter analyzes the reduction of form factors in the heavy quark limit, and there are two IWFs $\xi_{E},\xi_{F}v_{\alpha}$ for $B_{c}$ to P-wave charmonium. While we only need IWF $\xi_{11}$ for these processes in the leading order. Ref. Wang and Zhu (2018) does not further describe the used IWFs. The difference needs further examination. Note that the above results are not confined to the processes of $B_{c}$ to charmonium, but hold true for any possible process where the initial and final mesons are corresponding $J^{PC}$ states. In the next section, we will give the numerical results and discussions on the specific processes.

V Results and Discussions

The parameters used in this paper: $\Gamma_{B_{c}}=1.298\times 10^{-12}~{}\mathrm{GeV},G_{F}=1.166\times 10^{-5}~{}\mathrm{GeV^{-2}},m_{b}=4.96~{}\mathrm{GeV},m_{c}=1.62~{}\mathrm{GeV}$, $M_{h_{c}(2P)}=3.887~{}\mathrm{GeV},M_{\chi_{c0}(2P)}=3.862~{}\mathrm{GeV},M_{\chi_{c1}(2P)}=3.872~{}\mathrm{GeV},M_{\chi_{c2}(2P)}=3.927~{}\mathrm{GeV}$.

After solving the corresponding full Salpeter equations, the numerical wave functions for different mesons are obtained and shown in figures 2-3. When $|\vec{q}\>|$ is large, the wave functions will decrease rapidly. So the weak-binding approximation Eq. (17) can be taken here, and the error from large $|\vec{q}\>|$ will be suppressed by the wave function. The numerical values of two dominate wave functions are almost equivalent for each meson, so the approximation that the four overlapping integrals in Eq. (18) are replaced by their average is reasonable. For $1^{--}$ or $2^{++}$ state, there are two other minor wave functions $g_{3},g_{4}$, and $g_{3}\approx-g_{4}$. Taking the approximation $g_{3}=-g_{4}$ and weak-binding approximation Eq. (17), these two minor wave functions $g_{3},g_{4}$ only appear in the $\mathcal{O}(q^{4})$ or higher order in the $1^{--}$ state BS wave function. Within the precision $\mathcal{O}(q^{3})$ of this study for process $0^{-}\to 1^{--}$, these two minor wave functions $g_{3},g_{4}$ disappear. It is same for $2^{++}$ state wave function. These are consistent with Eq. (27) and (28), in which there is only one radial wave function.

The behaviors of the Isgur-Wise function and high-order corrections, i.e., the overlapping integrals of the wave functions of the initial and final bound states, are computed numerically and plotted in Fig. 4-5, where $\omega=v\cdot v_{f}=\frac{P\cdot P_{f}}{MM_{f}}$. These IWFs can be classified into four categories according to the configurations $nL$ of initial and final states. They belong to the modes $1S\to 1S$, $1S\to 1P$, $1S\to 2S$ and $1S\to 2P$ respectively. With the same configurations of initial and final states, for example, in the processes $B_{c}\to\eta_{c}$ and $B_{c}\to J/\psi$, the behaviors of IWFs are virtually identical except the ranges of $\omega=v\cdot v_{f}$ slightly differ. Because these decay processes are just related by a rotation of the heavy-quark spin or the meson spin, and this rotation is a symmetry transformation in the infinite-mass limit. Note that the infinite-mass limit is not used in this paper, but this spin-symmetry reflected in the results automatically, as Fig. 4-5 shows. This indicates that spin-symmetry still mantains though the initial and final states are both the double-heavy mesons. When the configuration of initial or final state changes, for example, the final $\eta_{c}$ turns into $\eta_{c}(2S)$, the behaviors of IWFs become significantly different from before. Next we will discuss these four modes one by one.

The mode $1S\to 1S$ has been extensively studied in HQET. $B_{c}$ and $\eta_{c}$ are related by the replacement $b\to c$, while $\eta_{c}$ and $J/\psi$ are related by the transformation $c^{\Uparrow}\to c^{\Downarrow}$ here. These two rotations (flavour and spin rotations) are symmetry transformations in the infinite-mass limit. So the radial wave functions of these mesons will be identical in this limit, and the corresponding IWF whose general form is the overlapping integral of the initial and final wave functions will be the same as the normalization formula at zero recoil. It is very natural that $\xi(1)=1$ in HQET. In this paper we solve the full Salpeter equations without the infinite-mass limit, and the normalized radial wave function is approximated to $2\sqrt{M}f$ for $1S$ state. The normalized wave functions have little difference for $\eta_{c}$ and $J/\psi$ (two dominate wave functions), which is consistent with Eq. (27). But the discrepancy between $B_{c}$ and the former two is in the order of 30% (peak value), as Fig. 2(a)-2(c) shows. This indicates that in the double-heavy system the spin-symmetry keeps, while the flavour-symmetry breaks. The masses of quark and antiquark are in the same order of magnitude, and therefore the change of flavour will lead to a great impact. Although the behaviors of IWF $\xi_{00}$ are the same as $\xi$ in HQET, they are not strict unity at zero recoil in this paper, as Fig. 4(a)-4(b) shows. The relativistic correction reflected in IWF $\xi_{00}$ is around 10% at zero recoil. In the mode $1S\to 1S$, it is convenient to fit the IWF as

[TABLE]

where $\rho^{2}$ is the slope parameter and $c$ is the curvature parameter which characterizes the shape of the IWF. The slope and curvature by fitting are 2.25 and 1.74 respectively in $B_{c}\to\eta_{c}$, and they are 2.38 and 1.98 respectively in $B_{c}\to J/\psi$. The result is agree with the rule that the slope is bigger as the (reduced) mass is heavier Das et al. (2016); Hassanabadi et al. (2014b). The other functions $\xi_{qx}$ are the relativistic corrections to IWF $\xi_{00}$. The more $\vec{q}\>^{\prime}$ the correction function contains, the less contribution it makes. We may call the correction function with one relative momentum $\vec{q}\>^{\prime}$ as the first order correction, the correction function with two $\vec{q}\>^{\prime}$ as the second order correction, and so on. The values of $\xi_{11}$ is about 1/20 of $\xi_{00}$. Because the decay width is proportional to modular square of amplitude, the first order correction may reach 1/10 of the leading order result. It is important for accurate calculation. Our previous study shows that the higher order relativistic corrections also have considerable contributions, and the total relativistic correction can reach around 20% at the level of decay width Geng et al. (2019).

In the mode $1S\to 1P$, the configuration of initial state is $1S$, while the configuration of final state is $1P$. Their orbital angular momenta are different, so the symmetry transformations exist only between the final states, i.e., spin rotations. $\chi_{c0}$, $\chi_{c1}$ and $\chi_{c2}$ are spin triplet states that are related by the rotation of total spin component (the component of total spin in the direction of orbital angular momentum), while $h_{c}$ and the former three are related by the transformation $c^{\Uparrow}\to c^{\Downarrow}$. These two spin rotations are symmetry transformations in the infinite-mass limit, so the normalized radial wave functions of these mesons will be identical. In this paper we study these mesons without the infinite-mass limit, and their normalized radial wave functions are approximated to $\frac{2|\vec{q}\>|h}{\sqrt{3M}},\frac{2|\vec{q}\>|\phi}{\sqrt{M}},\frac{2\sqrt{2}|\vec{q}\>|\psi}{\sqrt{3M}}$ and $\frac{2\sqrt{M}|\vec{q}\>|\zeta}{\sqrt{3}}$ respectively. Their numerical results are almost the same, as Figs. 2(d)-2(f) shows, which is consistent with Eq. (28). This indicates that the spin-symmetry keeps in the P-wave charmonium though the quark and anti-quark have the same masses. Because P-wave function contains a $\vec{q}\>$, $\xi_{00}$ disappears, and IWF is $\xi_{11}$ which behavior is obvious different from $\xi_{00}$. Due to the presence of $\cos\theta$, see Eq. (20), the IWF $\xi_{11}$ is zero at zero recoil. And it is enhanced kinematically, as Figs. 4(c)-4(e) shows. This behavior is agree with Ref. Chang et al. (2002). There is a kinematically suppressed factor $1/|\vec{v}_{f}|$ in the form factors Eq. (29), and therefore the behaviors of the leading order form factors are not purely dependent on $\xi_{11}$. The $\xi_{20}$ and $\xi_{22}$ are comparable to the leading order $\xi_{11}$, especially at zero recoil. $\xi_{22}$ is smaller than $\xi_{20}$ due to the factor $\cos^{2}\theta$. They decrease slowly when the momentum recoil increases and therefore the relativistic corrections may be comparable to the nonrelativistic results in this mode. Although the other correction functions seem to be very small, they are still important for accurate calculation, just as the mode $1S\to 1S$. For the final states as $h_{c}$, $\chi_{c0}$, $\chi_{c1}$ and $\chi_{c2}$, the total relativistic corrections are 50%, 64%, 34% and 14% at the level of decay width respectivelyGeng et al. (2019). The total relativistic correction of $B_{c}\to\chi_{c2}$ is unusually small, because the different orders corrections cancel each other out. This can be seen in the following analysis of form factors.

In the mode $1S\to 2S$, the configuration of initial state is different from the final state. Similarly, the only symmetry transformation is the spin rotation $c^{\Uparrow}\to c^{\Downarrow}$ which relates $\eta_{c}(2S)$ with $\psi(2S)$. Their normalized wave functions are almost the same, as Fig. 3(a)-3(b) shows, which is consistent with Eq. (27). The numerical results of IWFs in this mode is negative. Though the negative and positive of IWFs dose not affect the result of width, but it indicates the negative parts of $2S$-wave functions play a primary role. The overlapping integral of wave functions $\int\mathrm{d}\vec{q}=\int\vec{q}\>^{2}\sin\theta\mathrm{d}|\vec{q}\>|\mathrm{d}\theta\mathrm{d}\phi$ contains a factor $\vec{q}\>^{2}$. It is suppressed when $|\vec{q}\>|<1$ while is enhanced when $|\vec{q}\>|>1$. The negative parts of $2S$-wave functions are mainly in the range of $|\vec{q}\>|>1$, so the negative parts play a primary role in the overlapping integral. The IWF $\xi_{00}$ is increasing together with the momentum recoil, as Figs. 5(a)-5(b) shows, and the leading order form factors have the same behaviors due to Eq. (24) and (26). The behaviors of the other correction functions in this mode are similar to $1S\to 1S$, but they make more contributions here. For example, $\xi_{11}$ and $\xi_{20}$ are about one fifth and one eighth of $\xi_{00}$ at the maximum recoil, respectively. So the relativistic corrections become greater, and our previous study shows they are about 19%–28% larger than those in the mode $1S\to 1S$ Geng et al. (2019).

Comparing to the mode $1S\to 1P$, the analysis about the symmetry and normalized wave functions is the same in $1S\to 2P$, as Figs. 5(c)-5(e) shows. The IWF $\xi_{11}$ is also zero at zero recoil, but increases more slowly. The $\xi_{20}$ and $\xi_{22}$ are comparable to IWF $\xi_{11}$, and are no more decreasing but increasing as the momentum recoil is increasing. So the relativistic corrections may be more significant in this mode. They are about 10%–16% larger than those in the mode $1S\to 1P$ Geng et al. (2019).

The above analysis of relativistic corrections is qualitative, because the kinematic factors multiplied by IWFs are different and complex. In order to discuss these relativistic corrections precisely, the form factors in these processes are calculated by different order corrections in turn. Their numerical results are compared with those calculated by instantaneous Bethe-Salpeter method directly, as Figs. 6-13 shows. In these Figs, $t\equiv(P-P_{f})^{2}$ is the momentum transfer, and $t_{m}$ is the maximum of $t$, so $t_{m}-t=2MM_{f}(v\cdot v^{\prime}-1)$. The circle-solid line (BSE) denotes the form factor calculated by instantaneous Bethe-Salpeter method directly, and we regard it as the relatively precise result because this method is almost covariant; the solid line denotes the leading order (LO) of form factor calculated only by IWF; the dash line denotes the result with IWF and first order (1st) correction; the dot-dash line denotes the result with IWF, the first and second order (2nd) corrections; the dot line denotes the result with IWF, the first, second and third order (3rd) corrections.

For the process $B_{c}\to\eta_{c}$, as Figs. 6(a)–6(b) shows, there is some gap between the Leading order $S_{+}$ and the result from BSE. The $S_{+}$ with 1st correction is close to BSE. When the 3rd correction is taken into account, the result becomes very accurate. The difference between the LO $S_{-}$ and BSE is slightly larger, but due to the small contribution of $S_{-}$ to the decay width, the nonrelativistic results may be approximate. Though the high order corrections do not make $S_{-}$ and BSE exactly the same, they become closer. For the process $B_{c}\to J/\psi$, as Figs. 6(c)–6(f) shows, the form factor $t_{3}$ makes main contribution to the decay width. The LO $t_{3}$ has a little gap with BSE, and the result with high order corrections is more accurate. The LO $t_{1}$ is zero that is agree with HQET, see Eq. (26), but far from BSE. The LO $t_{2}$ and $t_{4}$ are different from BSE similarly. The high order corrections bring them closer to BSE, of which 1st correction is the most important one. In the mode $1S\to 1S$, the nonrelativistic results may be approximate, but high order corrections can make the result more precise. Note that, the accurate result of $t_{1}$ cannot be obtained by correcting $\xi_{00}$, as HQET did, so we need to introduce new high-order correction functions.

For the process $B_{c}\to h_{c}$, as Figs. 7(a)–7(d) shows, the form factor $t_{2}$ makes the main contribution to the decay width. The LO $t_{2}$ is slightly different from BSE. The 1st, 2nd and 3rd corrections are not small but almost cancel each other out. These corrections make $t_{2}$ closer to BSE. The LO $t_{3}$ and $t_{4}$ is zero which can be used to examine our method. The $t_{1}$, $t_{2}$ and $t_{4}$ with 1st correction are still a lot different from BSE, and therefore the higher order corrections are necessary. For the process $B_{c}\to\chi_{c0}$, as Figs. 7(e)–7(f) shows, the difference between the LO $S_{+}$ and BSE is large. At least the 1st correction needs to be considered in order to reach an approximate result, though the $S_{-}$ with 1st correction is not accurate enough. The higher order corrections can make the results more accurate. For the process $B_{c}\to\chi_{c1}$, as Figs. 8(a)–8(d) shows, the form factor $t_{2}$ makes main contribution to the decay width. The corrections from 1st IWF are great except $t_{4}$. The higher order corrections are small but still important for accurate calculation. For the process $B_{c}\to\chi_{c2}$, as Figs. 9(a)–9(d) shows, the form factor $t_{3}$ makes main contribution to the decay width. The LO $t_{3}$ is close to BSE, and the high order corrections almost cancel each other out. It leads to the unusually small result of the total relativistic correction. The 1st correction makes $t_{1}$, $t_{2}$ and $t_{4}$ closer to BSE, but makes the main form factor $t_{3}$ farther from BSE. The result may be more imprecise if only the IWF and 1st correction are considered, so the higher-order corrections are very important. At zero recoil, the IWF $\xi_{11}$ is zero and the kinematical factor $(v\cdot v_{f})/|\vec{v}_{f}|$ will lead to the divergence, see Eq. (29). However most LO form factors are limited values at zero recoil. In general, the relativistic corrections are large, and the 1st corrections can only abtain the approximate results in the mode $1S\to 1P$.

In the modes $1S\to 2S$ and $1S\to 2P$, the form factors are no longer kinematically depressed but a little kinematically enhanced. The relativistic corrections are similar to those discussed above, but are greater and more complicated, as Figs. 10-13 shows. Generally, there are big gaps between the leading order form factors and those from BSE directly. The newly introduced high-order correction functions make significant contributions in these relativistic corrections.

VI Conclusion

In this work, we extract the Isgur-Wise functions in the framework of instantaneous Bethe-Salpeter equation. The Isgur-Wise function is the overlapping integrals of the wave functions for the initial state and the final state. The overlapping integrals which are with the relative momentum $\vec{q}\>^{\prime}$ being inserted are the relativistic corrections ($1/m_{q}$ corrections) to the Isgur-Wise function, and the number of $\vec{q}\>^{\prime}$ contained in the function corresponds to the order of the correction. We choose the semileptonic $B_{c}$ decays to charmonium to calculate the numerical results of Isgur-Wise functions and form factors, where the final states include $1S$, $1P$, $2S$ and $2P$. In the mode $1S\to 1S$, the Isgur-Wise function is coincident with the heavy quark effective theory, but it is not strict unity at zero recoil due to the relativistic correction. Another part of relativistic correction comes from the higher order correction functions, which make the results more accurate. The Isgur-Wise function can be generalized to other modes, including but not limited to $1S\to 1P$, $1S\to 2S$ and $1S\to 2P$ which this paper studies. The behavior of Isgur-Wise function almost exclusively depends on the configurations $nL$ of initial and final states, so they can still simplify the calculations of form factors. Some newly introduced higher order corrections provide a great relativistic corrections, even though the initial and final states are both the double-heavy mesons. These corrections can not be obtained by correcting the Isgur-Wise function, as HQET would do, so the higher order correction functions are necessary for accurate calculations.

VII Acknowledgments

This work was supported in part by the National Natural Science Foundation of China (NSFC) under Grant No. 11405037, No. 11575048 and No. 11505039, and also in part by PIRS of HIT No.B201506.

Appendix A Equation and solution for heavy mesons

BS equation for a quark-antiquark bound state generally is written as Chang et al. (2015)

[TABLE]

where $p_{1},p_{2};m_{1},m_{2}$ are the momenta and masses of the quark and antiquark, respectively; $\chi_{P}(q)$ is the BS wave function with the total momentum $P$ and relative momentum $q$; $V(P,k,q)$ is the kernel between the quark-antiquark in the bound state. $P$ and $q$ are defined as

[TABLE]

We divide the relative momentum $q$ into two parts, $q_{P_{||}}$ and $q_{P_{\perp}}$, a parallel part and an orthogonal one to $P$, respectively

[TABLE]

where $q_{P_{||}}^{\mu}\equiv(P\cdot q/M^{2})P^{\mu},~{}q_{P_{\perp}}^{\mu}\equiv q^{\mu}-q_{P_{||}}^{\mu}$, and $M$ is the mass of the relevant meson. Correspondingly, we have two Lorentz-invariant variables

[TABLE]

If we introduce two notations as below

[TABLE]

Then the BS equation can take the form as follow

[TABLE]

The propagators of the relevant particles with masses $m_{1}$ and $m_{2}$ can be decomposed as

[TABLE]

with

[TABLE]

where $i=1,2$ for quark and antiquark, respectively, and $J(i)=(-1)^{i+1}$.

Then the instantaneous Bethe-Salpeter equation can be decomposed into the coupled equations

[TABLE]

The instantaneous kernel has the following form

[TABLE]

especially when the two constituents of meson are very heavy. The kernel we used contains a linear scalar interaction for color-confinement, a vector interaction for one-gluon exchange and a constant $V_{0}$ which as a ‘zero-point’, i.e.

[TABLE]

where $\lambda$ is the ‘string constant’, $\alpha_{s}(r)$ is the running coupling constant. In order to avoid the infrared divergence, a factor $e^{-\alpha r}$ is introduced, i.e.

[TABLE]

In momentum space the kernel reads:

[TABLE]

where

[TABLE]

The fitted parameters are $a=e=2.7183$, $\alpha=0.06$ GeV, $\lambda=0.21$ ${\rm GeV}^{2}$, $\Lambda_{QCD}=0.27$ GeV; $V_{0}$ is fixed by fitting the mass of the ground state. With these parameters, the mass spectrums, decay constants and some branching fractions of the double heavy mesons, including $B_{c}$, charmonium and bottomium, can be obtained. These results are in good accord with the experimental data. More details can be referred to in the literatures Chang and Wang (2010); Fu et al. (2011).

The instantaneous Bethe-Salpeter wave function for $2^{++}$ states mesons have the general formWang (2009)

[TABLE]

with

[TABLE]

The wave function corresponding to the positive projection has the form

[TABLE]

where

[TABLE]

If the masses of the quark and antiquark are equal, the normalization condition reads as

[TABLE]