L(\Gamma,\mathcal{S},s)=-\delta(m_{\alpha},m_{\text{max}})\left[|\mathcal{A}|\log\left(1+\frac{1}{m_{\textrm{max}}}\right)+\log\left(1+\frac{2}{m_{\textrm{max}}}\right)\right]-\delta(m_{\alpha},0)\left(\log|\Gamma|+\log
m_{\text{max}}\right)$ (87)

The update of the encoding of the supernode labels,
$L(\mathcal{V}|H,\mathcal{S})$ (cf. Eq.7), is, again, affected by the graph
size, the growth of the supernode number and the incremental increase of a
graphlet occurrence. Denoting by
$M\equiv\sum_{\alpha^{\prime}}m_{\alpha^{\prime}}=|\mathcal{S}|$ the number of
supernodes,

$\displaystyle\Delta
L(\mathcal{V},s|H,\mathcal{S})=\log\left(\frac{m_{\alpha}+1}{M+1}\right)+\log\binom{N}{M}-\log\binom{N-n+1}{M+1}$
(88)

The update of reconstruction cost from $H$ to $G$ depends on the reduced graph
size, the associated graphlet orientation number and the subgraph’s
neighborhood:

$\displaystyle\Delta L(G,s|H,\mathcal{V},\mathcal{S},\Gamma)$
$\displaystyle=\log\frac{(N-n+1)!}{N!}+\frac{n_{\alpha}!}{|\mathrm{Aut}(\alpha)|}$
(89)
$\displaystyle+\sum_{j\in\mathcal{N}(H)\backslash\mathcal{V}}\log\binom{n}{\sum_{i\in\nu}A_{ij}}\binom{n}{\sum_{i\in\nu
s}A_{ji}}+\sum_{j_{s^{\prime}}\in\mathcal{V}}\log\binom{nn_{j_{s}^{\prime}}}{\sum_{i\in\nu}A_{ij_{s^{\prime}}}}\binom{nn_{j_{s}^{\prime}}}{\sum_{i\in\nu}A_{j_{s^{\prime}}i}},$

where $n_{j_{s}^{\prime}}$ is the subgraph size relative to the supernode
$j_{s^{\prime}}\in\mathcal{V}$, replacing the subgraph $s^{\prime}$. The two
sums represent the encoding of the nodes’ neighbours within $s$, i.e., how to
distribute the multiedge among the nodes that would be deleted. The first sum
correponds to regular node neighbors, while the second sum corresponds to
supernode neighbors.

Finally, the complete putative codelength difference is:

$\Delta L(G,\theta,s)=\Delta L_{\phi}(H,s)+\Delta L(\phi,s)+\Delta
L(G,s|H,\mathcal{V},\mathcal{S},\Gamma)+\Delta
L(\mathcal{V},s|H,\mathcal{S})+\Delta L(\Gamma,\mathcal{S},s)$ (90)

##  Supplementary Note S3 Entropy of simple the reciprocal configuration model

To derive the entropy of the reciprocal configuration model for simple graphs,
we follow the approach developed in [69]. Compared to the directed
configuration model, the reciprocal version has the mutual degree sequence,
i.e., the number of mutual stubs (corresponding to reciprocal edges) per node,
as an additional set of parameters. We let $\mathbf{u}$ denote a vector of
size $N$ filled with ones, and we define the microcanonical partition function
as a sum over a product of Dirac delta functions which define the constrained
parameter values of the model,

$\Omega{(\bm{\kappa}^{m},\bm{\kappa}^{+},\bm{\kappa}^{-})}=\sum_{\mathbf{A}}\delta\left(\bm{\kappa}^{m}-v_{D}\left(\mathbf{AA}^{T}\right)\right)\delta\left(\bm{\kappa}^{+}-\mathbf{Au}+v_{D}\left(\mathbf{AA}^{T}\right)\right)\delta\left(\bm{\kappa}^{-}-\mathbf{A}^{T}\mathbf{u}+v_{D}\left(\mathbf{AA}^{T}\right)\right),$
(91)

where $v_{D}$ is a function that maps the diagonal elements of a $N\times N$
matrix to a $N$-dimensional vector. The Dirac delta functions can be expanded
in terms of Fourier integrals to obtain:

$\displaystyle\Omega{(\bm{\kappa}^{m},\bm{\kappa}^{+},\bm{\kappa}^{-})}$
$\displaystyle=\int\frac{\text{d}\bm{\lambda}^{+}\text{d}\bm{\lambda}^{-}\text{d}\bm{\mu}}{(2\pi)^{3N}}\
e^{{\bm{\lambda}^{+}}^{T}\bm{\kappa}^{+}+{\bm{\lambda}^{-}}^{T}\bm{\kappa}^{-}+\bm{\mu}^{T}\bm{\kappa}^{m}}\sum_{\mathbf{A}}e^{-{\bm{\lambda}^{+}}^{T}\mathbf{Au}-{\bm{\lambda}^{-}}^{T}\mathbf{A}^{T}\mathbf{u}-{\bm{\mu}}^{T}\mathrm{diag}\left(\mathbf{AA}^{T}\right)},$
(92)

where $\bm{\lambda}^{+},\bm{\lambda}^{-},\bm{\mu}$ are identified as vectors
of Lagrange multipliers.

Compared to the classical configuration model, symmetric pairs of elements of
the adjacency matrix are not independent and need to be considered
simultaneously, such that

$\displaystyle\sum_{\mathbf{A}}e^{-{\bm{\lambda}^{+}}^{T}\mathbf{Au}-{\bm{\lambda}^{-}}^{T}\mathbf{A}^{T}\mathbf{u}-{\bm{\mu}}^{T}\mathrm{diag}\left(\mathbf{AA}^{T}\right)}$
$\displaystyle=\sum_{\\{A_{ij},A_{ji}\\}=\\{\\{0,0\\},\\{1,1\\},\\{0,1\\},\\{1,0\\}\\}}e^{-\sum_{ij}(\mu_{i}-\lambda_{i}^{+}-\lambda_{i}^{-})A_{ij}A_{ji}}e^{-\sum_{ij}(\lambda_{i}^{+}+\lambda_{j}^{-})A_{ij}}$
$\displaystyle=\prod_{i<j}\left(1+e^{-\mu_{i}-\mu_{j}}+e^{-\lambda_{i}^{+}-\lambda_{j}^{-}}+e^{-\lambda_{i}^{-}-\lambda_{j}^{+}}\right).$
(93)

For the sake of readability, we set $s_{ij}\equiv
e^{-\mu_{i}-\mu_{j}}+e^{-\lambda_{i}^{+}-\lambda_{j}^{-}}+e^{-\lambda_{j}^{+}-\lambda_{i}^{-}}$.
To estimate $\Omega{(\bm{\kappa}^{m},\bm{\kappa}^{+},\bm{\kappa}^{-})}$, we
apply a Laplace approximation to its Fourier integral form, and we thus we
seek the maximisation of the following quantity

$NQ(\bm{\mu},\bm{\lambda}^{+},\bm{\lambda}^{-}|\bm{\kappa}^{m},\bm{\kappa}^{+},\bm{\kappa}^{-})={\bm{\mu}}^{T}\bm{\kappa}^{m}+{\bm{\lambda}^{+}}^{T}\bm{\kappa}^{+}+{\bm{\lambda}^{-}}^{T}\bm{\kappa}^{-}+\sum_{i<j}\ln(1+s_{ij}).$
(94)

The saddle point equations to be solved are thus

$\displaystyle\frac{\partial Q}{\partial\mu_{i}}=0$
$\displaystyle\Leftrightarrow\kappa^{m}_{i}=e^{-\mu_{i}}\sum_{j\neq
i}\frac{e^{-\mu_{j}}}{1+s_{ij}},$ (95) $\displaystyle\frac{\partial
Q}{\partial\lambda_{i}^{+}}=0$
$\displaystyle\Leftrightarrow\kappa^{+}_{i}=e^{-\lambda_{i}^{+}}\sum_{j\neq
i}\frac{e^{-\lambda_{j}^{-}}}{1+s_{ij}},$ (96) $\displaystyle\frac{\partial
Q}{\partial\lambda_{i}^{-}}=0$
$\displaystyle\Leftrightarrow\kappa^{-}_{i}=e^{-\lambda_{i}^{-}}\sum_{j\neq
i}\frac{e^{-\lambda_{j}^{+}}}{1+s_{ij}}.$ (97)

Here, we are only interested in the sparse graph approximation. The
computation of the Hessian would be compulsory in a rigorous calculation.
However, when actually evaluated, it only leads to sums of log terms in the
degree sequences elements, which are negligible compared to the sums of log
factorial terms in the degree sequence elements (when the graph is of finite
size). This is the same observation that is found in [69], while never
explicitly justified there. For our results to be consistent with the standard
form of other simple network entropy expressions, we also choose to neglect
the contribution of the Hessian to the sparse graph approximation of the
microcanonical partition function.

At large $N$, it is reasonable to assume $s_{ij}=o(1)$ for the right hand
terms of the saddle point equations to be finite. This leads to the simplified
equations

