# Effect of Contact Angle on Capillary Pressure and Liquid Recovery from
Angular Pore Channels

Afshin Davarpanah111Correspondence to<EMAIL_ADDRESS>Simon Cox,
Department of Mathematics, Aberystwyth University,
Aberystwyth SY23 3BZ, United Kingdom

###### Abstract

We predict the volume of liquid recovered from different-shaped prismatic
channels following gas displacement. This recovery factor depends strongly on
the contact angle at which the gas-liquid interfaces meet the walls of the
channel. We show that scaling the contact angle by a critical value, dependent
on the channel shape, leads to a master curve for the recovery factor. Our
numerical simulations with the Surface Evolver are supported by theoretical
calculations for the capillary entry pressure in a channel with a triangular
or rectangular cross-section and for the interface shape in a “scalloped”
triangular channel formed between three parallel cylinders.

Keywords: Interface shape; Capillary pressure; Porous media; Oil Recovery;
Contact angle.

##  1 Introduction

The phenomenon of immiscible displacement is a crucial factor in enhanced oil
recovery from porous media. In pore-scale networks, fluid interfaces generate
capillary forces that can trap hydrocarbons and hence reduce the recovery
factor [1, 2]. The significant factors in determining capillary behaviour in
porous media are the pore structure itself and the wettability.

Wettability is defined as a fluid’s relative tendency to adhere to the solid
phase. It is measured by the contact angle at which a fluid interface meets a
solid surface at mechanical equilibrium. For contact angles below
$\frac{\pi}{2}$, the surface is known as wetting, while for contact angle
greater than $\frac{\pi}{2}$, it is called non-wetting.

The contact angle is directly related to the balance between the surface
tension $\gamma$ of the liquid-gas interface and the difference in surface
tension, which we denote $\gamma_{w}$, between the wetting film on the solid
walls and the liquid-solid interface inside the liquid. If we denote the
contact angle by $\theta$, then

$\cos\theta=\frac{\gamma_{w}}{\gamma}$ (1)

When the wall tension $\gamma_{w}$ is equal to $\gamma$ the contact angle is
zero, often referred to as the case of perfect wetting. As $\gamma_{w}$
decreases the contact angle increases. We assume here that the contact angle
is fixed, and does not depend on the direction of motion of the contact line
(i.e. there is no hysteresis between imbibition and drainage). Some authors
[3, 4] include a disjoining pressure in the thin wetting film; we neglect such
a contribution here.

The usual idealisations of a porous medium consists of packings of spheres [5]
or cylinders [6]. To ensure that modeling is capturing processes in real
porous materials, it is also important to consider angular pore structures. In
channels with a polygonal cross-section, for example, liquid preferentially
collects in the sharp corners, and it is important to determine the shape and
size of the corner meniscus as the contact angle varies [7, 11]. The capillary
pressure ($P_{c}$) is directly proportional to the mean curvature $C$ of the
meniscus, as determined by the Laplace-Young equation [12], $P_{c}=\gamma C$.

Ma et al. [7] investigated the relationship between liquid saturation, corner
angle, and contact angle on the meniscus curvature within straight channels
with regular polygonal cross-sections. They determined the meniscus curvature
for various $n-$sided channels (for example $n=3$ for an equilateral triangle
and $n=4$ for a square) with different contact angles. In terms of the half-
corner angle $\alpha=\pi(n-2)/2n$ (so $\alpha=\frac{\pi}{6}$ for a channel
with an equilateral triangular cross-section and $\alpha=\frac{\pi}{4}$ for a
channel with a square cross-section), they give the following expression for
the capillary pressure when the fluid in each of the $n$ corners occupies an
area $A$ of a channel with cross-sectional area $A_{T}$:

$P_{c}(\theta,\alpha,A)=\frac{\gamma}{R}\sqrt{\frac{\tan\alpha\;A_{T}\;F(\theta,\alpha)}{nA}},$
(2)

where $R=\frac{1}{2}W\tan\alpha$ is the radius of the in-circle of the cross-
section and the geometric function $F$ is

$F(\theta,\alpha)=\frac{\cos\theta\cos(\theta+\alpha)}{\sin\alpha}-\left(\frac{\pi}{2}-(\alpha+\theta)\right).$
(3)

For example, in a channel with a cross-section in the form of an equilateral
triangle (the case $n=3$, $\alpha=\frac{\pi}{6}$ shown in Figure 1(a)) with
side length $W$, $R=\sqrt{3}W/6$ and area $\sqrt{3}W^{2}/4$, the capillary
pressure is

$P_{c}(\theta,\frac{\pi}{6},A)=\gamma\sqrt{\frac{F(\theta,\frac{\pi}{6})}{A}}.$
(4)

As the liquid area increases the capillary pressure drops and, since
$F(\theta,\alpha)$ is a decreasing function of $\theta$ in the interval
$0\leq\theta\leq\frac{\pi}{3}$, the capillary pressure also drops as the
contact angle increases. When $\theta=\frac{\pi}{3}$
$F(\frac{\pi}{3},\frac{\pi}{6})=0$, the capillary pressure is zero and the
interface becomes planar.

We consider a straight 3D channel, to represent a pore body, with a fixed
cross-section that is initially filled with wetting fluid (“liquid”). A
second, non-wetting fluid (“gas”), is gradually introduced at one end of the
channel. The gas must overcome the capillary entry pressure (or threshold
pressure), denoted $P_{c}^{e}$, before the wetting phase is displaced from the
other end of the channel to be recovered. The interface between the phases
then moves downstream with a roughly hemispherical form, leaving behind thin
triangular regions of liquid in the corners, as shown in Figure 1(d)-(f) and
in the animation (in the case of a channel with rectangular cross-section) in
the Supplementary Information.

It is therefore necessary to determine the capillary entry pressure. For a
channel with the same regular polygonal cross-section as above, Ma et al. [7]
give

$P_{c}^{e}=\frac{\gamma}{R}\left(\cos\theta+\sqrt{\tan\alpha\left(\cos\theta\sin\theta+\frac{\pi}{2}-(\alpha+\theta)\right)}\right).$
(5)

The radius of the in-circle $R$ depends on $n$, via $\alpha$, as demonstrated
above. Taking the example of an equilateral triangle again, we have

$P_{c}^{e}=2\sqrt{3}\frac{\gamma}{W}\left(\cos\theta+\sqrt{\frac{1}{\sqrt{3}}\left(\cos\theta\sin\theta+\frac{\pi}{3}-\theta\right)}\right).$
(6)

Expressions for $P_{c}$ for channels with various shapes can be found in [8,
9, 10].

We can then determine the area of liquid remaining in a channel after the
drying front has moved far downstream by evaluating the shape of the liquid
regions when the capillary pressure of the liquid is equal to the capillary
entry pressure. In the case of an equilateral triangular cross-section, for
any contact angle, we equate eq. (6) with eq. (4) to find $A(\theta)$.

We describe the drainage of channels with non-uniform or non-polygonal cross-
sections, such as a rectangle or the gap between three cylinders, with further
analysis and three-dimensional (3D) numerical simulations. We first describe
the geometry of the channels under consideration in § 2, and outline our
numerical method. In § 3.1 we validate our simulations by comparing with
predictions for the capillary pressure, which requires that we derive an
expression for the capillary pressure in the gap between three cylinders. In §
3.2 we determine the maximum possible liquid volume in each channel before
snap-off occurs and the channel rapidly fills with liquid. And then in § 3.3
we calculate the capillary entry pressure in these non-polygonal channels,
which allows us to determine the liquid recovery factor as a function of
contact angle.

##  2 Channel Geometry and Numerical Method

We consider channels with a cross-section in one of three shapes shown in
Figure 1:

  * •

an equilateral triangular with side-length $W$; this shape has no free
parameters and allows us to validate our numerical method against known
results such as eq. (4).

  * •

a rectangle with width $W$ and height $H$; we can vary the aspect ratio $H/W$
to determine how the recovery factor depends on this ratio.

  * •

the gap between three touching, parallel, circular cylinders [10] of radius
$W$; we refer to this shape as a scalloped triangle, also considered
previously in the context of capillary rise [13], which provides a stringent
test for both theory and numerics because of the cusps in the corners of the
channel.

To determine the shape and capillary pressure of the liquid interfaces in the
different channels, we perform simulations with the Surface Evolver [14], a
finite element software designed to predict the static shapes of fluid
interfaces. By directly determining interface shape, without the need for a
discretisation of space, Evolver offers a fast and efficient alternative to
numerical methods such as level sets [15] and 3D finite element fluid flow
calculations with interface tracking [11].

Surface Evolver provides a surface discretation (into small triangles) and
routines for minimizing energy; in this case the energy is the surface area of
the interfaces multiplied by their surface tension, and the interface shapes
are constrained by the enclosed volumes. The contact angle can be accurately
described and quantities associated with the surfaces, such as area and
capillary pressure, determined to arbitrary accuracy by repeatedly refining
the triangulation into smaller triangles and performing further gradient
descent iterations. By slowly incrementing the values of the constraints, such
as the enclosed volume, we are able to perform quasi-static simulations of
drainage and imbibition, for example.

For 2D calculations in the cross-section of a channel there is a single
movable interface in each corner of the channel; for a channel size of $W=1$
we refine this interface into short straight segments using three levels of
refinement. In 3D we commence our calculations with the interface roughly
half-way along the channel and increase or decrease the gas volume to explore
how the capillary pressure varies during steady-state drainage. Without loss
of generality we set the interfacial tension to $\gamma=1$ and choose a value
of the wall tension $\gamma_{w}$ to set the desired contact angle via eq. (1).
In each case we converge using gradient descent and Hessian steps.

The capillary pressures that we report are in units of $\gamma/W$, i.e. the
Laplace pressure. For a dimensional tension $\gamma=30\times 10^{-3}{\rm N/m}$
and pore scale $W=1{\rm\mu m}$, we have that $P_{c}=1$ here corresponds to a
capillary pressure of $30\times 10^{3}{\rm Pa}$. Liquid areas are scaled by
the total cross-sectional area of the channel, and liquid volumes by the total
channel volume.

Figure 1: Channel shapes considered here, with contact angle $10^{\circ}$ in
each case. (a) Channel with a cross-section in the form of an equilateral
triangle, showing the width $W$, contact angle $\theta$ and three regions of
liquid in the corners. The height $H=\sqrt{3}W/2$, and the corner half-angle
$\alpha=\pi/6$. (b) Rectangular channel with width $W$ and height $H$. (c)
Scalloped triangle between three cylinders of radius $W$. (d) 3D channel
corresponding to the cross-section in (a), with length $L$. The interface
separating gas from liquid has penetrated almost half-way along the channel,
moving from left to right, with capillary pressure equal to the capillary
entry pressure. The wetting films on the channel walls behind the moving front
are not shown. (e) Equivalent image for a rectangular channel with $H=2W$. (f)
Equivalent image for the channel with scalloped triangular cross-section.