$\displaystyle\kappa^{m}_{i}$ $\displaystyle=K_{m}e^{-\mu_{i}},$
$\displaystyle K_{m}=\sum_{j}e^{-\mu_{j}}$ (98) $\displaystyle\kappa^{+}_{i}$
$\displaystyle=K_{-}e^{-\lambda^{+}_{i}},$ $\displaystyle
K_{-}=\sum_{j}e^{-\lambda_{j}^{-}}$ (99) $\displaystyle\kappa^{-}_{i}$
$\displaystyle=K_{+}e^{-\lambda^{-}_{i}},$ $\displaystyle
K_{+}=\sum_{j}e^{-\lambda_{j}^{+}}$ (100)

The constants $K_{m},K_{+},K_{-}$ are determined by global structural
constraints, which are the number of asymmetrically connected (or directed)
pairs of nodes and the number of mutual edges:

$\displaystyle 2E_{m}$
$\displaystyle=\sum_{i}\kappa^{m}_{i}=K_{m}^{2}\Leftrightarrow
K_{m}=\sqrt{2E_{m}},$ (101) $\displaystyle E_{d}$
$\displaystyle=\sum_{i}\kappa^{+}_{i}=K_{+}K_{-}\Leftrightarrow
K_{+}=K_{-}=\sqrt{E_{d}}.$ (102)

All is now set for a second-order estimation of the microcanonical partition
function in $s_{ij}$. We have

$\displaystyle\bm{\mu}^{T}\bm{\kappa}^{m}$
$\displaystyle=-\sum_{i}\kappa^{m}_{i}\ln\kappa^{m}_{i}+E_{m}\ln(2E_{m})\approx\ln\frac{(2E_{m})!!}{\prod_{i}\kappa^{m}_{i}!}-E_{m},$
(103)
$\displaystyle{\bm{\lambda}^{+}}^{T}\bm{\kappa}^{+}+{\bm{\lambda}^{-}}^{T}\bm{\kappa}^{-}$
$\displaystyle=-\sum_{i}\kappa^{+}_{i}\ln\kappa^{+}_{i}-\sum_{i}\kappa^{-}_{i}\ln\kappa^{-}_{i}+E_{d}\ln
E_{d}\approx\ln\frac{E_{d}!}{\prod_{i}\kappa^{+}_{i}!\kappa^{-}_{i}!}-E_{d},,$
(104) $\displaystyle s_{ij}$
$\displaystyle=\frac{\kappa^{m}_{i}\kappa^{m}_{j}}{2E_{m}}+\frac{\kappa^{+}_{i}\kappa^{-}_{j}+\kappa^{-}_{i}\kappa^{+}_{j}}{E_{d}},$
(105)

and

$\displaystyle\frac{1}{2}\sum_{ij}\left(s_{ij}-\frac{s_{ij}^{2}}{2}\right)$
$\displaystyle=E_{m}+E_{d}-\frac{1}{2}\underbrace{\left(\frac{1}{2}\frac{\langle{\kappa^{m}_{i}}^{2}\rangle^{2}}{\langle\kappa^{m}_{i}\rangle^{2}}+\frac{\langle{\kappa^{+}_{i}}^{2}\rangle\langle{\kappa^{-}_{i}}^{2}\rangle}{\langle\kappa^{+}_{i}\rangle\langle\kappa^{-}_{i}\rangle}+\frac{\langle\kappa^{+}_{i}\kappa^{-}_{i}\rangle^{2}}{\langle\kappa^{+}_{i}\rangle\langle\kappa^{-}_{i}\rangle}+\frac{\langle\kappa^{m}_{i}\kappa^{+}_{i}\rangle\langle\kappa^{m}_{i}\kappa^{-}_{i}\rangle}{\langle\kappa^{m}_{i}\rangle\langle\kappa^{+}_{i}\rangle}\right)}_{\Psi{(\bm{\kappa}^{m},\bm{\kappa}^{+},\bm{\kappa}^{-})}}.$
(106)

Putting it all together we obtain

$L_{(\bm{\kappa}^{m},\bm{\kappa}^{+},\bm{\kappa}^{-})}(G)=\log\frac{(2E_{m})!!}{\prod_{i}\kappa^{m}_{i}!}+\log\frac{E_{d}!}{\prod_{i}\kappa^{+}_{i}!\kappa^{-}_{i}!}-\frac{1}{2\ln
2}\Psi{(\bm{\kappa}^{m},\bm{\kappa}^{+},\bm{\kappa}^{-})}.$ (107)

Similarly to the sparse approximation of entropy of the simple graph
configuration model [69], we see that the entropy for the simple graph
reciprocal configuration model amounts to the multigraph codelength from which
is subtracted a functional cut-off that depends on the statistics of the
degree sequences.

##  Supplementary Note S4 Motif mining based on hypothesis testing

Given a graphlet $\alpha$, hypothesis-based motif inference qualifies $\alpha$
as a network motif in an empirical network $G$ if its frequency
$f_{\alpha}(G)$ in $G$ is significantly greater than in an ensemble of random
networks $\mathcal{G}_{\theta}$ sampled from a null model $P_{\theta}$.

For uniformly sampling simple random networks, we use the shuffling algorithms
described in Randomized networks in the Methods section. When the edge
swapping procedures are ergodic and unbiased, they are guaranteed to uniformly
sample the corresponding ensembles of random networks after a large enough
number of swaps [28]. However, the mixing time, i.e., the number of swaps
needed for the generated networks to be practically independent, is not known
in general [28]. To ensure that correlations between randomized networks are
not likely to influence results (and thus favor hypothesis-testing based
methods as much as possible), we perform $100E$ successful edge swaps to
generate each random network. This does not guarantee an absence of
correlations, but we note that the number of swaps is larger than what is
typically prescribed in the literature (for reference 0.2E edge-swaps were
used to generate each random network in [16], 3E in [101] and 6E in [60]).

We utilize the typical normality assumption of the graphlet frequencies under
the null and employ as test statistic the $Z$-score given by

$\displaystyle Z_{\alpha,\theta}(G)$
$\displaystyle=\frac{f_{\alpha}(G)-\mu_{\alpha,\theta}}{{\sigma}_{\alpha,\theta}},$
(108)

where

$\displaystyle\mu_{\alpha,\theta}$
$\displaystyle=\frac{1}{|\mathcal{G}_{\theta}|}\sum_{G^{\prime}\in\mathcal{G}_{\theta}}f_{\alpha}(G^{\prime})$
(109)

and

$\displaystyle{\sigma}^{2}_{\alpha,\theta}$
$\displaystyle=\frac{1}{|\mathcal{G}_{\theta}|-1}\sum_{G^{\prime}\in\mathcal{G}_{\theta}}\left[f_{\alpha}^{2}(G^{\prime})-{\mu}^{2}_{\alpha,\theta}\right].$
(110)

In all experiments, the size of $\mathcal{G}_{\theta}$ is set to 100 and the
significance threshold (nominal alpha-level) is fixed at 0.01. To correct for
multiple testing (one test for each graphlet), we employ a Bonferroni
correction, which multiplies the raw $p$-values obtained directly from the
$Z$-scores by $|\Gamma|\approx 10^{4}$. As displayed in Fig. 3A–D, depending
on the choice of the null, a considerable number of motifs can be falsely
detected. A similar effect can also be seen in empirical data, where the
number of motifs found varies enormously with the choice of null model
(Supplementary Fig. S8), even though we corrected for multiple testing with
the maximally conservative Bonferroni correction. Supplementary Fig. S8 also
demonstrates that the motifs found vary significantly depending on the null,
and that the smallest number of motifs is not necessarily found under the most
restricted null hypothesis.

##  Supplementary Note S5 Measures of graphlet topology

The _density_ $\rho$ measures the ratio fraction of node pairs in a simple
graph $G=(N,E)$ that are connected by an edge [81],

$\rho=\frac{E}{N(N-1)}.$ (111)

A simple graph is said to be _sparse_ if its density is close to zero, and
_dense_ if its density is close to one.

The _reciprocity_ $r$ measures the fraction of edges in a graph $G$ that are
reciprocated [81],

$r=\frac{1}{E}\sum_{ij}A_{ij}A_{ji},$ (112)

where $A_{ij}\in\\{0,1\\}$ are the entries of the adjacency matrix for $G$.

The _number of cycles_ of a simple graph $G$ is the number of closed paths in
$G$ where no node appears twice [81].

The _graph polynomial root_ GPR is a measure of the symmetry of a graph. It is
related to the so-called orbit-polynomial [80] $\Pi_{G}(z)$ and allows a
ranking of graphs based on the distribution of their orbit sizes. Let
$c_{o_{l}}$ be the number of orbits of size $o_{l}$, where
$l\in\\{1,2,\ldots,L\\}$ and $L$ is the number of different orbit sizes. The
graph polynomial is then defined as

$\Pi_{G}(z)=\sum_{l=1}^{L}c_{o_{l}}z^{o_{l}}.$ (113)

The GPR, denoted $z^{*}$, is the unique solution of the following equation

$\Pi_{G}(z^{*})=1,$ (114)

which can be solved numerically. Orbit sizes are determined using McKay’s
nauty algorithm [102]. A strong degree of symmetry is affiliated with a high
GPR, while an asymmetric structure corresponds to a low GPR. Supplementary
Figure S9 shows the values of the GPR of all 9 364 three- to five-node
graphlets.


# The ${\Upsilon}(1S)$ ${\to}$ $B_{c}{\rho}$ decay with perturbative QCD
approach

Junfeng Sun  Institute of Particle and Nuclear Physics, Henan Normal
University, Xinxiang 453007, China  Yueling Yang  Institute of Particle and
Nuclear Physics, Henan Normal University, Xinxiang 453007, China  Qingxia Li
Institute of Particle and Nuclear Physics, Henan Normal University, Xinxiang
453007, China  Gongru Lu  Institute of Particle and Nuclear Physics, Henan
Normal University, Xinxiang 453007, China  Jinshu Huang  College of Physics
and Electronic Engineering, Nanyang Normal University, Nanyang 473061, China
Qin Chang  Institute of Particle and Nuclear Physics, Henan Normal University,
Xinxiang 453007, China

###### Abstract

With the potential prospects of the ${\Upsilon}(1S)$ data samples at the
running LHC and upcoming SuperKEKB, the ${\Upsilon}(1S)$ ${\to}$ $B_{c}{\rho}$
weak decay is studied with the pQCD approach. It is found that (1) the lion’s
share of branching ratio comes from the longitudinal polarization helicity
amplitudes; (2) branching ratio for the ${\Upsilon}(1S)$ ${\to}$ $B_{c}{\rho}$
decay can reach up to ${\cal O}(10^{-9})$, which might be hopefully
measurable.

###### pacs:

13.25.Gv 12.39.St 14.40.Pq

##  I Introduction

The ${\Upsilon}(1S)$ meson consists of the bottom quark and antiquark pair
$b\bar{b}$, carries the definitely established quantum numbers of
$I^{G}J^{PC}$ $=$ $0^{-}1^{--}$ pdg , and lies below the kinematic $B\bar{B}$
threshold. The ${\Upsilon}(1S)$ meson decay mainly through the strong
interaction, the electromagnetic interaction and radiative transition.
Besides, the ${\Upsilon}(1S)$ meson can also decay via the weak interactions
within the standard model. More than $10^{8}$ ${\Upsilon}(1S)$ data samples
have been accumulated at Belle epjc74 . More and more upsilon data samples
with high precision are promisingly expected at the running LHC and the
forthcoming SuperKEKB. Although the branching ratio for the ${\Upsilon}(1S)$
weak decay is tiny, it seems to exist a realistic possibility to search for
the signals of the ${\Upsilon}(1S)$ weak decay at future experiments. In this
paper, we will study the ${\Upsilon}(1S)$ ${\to}$ $B_{c}{\rho}$ weak decay
with the perturbative QCD (pQCD) approach pqcd1 ; pqcd2 ; pqcd3 .

Experimentally, there is no report on the ${\Upsilon}(1S)$ ${\to}$
$B_{c}{\rho}$ weak decay so far. The signals for the ${\Upsilon}(1S)$ ${\to}$
$B_{c}{\rho}$ weak decay should, in principle, be easily identified, due to
the facts that the final states have different electric charges, have definite
momentum and energy, and are back-to-back in the rest frame of the
${\Upsilon}(1S)$ meson. In addition, the identification of a single flavored
$B_{c}$ meson could be used to effectively enhance signal-to-background ratio.
Another important and fashionable motivation is that evidences of an
abnormally large branching ratio for the ${\Upsilon}(1S)$ weak decay might be
a hint of new physics.

Theoretically, the ${\Upsilon}(1S)$ ${\to}$ $B_{c}{\rho}$ weak decay belongs
to the external $W$ emission topography, and is favored by the Cabibbo-
Kabayashi-Maskawa (CKM) matrix elements ${|}V_{cb}V_{ud}^{\ast}{|}$. So it
should have relatively large branching ratio among the ${\Upsilon}(1S)$ weak
decays, which has been studied with the naive factorization (NF) approximation
ijma14 ; adv2013 . Recently, some attractive methods have been developed, such
as the pQCD approach pqcd1 ; pqcd2 ; pqcd3 , the QCD factorization approach
qcdf1 ; qcdf2 ; qcdf3 , soft and collinear effective theory scet1 ; scet2 ;
scet3 ; scet4 , and applied widely to accommodate measurements on the $B$
meson weak decays. The ${\Upsilon}(1S)$ ${\to}$ $B_{c}{\rho}$ decay permit one
to cross check parameters obtained from the $B$ meson decay, to test the
practical applicability of various phenomenological models in the vector meson
weak decays, and to further explore the underlying dynamical mechanism of the
heavy quark weak decay. In addition, as it is well known, the $B_{c}$ meson
carries two explicit heavy flavors and has extremely abundant decay modes, but
its hadronic production is suppressed compared with that for hidden-flavor
quarkonia and heavy-light mesons, due to higher order in QCD coupling
constants ${\alpha}_{s}$ and the presence of additional heavy quarks bc1 ; bc2
. The ${\Upsilon}(1S)$ ${\to}$ $B_{c}{\rho}$ decay offers another platform to
study the $B_{c}$ meson production at high energy colliders.

This paper is organized as follows. In section II, we present the theoretical
framework and the amplitudes for the ${\Upsilon}(1S)$ ${\to}$ $B_{c}{\rho}$
decay with the pQCD approach. Section III is devoted to numerical results and
discussion. The last section is our summary.

##  II theoretical framework

###  II.1 The effective Hamiltonian

The effective Hamiltonian responsible for the ${\Upsilon}(1S)$ ${\to}$
$B_{c}{\rho}$ weak decay is 9512380

${\cal H}_{\rm eff}\ =\
\frac{G_{F}}{\sqrt{2}}\,V_{cb}V_{ud}^{\ast}\,\Big{\\{}C_{1}({\mu})\,Q_{1}({\mu})+C_{2}({\mu})\,Q_{2}({\mu})\Big{\\}}+{\rm
H.c.},$ (1)

where $G_{F}$ ${\simeq}$ $1.166{\times}10^{-5}\,{\rm GeV}^{-2}$ pdg  is the
Fermi coupling constant; the CKM factor is written as a power series in the
Wolfenstein parameter ${\lambda}$ ${\simeq}$ $0.2$ pdg ,

$V_{cb}V_{ud}^{\ast}\ =\
A{\lambda}^{2}-\frac{1}{2}A{\lambda}^{4}-\frac{1}{8}A{\lambda}^{6}+{\cal
O}({\lambda}^{8}).$ (2)

The local operators are defined as follows:

$\displaystyle Q_{1}$ $\displaystyle=$
$\displaystyle[\bar{c}_{\alpha}{\gamma}_{\mu}(1-{\gamma}_{5})b_{\alpha}][\bar{q}_{\beta}{\gamma}^{\mu}(1-{\gamma}_{5})u_{\beta}],$
(3) $\displaystyle Q_{2}$ $\displaystyle=$
$\displaystyle[\bar{c}_{\alpha}{\gamma}_{\mu}(1-{\gamma}_{5})b_{\beta}][\bar{q}_{\beta}{\gamma}^{\mu}(1-{\gamma}_{5})u_{\alpha}],$
(4)

where ${\alpha}$ and ${\beta}$ are color indices.

From Eq.(1), it is clearly seen that only the tree operators contribute to the
concerned process, and there is no pollution from penguin and annihilation
contributions. As it is well known, degrees of freedom with mass scales above
${\mu}$ are integrated out into the Wilson coefficients $C_{1,2}(\mu)$
typically using the renormalization group assisted perturbation theory. The
physical contributions below the scale of ${\mu}$ are included in the hadronic
matrix elements (HME) where the local operators sandwiched between initial and
final hadron states. The most complicated part is the treatment on HME, where
the perturbative and nonperturbative effects entangle with each other. To
obtain the decay amplitudes, the remaining work is to calculate HME properly.

###  II.2 Hadronic matrix elements

With the Lepage-Brodsky approach for exclusive processes prd22 , HME could be
expressed as the convolution of hard scattering subamplitudes containing
perturbative contributions with the universal wave functions reflecting the
nonperturbative contributions. To eliminate the endpoint singularities
appearing in the collinear factorization approximation, the pQCD approach
suggests pqcd1 ; pqcd2 ; pqcd3  retaining the transverse momentum of quarks
and introducing the Sudakov factor. Finally, the decay amplitudes could be
factorized into three parts pqcd2 ; pqcd3 : the hard effects enclosed by the
Wilson coefficients $C_{i}$, the heavy quark decay subamplitudes ${\cal H}$,
and the universal wave functions ${\Phi}$,

${\int}dk\,C_{i}(t)\,{\cal H}(t,k)\,{\Phi}(k)\,e^{-S},$ (5)

where $t$ is a typical scale, $k$ is the momentum of the valence quarks, and
the Sudakov factor $e^{-S}$ can effectively suppress the long-distance
contributions and make the hard scattering more perturbative.

###  II.3 Kinematic variables

The light cone kinematic variables in the ${\Upsilon}(1S)$ rest frame are
defined as follows:

$p_{{\Upsilon}}\,=\,p_{1}\,=\,\frac{m_{1}}{\sqrt{2}}(1,1,0),$ (6)
$p_{B_{c}}\,=\,p_{2}\,=\,(p_{2}^{+},p_{2}^{-},0),$ (7)
$p_{\rho}\,=\,p_{3}\,=\,(p_{3}^{-},p_{3}^{+},0),$ (8)
$k_{i}\,=\,x_{i}\,p_{i}+(0,0,\vec{k}_{i{\perp}}),$ (9)
${\epsilon}_{i}^{\parallel}\,=\,\frac{p_{i}}{m_{i}}-\frac{m_{i}}{p_{i}{\cdot}n_{+}}n_{+},$
(10) ${\epsilon}_{i}^{\perp}\,=\,(0,0,\vec{1}),$ (11) $n_{+}=(1,0,0),$ (12)
$p_{i}^{\pm}\,=\,(E_{i}\,{\pm}\,p)/\sqrt{2},$ (13)
$s\,=\,2\,p_{2}{\cdot}p_{3},$ (14) $t\,=\,2\,p_{1}{\cdot}p_{2}\,=\
2\,m_{1}\,E_{2},$ (15) $u\,=\,2\,p_{1}{\cdot}p_{3}\,=\ 2\,m_{1}\,E_{3},$ (16)
$p=\frac{\sqrt{[m_{1}^{2}-(m_{2}+m_{3})^{2}]\,[m_{1}^{2}-(m_{2}-m_{3})^{2}]}}{2\,m_{1}},$
(17)

where $x_{i}$ and $\vec{k}_{i{\perp}}$ are the longitudinal momentum fraction
and transverse momentum of the valence quark, respectively;
${\epsilon}_{i}^{\parallel}$ and ${\epsilon}_{i}^{\perp}$ are the longitudinal
and transverse polarization vectors, respectively, satisfying with the
relations ${\epsilon}_{i}^{2}$ $=$ $-1$ and ${\epsilon}_{i}{\cdot}p_{i}$ $=$
$0$; the subscript $i$ $=$ $1$, $2$, $3$ on variables ($p_{i}$, $E_{i}$,
$m_{i}$ and ${\epsilon}_{i}^{{\parallel},{\perp}}$) correspond to the
${\Upsilon}(1S)$, $B_{c}$ and ${\rho}$ mesons, respectively; $n_{+}$ is the
null vector; $s$, $t$ and $u$ are the Lorentz-invariant variables; $p$ is the
common momentum of final states. The notation of momentum is displayed in
Fig.1(a).

###  II.4 Wave functions

With the notation in prd65 ; jhep0605 , the definitions of the diquark
operator HME are

${\langle}0{|}b_{i}(z)\bar{b}_{j}(0){|}{\Upsilon}(p_{1},{\epsilon}_{1}^{\parallel}){\rangle}\,=\,\frac{f_{\Upsilon}}{4}{\int}d^{4}k_{1}\,e^{-ik_{1}{\cdot}z}\Big{\\{}\\!\\!\not{\epsilon}_{1}^{{\parallel}}\Big{[}m_{1}\,{\Phi}_{\Upsilon}^{v}(k_{1})-\\!\\!\not{p}_{1}\,{\Phi}_{\Upsilon}^{t}(k_{1})\Big{]}\Big{\\}}_{ji},$
(18)
${\langle}0{|}b_{i}(z)\bar{b}_{j}(0){|}{\Upsilon}(p_{1},{\epsilon}_{1}^{{\perp}}){\rangle}\,=\,\frac{f_{\Upsilon}}{4}{\int}d^{4}k_{1}\,e^{-ik_{1}{\cdot}z}\Big{\\{}\\!\\!\not{\epsilon}_{1}^{{\perp}}\Big{[}m_{1}\,{\Phi}_{\Upsilon}^{V}(k_{1})-\\!\\!\not{p}_{1}\,{\Phi}_{\Upsilon}^{T}(k_{1})\Big{]}\Big{\\}}_{ji},$
(19)
${\langle}B_{c}(p_{2}){|}\bar{c}_{i}(z)b_{j}(0){|}0{\rangle}\,=\,\frac{i}{4}f_{B_{c}}{\int}dx_{2}\,e^{ix_{2}p_{2}{\cdot}z}\,\Big{\\{}{\gamma}_{5}\Big{[}\\!\\!\not{p}_{2}+m_{2}\Big{]}{\phi}_{B_{c}}(x_{2})\Big{\\}}_{ji},$
(20)
${\langle}{\rho}(p_{3},{\epsilon}_{3}^{{\parallel}}){|}u_{i}(0)\bar{d}_{j}(z){|}0{\rangle}\,=\,\frac{1}{4}{\int}_{0}^{1}dk_{3}\,e^{ik_{3}{\cdot}z}\Big{\\{}\\!\\!\not{\epsilon}_{3}^{{\parallel}}\,m_{3}\,{\Phi}_{\rho}^{v}(k_{3})+\\!\\!\not{\epsilon}_{3}^{{\parallel}}\\!\\!\not{p}_{3}\,{\Phi}_{\rho}^{t}(k_{3})+m_{3}\,{\Phi}_{\rho}^{s}(k_{3})\Big{\\}}_{ji},$
(21)
$\displaystyle{\langle}{\rho}(p_{3},{\epsilon}_{3}^{{\perp}}){|}u_{i}(0)\bar{d}_{j}(z){|}0{\rangle}\
=\
\frac{1}{4}{\int}_{0}^{1}dk_{3}\,e^{ik_{3}{\cdot}z}\Big{\\{}\\!\not{\epsilon}_{3}^{{\perp}}m_{3}\,{\Phi}_{\rho}^{V}(k_{3})$
(22)
$\displaystyle+\\!\not{\epsilon}_{3}^{{\perp}}\\!\not{p}_{3}\,{\Phi}_{\rho}^{T}(k_{3})+\frac{i\,m_{3}}{p_{3}{\cdot}n_{+}}{\varepsilon}_{{\mu}{\nu}{\alpha}{\beta}}\,{\gamma}_{5}\,{\gamma}^{\mu}\,{\epsilon}_{3}^{{\perp},{\nu}}\,p_{3}^{\alpha}\,n_{+}^{\beta}\,{\Phi}_{\rho}^{A}(k_{3})\Big{\\}}_{ji},$

where $f_{\Upsilon}$ and $f_{B_{c}}$ are decay constants; the definitions of
wave functions ${\Phi}_{\rho}^{v,t,s}$ and ${\Phi}_{\rho}^{V,T,A}$ can be
found in Ref. prd65 ; jhep0605 . In fact, for the ${\rho}$ meson, only three
wave functions ${\Phi}_{\rho}^{v}$ and ${\Phi}_{\rho}^{V,A}$ are involved in
the decay amplitudes (see Appendix A). The twist-2 distribution amplitude for
the longitudinal polarization ${\rho}$ meson is prd65 ; jhep0605 :

${\phi}_{\rho}^{v}(x)\,=\,f_{\rho}\,6\,x\,\bar{x}\sum\limits_{i=0}a^{\parallel}_{2i}\,C_{2i}^{3/2}(t),$
(23)

where $f_{\rho}$ is the decay constant; $\bar{x}$ $=$ $1$ $-$ $x$; $t$ $=$
$\bar{x}$ $-$ $x$; $a^{\parallel}_{i}$ and $C_{i}^{3/2}(t)$ are the Gegenbauer
moment and polynomial, respectively; $a^{\parallel}_{i}$ $=$ $0$ for odd $i$
due to the $G$-parity invariance of the ${\rho}$ distribution amplitudes. As
to the twist-3 distribution amplitudes of the transverse polarization ${\rho}$
meson, for simplicity, we will take their asymptotic forms prd65 ; jhep0605 :

${\phi}_{\rho}^{V}(x)\,=\,f_{\rho}\,\frac{3}{4}\,(1+t^{2}),$ (24)
${\phi}_{\rho}^{A}(x)\,=\,f_{\rho}\,\frac{3}{2}\,(-t).$ (25)

Because of $m_{{\Upsilon}(1S)}$ ${\simeq}$ $2m_{b}$ and $m_{B_{c}}$ ${\simeq}$
$m_{b}$ $+$ $m_{c}$, both ${\Upsilon}(1S)$ and $B_{c}$ systems are nearly
nonrelativistic. Nonrelativistic quantum chromodynamics (NRQCD) prd46 ; prd51
; rmp77  and Schrödinger equation can be used to describe their spectrum. The
eigenfunction of the time-independent Schrödinger equation with scalar
harmonic oscillator potential corresponding to the quantum numbers $nL$ $=$
$1S$ is written as

${\phi}(\vec{k})\ {\sim}\ e^{-\vec{k}^{2}/2{\beta}^{2}},$ (26)

where parameter ${\beta}$ determines the average transverse momentum, i.e.,
${\langle}1S{|}\vec{k}^{2}_{\perp}{|}1S{\rangle}$ $=$ ${\beta}^{2}$. Employing
the Brodsky-Huang-Lepage ansatz huang ; xiao  which has been used to structure
wave functions for light and heavy mesons wu2 ,

$\vec{k}^{2}\ {\to}\
\frac{1}{4}\sum\limits_{i}\frac{\vec{k}_{i\perp}^{2}+m_{q_{i}}^{2}}{x_{i}},$
(27)

where $x_{i}$, $\vec{k}_{i\perp}$, $m_{q_{i}}$ are the longitudinal momentum
fraction, transverse momentum, mass of the valence quarks in hadrons,
respectively, with the relations ${\sum}x_{i}$ $=$ $1$ and
$\sum\vec{k}_{i\perp}$ $=$ $0$, then integrating out $\vec{k}_{i\perp}$ and
combining with their asymptotic forms, one can obtain prd65 ; fbb