##  3 Results

###  3.1 Dependence of capillary pressure on contact angle and liquid volume

Figure 2 shows how the capillary pressure varies in the different channels
considered here. As the liquid areas increase, the capillary pressures
decrease.

Increasing the contact angle also causes the capillary pressure to decrease,
leading to a reduction in the amount of liquid left in the channel. At a
critical contact angle $\theta_{0}$, which depends on the shape of the
channel, the capillary pressure reaches zero. For a channel with a regular
polygonal cross-section the critical angle is
$\theta_{0}=\frac{\pi}{2}-\alpha=\frac{\pi}{n}$, so for a triangular channel
(Figure 2(a)) we have $\theta_{0}=\frac{\pi}{3}$. A rectangular channel takes
the same value as a square channel, $\theta_{0}=\frac{\pi}{4}$ [11], while for
a scalloped triangular channel, which is not polygonal, it is
$\theta_{0}=\frac{\pi}{2}$.

(a)

(b)

(c)

Figure 2: Capillary pressure versus total liquid area in a 2D cross-section
through different channels. Large discs denote the capillary entry pressure;
small squares indicate the maximum liquid area at which two (or more) liquid
regions merge, from eqs. (10) and (11). (a) An equilateral triangular channel,
showing the capillary entry pressure from eq. (5) superimposed on the lines
given by eq. (2), for various contact angles. Simulation results are shown as
points, validating the simulation technique. The liquid area is scaled by the
area of the triangle, $A_{T}=\sqrt{3}W^{2}/4$. (b) Channels with rectangular
cross-section, for contact angle $\theta=\frac{\pi}{18}$. Eq. (2) with $n=4$
and $\alpha=\frac{\pi}{4}$ gives the lines, the capillary entry pressure is
determined by eq. (15), and the liquid areas are scaled by $A_{T}=HW$. (c)
Channels with scalloped triangle cross-section for various contact angles,
showing the capillary pressure from eqs. (7) and (8) as lines, compared with
the result of numerical simulation in a 2D cross-section in the case
$\theta=\pi/3$, and the capillary entry pressure from simulations in a 3D
channel. The total channel area is
$A_{T}=\left(\sqrt{3}-\frac{\pi}{2}\right)W^{2}$.

Compared to the triangular channel, the pressures are lower for given liquid
area in a square channel (Figure 2(a) and (b)) because the curvature of the
interface required to meet the walls at the required contact angle is less
when the angle between the walls is larger. In the scalloped triangular
channel (Figure 2(c)) the capillary pressures are therefore even higher.

The capillary pressure of the small regions in the corners of the rectangular
channel are identical to those that are found in a square channel up until the
point at which they meet (i.e. with capillary pressure given by eq. (2) with
$n=4$, $\alpha=\pi/4$). For channels with different aspect ratios but the same
total area ($HW=1$), Figure 2 shows that the data collapse onto one curve,
with an increase in the capillary entry pressure as the aspect ratio $H/W$
increases. However, in the rectangular channels with higher aspect ratio the
liquid regions meet when their area is smaller. Increasing the channel size
(the figure shows the case $HW=2$) reduces the capillary pressures, although
the relative liquid area at which the capillary pressure occurs depends only
on the channel aspect ratio.

(a)

(b)

(c)

Figure 3: (a) From left to right, increasing liquid areas in the scalloped
triangular channel, with $3A/A_{T}=0.2,0.3,0.4,0.5,0.58$. The right-most image
corresponds to the point at which the liquid regions merge. (b) Geometry of
the liquid in one corner of a cross-section through the scalloped triangular
channel. The two arcs enclosing the liquid meet at the contact angle $\theta$,
and this point on the channel wall subtends an angle $\Psi$ at the centre of
the cylinder. The inner interface has radius of curvature $r$, which sets the
capillary pressure $P_{c}$. (c) Solution of eq. (16) for the value of $\Psi$
at which the capillary entry pressure is reached as the contact angle varies.

For the scalloped triangular channel, we determine the dependence of the
liquid area and capillary pressure on the contact angle in terms of a
parameter $\Psi$. This is the angle subtended by the point at which the
meniscus meets the wall of the channel, measured from the corner of the
triangle (see figure 3(b)), and therefore taking values in the range
$[0:\pi/6]$. We find

$A(\theta,\Psi)=W^{2}\left[\tan\Psi-\Psi+\left(\frac{1-\cos\Psi}{\cos(\theta+\Psi)}\right)^{2}\left(\frac{\sin\theta\cos(\theta+\Psi)}{\cos\Psi}+\theta+\Psi-\frac{\pi}{2}\right)\right]$
(7)

and

$P_{c}(\theta,\Psi)=\frac{\gamma}{W}\frac{\cos(\theta+\Psi)}{1-\cos(\Psi)}.$
(8)

It is then possible to invert eq. (8) to give $\Psi$ in terms of $p_{c}$, and
hence write $A$ in terms of $p_{c}$. The result is shown in figure 2, and this
appears to agree with the calculations in [10] (Table 3).

###  3.2 Merging of liquid regions

We now determine the point at which the capillary pressure curves in Figure 2
stop, often well before the point at which liquid fills the channel
($nA=A_{T}$), because two (or more) liquid regions have come into contact, or
fused [7]. Figure 3(a) shows the point just before this happens, in the case
of the scalloped triangle. In other words, we seek the point at which the
wetted length on the channel wall is equal to half the length of the wall.

For a polygonal channel, the wetted length can be written in terms of each
liquid area $A$ or the capillary pressure [7, 16] as

$L=\frac{\cos(\theta+\alpha)}{\sin\alpha}\sqrt{\frac{A}{F(\theta,\alpha)}}=\frac{\cos(\theta+\alpha)}{\sin\alpha}\frac{\gamma}{P_{c}},$
(9)

where $F$ is defined in eq. (3) and $A$ is the area of liquid in each corner,
as before. When $L$ becomes equal to $W/2$, we find the critical capillary
pressure and liquid areas:

$P_{c}^{\rm
crit}=\frac{2\gamma}{W}\frac{\cos(\theta+\alpha)}{\sin\alpha},\qquad\quad
A^{\rm crit}=F(\theta,\alpha)\left(\frac{\gamma}{P_{c}^{\rm
crit}}\right)^{2}=\frac{W^{2}}{4}\frac{\sin\alpha\;F(\theta,\alpha)}{\cos(\theta+\alpha)}.$
(10)

These are marked in Figure 2(a) for the equilateral triangle
($W=1,\alpha=\pi/6)$. For the rectangular channel, merging occurs first on the
shortest side. Without loss of generality we can take this to be the width
$W$, and then eq. (10) applies – see Figure 2(b).

For the scalloped triangle, eq. (10) does not apply. Instead we determine the
critical capillary pressure and liquid area from eqs. (7) and (8) with
$\Psi=\pi/6$:

$P_{c}^{\rm
crit}=\frac{2\gamma}{W}\frac{\cos\left(\theta+\frac{\pi}{6}\right)}{(2-\sqrt{3})},\qquad\quad
A^{\rm
crit}=W^{2}\left[\frac{1}{\sqrt{3}}-\frac{\pi}{6}+\frac{(7-4\sqrt{3})}{4\cos^{2}\left(\theta+\frac{\pi}{6}\right)}\left(\frac{2}{\sqrt{3}}\sin\theta\cos\left(\theta+\frac{\pi}{6}\right)+\theta-\frac{\pi}{3}\right)\right].$
(11)

This is shown in Figure 2(c), again in good agreement with the data.

Following the merging of two liquid regions, the capillary pressure would
increase as the total liquid area increased. This is an unstable situation,
leading to what is known as snap-off [17] and the channel filling with liquid.
We do not consider it further here, but it is important to determine the
liquid volumes at which it occurs.

###  3.3 Capillary Entry Pressure and Liquid Recovery Factor

The data of Figure 2 can be summarised by calculating the liquid recovery
factor (RF), that is, the relative volume of liquid (wetting phase) removed
from the channel:

$RF=1-\frac{nAL}{A_{T}L}=1-\frac{nA}{A_{T}},$ (12)

where $A$ is measured at the capillary entry pressure, $P_{c}^{e}$,
corresponding to the area of liquid remaining in the triangular regions of
liquid in the corners of the channel after the roughly hemispherical front has
passed through the channel (Figure 2).

For the triangular channel we have

$RF_{\rm tri}=2-\frac{4\sqrt{3}}{W}\frac{\gamma\cos\theta}{P_{c}^{e}},$ (13)

where $P_{c}^{e}$ is given by eq. (6).

For a rectangular channel we have

$RF_{\rm
rect}=2-\left(\frac{1}{W}+\frac{1}{H}\right)\frac{2\gamma\cos\theta}{P_{c}^{e}},$
(14)

with the capillary entry pressure (from [9] (eq. 51), re-cast in our notation)
given by:

$P_{c}^{e}=\gamma\left(\left(\frac{1}{W}+\frac{1}{H}\right)\cos\theta+\sqrt{\left(\frac{1}{W}+\frac{1}{H}\right)^{2}\cos^{2}\theta-\frac{4}{WH}\left[\cos^{2}\theta-\cos\theta\sin\theta-\left(\frac{\pi}{4}-\theta\right)\right]}\right).$
(15)

This is derived using what Mason and Morrow [18] call the “MS-P” method, after
Mayer and Stowe [19] and Princen [13], and improves upon the formula given by
Lenormand et al. [8] for a rectangle with zero contact angle, which omits the
term proportional to $4/WH$ (and sets $\cos\theta=1$).

For the scalloped triangle, the expressions required to derive the dependence
of the capillary entry pressure on the radius $W$ and contact angle $\theta$
are given by Mayer and Stowe [19] and, in the slightly more general case where
the circles need not touch, by Princen [13]. These formulae do not admit a
closed form solution, but require numerical solution. In our notation, the
MS-P method leads to the following nonlinear equation for the critical value
of the parameter $\Psi$ (which we denote $\hat{\Psi}$, a function of the
contact angle $\theta$) at which the capillary entry pressure occurs:

$\left(\frac{1}{\sqrt{3}}-\frac{\pi}{6}+\hat{\Psi}-\tan\hat{\Psi}\right)=\frac{(1-\cos\hat{\Psi})^{2}}{\cos^{2}(\theta+\hat{\Psi})}\left(\frac{\sin\theta\cos(\theta+\hat{\Psi})}{\cos\hat{\Psi}}+\frac{\pi}{2}-\theta-\hat{\Psi}\right)+\frac{2(1-\cos\hat{\Psi})\cos\theta}{\cos(\theta+\hat{\Psi})}\left(\frac{\pi}{6}-\hat{\Psi}\right),$
(16)

and then the capillary entry pressure is given by evaluating eq. (8) with
$\Psi=\hat{\Psi}$, i.e. $P_{c}^{e}=P_{c}(\theta,\hat{\Psi})$. We solve eq.
(16) using a numerical root-finding algorithm, the result of which is shown in
Figure 3(c), to give the values of $P_{c}^{e}$ shown in figure 2(c). The
recovery factor is then

$RF_{\rm sc-
tri}=\left(\frac{6}{\sqrt{3}-\frac{1}{2}\pi}\right)\frac{\gamma}{WP_{c}^{e}}\left[\frac{\gamma}{WP_{c}^{e}}\left(\frac{\pi}{2}-\theta-\Psi\right)+\cos\theta\left(\frac{\pi}{6}-\Psi\right)\right].$
(17)

Figure 4 shows that the recovery factor depends strongly on the contact angle
for all three channel shapes, reaching a value of one (i.e. complete emptying
of liquid from the channel) when the critical contact angle $\theta_{0}$ is
reached. Our analytic theory is in excellent agreement with the simulations,
although in the case of the scalloped triangle we need a straightforward
numerical solution of eq. (16), or reference to Figure 3(c), to completely
determine $RF$ for given contact angle $\theta$.

The volume of liquid remaining in the channel in the case of zero contact
angle is greatest for the channel with the scalloped triangular section and
smallest for a rectangular channel with large aspect ratio. Scaling the
recovery factor by its value for zero contact angle, as shown in Figure 4(b),
leads to a master curve of the form

$RF(\theta)=RF(0)\;\;g(\theta/\theta_{0})$ (18)

in which the contact angle is scaled by its critical value, which implicitly
describes the shape of the cross-section of the channel. We have not
determined a precise form for the function $g$, which is roughly sigmoidal.

(a)

(b)

Figure 4: Liquid recovery factor for different contact angles in different
shaped channels. The rectangular channel has aspect ratio $H/W=2$. (a)
Recovery factor versus contact angle. (b) Recovery factor scaled by its value
for contact angle $\theta=0$ versus relative contact angle, defined as the
contact angle relative to the critical contact angle at which 100% recovery is
achieved. Numerical simulations are shown with dots. Lines are from eqs. (13),
(14) and (17).

##  4 Discussion and Conclusions

In a cylindrical channel the recovery factor is equal to one, irrespective of
the contact angle. The calculations that we have described here for channels
with corners give a more realistic idea of what is likely to happen in
geological pores. The channel with a cross-section in the form of a scalloped
triangle is probably one of the more extreme cases likely to be encountered,
with its sharp cusp-like corners, and in the case of perfect wetting (zero
contact angle) the recovery factor falls to around 73%.

Our theoretical predictions have largely captured the dependence of RF on
contact angle for the channel shapes that we consider here. The Surface
Evolver simulations provide useful validation and visualisation of the
results, and for more complicated channel shapes are likely to be
indispensable.

The collapse of the data in Figure 4(b) suggests that for any channel shape
for which an estimate of the recovery factor in the case of perfect wetting is
available, and if it possible to estimate the critical angle $\theta_{0}$,
then the recovery factor for different contact angles follows a similar curve
for all channel shapes, and an approximate prediction can be obtained from
this almost universal curve.

An even greater level of realism will be found by considering channels with a
cross-section that varies along its length. This of course brings new
theoretical and computational challenges. To make progress in their analysis
of flow between packed spheres, Singh et al [5] assumed that the interface
itself is spherical, while Hilden and Trumble [20] used Surface Evolver to
accurately calculate the liquid shape within a planar array of spheres. More
general situations will surely require more ambitious numerical simulations
than those described here, but we expect that the procedure that we have
outlined will be required.

## Acknowledgments

We thank K. Brakke for the development and support of the Surface Evolver. An
AberDOC PhD scholarship from Aberystwyth University is gratefully
acknowledged.

## References

  * [1] Dullien, F.A. Porous media: fluid transport and pore structure. Academic press; 2012.
  * [2] de Castro, A.R., Oostrom, M., Shokri, N.. Effects of shear-thinning fluids on residual oil formation in microfluidic pore networks. J. Colloid Interface Sci. 2016; 472:34-43.
  * [3] Wong, H., Morris, S. and Radke, C.J., Three-Dimensional Menisci in Polygonal Capillaries. J. Colloid Interface Sci. 1992; 148:317–336.
  * [4] Kagan, M. and Pinczewski, W.V. Menisci in a Diamond-Shaped Capillary. J. Colloid Interface Sci. 2000; 230:452–454.
  * [5] Singh K., Scholl H., Brinkmann M., Di Michiel M., Scheel M., Herminghaus S., Seemann R. The role of local instabilities in fluid invasion into permeable media. Scientific reports. 2017; 7:444.
  * [6] Sadiq T.A., Advani S.G., Parnas R.S. Experimental investigation of transverse flow through aligned cylinders. International Journal of Multiphase Flow. 1995; 21:755-74.
  * [7] Ma, S., Mason, G., Morrow, N.R. Effect of contact angle on drainage and imbibition in regular polygonal tubes. Coll. Surf. A: Physicochem. Eng. Aspects. 1996; 117:273-91.
  * [8] Lenormand, R., Zarcone, C. and Sarr, A. Mechanisms of the displacement of one fluid by another in a network of capillary ducts. Journal of Fluid Mechanics. 1983; 135:337-53.
  * [9] Lago, M. and Araujo, M. Threshold Pressure in Capillaries with Polygonal Cross Section. J. Coll. Interf. Sci., 2001; 243:219–226.
  * [10] Lago, M. and Araujo, M. Threshold capillary pressure in capillaries with curved sides. Physica A, 2003; 319:175–187.
  * [11] Rabbani, H.S., Joekar-Niasar, V., Shokri, N. Effects of intermediate wettability on entry capillary pressure in angular pores. J. Colloid Interface Sci. 2016; 473:34–43.
  * [12] De Gennes, P.G., Brochard-Wyart, F. and Quéré, D., 2004. Capillarity and wetting phenomena: drops, bubbles, pearls, waves, Springer, New York.
  * [13] Princen, H.M. Capillary phenomena in assemblies of parallel cylinders: II. Capillary rise in systems with more than two cylinders. J. Colloid Interface Sci., 1969; 30:359–371.
  * [14] Brakke K.A. The Surface Evolver. Exp. Math. 1992;1:141-65.
  * [15] Jettestuen, E., Helland, J. O., and Prodanović, M., A level set method for simulating capillary-controlled displacements at the pore scale with nonzero contact angles, Water Resour. Res.; 49: 4645-4661.
  * [16] Helland, J.O. and Skjæveland, S.M., Physically-Based Capillary Pressure Correlation For Mixed-Wet Reservoirs From A Bundle-Of-Tubes Model. SPE Journal. 2006; 11:171–180.
  * [17] Rossen, W.R. A critical review of Roof snap-off as a mechanism of steady-state foam generation in homogeneous porous media. Coll. Surf. A: Physicochem. Eng. Aspects, 2003; 225:1-24.
  * [18] Mason, G. and Morrow, N.R. Meniscus Curvatures in Capillaries of Uniform Cross-section. J. Chem. Soc., Faraday Trans. 1. 1984; 80:2375–2393.
  * [19] Mayer, R.P. and Stowe, R.A. Mercury porosimetry – breakthrough pressure for penetration between packed spheres. J. Coll. Sci., 1965; 20:893–911.
  * [20] Hilden, J.L. and Trumble, K.P. Numerical analysis of capillarity in packed spheres: Planar hexagonal-packed spheres. J. Coll. Interf. Sci., 2003; 267:463–474.


# Study of $J/{\psi}$ ${\to}$ $D_{s,d}V$ decays with perturbative QCD approach

Yueling Yang  Institute of Particle and Nuclear Physics, Henan Normal
University, Xinxiang 453007, China  Junfeng Sun  Institute of Particle and
Nuclear Physics, Henan Normal University, Xinxiang 453007, China  Jie Gao
Institute of Particle and Nuclear Physics, Henan Normal University, Xinxiang
453007, China  Qin Chang  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
Gongru Lu  Institute of Particle and Nuclear Physics, Henan Normal University,
Xinxiang 453007, China

###### Abstract

Inspired by the recent measurements on two-body nonleptonic $J/{\psi}$ weak
decay at BESIII, the charm-changing $J/{\psi}$ ${\to}$ $D_{s,d}V$ weak decays
are studied with perturbative QCD approach, where $V$ denotes ${\rho}$ and
$K^{\ast}$ vector mesons. It is found that branching ratio for $J/{\psi}$
${\to}$ $D_{s}{\rho}$ decay can reach up to ${\cal O}(10^{-9})$, which is
within the potential measurement capability of the future high-luminosity
experiments.

$J/{\psi}$ meson; weak decay; branching ratio; perturbative QCD

###### pacs:

13.25.Gv 12.39.St 14.40.Pq 14.65.Dw

##  I Introduction

The $J/{\psi}$ particle is bound state of $c\bar{c}$ pair with given quantum
numbers $I^{G}J^{PC}$ $=$ $0^{-}1^{--}$ pdg . Since its discovery in 1974 bnl
; slac , the $J/{\psi}$ meson is always a hot and active topic for particle
physicists. The $c\bar{c}$ pair of the $J/{\psi}$ meson annihilate mainly into
gluons, which provides a valuable resource to explore the properties of the
quark-gluon coupling and the invisible gluons, to search for various glueballs
and possible exotic hadrons. There are two hierarchies in the $J/{\psi}$ meson
and other heavy quarkonium, one is dynamical energy scales responsible for
production and decay interactions of particles, and the other is relative
velocity of $c$ quark111According to the power counting rules of
nonrelativistic quantum chromodynamics (NRQCD) prd46 ; prd51 ; rmp77 , there
are several distinct energy scales in charmonium, for example, typical three-
momentum $m_{c}v$ and kinetic energy $m_{c}v^{2}/2$, where $v$ ${\ll}$ $1$ is
the typical relative velocity of heavy quark. Those energy scales satisfy a
hierarchy relation $m_{c}$ ${\gg}$ $m_{c}v$ ${\gg}$ $m_{c}v^{2}$.. The
$J/{\psi}$ meson plays a prominent role in investigation of QCD dynamical.

A conspicuous property of the $J/{\psi}$ meson is its narrow decay width, only
about 30 ppm222ppm means percent per million, i.e. $10^{-6}$. of its mass. The
$J/{\psi}$ meson lies below the kinematic $D\bar{D}$ threshold. Its hadronic
decay into light hadrons violates the phenomenological Okubo-Zweig-Iizuka
rules o ; z ; i . Besides the decay dominated by the strong and
electromagnetic interactions, the $J/{\psi}$ can also decay via the weak
interaction within the standard model. In this paper, we will study the
$J/{\psi}$ ${\to}$ $D_{s,d}{\rho}$, $D_{s,d}K^{\ast}$ weak decays with
perturbative QCD (pQCD) approach pqcd1 ; pqcd2 ; pqcd3 .

Experimentally, thanks to the good performance of CLEO-c, BES, LHCb,
B-factories, and so on, plenty of $J/{\psi}$ data samples have been
accumulated. Recently, the $J/{\psi}$ ${\to}$ $D_{s}{\rho}$, $D_{u}K^{\ast}$
weak decays have been searched for at BESIII using part of the available
$J/{\psi}$ samples prd89.071101 . It is eagerly expected to have about
$10^{10}$ $J/{\psi}$ samples at BESIII per year with the designed luminosity
cpc36 , and over $10^{10}$ prompt $J/{\psi}$ samples at LHCb per $fb^{-1}$
data epjc71 , which offers opportunities to discover phenomena that have been
previously overlooked because of statistical limitations. So a careful
scrutiny of $J/{\psi}$ weak decays at high-luminosity dedicated experiments
may be possible in the future. In particular, the “flavor tag” of a single
charged $D$ meson from $J/{\psi}$ decay will precisely identify potential
signal from massive background. In addition, an abnormal large production rate
of single $D$ meson from $J/{\psi}$ decay would be a hint of new physics.

Theoretically, the $J/{\psi}$ ${\to}$ $D_{q}V$ decay is, in fact, induced by
$c$ ${\to}$ $q$ $+$ $W^{+}$ transition at quark level, where $q$ $=$ $s$ and
$d$, the virtual $W^{+}$ boson materializes into a pair of quarks which then
hadronizes into a vector meson $V$ $=$ ${\rho}$ and $K^{\ast}$. As it is well
known, there must be the participation of strong interaction in nonleptonic
$J/{\psi}$ weak decay, and $c$ quark mass is between perturbative and
nonperturbative domain. In recent years, some QCD-inspired methods, such as
pQCD approach pqcd1 ; pqcd2 ; pqcd3 , QCD factorization approach qcdf1 ; qcdf2
; qcdf3 , soft and collinear effective theory scet1 ; scet2 ; scet3 ; scet4 ,
have been fully formulated to explain nonleptonic $B$ decays. The $J/{\psi}$
${\to}$ $D_{q}V$ decays have been investigated based on collinear
approximation plb252 ; ijmpa14 ; ahep2013 ; ijmpa30 . In this paper, the
$J/{\psi}$ ${\to}$ $D_{s,d}V$ decays will be restudied based on $k_{T}$
factorization. It is expected to glean new insights into factorization
mechanism, nonperturbative dynamics, final state interactions, and so on, from
nonleptonic $J/{\psi}$ weak decay.

This paper is organized as follows. The theoretical framework and amplitudes
for $J/{\psi}$ ${\to}$ $D_{s,d}V$ decays are given in section II, followed by
numerical results and discussion in section III. Finally, we summarize in the
last section.

##  II theoretical framework

###  II.1 The effective Hamiltonian

Theoretically, one usually uses the effective Hamiltonian to describe hadron
weak decay, where hard contributions can be decently factorized based on
operator product expansion and the renormalization group (RG) method. The
effective Hamiltonian responsible for $J/{\psi}$ ${\to}$ $D_{s,d}V$ decay
could be written as 9512380 ,

${\cal H}_{\rm eff}\ =\
\frac{G_{F}}{\sqrt{2}}\,\sum\limits_{q_{1},q_{2}}V_{cq_{1}}V_{uq_{2}}^{\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}\,\text{GeV}^{-2}$ pdg  is Fermi
constant; $q_{1,2}$ $=$ $d$ and $s$.

The Cabibbo-Kobayashi-Maskawa (CKM) factors are written as

$\begin{array}[]{lcl}V_{cs}V_{ud}^{\ast}\ =\
1-{\lambda}^{2}-\frac{1}{2}A^{2}{\lambda}^{4}+{\cal
O}({\lambda}^{6}),&{}{}\hfil&\text{for}~{}J/{\psi}{\to}D_{s}{\rho}~{}\text{decay}\\\
V_{cs}V_{us}^{\ast}\ =\
{\lambda}-\frac{1}{2}{\lambda}^{3}-\frac{1}{8}{\lambda}^{5}(1+4A^{2})+{\cal
O}({\lambda}^{6}),&{}{}\hfil&\text{for}~{}J/{\psi}{\to}D_{s}K^{\ast}~{}\text{decay}\\\
V_{cd}V_{ud}^{\ast}\ =\
-V_{cs}V_{us}^{\ast}-A^{2}{\lambda}^{5}({\varrho}+i{\eta})+{\cal
O}({\lambda}^{6}),&{}{}\hfil&\text{for}~{}J/{\psi}{\to}D_{d}{\rho}~{}\text{decay}\\\
V_{cd}V_{us}^{\ast}\ =\ -{\lambda}^{2}+{\cal
O}({\lambda}^{6}),&{}{}\hfil&\text{for}~{}J/{\psi}{\to}D_{d}K^{\ast}~{}\text{decay}\end{array}$
(2)

where $A$, ${\lambda}$, ${\varrho}$, ${\eta}$ are the Wolfenstein parameters;
${\lambda}$ $=$ ${\sin}{\theta}_{c}$ ${\simeq}$ $0.2$ pdg  and ${\theta}_{c}$
is the Cabibbo angle. It is clearly seen that the $J/{\psi}$ ${\to}$
$D_{s}{\rho}$ decay is favored by the CKM factor $V_{cs}V_{ud}^{\ast}$.

The Wilson coefficients $C_{1,2}(\mu)$ summarize the physical contributions
above the scales of ${\mu}$. They are calculated at scale of the $W$ boson
mass ${\mu}$ ${\sim}$ ${\cal O}(m_{W})$ with perturbation theory, and then
evolved to scale of the $c$ quark mass ${\mu}$ ${\sim}$ ${\cal O}(m_{c})$ with
RG evolution function,

$\vec{C}({\mu})=U_{4}({\mu},m_{b})U_{5}(m_{b},m_{W})\vec{C}(m_{W}),$ (3)

where $U_{f}({\mu}_{j},{\mu}_{i})$ is RG evolution matrix 9512380 . The Wilson
coefficients are independent of a particular process in the same role of
universal gauge couplings. They have properly been evaluated to the next-to-
leading order.

Generally, the penguin contributions induced by flavor changing neutral
current transitions are proportional to small Wilson coefficients relative to
tree contributions. Besides, for $c$ quark decay, the penguin contributions
are also severely suppressed by the CKM factors $V_{cd}V_{ud}^{\ast}$ $+$
$V_{cs}V_{us}^{\ast}$ $=$ $-V_{cb}V_{ub}^{\ast}$ ${\sim}$ ${\cal
O}({\lambda}^{5})$. Hence, only the tree operators related to $W$ emission
contributions are considered here. The expressions of tree operators are

$\displaystyle Q_{1}$ $\displaystyle=$
$\displaystyle[\bar{q}_{1,{\alpha}}{\gamma}_{\mu}(1-{\gamma}_{5})c_{\alpha}][\bar{u}_{\beta}{\gamma}^{\mu}(1-{\gamma}_{5})q_{2,{\beta}}],$
(4) $\displaystyle Q_{2}$ $\displaystyle=$
$\displaystyle[\bar{q}_{1,{\alpha}}{\gamma}_{\mu}(1-{\gamma}_{5})c_{\beta}][\bar{u}_{\beta}{\gamma}^{\mu}(1-{\gamma}_{5})q_{2,{\alpha}}],$
(5)

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

The physical contributions below scales of ${\mu}$ are included in hadronic
matrix elements (HME). Because of the participation of the strong interaction,
the entanglement perturbative and nonperturbative effects, the inadequate
comprehension of hadronization mechanism and low energy QCD behavior, HME is
the most complicated and intractable part. To get the amplitude, one has to
face directly the HME calculation.

###  II.2 Hadronic matrix elements

Phenomenologically, the simplest approximation is that HME is parameterized
into the production of transition form factors and decay constant based on
naive factorization (NF) scheme nf . The NF treatment on HME deprives any
physical mechanism that could illustrate strong phases and rescattering among
participating hadrons, and loses the ${\mu}$ dependence of HME which must
exist to cancel that of Wilson coefficients. So the Lepage-Brodsky hard
scattering approach prd22  is usually used, and HME is generally expressed as
the convolution of hard scattering kernel with distribution amplitudes (DAs),
where DAs reflect nonperturbative contributions but are universal. The hard
part is, in principle, perturbatively calculable as a power of series of
coupling ${\alpha}_{s}$. To suppress the soft contributions and avoid the
problem of the endpoint singularity from collinear assumption qcdf1 ; qcdf2 ;
qcdf3 , the transverse momentum of quarks are retained explicitly and the
Sudakov factors are introduced for each of meson wave functions in evaluation
of potentially infrared contributions with pQCD approach pqcd1 ; pqcd2 ; pqcd3
. Finally, a decay amplitude could be written as a convolution integral of
three parts pqcd1 ; pqcd2 ; pqcd3 : the hard effects enclosed by the Wilson
coefficients $C_{i}$, the rescattering kernel amplitudes ${\cal H}$, and
process-independent wave functions ${\Phi}$,

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

where $k$ is the momentum of valence quarks, $t$ is a typical scale and
$e^{-S}$ is a Sudakov factor.

###  II.3 Kinematic variables

In the center-of-mass frame of $J/{\psi}$ meson, the light-cone kinematic
variables are defined as follows.

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

where the subscript $i$ $=$ $1$, $2$, $3$ on variables, including polarization
vector ${\epsilon}_{i}$, four dimensional momentum $p_{i}$, energy $E_{i}$ and
mass $m_{i}$, correspond to initial $J/{\psi}$ meson, recoiled $D$ meson,
emitted vector meson $V$ $=$ ${\rho}$ and $K^{\ast}$, respectively; $x_{i}$
and $k_{i}$ ($\vec{k}_{iT}$) denote the longitudinal momentum fraction and
(transverse) momentum of valence quarks, respectively; $n_{+}$ is the plus
null vector; $s$, $t$ and $u$ are Lorentz transformation scalars; $p$ is the
common momentum of final states. The kinematic variables are displayed in
Fig.2(a).