${\phi}_{B_{c}}(x)=A\,x\bar{x}\,{\exp}\Big{\\{}-\frac{\bar{x}\,m_{c}^{2}+x\,m_{b}^{2}}{8\,{\beta}_{2}^{2}\,x\,\bar{x}}\Big{\\}},$
(28)
${\phi}_{\Upsilon}^{v}(x)={\phi}_{\Upsilon}^{T}(x)=B\,x\bar{x}\,{\exp}\Big{\\{}-\frac{m_{b}^{2}}{8\,{\beta}_{1}^{2}\,x\,\bar{x}}\Big{\\}},$
(29)
${\phi}_{\Upsilon}^{t}(x)=C\,t^{2}\,{\exp}\Big{\\{}-\frac{m_{b}^{2}}{8\,{\beta}_{1}^{2}\,x\,\bar{x}}\Big{\\}},$
(30)
${\phi}_{\Upsilon}^{V}(x)=D\,(1+t^{2})\,{\exp}\Big{\\{}-\frac{m_{b}^{2}}{8\,{\beta}_{1}^{2}\,x\,\bar{x}}\Big{\\}},$
(31)

where the exponential function represents the transverse momentum distribution
and can suppress the end-point singularity; ${\beta}_{i}$ ${\simeq}$
${\xi}_{i}{\alpha}_{s}({\xi}_{i})$ with ${\xi}_{i}$ $=$ $m_{i}/2$ based on the
NRQCD power counting rules prd46 ; parameters $A$, $B$, $C$, $D$ are the
normalization coefficients satisfying the conditions

${\int}_{0}^{1}dx\,{\phi}_{B_{c}}(x)=1,\quad{\int}_{0}^{1}dx\,{\phi}_{\Upsilon}^{v,t}(x)={\int}_{0}^{1}dx\,{\phi}_{\Upsilon}^{V,T}(x)=1.$
(32)

The shape lines for the normalized distribution amplitudes of
${\phi}_{B_{c}}(x)$ and ${\phi}_{\Upsilon}^{v,t,V,T}(x)$ have been displayed
in Fig.1 of Ref.sun16 , from which one can see that Eqs.(28)-(31) reflect
generally the feature that valence quarks of hadrons share momentum fractions
according to their masses.

Figure 1: Feynman diagrams for the ${\Upsilon}$ ${\to}$ $B_{c}{\rho}$ decay
with the pQCD approach, where (a) and (b) are factorizable emission diagrams,
(c) and (d) are nonfactorizable emission diagrams.

###  II.5 Decay amplitudes

The Feynman diagrams for the ${\Upsilon}(1S)$ ${\to}$ $B_{c}{\rho}$ decay are
shown in Fig.1, including factorizable emission topologies (a) and (b) where
gluon connects to the quarks in the same meson, and nonfactorizable emission
topologies (c) and (d) where gluon attaches to the quarks in two different
mesons.

The amplitude for the ${\Upsilon}(1S)$ ${\to}$ $B_{c}{\rho}$ decay is defined
as below prd66 ,

${\cal A}({\Upsilon}(1S){\to}B_{c}{\rho})\ =\ {\cal
A}_{L}({\epsilon}_{1}^{{\parallel}},{\epsilon}_{3}^{{\parallel}})+{\cal
A}_{N}({\epsilon}_{1}^{{\perp}},{\epsilon}_{3}^{{\perp}})+i\,{\cal
A}_{T}\,{\varepsilon}_{{\mu}{\nu}{\alpha}{\beta}}\,{\epsilon}_{1}^{{\mu}}\,{\epsilon}_{3}^{{\nu}}\,p_{1}^{\alpha}\,p_{3}^{\beta},$
(33)

which is conventionally written as the helicity amplitudes prd66 ,

${\cal A}_{0}\ =\ -C_{\cal A}\,\sum\limits_{i}{\cal
A}_{L}^{i}({\epsilon}_{1}^{{\parallel}},{\epsilon}_{3}^{{\parallel}}),$ (34)
${\cal A}_{\parallel}\ =\ \sqrt{2}\,C_{\cal A}\sum\limits_{i}{\cal
A}_{N}^{i}({\epsilon}_{1}^{{\perp}},{\epsilon}_{3}^{{\perp}}),$ (35) ${\cal
A}_{\perp}\ =\ \sqrt{2}\,C_{\cal A}\,m_{1}\,p\sum\limits_{i}{\cal A}_{T}^{i},$
(36) $C_{\cal A}\ =\
i\frac{G_{F}}{\sqrt{2}}\,\frac{C_{F}}{N}\,{\pi}\,f_{\Upsilon}\,f_{B_{c}}\,V_{cb}V_{ud}^{\ast},$
(37)

where $C_{F}$ $=$ $4/3$ and the color number $N$ $=$ $3$; the superscript $i$
on ${\cal A}_{L,N,T}^{i}$ corresponds to the indices of Fig.1. The explicit
expressions of building blocks ${\cal A}_{L,N,T}^{i}$ are collected in
Appendix A.

##  III Numerical results and discussion

In the rest frame of the ${\Upsilon}(1S)$ meson, branching ratio (${\cal
B}r$), polarization fractions ($f_{0,{\parallel},{\perp}}$) and relative phase
between helicity amplitudes (${\phi}_{{\parallel},{\perp}}$) for the
${\Upsilon}(1S)$ ${\to}$ $B_{c}{\rho}$ weak decay are defined as

${\cal B}r\ =\
\frac{1}{12{\pi}}\,\frac{p}{m_{{\Upsilon}}^{2}{\Gamma}_{{\Upsilon}}}\,\Big{\\{}{|}{\cal
A}_{0}{|}^{2}+{|}{\cal A}_{\parallel}{|}^{2}+{|}{\cal
A}_{\perp}{|}^{2}\Big{\\}},$ (38) $f_{0,{\parallel},{\perp}}\ =\
\frac{{|}{\cal A}_{0,{\parallel},{\perp}}{|}^{2}}{{|}{\cal
A}_{0}{|}^{2}+{|}{\cal A}_{\parallel}{|}^{2}+{|}{\cal A}_{\perp}{|}^{2}},$
(39) ${\phi}_{{\parallel},{\perp}}\ =\ {\arg}({\cal
A}_{{\parallel},{\perp}}/{\cal A}_{0}),$ (40)

where mass $m_{{\Upsilon}}$ $=$ $9460.30{\pm}0.26$ MeV and decay width
${\Gamma}_{\Upsilon}$ $=$ $54.02{\pm}1.25$ keV pdg .

The values of other input parameters are listed as follows. If not specified
explicitly, we will take their central values as default inputs.

(1) Wolfenstein parameters pdg : $A$ $=$ $0.814^{+0.023}_{-0.024}$ and
${\lambda}$ $=$ $0.22537{\pm}0.00061$.

(2) Masses of quarks pdg : $m_{c}$ $=$ $1.67{\pm}0.07$ GeV and $m_{b}$ $=$
$4.78{\pm}0.06$ GeV.

(3) Gegenbauer moments111$a^{\parallel}_{0}$ $=$ $1$ is due to the
normalization condition ${\int}_{0}^{1}{\phi}_{\rho}^{v}(x)dx$ $=$ $1$. More
discussion on the ${\rho}$ wave functions and Gegenbauer moments
$a^{\parallel}_{2}$ can be found in the recent references, such as Ref.wu3 .
$a^{\parallel}_{0}$ $=$ $1$ and $a^{\parallel}_{2}$ $=$ $0.15{\pm}0.07$ for
twist-2 distribution amplitudes of the ${\rho}$ meson jhep0605 .

(4) Decay constants: $f_{{\Upsilon}}$ $=$ $(676.4{\pm}10.7)$ MeV fbb ,
$f_{B_{c}}$ $=$ $489{\pm}5$ MeV fbc , $f_{\rho}$ $=$ $216{\pm}3$ MeV jhep0605
.

Our numerical results are presented as follows:

${\cal
B}r\,=\,(8.34^{+0.47+1.35+0.40+1.44}_{-0.69-0.88-0.40-1.26}){\times}10^{-9},$
(41) $f_{0}\,=\,(82.2^{+0.0+1.1+0.0}_{-0.7-1.3-0.0})\%,$ (42)
$f_{{\parallel}}\,=\,(15.0^{+0.6+1.0+0.0}_{-0.0-0.8-0.0})\%,$ (43)
$f_{\perp}\,=\,(2.8^{+0.1+0.3+0.0}_{-0.0-0.3-0.0})\%,$ (44)
${\phi}_{{\parallel}}\,{\simeq}\,0,\qquad{\phi}_{{\perp}}\,{\simeq}\,{\pi},$
(45)

where the first uncertainty comes from the choice of the typical scale
$(1{\pm}0.1)t_{i}$, and the expression $t_{i}$ is given in Eq.(71) and
Eq.(72); the second uncertainty is from masses $m_{b}$ and $m_{c}$; the third
uncertainty is from hadronic parameters including decay constants and
Gegenbauer moments; and the fourth uncertainty of branching ratio comes from
the CKM parameters. The following are some comments.