###  II.4 Wave functions

Taking the convention of Ref. prd65 ; jhep0703 , HME of the diquark operators
squeezed between the vacuum and meson state is defined as below.

${\langle}0{|}c_{i}(z)\bar{c}_{j}(0){|}{\psi}(p_{1},{\epsilon}_{1}^{\parallel}){\rangle}\,=\,\frac{f_{\psi}}{4}{\int}d^{4}k_{1}\,e^{-ik_{1}{\cdot}z}\Big{\\{}\\!\\!\not{\epsilon}_{1}^{\parallel}\Big{[}m_{1}\,{\phi}_{\psi}^{v}(k_{1})-\\!\\!\not{p}_{1}\,{\phi}_{\psi}^{t}(k_{1})\Big{]}\Big{\\}}_{ji},$
(19)
${\langle}0{|}c_{i}(z)\bar{c}_{j}(0){|}{\psi}(p_{1},{\epsilon}_{1}^{\perp}){\rangle}\,=\,\frac{f_{\psi}}{4}{\int}d^{4}k_{1}\,e^{-ik_{1}{\cdot}z}\Big{\\{}\\!\\!\not{\epsilon}_{1}^{\perp}\Big{[}m_{1}\,{\phi}_{\psi}^{V}(k_{1})-\\!\\!\not{p}_{1}\,{\phi}_{\psi}^{T}(k_{1})\Big{]}\Big{\\}}_{ji},$
(20)
${\langle}D_{q}(p_{2}){|}\bar{q}_{i}(z)c_{j}(0){|}0{\rangle}=\frac{if_{D_{q}}}{4}{\int}d^{4}k_{2}\,e^{ik_{2}{\cdot}z}\,\Big{\\{}{\gamma}_{5}\Big{[}\\!\\!\not{p}_{2}\,{\Phi}_{D}^{a}(k_{2})+m_{2}\,{\Phi}_{D}^{p}(k_{2})\Big{]}\Big{\\}}_{ji},$
(21)
$\displaystyle{\langle}V(p_{3},{\epsilon}_{3}^{\parallel}){|}u_{i}(z)\bar{q}_{j}(0){|}0{\rangle}$
(22) $\displaystyle=$
$\displaystyle\frac{1}{4}{\int}_{0}^{1}dk_{3}\,e^{ik_{3}{\cdot}z}\Big{\\{}\\!\\!\not{\epsilon}_{3}^{\parallel}\,m_{3}\,{\Phi}_{V}^{v}(k_{3})+\\!\\!\not{\epsilon}_{3}^{\parallel}\\!\\!\not{p}_{3}\,{\Phi}_{V}^{t}(k_{3})-m_{3}\,{\Phi}_{V}^{s}(k_{3})\Big{\\}}_{ji},$
$\displaystyle{\langle}V(p_{3},{\epsilon}_{3}^{{\perp}}){|}u_{i}(z)\bar{q}_{j}(0){|}0{\rangle}\
=\
\frac{1}{4}{\int}_{0}^{1}dk_{3}\,e^{ik_{3}{\cdot}z}\Big{\\{}\\!\\!\not{\epsilon}_{3}^{{\perp}}\,m_{3}\,{\Phi}_{V}^{V}(k_{3})$
$\displaystyle+\\!\\!\not{\epsilon}_{3}^{{\perp}}\\!\\!\not{p}_{3}\,{\Phi}_{V}^{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}_{V}^{A}(k_{3})\Big{\\}}_{ji},$
(23)

where $f_{\psi}$ and $f_{D_{q}}$ are decay constants; wave functions
${\Phi}_{{\psi},V}^{v,T}$ and ${\Phi}_{D}^{a}$ are twist-2; wave functions
${\Phi}_{{\psi},V}^{t,s,V,A}$ and ${\Phi}_{D}^{p}$ are twist-3. For wave
functions of light vector meson, only ${\Phi}_{V}^{v}$ and ${\Phi}_{V}^{V,A}$
are involved in decay amplitudes (see Appendix A). Their expressions are prd65
; jhep0703 :

${\phi}_{V}^{v}(x)=6\,x\,\bar{x}\,\Big{\\{}1+\sum\limits_{i=1}a_{i}^{V}\,C_{i}^{3/2}(t)\Big{\\}},$
(24) ${\phi}_{V}^{V}(x)\,=\,\frac{3}{4}\,(1+t^{2}),$ (25)
${\phi}_{V}^{A}(x)\,=\,\frac{3}{2}\,(-t).$ (26)

where $\bar{x}$ $=$ $1$ $-$ $x$ and $t$ $=$ $\bar{x}$ $-$ $x$; $a_{i}^{V}$ is
nonperturbative Gegenbauer moment and corresponds to Gegenbauer polynomial
$C_{i}^{3/2}(t)$.

$C_{1}^{3/2}(t)=3\,t,\quad
C_{2}^{3/2}(t)=\frac{3}{2}\,(5\,t^{2}-1),\quad{\cdots}$ (27)

With the relation of $m_{J/{\psi}}$ ${\simeq}$ $2m_{c}$ and $m_{D_{q}}$
${\simeq}$ $m_{c}$ $+$ $m_{q}$, it is suspected that the motion of valence
quarks in $J/{\psi}$ and $D_{q}$ mesons is nearly nonrelativistic. So their
spectrum can be described with time-independent Schrödinger equation. Suppose
the interaction between valence quarks is an isotropic harmonic oscillator
potential, the ground state eigenfunction with quantum numbers $nL$ $=$ $1S$
is expressed as:

${\phi}_{1S}(\vec{k})\ {\sim}\ e^{-\vec{k}^{2}/2{\omega}^{2}},$ (28)

where parameter ${\omega}$ determines the average transverse momentum,
${\langle}1S{|}k^{2}_{T}{|}1S{\rangle}$ $=$ ${\omega}^{2}$. By using the
transformation xiao ,

$\vec{k}^{2}\ {\to}\
\frac{1}{4}\Big{(}\frac{\vec{k}_{T}^{2}+m_{q_{1}}^{2}}{x_{1}}+\frac{\vec{k}_{T}^{2}+m_{q_{2}}^{2}}{x_{2}}\Big{)},$
(29)

then integrating out transverse momentum $k_{T}$ and combining with their
asymptotic forms, finally, DAs for $J/{\psi}$ and $D$ mesons are written as

${\phi}_{\psi}^{v}(x)={\phi}_{\psi}^{T}(x)=A\,x\,\bar{x}\,{\exp}\Big{\\{}-\frac{m_{c}^{2}}{8\,{\omega}_{1}^{2}\,x\,\bar{x}}\Big{\\}},$
(30)
${\phi}_{\psi}^{t}(x)=B\,t^{2}\,{\exp}\Big{\\{}-\frac{m_{c}^{2}}{8\,{\omega}_{1}^{2}\,x\,\bar{x}}\Big{\\}},$
(31)
${\phi}_{\psi}^{V}(x)=C\,(1+t^{2})\,{\exp}\Big{\\{}-\frac{m_{c}^{2}}{8\,{\omega}_{1}^{2}\,x\,\bar{x}}\Big{\\}},$
(32)
${\phi}_{D}^{a}(x)=D\,x\,\bar{x}\,{\exp}\Big{\\{}-\frac{\bar{x}\,m_{q}^{2}+x\,m_{c}^{2}}{8\,{\omega}_{2}^{2}\,x\,\bar{x}}\Big{\\}},$
(33)
${\phi}_{D}^{p}(x)=E\,{\exp}\Big{\\{}-\frac{\bar{x}\,m_{q}^{2}+x\,m_{c}^{2}}{8\,{\omega}_{2}^{2}\,x\,\bar{x}}\Big{\\}},$
(34)

where ${\omega}_{i}$ $=$ $m_{i}{\alpha}_{s}$ according to the NRQCD power
counting rules prd46 , coefficients of $A$, $B$, $C$, $D$, $E$ are determined
by the normalization conditions,

${\int}_{0}^{1}dx\,{\phi}_{\psi}^{v,t,V,T}(x)=1,$ (35)
${\int}_{0}^{1}dx\,{\phi}_{D}^{a,p}(x)=1.$ (36)

Here, it should be pointed out that there are many wave function models for
$D$ meson, for example, Eq.(30) in Ref. prd78lv . The preferred one in Ref.
prd78lv  is:

${\phi}_{D}(x,b)=6\,x\bar{x}\,\Big{\\{}1+C_{D}(1-2x)\Big{\\}}{\exp}\Big{\\{}-\frac{1}{2}w^{2}b^{2}\Big{\\}},$
(37)

where $C_{D}$ $=$ $0.4$ and $w$ $=$ $0.2$ GeV for $D_{s}$ meson; $C_{D}$ $=$
$0.5$ and $w$ $=$ $0.1$ GeV for $D_{d}$ meson. In addition, the same form of
Eq.(37), without a distinction between twist-2 and twist-3, is used in many
practical calculation.

Figure 1: The shape lines of DAs for $J/{\psi}$ meson in (a) and $D_{d,s}$
mesons in (b), where ${\phi}_{\psi}^{v,t,V}(x)$, ${\phi}_{D}^{a,p}(x)$ and
${\phi}_{D}(x,b)$ correspond to Eqs.(30,31,32), Eqs.(33,34) and Eq.(37),
respectively.

The shape lines of DAs for $J/{\psi}$ and $D_{s,d}$ mesons are displayed in
Fig.1. It is clearly seen that (1) DAs for $J/{\psi}$ meson is symmetric
versus $x$, and a broad peak of ${\phi}_{D}^{a,p}(x)$ appears at $x$ $<$ $0.5$
regions, which is basically in line with the picture that momentum fraction is
proportional to valence quark mass. (2) under the influence of exponential
functions, DAs of Eqs.(30—34) fall quickly down to zero at endpoint $x$,
$\bar{x}$ ${\to}$ $0$, which is bound to suppress soft contributions. (3) The
flavor symmetry breaking effects between $D_{d}$ and $D_{s}$ mesons, and
difference between twist-2 and twist-3 are obvious in Eqs.(33,34) rather than
Eq.(37). In this paper, we will use DAs of Eqs.(33,34) for $D$ meson.

###  II.5 Decay amplitudes