(1) Branching ratio for the ${\Upsilon}(1S)$ ${\to}$ $B_{c}{\rho}$ decay with
the pQCD approach is different from previous estimation ijma14 ; adv2013  with
the NF approximation. Many factors lead to these differences. For example, as
it is showed in Ref. adv2013 , the values of form factors for ${\Upsilon}(1S)$
${\to}$ $B_{c}$ transition are very sensitive to the choice of wave functions.
In addition, form factors written as the convolution integral of wave
functions in Ref. adv2013  are usually enhanced by one-gluon-exchange
scattering amplitudes with the pQCD approach. These discrepancy deserve much
dedicated study and should be carefully tested by the future experiments.

(2) Branching ratio for the ${\Upsilon}(1S)$ ${\to}$ $B_{c}{\rho}$ decay can
reach up to ${\cal O}(10^{-9})$, which might be measurable at the running LHC
and forthcoming SuperKEKB. For example, the ${\Upsilon}(1S)$ production cross
section in p-Pb collision is about a few ${\mu}b$ at LHCb jhep1407  and ALICE
plb740 . Over $10^{12}$ ${\Upsilon}(1S)$ data samples per $ab^{-1}$ data
collected at LHCb and ALICE are in principle available, corresponding to a few
thousands of the ${\Upsilon}(1S)$ ${\to}$ $B_{c}{\rho}$ events.

(3) There is a hierarchical pattern among the longitudinal $f_{0}$, parallel
$f_{{\parallel}}$, and perpendicular $f_{{\perp}}$ polarization fractions,
i.e.,

$f_{0}:f_{{\parallel}}:f_{{\perp}}\ {\simeq}\
1:\frac{p}{\sqrt{2}m_{{\Upsilon}(1S)}}:\frac{p^{2}}{2m_{{\Upsilon}(1S)}^{2}},$
(46)

where $p$ is the common momentum of final state in the rest frame of the
${\Upsilon}(1S)$ meson. The relation Eq.(46) is basically agree with previous
estimation adv2013 . It means that the contributions to branching ratio for
the ${\Upsilon}(1S)$ ${\to}$ $B_{c}{\rho}$ decay mainly come from the
longitudinal polarization fractions, because of $f_{0}$ $>$ $f_{{\parallel}}$
$>$ $f_{{\perp}}$.

(4) The relative phase ${\phi}_{{\parallel}}$ is close to zero. The reason is
that the factorizable contributions from diagrams Fig.1(a,b) is real and
proportional to the large coefficient $a_{1}$, while the nonfactorizable
contributions from diagrams Fig.1(c,d) is suppressed by the color factor and
proportional to the small Wilson coefficient $C_{2}$, and the strong phases
arise only from the nonfactorizable contributions, which is consistent with
the prediction of the QCD factorization approach qcdf1 ; qcdf2  where the
strong phase arising from nonfactorizable contributions is suppressed by color
and ${\alpha}_{s}$ for the $a_{1}$-dominated processes. The relative phases,
if they could be determined experimentally, will improve our understanding on
the strong interactions.

##  IV Summary

The ${\Upsilon}(1S)$ weak decay is allowable within the standard model. In
this paper, the ${\Upsilon}(1S)$ ${\to}$ $B_{c}{\rho}$ weak decays are studied
with the pQCD approach. It is found that with the nonrelativistic wave
functions for ${\Upsilon}(1S)$ and $B_{c}$ mesons, the longitudinal
polarization fraction is the largest one, and branching ratios for the
${\Upsilon}(1S)$ ${\to}$ $B_{c}{\rho}$ decay can reach up to ${\cal
O}(10^{-9})$, which might be detectable at the future experiments.

## Acknowledgments

We thank Professor Dongsheng Du (IHEP@CAS), Professor Caidian Lü (IHEP@CAS)
and Professor Yadong Yang (CCNU) for helpful discussion. The work is supported
by the National Natural Science Foundation of China (Grant Nos. 11547014,
11475055, U1332103 and 11275057).

##  Appendix A Building blocks of decay amplitudes

For the sake of simplicity, the amplitude for the ${\Upsilon}(1S)$ ${\to}$
$B_{c}{\rho}$ decay, Eq.(33), is decomposed into building blocks ${\cal
A}^{i}_{L,N,T}$, where the superscript $i$ corresponds to the indices of
Fig.1. With the pQCD master formula Eq.(5), the explicit expressions of ${\cal
A}^{i}_{L,N,T}$ are written as follows:

$\displaystyle{\cal A}_{L}^{a}$ $\displaystyle=$
$\displaystyle{\int}_{0}^{1}dx_{1}{\int}_{0}^{1}dx_{2}{\int}_{0}^{\infty}b_{1}db_{1}{\int}_{0}^{\infty}b_{2}db_{2}\,{\phi}_{\Upsilon}^{v}(x_{1})$
(47)
$\displaystyle{\phi}_{B_{c}}(x_{2})\,E_{f}(t_{a})\,{\alpha}_{s}(t_{a})\,a_{1}(t_{a})\,H_{f}({\alpha}_{e},{\beta}_{a},b_{1},b_{2})$
$\displaystyle\Big{\\{}m_{1}^{2}\,s+m_{2}\,m_{b}\,u-(4\,m_{1}^{2}\,p^{2}+m_{2}^{2}\,u)\,\bar{x}_{2}\Big{\\}},$
$\displaystyle{\cal A}_{N}^{a}$ $\displaystyle=$ $\displaystyle
m_{1}\,m_{3}{\int}_{0}^{1}dx_{1}{\int}_{0}^{1}dx_{2}{\int}_{0}^{\infty}b_{1}db_{1}{\int}_{0}^{\infty}b_{2}db_{2}\,{\phi}_{\Upsilon}^{V}(x_{1})$
(48)
$\displaystyle{\phi}_{B_{c}}(x_{2})\,E_{f}(t_{a})\,{\alpha}_{s}(t_{a})\,a_{1}(t_{a})\,H_{f}({\alpha}_{e},{\beta}_{a},b_{1},b_{2})\,$
$\displaystyle\Big{\\{}2\,m_{2}^{2}\,\bar{x}_{2}-2\,m_{2}\,m_{b}-t\Big{\\}},$
$\displaystyle{\cal A}_{T}^{a}$ $\displaystyle=$ $\displaystyle
2\,m_{1}\,m_{3}\,{\int}_{0}^{1}dx_{1}{\int}_{0}^{1}dx_{2}{\int}_{0}^{\infty}b_{1}db_{1}{\int}_{0}^{\infty}b_{2}db_{2}\,{\phi}_{\Upsilon}^{V}(x_{1})$
(49)
$\displaystyle{\phi}_{B_{c}}(x_{2})\,E_{f}(t_{a})\,{\alpha}_{s}(t_{a})\,a_{1}(t_{a})\,H_{f}({\alpha}_{e},{\beta}_{a},b_{1},b_{2}),$
$\displaystyle{\cal A}_{L}^{b}$ $\displaystyle=$
$\displaystyle{\int}_{0}^{1}dx_{1}{\int}_{0}^{1}dx_{2}{\int}_{0}^{\infty}b_{1}db_{1}{\int}_{0}^{\infty}b_{2}db_{2}$
(50)
$\displaystyle{\phi}_{B_{c}}(x_{2})\,E_{f}(t_{b})\,{\alpha}_{s}(t_{b})\,a_{1}(t_{b})H_{f}({\alpha}_{e},{\beta}_{b},b_{2},b_{1})$
$\displaystyle\Big{\\{}{\phi}_{\Upsilon}^{v}(x_{1})\,\Big{[}m_{1}^{2}\,(s-4\,p^{2})\,\bar{x}_{1}+2\,m_{2}\,m_{c}\,u-m_{2}^{2}\,u\Big{]}$
$\displaystyle+{\phi}_{\Upsilon}^{t}(x_{1})\,m_{1}\,\Big{[}s\,(2\,m_{2}-m_{c})-2\,m_{2}\,u\,\bar{x}_{1}\Big{]}\Big{\\}},$
$\displaystyle{\cal A}_{N}^{b}$ $\displaystyle=$ $\displaystyle
m_{3}{\int}_{0}^{1}dx_{1}{\int}_{0}^{1}dx_{2}{\int}_{0}^{\infty}b_{1}db_{1}{\int}_{0}^{\infty}b_{2}db_{2}$
(51)
$\displaystyle{\phi}_{B_{c}}(x_{2})\,E_{f}(t_{b})\,{\alpha}_{s}(t_{b})\,a_{1}(t_{b})\,H_{f}({\alpha}_{e},{\beta}_{b},b_{2},b_{1})$
$\displaystyle\Big{\\{}{\phi}_{\Upsilon}^{V}(x_{1})\,m_{1}\,\Big{[}2\,m_{2}^{2}-4\,m_{2}\,m_{c}-t\,\bar{x}_{1}\Big{]}$
$\displaystyle+{\phi}_{\Upsilon}^{T}(x_{1})\,\Big{[}t\,(m_{c}-2\,m_{2})+4\,m_{1}^{2}\,m_{2}\,\bar{x}_{1}\Big{]}\Big{\\}},$
$\displaystyle{\cal A}_{T}^{b}$ $\displaystyle=$
$\displaystyle-2\,m_{3}\,{\int}_{0}^{1}dx_{1}{\int}_{0}^{1}dx_{2}{\int}_{0}^{\infty}b_{1}db_{1}{\int}_{0}^{\infty}b_{2}db_{2}$
(52)
$\displaystyle{\phi}_{B_{c}}(x_{2})\,E_{f}(t_{b})\,{\alpha}_{s}(t_{b})\,a_{1}(t_{b})\,H_{f}({\alpha}_{e},{\beta}_{b},b_{2},b_{1})$
$\displaystyle\Big{\\{}{\phi}_{\Upsilon}^{V}(x_{1})\,m_{1}\,\bar{x}_{1}+{\phi}_{\Upsilon}^{T}(x_{1})\,(m_{c}-2\,m_{2})\Big{\\}},$
$\displaystyle{\cal A}_{L}^{c}$ $\displaystyle=$
$\displaystyle\frac{1}{N_{c}}{\int}_{0}^{1}dx_{1}{\int}_{0}^{1}dx_{2}{\int}_{0}^{1}dx_{3}{\int}_{0}^{\infty}db_{1}{\int}_{0}^{\infty}b_{2}db_{2}{\int}_{0}^{\infty}b_{3}db_{3}$
(53)
$\displaystyle{\phi}_{B_{c}}(x_{2})\,{\phi}_{\rho}^{v}(x_{3})\,E_{n}(t_{c})\,{\alpha}_{s}(t_{c})\,C_{2}(t_{c})\,H_{n}({\alpha}_{e},{\beta}_{c},b_{2},b_{3})$
$\displaystyle{\delta}(b_{1}-b_{2})\,\Big{\\{}{\phi}_{\Upsilon}^{v}(x_{1})\,u\,\Big{[}t\,x_{1}-2\,m_{2}^{2}\,x_{2}-s\,\bar{x}_{3}\Big{]}$
$\displaystyle+{\phi}_{\Upsilon}^{t}(x_{1})\,m_{1}\,m_{2}\,\Big{[}s\,x_{2}+2\,m_{3}^{2}\,\bar{x}_{3}-u\,x_{1}\Big{]}\Big{\\}},$
$\displaystyle{\cal A}_{N}^{c}$ $\displaystyle=$
$\displaystyle\frac{m_{3}}{N_{c}}{\int}_{0}^{1}dx_{1}{\int}_{0}^{1}dx_{2}{\int}_{0}^{1}dx_{3}{\int}_{0}^{\infty}db_{1}{\int}_{0}^{\infty}b_{2}db_{2}{\int}_{0}^{\infty}b_{3}db_{3}$
(54)
$\displaystyle{\phi}_{B_{c}}(x_{2})\,E_{n}(t_{c})\,{\alpha}_{s}(t_{c})\,C_{2}(t_{c})\,H_{n}({\alpha}_{e},{\beta}_{c},b_{2},b_{3})$
$\displaystyle\Big{\\{}{\phi}_{\Upsilon}^{V}(x_{1})\,{\phi}_{\rho}^{V}(x_{3})\,m_{1}\,\Big{[}2\,s\,\bar{x}_{3}+4\,m_{2}^{2}\,x_{2}\,-2\,t\,x_{1}\Big{]}$
$\displaystyle+{\phi}_{\Upsilon}^{T}(x_{1})\,{\phi}_{\rho}^{V}(x_{3})\,m_{2}\,\Big{[}2\,m_{1}^{2}\,x_{1}-t\,x_{2}-u\,\bar{x}_{3}\Big{]}$
$\displaystyle+{\phi}_{\Upsilon}^{T}(x_{1})\,{\phi}_{\rho}^{A}(x_{3})\,2\,m_{1}\,m_{2}\,p\,(x_{2}-\bar{x}_{3})\Big{\\}}\,{\delta}(b_{1}-b_{2}),$
$\displaystyle{\cal A}_{T}^{c}$ $\displaystyle=$
$\displaystyle\frac{m_{3}}{N_{c}\,p}{\int}_{0}^{1}dx_{1}{\int}_{0}^{1}dx_{2}{\int}_{0}^{1}dx_{3}{\int}_{0}^{\infty}db_{1}{\int}_{0}^{\infty}b_{2}db_{2}{\int}_{0}^{\infty}b_{3}db_{3}\,$
(55)
$\displaystyle{\phi}_{B_{c}}(x_{2})\,E_{n}(t_{c})\,{\alpha}_{s}(t_{c})\,C_{2}(t_{c})\,H_{n}({\alpha}_{e},{\beta}_{c},b_{2},b_{3})$
$\displaystyle\Big{\\{}{\phi}_{\Upsilon}^{V}(x_{1})\,{\phi}_{\rho}^{A}(x_{3})\,\Big{[}2\,s\,\bar{x}_{3}+4\,m_{2}^{2}\,x_{2}\,-2\,t\,x_{1}\Big{]}$
$\displaystyle+{\phi}_{\Upsilon}^{T}(x_{1})\,{\phi}_{\rho}^{A}(x_{3})\,r_{2}\,\Big{[}2\,m_{1}^{2}\,x_{1}-t\,x_{2}-u\,\bar{x}_{3}\Big{]}$
$\displaystyle+2\,m_{2}\,p\,{\phi}_{\Upsilon}^{T}(x_{1})\,{\phi}_{\rho}^{V}(x_{3})\,(x_{2}-\bar{x}_{3})\Big{\\}}\,{\delta}(b_{1}-b_{2}),$
$\displaystyle{\cal A}_{L}^{d}$ $\displaystyle=$
$\displaystyle\frac{1}{N_{c}}{\int}_{0}^{1}dx_{1}{\int}_{0}^{1}dx_{2}{\int}_{0}^{1}dx_{3}{\int}_{0}^{\infty}db_{1}{\int}_{0}^{\infty}b_{2}db_{2}{\int}_{0}^{\infty}b_{3}db_{3}$
(56)
$\displaystyle{\phi}_{B_{c}}(x_{2})\,{\phi}_{\rho}^{v}(x_{3})\,E_{n}(t_{d})\,{\alpha}_{s}(t_{d})\,C_{2}(t_{d})\,H_{n}({\alpha}_{e},{\beta}_{d},b_{2},b_{3})$
$\displaystyle{\delta}(b_{1}-b_{2})\,\Big{\\{}{\phi}_{\Upsilon}^{t}(x_{1})\,m_{1}\,m_{2}\,\Big{[}s\,x_{2}+2\,m_{3}^{2}\,x_{3}-u\,x_{1}\Big{]}$
$\displaystyle+{\phi}_{\Upsilon}^{v}(x_{1})\,4\,m_{1}^{2}\,p^{2}\,(x_{3}-x_{2})\Big{\\}},$
$\displaystyle{\cal A}_{N}^{d}$ $\displaystyle=$
$\displaystyle\frac{m_{3}}{N_{c}}{\int}_{0}^{1}dx_{1}{\int}_{0}^{1}dx_{2}{\int}_{0}^{1}dx_{3}{\int}_{0}^{\infty}db_{1}{\int}_{0}^{\infty}b_{2}db_{2}{\int}_{0}^{\infty}b_{3}db_{3}$
(57)
$\displaystyle{\phi}_{\Upsilon}^{T}(x_{1})\,{\phi}_{B_{c}}(x_{2})\,E_{n}(t_{d})\,{\alpha}_{s}(t_{d})\,C_{2}(t_{d})\,H_{n}({\alpha}_{e},{\beta}_{d},b_{2},b_{3})$
$\displaystyle{\delta}(b_{1}-b_{2})\,\Big{\\{}{\phi}_{\rho}^{V}(x_{3})\,m_{2}\,\Big{[}2\,m_{1}^{2}\,x_{1}-t\,x_{2}\,-u\,x_{3}\Big{]}$
$\displaystyle+2\,m_{1}\,m_{2}\,p\,{\phi}_{\rho}^{A}(x_{3})\,(x_{2}-x_{3})\Big{\\}},$
$\displaystyle{\cal A}_{T}^{d}$ $\displaystyle=$
$\displaystyle\frac{m_{3}}{N_{c}\,p}{\int}_{0}^{1}dx_{1}{\int}_{0}^{1}dx_{2}{\int}_{0}^{1}dx_{3}{\int}_{0}^{\infty}db_{1}{\int}_{0}^{\infty}b_{2}db_{2}{\int}_{0}^{\infty}b_{3}db_{3}$
(58)
$\displaystyle{\phi}_{\Upsilon}^{T}(x_{1})\,{\phi}_{B_{c}}(x_{2})\,E_{n}(t_{d})\,{\alpha}_{s}(t_{d})\,C_{2}(t_{d})\,H_{n}({\alpha}_{e},{\beta}_{d},b_{2},b_{3})$
$\displaystyle{\delta}(b_{1}-b_{2})\,\Big{\\{}{\phi}_{\rho}^{A}(x_{3})\,r_{2}\,\Big{[}2\,m_{1}^{2}\,x_{1}-t\,x_{2}-u\,x_{3}\Big{]}$
$\displaystyle+2\,m_{2}\,p\,{\phi}_{\rho}^{V}(x_{3})\,(x_{2}-x_{3})\Big{\\}},$

where $\bar{x}_{i}$ $=$ $1$ $-$ $x_{i}$; variable $x_{i}$ and $b_{i}$ are the
longitudinal momentum fraction and the conjugate variable of the transverse
momentum $k_{i{\perp}}$ of the valence quark, respectively; ${\alpha}_{s}$ is
the QCD coupling; $a_{1}$ $=$ $C_{1}$ $+$ $C_{2}/N$; $C_{1,2}$ are the Wilson
coefficients.