The Feynman diagrams for $J/{\psi}$ ${\to}$ $D_{s}{\rho}$ decay are shown in
Fig.2, including factorizable emission topologies (a) and (b) where gluon
connects $J/{\psi}$ with $D_{s}$ meson, and nonfactorizable emission
topologies (c) and (d) where gluon couples the spectator quark with emitted
${\rho}$ meson.

Figure 2: Feynman diagrams for $J/{\psi}$ ${\to}$ $D_{s}{\rho}$ decay,
including factorizable diagrams (a) and (b), and nonfactorizable diagrams (c)
and (d).

The amplitude for $J/{\psi}$ ${\to}$ $D_{q}V$ decay is written as prd66

${\cal A}(J/{\psi}{\to}D_{q}V)={\cal
A}_{L}({\epsilon}_{1}^{{\parallel}},{\epsilon}_{3}^{{\parallel}})+{\cal
A}_{N}({\epsilon}_{1}^{{\perp}}{\cdot}{\epsilon}_{3}^{{\perp}})+i\,{\cal
A}_{T}\,{\varepsilon}_{{\mu}{\nu}{\alpha}{\beta}}\,{\epsilon}_{1}^{{\mu}}\,{\epsilon}_{3}^{{\nu}}\,p_{1}^{\alpha}\,p_{3}^{\beta},$
(38)

which is conventionally written as helicity amplitudes prd66 ,

${\cal A}_{0}\ =\ -{\cal F}\,\sum\limits_{i}{\cal
A}_{i,L}({\epsilon}_{1}^{{\parallel}},{\epsilon}_{3}^{{\parallel}}),$ (39)
${\cal A}_{\parallel}\ =\ \sqrt{2}\,{\cal F}\sum\limits_{i}{\cal A}_{i,N},$
(40) ${\cal A}_{\perp}\ =\ \sqrt{2}\,{\cal F}\,m_{1}\,p\sum\limits_{i}{\cal
A}_{i,T},$ (41) ${\cal F}\ =\
i\frac{G_{F}}{\sqrt{2}}\,\frac{C_{F}}{N_{c}}\,{\pi}\,f_{\psi}\,f_{D_{q}}\,f_{V}\,V_{cq_{1}}V_{uq_{2}}^{\ast},$
(42)

where the color number $N_{c}$ $=$ $3$ and color factor $C_{F}$ $=$
$(N_{c}^{2}-1)/2N_{c}$; the subscript $i$ on ${\cal A}_{i,j}$ corresponds to
indices of Fig.2. The expressions of building blocks ${\cal A}_{i,j}$ can be
found in Appendix A. Our results show that (1) factorizable contributions
[Fig.2 (a) and (b)] are color-favored, i.e., $a_{1}$-dominated; (2)
nonfactorizable contributions [Fig.2 (c) and (d)] are proportion to small
Wilson coefficient $C_{2}$ and suppressed by color factor $1/N_{c}$.

##  III Numerical results and discussion

In the rest frame of $J/{\psi}$ meson, branching ratio is defined as

${\cal B}r\
=\frac{1}{12{\pi}}\,\frac{p}{m_{\psi}^{2}{\Gamma}_{\psi}}\,\Big{\\{}{|}{\cal
A}_{0}{|}^{2}+{|}{\cal A}_{\parallel}{|}^{2}+{|}{\cal
A}_{\perp}{|}^{2}\Big{\\}}.$ (43)

Table 1: The numerical values of input parameters. CKM parameters333The
relations between CKM parameters (${\varrho}$, ${\eta}$) and ($\bar{\varrho}$,
$\bar{\eta}$) are pdg : $({\varrho}+i{\eta})$ $=$
$\displaystyle\frac{\sqrt{1-A^{2}{\lambda}^{4}}(\bar{\varrho}+i\bar{\eta})}{\sqrt{1-{\lambda}^{2}}[1-A^{2}{\lambda}^{4}(\bar{\varrho}+i\bar{\eta})]}$.
pdg
---
$A$ $=$ $0.814^{+0.023}_{-0.024}$, ${\lambda}$ $=$ $0.22537{\pm}0.00061$,
$\bar{\varrho}$ $=$ $0.117{\pm}0.021$, $\bar{\eta}$ $=$ $0.353{\pm}0.013$,
mass and decay constants
$m_{\psi}$ $=$ $3096.916{\pm}0.011$ MeV pdg , |  $m_{D_{s}}$ $=$ $1968.30{\pm}0.11$ MeV pdg , |  $m_{D_{d}}$ $=$ $1869.61{\pm}0.10$ MeV pdg ,
$f_{\psi}$ $=$ $395.1{\pm}5.0$ MeV444The decay constant $f_{\psi}$ can be obtained from experimental branching ratios for electromagnetic $J/{\psi}$ decay into charged lepton pairs through the formula  ${\cal B}r(J/{\psi}{\to}{\ell}^{+}{\ell}^{-})=\frac{16{\pi}}{27}f_{\psi}^{2}\frac{{\alpha}_{\rm QED}^{2}}{m_{\psi}\,{\Gamma}_{\psi}}\sqrt{1-4\frac{m_{\ell}^{2}}{m_{\psi}^{2}}}\Big{\\{}1+2\frac{m_{\ell}^{2}}{m_{\psi}^{2}}\Big{\\}},$ (44) where ${\alpha}_{\rm QED}$ is the fine-structure constant, ${\ell}$ $=$ $e$ and ${\mu}$. One can get pdg  $f_{\psi}$ $=$ $395.4{\pm}7.0$ MeV with ${\cal B}r(J/{\psi}{\to}e^{+}e^{-})$ $=$ $(5.971{\pm}0.032)\%$, and $f_{\psi}$ $=$ $394.8{\pm}7.1$ MeV with ${\cal B}r(J/{\psi}{\to}{\mu}^{+}{\mu}^{-})$ $=$ $(5.961{\pm}0.033)\%$, respectively, where the errors arise from mass $m_{\psi}$, decay width ${\Gamma}_{\psi}$ and branching ratios. The weighted average is $f_{\psi}$ $=$ $395.1{\pm}5.0$ MeV., |  $f_{D_{s}}$ $=$ $257.5{\pm}4.6$ MeV pdg , |  $f_{D_{d}}$ $=$ $204.6{\pm}5.0$ MeV pdg ,
$m_{K^{\ast}}$ $=$ $891.66{\pm}0.26$ MeV pdg , |  $m_{\rho}$ $=$ $775.26{\pm}0.25$ MeV pdg , |  $m_{c}$ $=$ $1.67{\pm}0.07$ GeV pdg ,
$f_{K^{\ast}}$ $=$ $220{\pm}5$ MeV jhep0703 , |  $f_{\rho}$ $=$ $216{\pm}3$ MeV jhep0703 , |  ${\Gamma}_{\psi}$ $=$ $92.9{\pm}2.8$ keV pdg ,
$m_{s}$ ${\approx}$ $510$ MeV uds , |  $m_{d}$ ${\approx}$ $310$ MeV uds , |
Gegenbauer moments at ${\mu}$ $=$ 1 GeV jhep0703
$a_{1}^{K^{\ast}}$ $=$ $-0.03{\pm}0.02$, $a_{2}^{K^{\ast}}$ $=$
$0.11{\pm}0.09$, $a_{1}^{\rho}$ $=$ $0$, $a_{2}^{\rho}$ $=$ $0.15{\pm}0.07$.
Table 2: Branching ratios for $J/{\psi}$ ${\to}$ $DV$ decays, where uncertainties of our results come from scale $(1{\pm}0.1)t_{i}$, quark mass $m_{c}$, hadronic parameters and CKM parameters, respectively. Reference |  ijmpa14 555The updated results are listed in Table 4 of Ref. ahep2013 . | ahep2013  | ijmpa30  | epjc55  | this work
---|---|---|---|---|---
$10^{9}{\times}{\cal B}r(J/{\psi}{\to}D_{s}{\rho})$ | $2.54$ | $5.1$ | $2.2$ | $1.3$ | $3.33^{+0.97+0.47+0.17+0.002}_{-0.42-0.51-0.17-0.002}$
$10^{10}{\times}{\cal B}r(J/{\psi}{\to}D_{s}K^{\ast})$ | $1.48$ | $2.8$ | $1.2$ | $0.8$ | $1.86^{+0.57+0.28+0.12+0.010}_{-0.24-0.35-0.12-0.010}$
$10^{10}{\times}{\cal B}r(J/{\psi}{\to}D_{d}{\rho})$ | $1.54$ | $2.2$ | $1.1$ | $0.4$ | $1.32^{+0.37+0.14+0.08+0.007}_{-0.16-0.19-0.08-0.007}$
$10^{11}{\times}{\cal B}r(J/{\psi}{\to}D_{d}K^{\ast})$ | … | $1.3$ | $0.6$ | … | $0.80^{+0.23+0.05+0.06+0.009}_{-0.10-0.12-0.06-0.009}$

The values of input parameters are listed in Table 1, where if it is not
specified explicitly, their central values will be taken as the default
inputs. Our numerical results are presented in Table 2, where the first
uncertainty comes from the choice of the typical scale $(1{\pm}0.1)t_{i}$, and
expression of $t_{i}$ is given in Eq.(70) and Eq.(71); the second uncertainty
is from quark mass $m_{c}$; the third uncertainty is from hadronic parameters
including decay constants and Gegenbauer moments; and the fourth uncertainty
of branching ratio comes from CKM parameters. The following are some comments.

(1) As it is aforementioned, the $J/{\psi}$ decay modes considered here are
dominated by the color-favored factorizable contributions and insensitive to
nonfactorizable contributions. So, generally, branching ratio for a given
$J/{\psi}$ ${\to}$ $D_{s,d}V$ decay has the same order of magnitude even with
different phenomenological models.

(2) There is a clear hierarchical pattern among branching ratios, mainly
resulting from the hierarchical structure of CKM factors in Eq.(2), i.e.,

${\cal B}r(J/{\psi}{\to}D_{s}{\rho})\,{\gg}\,{\cal
B}r(J/{\psi}{\to}D_{s}K^{\ast})\,{\sim}\,{\cal
B}r(J/{\psi}{\to}D_{d}{\rho})\,{\gg}\,{\cal B}r(J/{\psi}{\to}D_{d}K^{\ast}).$
(45)

In addition, because nonfactorizable contributions are suppressed by both
small $C_{2}$ and color factor $1/N_{c}$, there is an approximate
relationship,

$\frac{{\cal B}r(J/{\psi}{\to}D_{s}K^{\ast})}{{\cal
B}r(J/{\psi}{\to}D_{s}{\rho})}\ {\approx}\ \frac{{\cal
B}r(J/{\psi}{\to}D_{d}K^{\ast})}{{\cal B}r(J/{\psi}{\to}D_{d}{\rho})}\
{\approx}\ {\lambda}^{2}\,\frac{f^{2}_{K^{\ast}}}{f^{2}_{\rho}}.$ (46)

Above all, the Cabibbo- and color-favored $J/{\psi}$ ${\to}$ $D_{s}{\rho}$
decay has branching ratio ${\sim}$ ${\cal O}(10^{-9})$, which is well within
the measurement capability of the future high-luminosity experiments, such as
super tau-charm factory, LHC and SuperKEKB.

Figure 3: Contributions to branching ratio ${\cal
B}r(J/{\psi}{\to}D_{s}{\rho})$ versus ${\alpha}_{s}/{\pi}$, where the numbers
over histogram denote the percentage of the corresponding contributions.

(3) Here, one might question the practicability of pQCD approach and the
feasibility of perturbative calculation because $c$ quark mass seems to be not
large enough. To clear this issue up or to check what percentage of
contributions come from perturbative domain, contributions to branching ratio
${\cal B}r(J/{\psi}{\to}D_{s}{\rho})$ from different ${\alpha}_{s}/{\pi}$
region are displayed in Fig.3. It is easily seen that about 80% contributions
come from ${\alpha}_{s}/{\pi}$ ${\leq}$ 0.4 regions, which implies that the
calculation with pQCD approach is valid. One of crucial reasons for the small
percentage in the region ${\alpha}_{s}/{\pi}$ ${\leq}$ $0.1$ is that the
absolute values of Wilson coefficients $C_{1,2}$, $a_{1}$ and coupling
${\alpha}_{s}$ decrease along with the increase of renormalization scale
${\mu}$. Of course, a perturbative calculation with pQCD approach is
influenced by many factors, such as Sudakov factors, the choice of scale $t$,
models of wave functions, etc., which deserve much attention but beyond the
scope of this paper.

(4) There are many uncertainties on branching ratios. The first uncertainty
from scale $t$ could be reduced by the inclusion of higher order corrections
to HME and an improved control on nonperturbative contributions. The second
uncertainty from wave function models or parameter $m_{c}$ will be greatly
lessened with the relative rate of branching ratios, for example, Eq.(46). The
third uncertainty is dominated by decay constants whose effects will be
weakened with the increasing precision of experimental measurements and/or
theoretical calculation using nonperturbative methods (such as lattice QCD and
so on). The uncertainty from CKM factor is small. Moreover, other factors,
such as the final state interactions which is important and necessary for $c$
quark decay, are not properly considered here, but deserve massive dedicated
study. Our results just provide an order of magnitude estimation on branching
ratio.

##  IV Summary

Within the standard model, the $J/{\psi}$ meson can decay via the weak
interaction, besides the strong and electromagnetic interactions. With
anticipation of copious $J/{\psi}$ data samples at the future high-luminosity
experiments and gradual improvement of particle identification techniques, we
investigated the charm-changing $J/{\psi}$ ${\to}$ $D_{s,d}{\rho}$,
$D_{s,d}K^{\ast}$ weak decays with pQCD approach. It is found that the
estimated branching ratio for the color- and CKM-favored $J/{\psi}$ ${\to}$
$D_{s}{\rho}$ decay can be up to ${\cal O}(10^{-9})$, which is very likely to
be measured in the future.

## Acknowledgments

We thank Professor Dongsheng Du (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, 11275057 and U1332103).

##  Appendix A Building blocks of decay amplitudes

The expressions of building blocks ${\cal A}_{i,j}$ are listed as follows,
where subscript $i$ corresponds to indices of Fig.2; and $j$ corresponds to
helicity amplitudes.