The function $H_{f,n}$ and Sudakov factor $E_{f,n}$ are defined as follows,
where the subscripts $f$ and $n$ correspond to factorizable and
nonfactorizable topologies, respectively.

$H_{f}({\alpha}_{e},{\beta},b_{i},b_{j})\ =\
K_{0}(\sqrt{-{\alpha}_{e}}b_{i})\Big{\\{}{\theta}(b_{i}-b_{j})K_{0}(\sqrt{-{\beta}}b_{i})I_{0}(\sqrt{-{\beta}}b_{j})+(b_{i}{\leftrightarrow}b_{j})\Big{\\}},$
(59) $\displaystyle H_{n}({\alpha}_{e},{\beta},b_{2},b_{3})$ $\displaystyle=$
$\displaystyle\Big{\\{}{\theta}(-{\beta})K_{0}(\sqrt{-{\beta}}b_{3})+\frac{{\pi}}{2}{\theta}({\beta})\Big{[}iJ_{0}(\sqrt{{\beta}}b_{3})-Y_{0}(\sqrt{{\beta}}b_{3})\Big{]}\Big{\\}}$
(60) $\displaystyle{\times}$
$\displaystyle\Big{\\{}{\theta}(b_{2}-b_{3})K_{0}(\sqrt{-{\alpha}_{e}}b_{2})I_{0}(\sqrt{-{\alpha}_{e}}b_{3})+(b_{2}{\leftrightarrow}b_{3})\Big{\\}},$
$E_{f}(w)\ =\ {\exp}\\{-S_{\Upsilon}(w)-S_{B_{c}}(w)\\},$ (61) $E_{n}(w)\ =\
{\exp}\\{-S_{\Upsilon}(w)-S_{B_{c}}(w)-S_{\rho}(w)\\},$ (62) $S_{\Upsilon}(w)\
=\
s(x_{1},p_{1}^{+},1/b_{1})+2{\int}_{1/b_{1}}^{w}\frac{d{\mu}}{\mu}{\gamma}_{q},$
(63) $S_{B_{c}}(w)\ =\
s(x_{2},p_{2}^{+},1/b_{2})+2{\int}_{1/b_{2}}^{w}\frac{d{\mu}}{\mu}{\gamma}_{q},$
(64) $S_{\rho}(w)\ =\
s(x_{3},p_{3}^{+},1/b_{3})+s(\bar{x}_{3},p_{3}^{+},1/b_{3})+2{\int}_{1/b_{3}}^{w}\frac{d{\mu}}{\mu}{\gamma}_{q},$
(65)

where $J_{0}$ and $Y_{0}$ ($I_{0}$ and $K_{0}$) are the (modified) Bessel
function of the first and second kind, respectively; ${\gamma}_{q}$ $=$
$-{\alpha}_{s}/{\pi}$ is the quark anomalous dimension; the expression of
$s(x,Q,1/b)$ can be found in the appendix of Ref.pqcd1 ; ${\alpha}_{e}$ is the
gluon virtuality; the subscript of the quark virtuality ${\beta}_{i}$
corresponds to the indices of Fig.1. The definitions of the particle
virtuality and typical scale $t_{i}$ are listed as follows:

$\displaystyle{\alpha}_{e}$ $\displaystyle=$
$\displaystyle\bar{x}_{1}^{2}m_{1}^{2}+\bar{x}_{2}^{2}m_{2}^{2}-\bar{x}_{1}\bar{x}_{2}t,$
(66) $\displaystyle{\beta}_{a}$ $\displaystyle=$ $\displaystyle
m_{1}^{2}-m_{b}^{2}+\bar{x}_{2}^{2}m_{2}^{2}-\bar{x}_{2}t,$ (67)
$\displaystyle{\beta}_{b}$ $\displaystyle=$ $\displaystyle
m_{2}^{2}-m_{c}^{2}+\bar{x}_{1}^{2}m_{1}^{2}-\bar{x}_{1}t,$ (68)
$\displaystyle{\beta}_{c}$ $\displaystyle=$ $\displaystyle
x_{1}^{2}m_{1}^{2}+x_{2}^{2}m_{2}^{2}+\bar{x}_{3}^{2}m_{3}^{2}$ (69)
$\displaystyle-$ $\displaystyle
x_{1}x_{2}t-x_{1}\bar{x}_{3}u+x_{2}\bar{x}_{3}s,$ $\displaystyle{\beta}_{d}$
$\displaystyle=$ $\displaystyle
x_{1}^{2}m_{1}^{2}+x_{2}^{2}m_{2}^{2}+x_{3}^{2}m_{3}^{2}$ (70)
$\displaystyle-$ $\displaystyle x_{1}x_{2}t-x_{1}x_{3}u+x_{2}x_{3}s,$
$\displaystyle t_{a(b)}$ $\displaystyle=$
$\displaystyle{\max}(\sqrt{-{\alpha}_{e}},\sqrt{-{\beta}_{a(b)}},1/b_{1},1/b_{2}),$
(71) $\displaystyle t_{c(d)}$ $\displaystyle=$
$\displaystyle{\max}(\sqrt{-{\alpha}_{e}},\sqrt{{|}{\beta}_{c(d)}{|}},1/b_{2},1/b_{3}).$
(72)

## References

  * (1) K. Olive et al. (Particle Data Group), Chin. Phys. C 38, 090001 (2014).
  * (2) Ed. A. Bevan et al., Eur. Phys. J. C 74, 3026 (2014).
  * (3) H. Li, Phys. Rev. D 52, 3958 (1995).
  * (4) C. Chang, H. Li, Phys. Rev. D 55, 5577 (1997).
  * (5) T. Yeh, H. Li, Phys. Rev. D 56, 1615 (1997).
  * (6) K. K. Sharma and R. C. Verma, Int. J. Mod. Phys. A 14, 937 (1999).
  * (7) R. Dhir, R. C. Verma and A. Sharma, Adv. High Energy Phys. 2013, 706543 (2013).
  * (8) M. Beneke et al., Phys. Rev. Lett. 83, 1914 (1999).
  * (9) M. Beneke et al., Nucl. Phys. B 591, 313 (2000).
  * (10) M. Beneke et al., Nucl. Phys. B 606, 245 (2001).
  * (11) C. Bauer et al., Phys. Rev. D 63, 114020 (2001).
  * (12) C. Bauer, D. Pirjol, I. Stewart, Phys. Rev. D 65, 054022 (2002).
  * (13) C. Bauer et al., Phys. Rev. D 66, 014017 (2002).
  * (14) M. Beneke et al., Nucl. Phys. B 643, 431 (2002).
  * (15) N. Brambilla et al. (QWG group), arXiv:hep-ph/0412158.
  * (16) C. Chang and X. Wu, Eur. Phys. J. C 38, 267 (2004).
  * (17) G. Buchalla, A. Buras, M. Lautenbacher, Rev. Mod. Phys. 68, 1125, (1996).
  * (18) G. Lepage, S. Brodsky, Phys. Rev. D 22, 2157 (1980).
  * (19) T. Kurimoto, H. Li, A. Sanda, Phys. Rev. D 65, 014007 (2001).
  * (20) P. Ball, G. Jones, JHEP, 0703, 069, (2007).
  * (21) G. Lepage et al., Phys. Rev. D 46, 4052 (1992).
  * (22) G. Bodwin, E. Braaten, G. Lepage, Phys. Rev. D 51, 1125 (1995).
  * (23) N. Brambilla et al., Rev. Mod. Phys. 77, 1423 (2005).
  * (24) S. J. Brodsky, T. Huang and G. P. Lepage, in Particles and Fields-2, Proceedings of the Banff Summer Institute, Banff, Alberta, 1981, edited by A. Z. Capri and A. N. Kamal (Plenum, New York, 1983), P143; T. Huang, in Proceedings of XXth International Conference on High Energy Physics, Madison, Wisconsin, 1980, edited by L. Durand and L. G. Pondrom, AIP Conf. Proc. No. 69 (AIP, New York, 1981), p1000.
  * (25) B. Xiao, X. Qin, B. Ma, Eur. Phys. J. A 15, 523 (2002).
  * (26) T. Huang, X. Wu and X. Wu, Phys. Rev. D 70, 053007 (2004); T. Huang, T. Zhong and X. Wu, Phys. Rev. D 88, 034013 (2013); T. Huang and F. Zuo, Eur. Phys. J. C 51, 833 (2007); X. Wu and T. Huang, Phys. Rev. D 84, 074011 (2011); Y. Sun et al., Eur. Phys. J. C 67, 117 (2010).
  * (27) J. Sun et al., Int. J. Mod. Phys. A 31, 1650061 (2016).
  * (28) J. Sun et al., Phys. Rev. D 92, 074028 (2015).
  * (29) C. Chen, Y. Keum, H. Li, Phys. Rev. D 66, 054013 (2002)
  * (30) H. Fu et al., Phys. Lett. B 738, 228 (2014); J. Phys. G 42, 055002 (2015).
  * (31) T. Chiu, T. Hsieh, C. Huang, K. Ogawa, Phys. Lett. B 651, 171 (2007).
  * (32) R. Aaij et al. (LHCb Collaboration), JHEP 1407, 094 (2014).
  * (33) B. Abelev et al. (ALICE Collaboration), Phys. Lett. B 740, 105 (2015).