$\displaystyle{\cal A}_{a,L}$ $\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}_{\psi}^{v}(x_{1})\,E_{a}(t_{a})$
(47) $\displaystyle{\times}$ $\displaystyle
H_{a}({\alpha},{\beta}_{a},b_{1},b_{2})\,{\alpha}_{s}(t_{a})\,a_{1}(t_{a})\,\Big{\\{}{\phi}_{D}^{p}(x_{2})\,m_{2}\,m_{c}\,u$
$\displaystyle+{\phi}_{D}^{a}(x_{2})\,\Big{[}m_{1}^{2}\,s-(4\,m_{1}^{2}\,p^{2}+m_{2}^{2}\,u)\,\bar{x}_{2}\Big{]}\Big{\\}},$
$\displaystyle{\cal A}_{a,N}$ $\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}_{\psi}^{V}(x_{1})$
(48) $\displaystyle{\times}$ $\displaystyle
E_{a}(t_{a})\,H_{a}({\alpha},{\beta}_{a},b_{1},b_{2})\,\Big{\\{}-2\,m_{2}\,m_{c}\,{\phi}_{D}^{p}(x_{2})$
$\displaystyle+{\phi}_{D}^{a}(x_{2})\,\Big{[}2\,m_{2}^{2}\,\bar{x}_{2}-t\Big{]}\Big{\\}}\,{\alpha}_{s}(t_{a})\,a_{1}(t_{a}),$
$\displaystyle{\cal A}_{a,T}$ $\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}_{\psi}^{V}(x_{1})$
(49) $\displaystyle{\times}$
$\displaystyle{\phi}_{D}^{a}(x_{2})\,E_{a}(t_{a})\,H_{a}({\alpha},{\beta}_{a},b_{1},b_{2})\,{\alpha}_{s}(t_{a})\,a_{1}(t_{a}),$
$\displaystyle{\cal A}_{b,L}$ $\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}\,H_{b}({\alpha},{\beta}_{b},b_{2},b_{1})$
(50) $\displaystyle{\times}$ $\displaystyle
E_{b}(t_{b})\,\Big{\\{}{\phi}_{\psi}^{v}(x_{1})\,{\phi}_{D}^{a}(x_{2})\,\Big{[}m_{1}^{2}\,(s-4\,p^{2})\,\bar{x}_{1}-m_{2}^{2}\,u\Big{]}$
$\displaystyle+2\,m_{1}\,m_{2}\,{\phi}_{\psi}^{t}(x_{1})\,{\phi}_{D}^{p}(x_{2})\,(s-u\,\bar{x}_{1})\Big{\\}}\,{\alpha}_{s}(t_{b})\,a_{1}(t_{b}),$
$\displaystyle{\cal A}_{b,N}$ $\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}\,H_{b}({\alpha},{\beta}_{b},b_{2},b_{1})$
(51) $\displaystyle{\times}$ $\displaystyle
E_{b}(t_{b})\,{\alpha}_{s}(t_{b})\,\Big{\\{}m_{1}\,m_{3}\,{\phi}_{\psi}^{V}(x_{1})\,{\phi}_{D}^{a}(x_{2})\,(2\,m_{2}^{2}-t\,\bar{x}_{1})$
$\displaystyle+2\,m_{2}\,m_{3}\,{\phi}_{\psi}^{T}(x_{1})\,{\phi}_{D}^{p}(x_{2})\,(2\,m_{1}^{2}\,\bar{x}_{1}-t)\Big{\\}}\,a_{1}(t_{b}),$
$\displaystyle{\cal A}_{b,T}$ $\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}\,H_{b}({\alpha}_{e},{\beta}_{b},b_{2},b_{1})\,E_{b}(t_{b})$
(52) $\displaystyle{\times}$
$\displaystyle{\alpha}_{s}(t_{b})\,a_{1}(t_{b})\,\Big{\\{}2\,m_{2}\,{\phi}_{\psi}^{T}(x_{1})\,{\phi}_{D}^{p}(x_{2})-m_{1}\,{\phi}_{\psi}^{V}(x_{1})\,{\phi}_{D}^{a}(x_{2})\,\bar{x}_{1}\Big{\\}},$
$\displaystyle{\cal A}_{c,L}$ $\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{\times}$
$\displaystyle{\phi}_{\rho}^{v}(x_{3})\,E_{c}(t_{c})\,H_{c}({\alpha},{\beta}_{c},b_{2},b_{3})\,{\alpha}_{s}(t_{c})\,{\delta}(b_{1}-b_{2})$
$\displaystyle{\times}$ $\displaystyle
C_{2}(t_{c})\,\Big{\\{}{\phi}_{\psi}^{v}(x_{1})\,{\phi}_{D}^{a}(x_{2})\,u\,\Big{[}t\,(\bar{x}_{2}-\bar{x}_{1})+s\,(x_{3}-\bar{x}_{2})\Big{]}$
$\displaystyle+{\phi}_{\psi}^{t}(x_{1})\,{\phi}_{D}^{p}(x_{2})\,m_{1}\,m_{2}\,\Big{[}u\,(\bar{x}_{1}-x_{3})+s\,(x_{3}-\bar{x}_{2})\Big{]}\Big{\\}},$
$\displaystyle{\cal A}_{c,N}$ $\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{\times}$ $\displaystyle
E_{c}(t_{c})\,H_{c}({\alpha},{\beta}_{c},b_{2},b_{3})\,{\alpha}_{s}(t_{c})\,C_{2}(t_{c})\,{\delta}(b_{1}-b_{2})$
$\displaystyle{\times}$
$\displaystyle\Big{\\{}{\phi}_{\psi}^{V}(x_{1})\,{\phi}_{D}^{a}(x_{2})\,{\phi}_{\rho}^{V}(x_{3})\,2\,m_{1}\,\Big{[}s\,(\bar{x}_{2}-x_{3})+t\,(\bar{x}_{1}-\bar{x}_{2})\Big{]}$
$\displaystyle+{\phi}_{\psi}^{T}(x_{1})\,{\phi}_{D}^{p}(x_{2})\,{\phi}_{\rho}^{V}(x_{3})\,m_{2}\,\Big{[}u\,(x_{3}-\bar{x}_{1})+t\,(\bar{x}_{2}-\bar{x}_{1})\Big{]}$
$\displaystyle+{\phi}_{\psi}^{T}(x_{1})\,{\phi}_{D}^{p}(x_{2})\,{\phi}_{\rho}^{A}(x_{3})\,2\,m_{1}\,m_{2}\,p\,(x_{3}-\bar{x}_{2})\Big{\\}},$
$\displaystyle{\cal A}_{c,T}$ $\displaystyle=$
$\displaystyle\frac{1}{N_{c}}\frac{m_{3}}{m_{1}\,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{\times}$ $\displaystyle
E_{c}(t_{c})\,H_{c}({\alpha},{\beta}_{c},b_{2},b_{3})\,{\alpha}_{s}(t_{c})\,C_{2}(t_{c})\,{\delta}(b_{1}-b_{2})$
$\displaystyle{\times}$
$\displaystyle\Big{\\{}{\phi}_{\psi}^{V}(x_{1})\,{\phi}_{D}^{a}(x_{2})\,{\phi}_{\rho}^{A}(x_{3})\,2\,m_{1}\,\Big{[}s\,(\bar{x}_{2}-x_{3})+t\,(\bar{x}_{1}-\bar{x}_{2})\Big{]}$
$\displaystyle+{\phi}_{\psi}^{T}(x_{1})\,{\phi}_{D}^{p}(x_{2})\,{\phi}_{\rho}^{A}(x_{3})\,m_{2}\,\Big{[}u\,(x_{3}-\bar{x}_{1})+t\,(\bar{x}_{2}-\bar{x}_{1})\Big{]}$
$\displaystyle+{\phi}_{\psi}^{T}(x_{1})\,{\phi}_{D}^{p}(x_{2})\,{\phi}_{\rho}^{V}(x_{3})\,2\,m_{1}\,m_{2}\,p\,(x_{3}-\bar{x}_{2})\Big{\\}},$
$\displaystyle{\cal A}_{d,L}$ $\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{\times}$
$\displaystyle{\phi}_{\rho}^{v}(x_{3})\,E_{d}(t_{d})\,H_{d}({\alpha},{\beta}_{d},b_{2},b_{3})\,{\alpha}_{s}(t_{d})\,{\delta}(b_{1}-b_{2})$
$\displaystyle{\times}$
$\displaystyle\Big{\\{}{\phi}_{\psi}^{t}(x_{1})\,{\phi}_{D}^{p}(x_{2})\,m_{1}\,m_{2}\,\Big{[}u\,(x_{3}-x_{1})+s\,(x_{2}-x_{3})\Big{]}$
$\displaystyle+{\phi}_{\psi}^{v}(x_{1})\,{\phi}_{D}^{a}(x_{2})\,4\,m_{1}^{2}\,p^{2}\,(x_{3}-x_{2})\Big{\\}}\,C_{2}(t_{d}),$
$\displaystyle{\cal A}_{d,N}$ $\displaystyle=$
$\displaystyle\frac{m_{2}\,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{\times}$
$\displaystyle{\phi}_{\psi}^{T}(x_{1})\,{\phi}_{D}^{p}(x_{2})\,E_{d}(t_{d})\,H_{d}({\alpha},{\beta}_{d},b_{2},b_{3})\,{\alpha}_{s}(t_{d})\,C_{2}(t_{d})\,{\delta}(b_{1}-b_{2})$
$\displaystyle{\times}$
$\displaystyle\Big{\\{}{\phi}_{\rho}^{V}(x_{3})\,\Big{[}2\,m_{1}^{2}\,x_{1}-t\,x_{2}-u\,x_{3}\Big{]}+2\,m_{1}\,p\,{\phi}_{\rho}^{A}(x_{3})\,(x_{2}-x_{3})\Big{\\}},$
$\displaystyle{\cal A}_{d,T}$ $\displaystyle=$
$\displaystyle\frac{1}{N_{c}}\,\frac{m_{2}\,m_{3}}{m_{1}\,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{\times}$
$\displaystyle{\phi}_{\psi}^{T}(x_{1})\,{\phi}_{D}^{p}(x_{2})\,E_{d}(t_{d})\,H_{d}({\alpha},{\beta}_{d},b_{2},b_{3})\,{\alpha}_{s}(t_{d})\,C_{2}(t_{d})\,{\delta}(b_{1}-b_{2})$
$\displaystyle{\times}$
$\displaystyle\Big{\\{}{\phi}_{\rho}^{A}(x_{3})\,\Big{[}2\,m_{1}^{2}\,x_{1}-t\,x_{2}-u\,x_{3}\Big{]}+2\,m_{1}\,p\,{\phi}_{\rho}^{V}(x_{3})\,(x_{2}-x_{3})\Big{\\}},$

where $b_{i}$ is the conjugate variable of the transverse momentum $k_{iT}$;
${\alpha}_{s}$ is the QCD running coupling; $a_{1}$ $=$ $C_{1}$ $+$
$C_{2}/N_{c}$; $C_{1,2}$ are the Wilson coefficients.

The hard scattering function $H_{i}$ and Sudakov factor $E_{i}$ are defined as
follows.

$H_{a(b)}({\alpha},{\beta},b_{i},b_{j})=K_{0}(b_{i}\sqrt{-{\alpha}})\Big{\\{}{\theta}(b_{i}-b_{j})K_{0}(b_{i}\sqrt{-{\beta}})I_{0}(b_{j}\sqrt{-{\beta}})+(b_{i}{\leftrightarrow}b_{j})\Big{\\}},$
(59) $\displaystyle H_{c(d)}({\alpha},{\beta},b_{2},b_{3})$ $\displaystyle=$
$\displaystyle\Big{\\{}{\theta}(-{\beta})K_{0}(b_{3}\sqrt{-{\beta}})+\frac{{\pi}}{2}{\theta}({\beta})\Big{[}iJ_{0}(b_{3}\sqrt{{\beta}})-Y_{0}(b_{3}\sqrt{{\beta}})\Big{]}\Big{\\}}$
(60) $\displaystyle{\times}$
$\displaystyle\Big{\\{}{\theta}(b_{2}-b_{3})K_{0}(b_{2}\sqrt{-{\alpha}})I_{0}(b_{3}\sqrt{-{\alpha}})+(b_{2}{\leftrightarrow}b_{3})\Big{\\}},$
$E_{i}(t)\ =\
\left\\{\begin{array}[]{ll}{\exp}\\{-S_{\psi}(t)-S_{D}(t)\\},&~{}~{}\text{ for
}i=a,b\\\ {\exp}\\{-S_{\psi}(t)-S_{D}(t)-S_{V}(t)\\},&~{}~{}\text{ for
}i=c,d\end{array}\right.$ (61) $S_{\psi}(t)\ =\
s(x_{1},p_{1}^{+},1/b_{1})+2{\int}_{1/b_{1}}^{t}\frac{d{\mu}}{\mu}{\gamma}_{q},$
(62) $S_{D}(t)\ =\
s(x_{2},p_{2}^{+},1/b_{2})+2{\int}_{1/b_{2}}^{t}\frac{d{\mu}}{\mu}{\gamma}_{q},$
(63) $S_{V}(t)\ =\
s(x_{3},p_{3}^{+},1/b_{3})+s(\bar{x}_{3},p_{3}^{+},1/b_{3})+2{\int}_{1/b_{3}}^{t}\frac{d{\mu}}{\mu}{\gamma}_{q},$
(64)

where $I_{0}$, $J_{0}$, $K_{0}$, $Y_{0}$ are Bessel functions; the expression
of $s(x,Q,1/b)$ can be found in Ref.pqcd1 ; ${\gamma}_{q}$ $=$
$-{\alpha}_{s}/{\pi}$ is the quark anomalous dimension; ${\alpha}$ and
${\beta}_{i}$ are gluon and quark virtuality, respectively, where subscript
$i$ on ${\beta}_{i}$ corresponds to indices of Fig.2.

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

## References

  * (1) K. Olive et al. (Particle Data Group), Chin. Phys. C 38, 090001 (2014).
  * (2) J. Aubert et al., Phys. Rev. Lett. 33, 1404 (1974).
  * (3) J. Augustin et al., Phys. Rev. Lett. 33, 1406 (1974).
  * (4) G. Lepage et al., Phys. Rev. D 46, 4052 (1992).
  * (5) G. Bodwin, E. Braaten, G. Lepage, Phys. Rev. D 51, 1125 (1995).
  * (6) N. Brambilla et al., Rev. Mod. Phys. 77, 1423 (2005).
  * (7) S. Okubo, Phys. Lett. 5, 165 (1963).
  * (8) G. Zweig, CERN-TH-401, 402, 412 (1964).
  * (9) J. Iizuka, Prog. Theor. Phys. Suppl. 37-38, 21 (1966).
  * (10) H. Li, Phys. Rev. D 52, 3958 (1995).
  * (11) C. Chang, H. Li, Phys. Rev. D 55, 5577 (1997).
  * (12) T. Yeh, H. Li, Phys. Rev. D 56, 1615 (1997).
  * (13) M. Ablikim et al. (BESIII Collaboration), Phys. Rev. D 89, 071101 (2014).
  * (14) H. Li, S. Zhu, Chin. Phys. C 36, 932 (2012).
  * (15) R. Aaij et al. (LHCb Collaboration), Eur. Phys. J. C 71, 1645 (2011).
  * (16) M. Beneke et al., Phys. Rev. Lett. 83, 1914 (1999).
  * (17) M. Beneke et al., Nucl. Phys. B 591, 313 (2000).
  * (18) M. Beneke et al., Nucl. Phys. B 606, 245 (2001).
  * (19) C. Bauer et al., Phys. Rev. D 63, 114020 (2001).
  * (20) C. Bauer, D. Pirjol, I. Stewart, Phys. Rev. D 65, 054022 (2002).
  * (21) C. Bauer et al., Phys. Rev. D 66, 014017 (2002).
  * (22) M. Beneke et al., Nucl. Phys. B 643, 431 (2002).
  * (23) R. Verma, A. Kamal and A. Czarnecki, Phys. Lett. B 252, 690 (1990).
  * (24) K. Sharma and R. Verma, Int. J. Mod. Phys. A 14, 937 (1999).
  * (25) R. Dhir, R. Verma and A. Sharma, Adv. High Energy Phys, 2013, 706543 (2013).
  * (26) J. Sun et al., Int. J. Mod. Phys. A 30, 1550094 (2015).
  * (27) G. Buchalla, A. Buras, M. Lautenbacher, Rev. Mod. Phys. 68, 1125, (1996).
  * (28) D. Fakirov and B. Stech, Nucl. Phys. B 133, 315 (1978).
  * (29) G. Lepage, S. Brodsky, Phys. Rev. D 22, 2157 (1980).
  * (30) T. Kurimoto, H. Li, A. Sanda, Phys. Rev. D 65, 014007 (2001).
  * (31) P. Ball and G. Jones, JHEP, 0703, 069, (2007).
  * (32) B. Xiao, X. Qin, B. Ma, Eur. Phys. J. A 15, 523 (2002).
  * (33) R. Li, C. Lü, H. Zou, Phys. Rev. D 78, 014018 (2008).
  * (34) C. Chen, Y. Keum, H. Li, Phys. Rev. D 66, 054013 (2002)
  * (35) A. Kamal, Particle Physics, Springer, 2014, p. 298.
  * (36) Y. Wang et al., Eur. Phys. J. C 55, 607 (2008).