# Self-sustained elastoinertial Tollmien-Schlichting waves

Ashwin Shekar1  Ryan M. McMullen2  Beverley J. McKeon2  Michael D. Graham1
<EMAIL_ADDRESS>1Department of Chemical and Biological Engineering,
University of Wisconsin-Madison, Madison, WI 53706, USA 2Graduate Aerospace
Laboratories, California Institute of Technology, Pasadena CA 91125, USA

###### Abstract

Direct simulations of two-dimensional plane channel flow of a viscoelastic
fluid at Reynolds number $\Rey=3000$ reveal the existence of a family of
attractors whose structure closely resembles the linear Tollmien-Schlichting
(TS) mode, and in particular exhibits strongly localized stress fluctuations
at the critical layer position of the TS mode. At the parameter values chosen,
this solution branch is not connected to the nonlinear TS solution branch
found for Newtonian flow, and thus represents a solution family that is
nonlinearly self-sustained by viscoelasticity. The ratio between stress and
velocity fluctuations is in quantitative agreement for the attractor and the
linear TS mode, and increases strongly with Weissenberg number, Wi. For the
latter, there is a transition in the scaling of this ratio as Wi increases,
and the Wi at which the nonlinear solution family comes into existence is just
above this transition. Finally, evidence indicates that this branch is
connected through an unstable solution branch to two-dimensional
elastoinertial turbulence (EIT). These results suggest that, in the parameter
range considered here, the bypass transition leading to EIT is mediated by
nonlinear amplification and self-sustenance of perturbations that excite the
Tollmien-Schlichting mode.

##  1 Introduction

Adding minute quantities (parts per million) of long-chain polymer additives
can dramatically change the turbulent flow of Newtonian fluids, the most
significant consequence being the considerable drop in friction factor, which
is commonly referred to as the Toms effect (Virk, 1975; White & Mungal, 2008;
Graham, 2014). Accompanying this overall change is a structural change to the
flow. At high levels of viscoelasticity, Samanta et al. (2013) and Dubief et
al. (2013) have shown that trains of weak spanwise-oriented flow structures
with inclined sheets of polymer stretch dominate the flow, denoting this
regime as elastoinertial turbulence (EIT). These sheets of polymer stretch
correspond with a layer near each wall of localized spanwise vortex motion. In
further contrast to the 3D structures that sustain Newtonian turbulence, Sid
et al. (2018) demonstrated that EIT is fundamentally 2D in nature by showing
that the structures sustaining 2D EIT in channel flow simulations are similar
to those in 3D. In a computational study of EIT in pipe flow, Lopez et al.
(2019) also observe nearly two-dimensional spanwise vortices in the near-wall
region, indicating the structural similarities between EIT in channels and
pipes.

Choueiri et al. (2018) experimentally studied the path to EIT in pipe flow by
varying polymer concentration at fixed Reynolds number, $\Rey$. For
sufficiently low $\Rey$, they observed an initial drop in friction factor
(i.e. modification of Newtonian turbulence) as concentration increased,
followed by re-laminarization and eventually by a reentrant transition to EIT,
where the flow had very different structure from Newtonian turbulence. These
observations point at two distinct types of turbulence in dilute polymer
solutions – one that is suppressed by viscoelasticity (Newtonian turbulence)
and one that is promoted (EIT).

Shekar et al. (2019) corroborated this observation of a reentrant transition
to EIT in simulations of channel flow with increasing Weissenberg number, Wi,
the ratio between the polymer relaxation time scale and the shear time scale.
They further showed that close to its inception, EIT exhibits localized
polymer stress fluctuations that bear strong resemblance to critical layer
structures predicted by linear analyses, i.e., sheetlike fluctuations
localized at wall-normal locations where the disturbance wavespeed equals the
base flow velocity. In particular, they demonstrated that the fluctuation
structure corresponding to the dominant spectral content strongly resembles
the viscoelastic extension of the linear Tollmien-Schlichting (TS) wave. This
is perhaps a surprising result, as the flow in the parameter regime considered
is linearly stable, and in Newtonian turbulence, the TS mode plays a very
limited role. Some light is shed on this issue through resolvent analysis,
i.e., determination of the response of the linearized dynamics to harmonic-in-
time disturbances, which shows that the linear TS mode becomes highly
amplified in the presence of viscoelasticity. This strong amplification
implies that even very weak disturbances may be sufficient to trigger the
nonlinear effects necessary to sustain EIT.

We note that similar structures have been observed by other researchers in
different contexts. Page & Zaki (2015) analyzed the evolution of vortical
perturbations in 2D viscoelastic simple shear flow. Their analysis reveals a
viscoelastic analogue of the Newtonian Orr mechanism. This “reverse-Orr”
mechanism generates tilted sheets of polymer stress fluctuations resembling
those seen at EIT and thus may play some role in this phenomenon.

Because prior work on elastoinertial turbulence reveals structures similar to
those seen in the linear Tollmien-Schliching mode, the present work focuses on
Tollmien-Schlichting waves, but in the fully nonlinear context of self-
sustained solutions in the channel flow geometry. (In the parameter regime
here, the laminar flow is always linearly stable.) After introducing the
formulation and computational methods, we show how the Newtonian nonlinear TS
wave branch is modified by viscoelasticity, resulting in its disappearance as
Wi increases. At still higher Wi, however, we demonstrate the onset of a new,
viscoelasticity-driven, nonlinear solution branch that strongly resembles the
linear Tollmien-Schlichting mode, and illustrate how it is related to the TS
mode of linear theory and to elastoinertial turbulence.

##  2 Formulation

This study focuses on two-dimensional pressure-driven channel flow with
constant mass flux. The $x$ and $y$ axes are aligned with the streamwise and
wall-normal directions, respectively. Lengths are scaled by the half channel
height $l$, so the dimensionless channel height $L_{y}=2$. The domain is
periodic in $x$ with length $L_{x}$. Velocity $\boldsymbol{v}$ is scaled with
the Newtonian laminar centerline velocity $U$; time $t$ with $l/U$, and
pressure $p$ with $\rho U^{2}$, where $\rho$ is the fluid density. The polymer
stress tensor $\boldsymbol{\tau}_{p}$ is related to the polymer conformation
tensor $\boldsymbol{\alpha}$ through the FENE-P constitutive relation, which
models each polymer molecule as a pair of beads connected by a nonlinear
spring with maximum extensibility $b$.

We solve the momentum, continuity and FENE-P equations:

$\displaystyle\frac{\partial\boldsymbol{v}}{\partial
t}+\boldsymbol{v}\cdot\boldsymbol{\nabla
v}=-\boldsymbol{\nabla}p+\frac{\beta}{Re}\nabla^{2}\boldsymbol{v}+\frac{\left(1-\beta\right)}{Re\textit{Wi}}\left(\boldsymbol{\nabla}\cdot\boldsymbol{\tau}_{\mathrm{p}}\right),$
(1) $\displaystyle\boldsymbol{\nabla}\cdot\boldsymbol{v}=0,$ (2)
$\displaystyle\boldsymbol{\tau}_{p}=\frac{\boldsymbol{\alpha}}{1-\frac{\mathrm{tr}(\boldsymbol{\alpha})}{b}}-\boldsymbol{I},$
(3) $\displaystyle\frac{\partial\boldsymbol{\alpha}}{\partial
t}+\boldsymbol{v}\cdot\boldsymbol{\nabla\alpha}-\boldsymbol{\alpha}\cdot\boldsymbol{\nabla
v}-\left(\boldsymbol{\alpha}\cdot\boldsymbol{\nabla
v}\right)^{\mathrm{T}}=\frac{-1}{\textit{Wi}}\boldsymbol{\tau}_{p}.$ (4)

Here $Re=\rho Ul/(\eta_{\mathrm{s}}+\eta_{\mathrm{p}})$, where $\eta_{s}$ and
$\eta_{p}$ are the solvent and polymer contributions to the zero-shear rate
viscosity. The viscosity ratio
$\beta=\eta_{\mathrm{s}}/(\eta_{\mathrm{s}}+\eta_{\mathrm{p}})$. We fix
$\beta=0.97$ and $b=6400$. Since $1-\beta$ is proportional to polymer
concentration and $b$ to the number of monomer units, this parameter set
corresponds to a dilute solution of a high molecular weight polymer. The
Weissenberg number $\textit{Wi}=\lambda U/l$, where $\lambda$ is the polymer
relaxation time.

For the nonlinear direct numerical simulations (DNS) described below, a finite
difference scheme and a fractional time step method are adopted for
integrating the Navier-Stokes equation. Second-order Adams-Bashforth and
Crank-Nicolson methods are used for convection and diffusion terms,
respectively. The FENE-P equation is discretized using a high resolution
central difference scheme (Kurganov & Tadmor, 2000; Vaithianathan et al.,
2006; Dallas et al., 2010) that guarantees positive definiteness of the
polymer conformation tensor without the need for any artificial diffusion. In
any case, the nonlinear solution branch on which we focus in this manuscript
displays weak fluctuations far from the limit of positive definiteness even at
the highest Wi of existence. A typical resolution for the following results is
$(N_{x},N_{y})$ = ${(79,402)}$. This resolution used was based on mesh
convergence results at $\textit{Wi}=45$. When the resolution was increased to
$(N_{x},N_{y})$ = ${(131,602)}$, the mean polymer stretch deviations from the
laminar base state change by less than 1 percent.

We also consider the linearized evolution of infinitesimal perturbations to
the laminar state with given streamwise wavenumber $k$. Two approaches are
used. The first is classical linear stability analysis, in which solutions of
the form $\boldsymbol{\hat{\phi}}(y)\exp\left[ik(x-ct)\right]$ are sought,
resulting in an eigenvalue problem for the complex wavespeed $c$ at a given
$k$. In the present study, $~{}\hat{}~{}$ always indicates deviation from the
laminar state. If all $c_{i}<0$, which is the case for all conditions
considered in the present study, the flow is linearly stable. A linearized
version of the DNS code was also developed using the numerical schemes
described above. Results were validated against linear stability analysis and
agreement to three significant digits was obtained for the value of $c$ for
the viscoelastic TS mode at the parameters of interest.

The second linear approach used here determines the linear response of the
laminar flow to external forcing with given wavenumber $k$ and frequency
$\omega$ using the resolvent operator of the linearized equations (Schmid,
2007; McKeon & Sharma, 2010). The norm used in the resolvent calculations is

$\lVert\boldsymbol{\hat{\phi}}\rVert^{2}_{\boldsymbol{A}}=\int_{-1}^{1}\left[\boldsymbol{\hat{v}}^{*}\boldsymbol{\hat{v}}+\mathrm{tr}\\!\left(\boldsymbol{A}^{-1}\boldsymbol{\hat{\alpha}}^{*}\boldsymbol{A}^{-1}\boldsymbol{\hat{\alpha}}\right)\right]\mathrm{d}y,$
(5)

where $\boldsymbol{A}$ is the conformation tensor in the laminar state. The
second term provides a measure of the conformation tensor perturbation
magnitude that is motivated by the non-Euclidean geometry of positive-definite
tensors (Hameduddin et al., 2019). For both the linear stability and resolvent
analyses, the equations are discretized with a Chebyshev pseudospectral method
using 401 Chebyshev polynomials. This number was arrived at by ensuring
convergence of the TS eigenvalue.

##  3 Results and discussion

###  3.1 Origin of Newtonian and viscoelastic nonlinear Tollmien-Schlichting
attractors

Figure 1: Structure of NNTSA at $\Rey=3000,\textit{Wi}=3$. Streamlines (blue)
in a reference frame moving at the wave speed $c=0.39$ are superimposed on
color contours of $\hat{\alpha}_{xx}$ where $\hat{}$ denotes deviations from
the laminar state. Blue dots indicate the locations of hyperbolic stagnation
points (in the traveling frame). Black contour lines of $\hat{v}$ are also
shown. For $\hat{v}$, dashed = negative and solid = positive.

In Newtonian flow, a family of nonlinear Tollmien-Schlichting waves bifurcates
subcritically from the laminar branch at $\Rey\approx 5772$ with $L_{x}\approx
2\pi/1.02\approx 6.15$. The lower limit of the parameter regime for which this
solution family exists is $\Rey\approx 2800,L_{x}\approx 2\pi/1.3\approx 4.83$
(Jiménez, 1990). Furthermore, in prior work on elastoinertial turbulence
(Shekar et al., 2019), as well as more recent simulations in long two-
dimensional domains, a strong peak in the spatial spectrum is found at
$L_{x}\approx 5$. Based on these observations, all of the results presented in
this study will be at $\Rey=3000,L_{x}=5$. (At $L_{x}=5$, Newtonian channel
flow is linearly stable at all $\Rey$.) In the Newtonian limit at these
parameters, there are upper and lower branch solutions (which merge in a
saddle-node bifurcation as $\Rey$ is lowered); the upper branch traveling wave
solution is linearly stable with respect to two-dimensional perturbations and
is thus easily computed via DNS. We call this solution branch, including its
viscoelastic extension, the Newtonian Nonlinear Tollmien-Schlichting attractor
(NNTSA). (The word “attractor” is chosen rather than “wave” because, depending
on parameters one can observe a pure traveling wave state or one with periodic
or nonperiodic modulations.)

On increasing $Wi$, the self-sustained nonlinear viscoelastic TS wave at
$Re=3000$ develops sheets of high polymer stretch that start out from near the
wall. These observations are evidence of the capability of nonlinear TS
critical layer mechanisms in generating sheets of polymer stretch. Figure 1
illustrates this point with a snapshot of $\hat{\alpha}_{xx}$ on the NNTSA
branch at $Wi=3$. The sheets originate in the nonlinear Kelvin cat’s eye
kinematics of TS waves at finite amplitude, as detailed in Shekar et al.
(2019). The NNTSA continues to display wall normal velocity fluctuations that
extend across the channel centerline – a signature of TS kinematics. Some of
the observations made in Shekar et al. (2019) are repeated here for
completeness, as they form the background for the new results of the present
study.

At the parameters chosen, the solution branch originating in the self-
sustained Newtonian TS wave bifurcates to a periodic orbit at
$\textit{Wi}\approx 3.5$ (cf. Lee & Zaki (2017)) before turning back into a
traveling wave and losing existence beyond $\textit{Wi}=3.875$, evidently in a
saddle-node bifurcation yielding a lower branch TS wave solution that becomes
the Newtonian solution as $\textit{Wi}\rightarrow 0$. Consistent with a
saddle-node bifurcation, if the solution at $\textit{Wi}=3.875$ is used as an
initial condition for a simulation at slightly higher Wi, the flow
laminarizes. This bifurcation scenario is shown on Fig. 2a in terms of average
wall shear rate vs. Wi. The unstable lower branch (dashed blue) was found
using edge tracking (Zammert & Eckhardt (2014)) between NNTSA and laminar
solutions at a given Wi. A bisection technique was used to arrive at
arbitrarily close initial conditions that are on either side of the edge. DNS
trajectories starting from such points stay on the edge for a while before
diverging to NNTSA and laminar.

As shown in Shekar et al. (2019), if Wi is large, sufficiently energetic
initial conditions lead to 2D EIT. Fig. 2(a) also shows the mean wall shear
rate for the EIT solution branch, which loses existence at finite amplitude
when $\textit{Wi}\lesssim 13$. The bifurcation underlying this transition is
presumably also of saddle-node form.

The central observation of the present paper arises from considering what
happens just below the onset of the EIT regime at $\textit{Wi}\approx 13$. We
do this by using a velocity and stress field from EIT at $\textit{Wi}=13$ as
an initial condition for a run at $\textit{Wi}=12$. This initial condition
persists as a slowly decaying form of EIT for hundreds of time units,
consistent with behavior just beyond a saddle-node bifurcation. As time
increases further, the structure continues to decay, but does not ultimately
reach the laminar state. Instead, it evolves to a nontrivial attractor state
that is very nearly a traveling wave, and in particular strongly resembles the
linear TS mode at these parameters. We call this new state the viscoelastic
nonlinear Tollmien-Schlichting attractor (VNTSA).

Figure 2: (a) Bifurcation diagram showing the evolution of the space and time-
averaged wall shear rate with $Wi$ for the 2D nonlinear NNTSA and EIT branches
at $Re=3000,L_{x}=5$. Point (a) corresponds to the structure shown in Figure
1. Labels ‘L’ and ‘VNTSA’ indicate that initial conditions starting at the
arrows evolve to laminar or VNTSA states, respectively. (b) Bifurcation
diagram of the $L_{\infty}$-norm of $\hat{\alpha}_{xx}$ with Wi for VNTSA and
EIT. Points (a)-(f) correspond to structures shown later in Figure 6. The
label ‘EIT’ indicates that initial conditions evolve to EIT. In both plots,
dashed lines correspond to the unstable solution branches obtained from edge
tracking.

###  3.2 Viscoelastic nonlinear Tollmien-Schlichting attractor

To elaborate on the relationship between the VNTSA and the linear TS mode, we
describe results at $Re=3000,\textit{Wi}=13,L_{x}=5$, i.e. close to the point
where the 2D EIT branch first comes into existence as shown in the bifurcation
diagram (Figure 2a). EIT and the VNTSA are coexisting attractors at these
parameter values. Figure 3 shows the evolution of the $L_{\infty}$ norm of
$\hat{\alpha}_{xx}$ starting from an initial condition consisting of the
laminar state plus some amplitude $\epsilon$ of the linear TS mode for this
parameter set. This mode, with $\epsilon=1$, is shown in Figure 4a. The
structure of the velocity field is virtually unchanged from the Newtonian case
and the polymer conformations are strongly localized to the critical layer
positions $y=\pm 0.82$. Sufficiently small perturbations, e.g. $\epsilon=1$,
decay to the laminar state, as they must since that state is linearly stable.
However, larger perturbations $\epsilon={\color[rgb]{0,0,0}5},10$ and 100,
where nonlinear mechanisms play a role, $\hat{\alpha}_{xx}$ settles to a
finite value corresponding to the VNTSA. The initial condition $\epsilon=5$
that starts from below the VNTSA in $||\hat{\alpha}_{xx}||_{\infty}$ shows an
initial growth phase before saturating onto the VNTSA, whereas $\epsilon=10$
and 100 relax onto the VNTSA from above.

Figure 3: Time evolution of the $L_{\infty}$-norm of $\hat{\alpha}_{xx}$ for
$Re=3000,\textit{Wi}=13,L_{x}=5$, starting from an initial condition of
laminar state + $\epsilon\,\times\,$TS-mode. Dashed lines correspond to
linearized runs starting from the same initial conditions for $\epsilon=10$
and $100$.

For comparison, the dashed lines on Figure 3 show the linearized evolution
starting from the same initial conditions; these all decay to laminar,
illustrating the role of nonlinearity in the transition to the VNTSA. This
state is robust: initial perturbation amplitudes over a wide range will evolve
to it. However, initial conditions with very large magnitudes (e.g. $\epsilon$
= 6000) evolve to EIT: as noted above, both EIT and VNTSA are attractors at
the chosen parameters (as is the laminar state).

We have also used initial conditions of the laminar state plus velocity
perturbations somewhat similar to the ones used by Page & Zaki (2015). These
perturbations satisfied incompressibility and were sinusoidal in nature with
streamwise and wall normal periods equal to the domain size. When appropriate
magnitudes are used, transient growth of polymer stress followed by a decay to
the VNTSA is observed.

In 3D channel flow simulations of Oldroyd-B fluids at $\Rey=3000$,
$\textit{Wi}=4$ to 6, $\beta=0.9$, Min et al. (2003) observe a transient
($\sim 500$ time units) before an abrupt jump to the fully developed state.
These observations were made starting from turbulent Newtonian initial
conditions. In the present work, we used an average of about 2500 TU of data,
and in some cases more than 5000 TU to calculate the statistics, and have
observed no similar transition in any simulations of VNTSA in the parameter
space considered here.

Figure 4b is a snapshot showing the typical fluctuation structure of the VNTSA
at $\textit{Wi}=13$. The streamwise conformation $\hat{\alpha}_{xx}$ has
tilted sheets highly localized near $y=\pm 0.82$ and contours of wall normal
velocity $\hat{v}$ span the entire channel. This structure bears strong
resemblance to the TS mode shown in Figure 4a. The VNTSA is thus a weakly
nonlinear self-sustaining state whose primary structure is the viscoelastic TS
mode. We elaborate in the following section on the linear TS mode and its
connections to the VNTSA.

Figure 4: (a) Structure of the linear TS mode at
$Re=3000,\textit{Wi}=13,L_{x}=5$. Magnitude of the eigenmode is arbitrary and
values shown here correspond to $\epsilon=1$. (b): Snapshot of the fluctuation
structure of the VNTSA at $Re=3000,\textit{Wi}=13$. (c) and (d): The
$kL_{x}/2\pi=1,2$ components respectively, of the snapshot shown in (b). Shown
are contour lines of $\hat{v}$ superimposed on color contours of
$\hat{\alpha}_{xx}$.

In the VNTSA state, the velocity fluctuations are very weak, and the mean wall
shear rate displays a very small change from laminar. This can be understood
on the grounds that changes of the mean wall shear rate correspond to
fluctuations with $k_{x}=0$, which arise only due to nonlinear interactions.
Since the primary velocity structure is very weak, the nonlinear effects will
be even weaker. To illustrate nonlinear effects, Figures 4c and 4d,
respectively, show the $kL_{x}/2\pi=1,2$, spatial Fourier components of the
snapshot shown in 4b. Figure 4c closely resembles the TS mode, with a slight
symmetry-breaking across the centerline $y=0$. The structure at
$kL_{x}/2\pi=2$ also displays polymer stress fluctuations localized around the
critical layer position, an observation that also holds for higher
wavenumbers.

Having established the structure of the flow on the VNTSA branch, we now
illustrate the bifurcation scenario of this solution branch by continuing in
Wi. The VNTSA branch loses existence at finite amplitude (i.e. in a saddle-
node bifurcation) for $\textit{Wi}\lesssim 6$, as we have confirmed both by
using the $\textit{Wi}=6$ solution as an initial condition for simulations at
lower Wi and by running simulations starting from the laminar state perturbed
by the TS mode with small $\epsilon$. For $Wi<6$ all these initial conditions
decay to laminar. On increasing Wi, the VNTSA branch seems to lose existence
beyond $\textit{Wi}\approx 49$, and initial conditions that land on the VNTSA
for $\textit{Wi}=49$ evolve to EIT at $\textit{Wi}=50$. Due to its weak nature
in relation to EIT, the bifurcation scenario associated with the two solution
branches is shown on Figure 2b in a log scale using the $L_{\infty}$ norm of
$\hat{\alpha}_{xx}$ as the amplitude measure. The EIT solutions in Figure 2a
are depicted again in Figure 2b as the upper stable branch (solid red) with
this measure.

Since the VNTSA and EIT are both stable states, we were also able to perform
edge tracking to find unstable solutions intermediate between these two
states. Five Wi values (13, 15, 20, 40 and 45) were studied. Repeated
bisections were performed until trajectories stayed on the edge for an average
of 300 TU. The red-dashed line on Figures 2a and b indicates solutions (all
time-dependent) on this intermediate branch. The magnitude of fluctuations
along this branch monotonically decreases on increasing Wi, displaying values
close to EIT at $\textit{Wi}=13$ and values to close to VNTSA at
$\textit{Wi}=45$. Furthermore, this branch also displays fluctuations that
resemble the viscoelastic TS mode, as illustrated below. A more detailed link
between VNTSA and EIT might be established through numerical continuations of
underlying traveling wave solutions; this is a topic of future endeavors.

To complete this discussion, we note that the bifurcation scenario we observe
implies the existence of an edge between the VNTSA and the laminar state,
which could in principle also be found using edge tracking. However, the weak
nature of the VNTSA implies that this edge would be even weaker, thus making
this a challenging task.

Figure 5a shows the fluctuation structure of the VNTSA at $\textit{Wi}=8$,
close to the point where it first comes into existence. The structure closely
resembles the TS mode and does not change appreciably with time. The flow is
almost a pure nonlinear traveling wave with some weak non-periodicity, as
indicated in the power density plot of the wall normal velocity at position
$(0,0.825)$ shown in Figure 5b. The spectrum is mainly composed of the
dominant TS mode frequency and its higher harmonics. The dynamics and
structures get more complicated as Wi increases. Figure 5c shows a typical
snapshot at $\textit{Wi}=20$, which clearly is more complex than a TS mode.
However, at this Wi, the VNTSA still intermittently displays clear TS-like
structures such as the snapshot in Figure 5d.

We now turn to the flow structures at various positions on the bifurcation
diagram, Figure 2b. Figures 6a, b and c are representative snapshots of VNTSA,
the intermediate branch and EIT respectively, at $\textit{Wi}=13$ i.e. near
the loss of existence of EIT. As detailed in the description of Figure 4 and
shown again in Figure 6a, VNTSA exhibits a weak structure that strongly
resembles the viscoelastic TS mode. At this low Wi, the intermediate branch
(Figure 6b) displays a much stronger fluctuation structure, on the same order
of magnitude as the structures seen at EIT (Figure 6c). Moreover, the
intermediate branch exhibits overlapping sheets of polymer stretch(especially
on the bottom side of the channel at the time instant shown) that strongly
resemble those seen at EIT.

Figures 6d, e and f show the structures at $\textit{Wi}=45$, close to the
point where the intermediate branch seems to turns around to merge with VNTSA.
The fluctuations on the intermediate branch (Figure 6e) decrease in magnitude
on increasing Wi, and at $\textit{Wi}=45$ become comparable to those on VNTSA
branch as shown in Figure 6d. Furthermore, structure on this branch goes from
displaying overlapping sheets of polymer stretch at low Wi to localized
striations at high Wi similar to those of VNTSA. In sharp contrast to the
localized structures just described, EIT (Figure 6f) continues to display
sheets of polymer stretch which get stronger on increasing Wi. All three
solution branches continue to intermittently display wall normal velocity
fluctuations that resemble the TS mode.

###  3.3 Linear analyses

In this section we elaborate on the linearized problem and its connection to
the attractors described above using linear stability and resolvent analyses.
The spectrum corresponding to disturbances with wavelengths equal to the DNS
box size, i.e. $k=2\pi/5$, has a least-stable eigenvalue at $c\approx
0.32-0.010i$, and the associated eigenfunction is the viscoelastic extension
of the TS mode. For low values of Wi, the mode is less stable than its
Newtonian counterpart, while for $\textit{Wi}\gtrsim 2$, it becomes more
stable with increasing elasticity; this non-monotonic behavior has been
reported by Zhang et al. (2013), who attribute it to viscoelastic modification
of the phase difference between $u$ and $v$. Nevertheless, over the range of
Wi considered here, the eigenvalue varies by less than 1% of the Newtonian
value. The linear stability of the laminar state in this range of Wi, and the
very weak dependence of $c$ on Wi continues up to at least $\Rey=6000$,
confirming the observations in Section 3.2 that finite amplitude disturbances
are required to trigger transition to EIT or the VNTSA. However, linear
instabilities not related to the TS mode have been found in other regions of
parameter space (Garg et al., 2018; Chaudhary et al., 2019), implying the
possibility of different attractor families in those regions.

Figure 5: (a): Fluctuation structure of the VNTSA and (b) power spectral
density (PSD) of $v$ at position $(0,0.825)$ at $\Rey=3000$, $\textit{Wi}=8$.
Frequencies corresponding to the TS mode (red-dashed) and its higher harmonics
(black-dashed) are also shown. (c) and (d): snapshots of the VNTSA structure
for $\Rey=3000$, $\textit{Wi}=20$. Contour plots follow the same format as in
Figure 4.

Figure 6: (a), (b) and (c): Snapshots of VNTSA, intermediate branch and EIT
respectively at $\textit{Wi}=13$. (d), (e) and (f): Snapshots of the same
solution branches at $\textit{Wi}=45$. Deviations from the laminar base state
shown as earlier. Color scales are intentionally centered about 0 for
comparison purposes.

A measure of the relative importance of the conformation tensor and velocity
disturbances is the ratio of the peak amplitudes of $\hat{\alpha}_{xx}$, (the
largest component of the conformation tensor), and $\hat{v}$. This ratio is
shown in Figure 7a. Two distinct regimes are apparent, with the transition
between the two occurring at $\textit{Wi}\approx 3.1$. The low Wi regime
scales as $\textit{Wi}^{2}$, which is the same scaling as in linear shear
flow. The amplitude ratio above the change in slope does not exhibit power law
scaling. The change in slope at $\textit{Wi}\approx 3.1$, can be understood by
examining the $\hat{\alpha}_{xx}$ mode shapes, the magnitudes of which are
plotted in Figure 7b for several values of Wi in the range shown in Figure 7a.
For small Wi, the disturbance is largest at the wall and decays rapidly away
from it. Therefore, the $\textit{Wi}^{2}$ scaling in this regime can be
explained by the fact that the leading-order approximation of the base flow
very near the wall is simple shear. As Wi increases, this value decreases,
while a new local maximum emerges and grows, becoming the global maximum just
above $\textit{Wi}=3$; the arrow in the figure indicates the profile where
this occurs. Upon further increase in Wi, the maximum gradually shifts away
from the wall, and the modes become increasingly localized around the location
of the critical layer $y_{c}$, at which the real part of the wavespeed equals
the base flow velocity. The critical layer for $\textit{Wi}=13$ is indicated
by the vertical dashed line. This suggests that a critical layer mechanism is
responsible for the change in scaling at large Wi, though at present we do not
understand the specific origin of this result. Interestingly, the Wi at which
the VNTSA comes into existence is only slightly larger than that at which the
transition to critical layer scaling occurs.

Figure 7: (a) Solid blue line: ratio of peak amplitudes of $\hat{\alpha}_{xx}$
and $\hat{v}$ for the linear TS mode as a function of Wi; circles: amplitude
ratio for the VNTSA. (b) Magnitude of $\hat{\alpha}_{xx}$ for the linear TS
mode for several values of Wi in the range $[1,\,20]$. Darker lines indicate
higher values of Wi. The thick red line shows the averaged magnitude of
$\hat{\alpha}_{xx}$ from the VNTSA for $\textit{Wi}=13$. For comparison, the
linear TS mode profile for the same Wi is shown in blue, and the vertical
dashed line marks the critical layer location $y_{c}=0.825$. The arrow
indicates the value of Wi corresponding to the arrow in (a). (c) First two
singular values of the resolvent operator for $k$ and $c$ corresponding to the
linear TS mode.

Also shown in Figure 7a is the amplitude ratio computed from the VNTSA for
several values of Wi. Excellent agreement between the linear and nonlinear
results quantitatively reinforces the TS-mode-like nature of the VNTSA.
Additionally, the profile of $|\hat{\alpha}_{xx}|$, averaged in the streamwise
direction and over many snapshots, for the VNTSA at $\textit{Wi}=13$ is shown
by the thick red line in Figure 7b, and the blue line highlights the linear
mode for the same Wi. The VNTSA profile exhibits the same localization, and
the location of the peak value is in close agreement with the critical layer
location.

Figure 7c shows the first two singular values of the resolvent operator for
$k$ and $c$ corresponding to the linear TS mode. Shekar et al. (2019) showed
that such modes are the most-amplified 2D disturbances in this parameter
regime and that the leading response mode is nearly identical to the TS
eigenmode; for this reason the resolvent modes are not plotted separately. The
substantial increase in the leading singular value with Wi indicates that this
amplification becomes much stronger with increasing elasticity, and
consequently that considerably smaller disturbances may be sufficient to
trigger self-sustaining nonlinear mechanisms. Further, the symmetry of the
flow geometry about $y=0$ means that resolvent modes typically come in pairs
having similar amplification, with one mode having a symmetric $\hat{v}$
response and the other having an antisymmetric $\hat{v}$ response. However,
the growing separation between the first and second singular values with
increasing Wi indicates that this pairing is broken by elasticity, and that
the symmetry exhibited by the TS mode is preferred in terms of linear
amplification.

###  3.4 Broader context: flow geometry and dimensionality

The results described above are limited to two-dimensional channel flow. Here
we describe how they may be viewed in the broader context of turbulent drag
reduction, first with regard to how they may relate to the pipe flow geometry
and then in the full context of three-dimensional turbulence.

As noted in the Introduction, elastoinertial turbulence with very similar
features has been observed in both channel and pipe flows. A natural question,
then, is to what extent the above channel flow results are relevant to the
pipe flow case. While it is true of course that there is no linear instability
of Newtonian pipe flow, the general structures of the linear stability problem
for the pipe and channel are very similar. The A, P, and S (wall, center,
strongly decaying) families of modes for the channel flow problem also arise
in the pipe, as illustrated in Figs. 4.18 and 4.19 of Drazin & Reid (2004).
(In channel flow, one of the wall modes, the TS mode, goes unstable.)
Additionally, critical layers are a strong source of linear amplification in
pipe flow, as they are in channel flow (McKeon & Sharma, 2010).

Turning to nonlinear behavior in Newtonian pipe flow, an important difference
from channel flow is the apparent absence of subcritical traveling wave
solutions that are analogous to the nonlinear Tollmien-Schlichting waves
(Patera & Orszag, 1981). That is, in pipe flow there is no nonlinear solution
branch that corresponds to the NNTSA described above. Given this observation,
one might wonder whether the results presented here for viscoelastic channel
flow are relevant for pipe flow.

To address this issue, we make the following remarks. The key observation of
the present work is that the TS mode is nonlinearly excited _by
viscoelasticity_. The resulting solution branch, the VNTSA, is not connected,
at least in this part of parameter space under consideration, to the Newtonian
branch, indicating that the nonlinear mechanism that sustains it is distinct
from the nonlinear mechanisms sustaining the Newtonian branch. So the absence
of a _Newtonian_ mechanism for nonlinear sustainment of Tollmien-Schlichting-
like traveling waves does not imply the absence of a _viscoelastic_ mechanism.
Furthermore, as noted in the Introduction, pipe flow simulations of EIT (Lopez
et al., 2019) display essentially two-dimensional velocity fluctuations
localized near the wall that are similar to those reported in channel flow.
These are precisely what would be expected in pipe flow based on the
observations we report here for channel flow, i.e. a critical layer mechanism
associated with excitation of a Tollmien-Schlichting-like mode. At the same
time, in more strongly viscoelastic regimes, Garg et al. (2018) and Chaudhary
et al. (2020) have found a center mode instability for pipe flow; for the
Oldroyd-B model with $\Rey=3500,\beta=0.9$, instability occurs when
$176.9<\textit{Wi}<4783.6$. A related instability might be present in the
channel flow problem. These results open up the possibility that other states
unrelated to the nonlinear excitation of a wall mode may also play a role at
EIT in both channels and pipes, especially at high Wi.

We now turn to the topic how the present results are related to the fully
three-dimensional context, beginning with a brief overview of Newtonian near-
wall turbulence.

Newtonian turbulence is of course, strongly three-dimensional, with the
dominant near-wall structure comprised of coherent wavy streamwise vortices.
In all of the canonical wall-bounded shear flow geometries (pipe flow, channel
flow, plane shear (Couette) flow), families of three-dimensional nonlinear
traveling wave solutions to the Navier-Stokes equations (NSE) have been
discovered (Waleffe, 1998, 2001, 2003; Wang et al., 2007; Hof et al., 2004;
Eckhardt et al., 2007, 2008; Duguet et al., 2008b; Wedin & Kerswell, 2004).
These solutions are often denoted “exact coherent states” (ECS), and their
predominant structure is very similar to that observed in wall turbulence: a
mean shear and wavy streamwise vortices. (In particular, they bear no
structural resemblance to Tollmien-Schlichting waves and do not arise from a
linear instability of the laminar state.) Related, but more complex states
have been found as well, that are not pure traveling waves but rather
“relative periodic orbits” that are time-periodic modulo a phase shift in one
of the translation-invariant spatial directions (Duguet et al., 2008a). In
minimal domains at Reynolds numbers near transition, the turbulent dynamics
have been found to be organized, at least in part, around these solutions
(see, e.g. Gibson et al. (2008); Kawahara et al. (2012); Park & Graham
(2015)).

It has long been known that one of the effects of viscoelasticity on wall
turbulence is to weaken and broaden the near-wall streamwise vortices (Kim et
al., 2007; White & Mungal, 2008). A number of studies have addressed this
observation by investigating the effect of viscoelasticity on ECS, which as
noted, capture this streamwise vortex structure (Stone et al., 2002; Stone &
Graham, 2003; Stone et al., 2004; Li et al., 2005, 2006; Li & Graham, 2007).
Indeed, the effect of viscoelasticity is to weaken these structures; the
polymer stresses directly counteract the streamwise vortices.

In particular, Li & Graham (2007) studied the bifurcation scenario for a
particular family of channel flow ECS (Waleffe (1998)) in a parameter regime
very close to that considered here. At $\Rey=1500$, this ECS family is
sufficiently weakened by viscoelasticity to lose existence at
$\textit{Wi}\approx 16$, somewhat above the value $\textit{Wi}\approx 10$
beyond which Newtonian turbulence cannot self-sustain in the DNS study of
Shekar et al. (2019). (This discrepancy is consistent with what we know about
transition in the Newtonian case: channel flow turbulence is self-sustaining
above $\Rey\approx 1000$, while ECS can exist in that case down to
$\Rey\approx 660$ (Shekar & Graham, 2018).) Extrapolating slightly from the
results of Li & Graham (2007), one can estimate that at $\Rey=3000$, this ECS
family loses existence at $\textit{Wi}\approx 25$.

Given that, in the present 2D work with $\Rey=3000$, EIT is found at
$\textit{Wi}\gtrsim 13$ and the VNTSA above at $\textit{Wi}\gtrsim 6$, there
would appear to be a regime $6\lesssim\textit{Wi}\lesssim 25$ in which both 2D
and 3D structures may exist and interact. Furthermore, existing results on the
effect of viscoelasticity on ECS are limited to one ECS family, and there
certainly may be others that can persist to higher Wi. Consistent with this
analysis, a number of studies have reported near-wall EIT-like spanwise
oriented structures, with 3D quasistreamwise vortices further away from the
wall (Dubief et al., 2013; Choueiri et al., 2018; Pereira et al., 2019a, b).
How the 2D and 3D structures interact is an important topic for future work.

##  4 Conclusion

This study focuses on two-dimensional plane channel flow of a very dilute
polymer solution at $\Rey=3000$. At sufficiently high Wi, elastoinertial
turbulence is observed in this parameter regime, and the focus of the present
work is to make progress toward understanding the structures and mechanisms
underlying the dynamics in this regime. We report here the existence of a new
attractor that is based on the viscoelastic linear Tollmien-Schlichting mode
and is nonlinearly sustained by viscoelastic stresses. We denote this as the
viscoelastic nonlinear Tollmien-Schlichting attractor (VNTSA). At the
parameters considered here, this solution branch is not connected to the
Newtonian branch of nonlinear self-sustained Tollmien-Schlichting waves; it
would be interesting to learn whether they become connected at higher $\Rey$.
In a domain of dimensionless length $5$, this solution comes into existence at
finite but very small amplitude when $\textit{Wi}\gtrsim 6$, increasing in
amplitude until $\textit{Wi}\approx 49$ where it loses existence again. At
higher Wi, initial conditions corresponding to this solution branch at lower
Wi evolve into elastoinertial turbulence. In general, we do not find pure
nonlinear traveling waves, but until Wi is large, the nonperiodic fluctuations
are very small. The connection of the VNTSA to the linear TS mode is
established via their strong structural similarities, including a quantitive
agreement between the relative magnitudes of the velocity and stress
fluctuations. The value of Wi at which the VNTSA comes into existence is close
to where the relative amplitude of the stress and velocity fluctuations for
the linear TS mode undergoes a change in scaling. Above this transition the
stress fluctuations become highly localized at the position of the critical
layer.

Taken together, these results suggest that, at least in the parameter range
considered here, the bypass transition leading to EIT is mediated by nonlinear
amplification and self-sustenance of perturbations that excite the Tollmien-
Schlichting mode. Gaining an understanding of the mechanism underlying this
phenomenon will shed light on the origin of elastoinertial turbulence.

## Acknowledgments

This work was supported by (UW) NSF CBET-1510291, AFOSR FA9550-18-1-0174, and
ONR N00014-18-1-2865, and (Caltech) ONR N00014-17-1-3022.

Declaration of interests: none.

## References

  * Chaudhary et al. (2019) Chaudhary, Indresh, Garg, Piyush, Shankar, V. & Subramanian, Ganesh 2019 Elasto-inertial wall mode instabilities in viscoelastic plane Poiseuille flow. Journal of Fluid Mechanics 881, 119–163.
  * Chaudhary et al. (2020) Chaudhary, Indresh, Garg, Piyush, Subramanian, Ganesh & Shankar, Viswanathan 2020 Linear instability of viscoelastic pipe flow. ArXiv , arXiv: 2003.09369v1.
  * Choueiri et al. (2018) Choueiri, George H, Lopez, Jose M & Hof, Björn 2018 Exceeding the Asymptotic Limit of Polymer Drag Reduction. Phys. Rev. Lett. 120 (12), 124501.
  * Dallas et al. (2010) Dallas, V, Vassilicos, J & Hewitt, G 2010 Strong polymer-turbulence interactions in viscoelastic turbulent channel flow. Phys. Rev. E 82 (6), 066303.
  * Drazin & Reid (2004) Drazin, P. G. & Reid, W. H. 2004 Hydrodynamic Stability, 2nd edn. Cambridge Mathematical Libraries . Cambridge University Press.
  * Dubief et al. (2013) Dubief, Yves, Terrapon, Vincent E & Soria, Julio 2013 On the mechanism of elasto-inertial turbulence. Phys. Fluids 25 (11), 110817.
  * Duguet et al. (2008a) Duguet, Yohann, Pringle, Chris C T & Kerswell, Rich R 2008a Relative periodic orbits in transitional pipe flow. Phys. Fluids 20 (11), 114102.
  * Duguet et al. (2008b) Duguet, Y, Willis, A P & Kerswell, R R 2008b Transition in pipe flow: the saddle structure on the boundary of turbulence. J. Fluid Mech. 613, 255–274.
  * Eckhardt et al. (2008) Eckhardt, Bruno, Faisst, Holger, Schmiegel, Armin & Schneider, Tobias M 2008 Dynamical systems and the transition to turbulence in linearly stable shear flows. Philosophical Transactions Of The Royal Society A-Mathematical Physical And Engineering Sciences 366 (1868), 1297–1315.
  * Eckhardt et al. (2007) Eckhardt, Bruno, Schneider, Tobias M, Hof, Bjorn & Westerweel, Jerry 2007 Turbulence transition in pipe flow. Annu Rev Fluid Mech 39, 447–468.
  * Garg et al. (2018) Garg, Piyush, Chaudhary, Indresh, Khalid, Mohammad, Shankar, V & Subramanian, Ganesh 2018 Viscoelastic pipe flow is linearly unstable. Phys. Rev. Lett. 121 (2), 024502\.
  * Gibson et al. (2008) Gibson, J F, Halcrow, J & Cvitanović, P 2008 Visualizing the geometry of state space in plane Couette flow. J. Fluid Mech. 611, 24.
  * Graham (2014) Graham, Michael D 2014 Drag reduction and the dynamics of turbulence in simple and complex fluids. Phys. Fluids 26, 101301\.
  * Hameduddin et al. (2019) Hameduddin, Ismail, Gayme, Dennice F & Zaki, Tamer A 2019 Perturbative expansions of the conformation tensor in viscoelastic flows. J. Fluid Mech. 858, 377–406.
  * Hof et al. (2004) Hof, Björn, van Doorne, Casimir WH, Westerweel, Jerry, Nieuwstadt, Frans TM, Faisst, Holger, Eckhardt, Bruno, Wedin, Håkan, Kerswell, Richard R & Waleffe, Fabian 2004 Experimental observation of nonlinear traveling waves in turbulent pipe flow. Science 305 (5690), 1594–1598.
  * Jiménez (1990) Jiménez, Javier 1990 Transition to turbulence in two-dimensional Poiseuille flow. J. Fluid Mech. 218, 265–297.
  * Kawahara et al. (2012) Kawahara, Genta, Uhlmann, Markus & van Veen, Lennaert 2012 The Significance of Simple Invariant Solutions in Turbulent Flows. Annu Rev Fluid Mech 44, 203–225.
  * Kim et al. (2007) Kim, Kyoungyoun, Li, Chang-F, Sureshkumar, R, Balachandar, S & Adrian, Ronald J 2007 Effects of polymer stresses on eddy structures in drag-reduced turbulent channel flow. J. Fluid Mech. 584, 281–299.
  * Kurganov & Tadmor (2000) Kurganov, Alexander & Tadmor, Eitan 2000 New high-resolution central schemes for nonlinear conservation laws and convection–diffusion equations. J. Comput. Phys. 160 (1), 241–282.
  * Lee & Zaki (2017) Lee, Sang Jin & Zaki, Tamer A 2017 Simulations of natural transition in viscoelastic channel flow. J. Fluid Mech. 820, 232–262.
  * Li & Graham (2007) Li, Wei & Graham, Michael D 2007 Polymer induced drag reduction in exact coherent structures of plane Poiseuille flow. Phys. Fluids 19 (8), 083101.
  * Li et al. (2005) Li, W, Stone, P A & Graham, Michael D 2005 Viscoelastic nonlinear traveling waves and drag reduction in plane Poiseuille flow. Fluid Mechanics and its Applications: Proceedings of the IUTAM Symposium on Laminar-Turbulent Transition and Finite Amplitute Solutions 77, 285–308.
  * Li et al. (2006) Li, Wei, Xi, Li & Graham, Michael D 2006 Nonlinear travelling waves as a framework for understanding turbulent drag reduction. J. Fluid Mech. 565, 353–362.
  * Lopez et al. (2019) Lopez, Jose M, Choueiri, George H & Hof, Björn 2019 Dynamics of viscoelastic pipe flow at low reynolds numbers in the maximum drag reduction limit. Journal of Fluid Mechanics 874, 699–719.
  * McKeon & Sharma (2010) McKeon, B J & Sharma, A S 2010 A critical-layer framework for turbulent pipe flow. J. Fluid Mech. 658, 336–382.
  * Min et al. (2003) Min, Taegee, Choi, Haecheon & Yoo, Jung Yul 2003 Maximum drag reduction in a turbulent channel flow by polymer additives. J. Fluid Mech. 492, 91–100.
  * Page & Zaki (2015) Page, Jacob & Zaki, Tamer A 2015 The dynamics of spanwise vorticity perturbations in homogeneous viscoelastic shear flow. J. Fluid Mech. 777, 327–363.
  * Park & Graham (2015) Park, Jae Sung & Graham, Michael D 2015 Exact coherent states and connections to turbulent dynamics in minimal channel flow. J. Fluid Mech. 782, 430–454.
  * Patera & Orszag (1981) Patera, Anthony T & Orszag, Steven A 1981 Finite-amplitude stability of axisymmetric pipe flow. J. Fluid Mech. 112, 467–474.
  * Pereira et al. (2019a) Pereira, Anselmo, Thompson, Roney L & Mompean, Gilmar 2019a Beyond the maximum drag reduction asymptote: the pseudo-laminar state. arXiv preprint arXiv:1911.00439 .
  * Pereira et al. (2019b) Pereira, Anselmo S, Thompson, Roney L & Mompean, Gilmar 2019b Common features between the newtonian laminar–turbulent transition and the viscoelastic drag-reducing turbulence. Journal of Fluid Mechanics 877, 405–428.
  * Samanta et al. (2013) Samanta, Devranjan, Dubief, Yves, Holzner, Markus, Schäfer, Christof, Morozov, Alexander N, Wagner, Christian & Hof, Björn 2013 Elasto-inertial turbulence. Proc. Nat. Acad. Sci. 110 (26), 10557–10562.
  * Schmid (2007) Schmid, P. J. 2007 Nonmodal stability theory. Annu. Rev. Fluid Mech. 39, 129–162.
  * Shekar & Graham (2018) Shekar, Ashwin & Graham, Michael D 2018 Exact coherent states with hairpin-like vortex structure in channel flow. J. Fluid Mech. 849, 76–89.
  * Shekar et al. (2019) Shekar, Ashwin, McMullen, Ryan M, Wang, Sung-Ning, McKeon, Beverley J & Graham, Michael D 2019 Critical-Layer Structures and Mechanisms in Elastoinertial Turbulence. Phys. Rev. Lett. 122 (12), 124503.
  * Sid et al. (2018) Sid, S, Terrapon, V E & Dubief, Y 2018 Two-dimensional dynamics of elasto-inertial turbulence and its role in polymer drag reduction. Phys. Rev. Fluids 3 (1), 011301\.
  * Stone & Graham (2003) Stone, PA & Graham, Michael D 2003 Polymer dynamics in a model of the turbulent buffer layer. Phys. Fluids 15 (5), 1247–1256.
  * Stone et al. (2002) Stone, PA, Waleffe, F & Graham, Michael D 2002 Toward a structural understanding of turbulent drag reduction: Nonlinear coherent states in viscoelastic shear flows. Phys. Rev. Lett. 89 (20), 208301.
  * Stone et al. (2004) Stone, Philip A, Roy, Anshuman, Larson, Ronald G, Waleffe, Fabian & Graham, Michael D 2004 Polymer drag reduction in exact coherent structures of plane shear flow. Phys. Fluids 16 (9), 3470–3482.
  * Vaithianathan et al. (2006) Vaithianathan, T, Robert, Ashish, Brasseur, James G & Collins, Lance R 2006 An improved algorithm for simulating three-dimensional, viscoelastic turbulence. J. Non-Newtonian Fluid Mech. 140 (1-3), 3–22.
  * Virk (1975) Virk, P S 1975 Drag reduction fundamentals. AIChE Journal 21 (4), 625–656.
  * Waleffe (1998) Waleffe, F 1998 Three-dimensional coherent states in plane shear flows. Phys. Rev. Lett. 81 (19), 4140–4143.
  * Waleffe (2001) Waleffe, F 2001 Exact coherent structures in channel flow. J. Fluid Mech. 435, 93–102.
  * Waleffe (2003) Waleffe, F 2003 Homotopy of exact coherent structures in plane shear flows. Phys. Fluids 15 (6), 1517–1534.
  * Wang et al. (2007) Wang, Jue, Gibson, John & Waleffe, Fabian 2007 Lower branch coherent states in shear flows: Transition and control. Phys. Rev. Lett. 98 (20), 204501.
  * Wedin & Kerswell (2004) Wedin, H & Kerswell, R R 2004 Exact coherent structures in pipe flow: travelling wave solutions. J. Fluid Mech. 508, 333–371.
  * White & Mungal (2008) White, Christopher M & Mungal, M Godfrey 2008 Mechanics and prediction of turbulent drag reduction with polymer additives. Annu Rev Fluid Mech 40, 235–256.
  * Zammert & Eckhardt (2014) Zammert, Stefan & Eckhardt, Bruno 2014 Streamwise and doubly-localised periodic orbits in plane Poiseuille flow. J. Fluid Mech. 761, 348–359.
  * Zhang et al. (2013) Zhang, Mengqi, Lashgari, Iman, Zaki, Tamer A & Brandt, Luca 2013 Linear stability analysis of channel flow of viscoelastic Oldroyd-B and FENE-P fluids. J. Fluid Mech. 737, 249–279.


# Semileptonic decays $B_{c}^{+}\to
D^{(*)}_{(s)}(l^{+}\nu_{l},l^{+}l^{-},\nu\bar{\nu})$ in the perturbative QCD
approach

Wen-Fei<EMAIL_ADDRESS>Xin Yu1, Cai-Dian Lü1, and Zhen-Jun Xiao2 1
Center for Future High Energy Physics, Institute of High Energy Physics,
Chinese Academy of Sciences, Beijing 100049, People’s Republic of China,
2 Department of Physics and Institute of Theoretical Physics,
Nanjing Normal University, Nanjing, Jiangsu 210023, People’s Republic of China

###### Abstract

In this paper we study the semileptonic decays of $B_{c}^{+}\to
D^{(*)}_{(s)}(l^{+}\nu_{l},l^{+}l^{-},\nu\bar{\nu})$ (here $l$ stands for $e$,
$\mu$, or $\tau$). After evaluating the $B_{c}^{+}\to(D_{(s)},D^{*}_{(s)})$
transition form factors $F_{0,+,T}(q^{2})$ and
$V(q^{2}),A_{0,1,2}(q^{2}),T_{1,2,3}(q^{2})$ by employing the perturbative QCD
factorization approach, we calculate the branching ratios for all these
semileptonic decays. Our predictions for the values of the $B_{c}^{+}\to
D_{(s)}$ and $B_{c}^{+}\to D^{*}_{(s)}$ transition form factors are consistent
with those obtained by using other methods. The branching ratios of the decay
modes with $\bar{\nu}\nu$ are almost an order of magnitude larger than the
corresponding decays with $l^{+}l^{-}$ after the summation over the three
neutrino generations. The branching ratios for the decays with $b\to d$
transitions are much smaller than those decays with the $b\to s$ transitions,
due to the Cabibbo-Kobayashi-Maskawa suppression. We define ratios $R_{D}$ and
$R_{D^{*}}$ for the branching ratios with the $\tau$ lepton versus $\mu$, $e$
lepton final states to cancel the uncertainties of the form factors, which
could possibly be tested in the near future.

###### pacs:

13.20.He, 12.38.Bx, 14.40.Nd

##  I Introduction

The $B_{c}$ meson is a pseudoscalar ground state of $b$ and $c$ quarks, and
thus the electromagnetic interaction cannot transform the $B_{c}$ meson into
other hadrons containing $b$ and $c$ quarks. The two difference of quark
flavors forbid its annihilation into gluons and being below the $B-D$
threshold makes the $B_{c}$ meson stable for strong interaction. The $B_{c}$
meson can only decay through weak interactions, so it is an ideal system to
study weak decays of heavy quarks. Either the heavy quark ($b$ or $c$) can
decay individually, which makes it different from the $B_{u,d}$ or $B_{s}$
meson. The phase space in the $c\to s$ transition is smaller than that in the
$b\to c$ transition, but the Cabibbo-Kobayashi-Maskawa (CKM) matrix element
$|V_{cs}|\sim 1$ is much larger than the CKM matrix element $|V_{cb}|\sim
0.04$. Thus the $c$-quark decays provide the dominant contribution (about
$70\%$) to the decay width of the $B_{c}$ meson PAN67-1559 . Because the mass
of a $B_{c}\bar{B}_{c}$ pair exceeds the threshold of $\Upsilon(4S)$, the
$B_{c}$ meson cannot be produced at the $B$ factories. So comparing with
$B_{u,d}$ or $B_{s}$ meson, the $B_{c}$ meson decays received much less
experimental attention in the past decades. However, at LHC experiments,
around $5\times 10^{10}$ $B_{c}$ events per year are expected pap-LHC events ;
PAN67-1559  due to the relatively large production cross section, which
provides a very good platform to study various $B_{c}$ meson decay modes.

Because there is only one hadronic final product, the $B_{c}$ meson
semileptonic decays among the abundant decay modes are relatively clean in the
theoretical treatment. These semileptonic decays provide good opportunities to
measure not only the CKM matrix elements, such as $|V_{cb}|,|V_{ub}|$, and
$|V_{cd}|$, but also the form factors of the $B_{c}$ to bottom and charmed
mesons transitions. The rare semileptonic decays governed by the flavor-
changing neutral currents are forbidden at tree level in the standard model
(SM). Those decays, which are very sensitive to the contributions of new
intermediate particles or interactions are especially interesting. There are
various approaches working on the semileptonic $B_{c}$ decays. In Ref.
arxiv0903.2234 , for example, Dhir and Verma presented a detailed analysis of
the exclusive semileptonic $B_{c}$ decays in the Bauer-Stech-Wirbel framework.
The authors of the Refs. epjcd4-18 ; prd65-094037 ; arxiv1006-4231  studied
the semileptonic $B_{c}$ decays in the relativistic and/or constituent quark
model. In Refs. M-Jamil-Aslam ; M-J-Aslam , $B_{c}\to D^{*}_{s}l^{+}l^{-}$
decays were studied in the SM with the fourth-generation and supersymmetric
models. The three point QCD sum rules approach was adopted to investigate the
$B_{c}^{+}\to D^{*+}_{(s)}l^{+}l^{-}$ in prd77-114024  and $B_{c}^{+}\to
D_{(s)}^{+}(l^{+}l^{-},\bar{\nu}\nu)$ in prd78-036005 .

In this paper, we will study the semileptonic decays of $B_{c}^{+}\to
D^{(*)}_{(s)}(l^{+}\nu_{l},l^{+}l^{-},\nu\bar{\nu})$ (here $l$ stands for
leptons $e,\mu,$ or $\tau$) the perturbative QCD (pQCD) approach pap-pQCD .
These semileptonic decays are governed by the form factors. At the maximum
recoil region, the final state meson is collinear with a large momentum. he
spectator $c$ quark in $B_{c}$ meson thus needs a hard gluon to kick it off
from almost zero momentum to a collinear state. However, when doing
integrations of momentum fractions of valence quarks, endpoint singularity
occurs. A natural way to kill this singularity is to pick up the neglected
transverse momentum in the collinear factorization. With the additional
transverse momentum cale $k_{T}$ pqcd , double logarithms appear in the
calculation. We have to use the renormalization group equation to perform the
resummation resulting in the so-called Sudakov form factors Suda-factor  and
make the perturbative calculation of the hard amplitudes (form factors)
infrared safe. The pQCD approach is widely adopted to calculate the transition
form factors of $B_{u,d}$ and $B_{s}$ mesonprd86-114025 ; prd87-097501 ; pQCD-
FF . Furthermore, various $B_{c}$ decay modes have also been studied in Refs.
cpc37-093102 ; PRD81-014022-EPJC45-711-EPJC60-107  in the pQCD approach.

The structure of this paper is as follows. After this Introduction, we collect
the distribution amplitudes of the $B_{c}$, $D^{(*)}$ and $D_{s}^{(*)}$ mesons
in Sec. II. Based on the $k_{\rm T}$ factorization formalism, we calculate and
present the expressions for the $B_{c}\to(D^{(*)},D_{s}^{(*)})$ transition
form factors in the large recoil regions in Section III. The numerical results
and relevant discussions are given in Sec. IV. And Sec. V contains a short
summary.

##  II Kinematics and the wave functions

Figure 1: The leading-order Feynman diagrams for the transition of
$B_{c}^{+}\to(D^{(*)},D_{s}^{(*)})$, where $M$ stands for a $D^{(*)}$ or
$D_{s}^{(*)}$ meson, and $\otimes$ is the weak vertex.

The lowest-order diagrams for $B_{c}^{+}\to(D^{(*)},D_{s}^{(*)})$ transitions
are displayed in Fig. 1, where $M$ stands for $D^{(*)}$ or $D_{s}^{(*)}$
meson, the $\otimes$ is the weak vertex for the leptonic pairs to come out. In
the rest frame of $B_{c}$ meson, with the $m_{B_{c}}$ standing for the mass of
the $B_{c}$ meson, and $m$ for the $D_{(s)}$ or $D^{*}_{(s)}$ mesons, the
momenta of $B_{c}$ and $D^{(*)}_{(s)}$ mesons are defined in the light-cone
coordinates as prd67-054028 ; cpc37-093102

$\displaystyle p_{1}=\frac{m_{B_{c}}}{\sqrt{2}}(1,1,0_{\bot}),\quad
p_{2}=\frac{m_{B_{c}}}{\sqrt{2}}(r\eta^{+},r\eta^{-},0_{\bot}),$ (1)

with $r=m/m_{B_{c}}$ and $\eta^{\pm}=\eta\pm\sqrt{\eta^{2}-1}$. As for the
$\eta$ in $\eta^{\pm}$, the expression

$\displaystyle\eta=\frac{1}{2r}\left[1+r^{2}-\frac{q^{2}}{m^{2}_{B_{c}}}\right]$
(2)

can be evaluated from $q^{2}=(p_{1}-p_{2})^{2}$ which is the invariant mass of
the lepton pairs. The momenta of the spectator quarks in the $B_{c}$ and
$D^{(*)}_{(s)}$ mesons are parameterized as

$\displaystyle
k_{1}=(x_{1}\frac{m_{B_{c}}}{\sqrt{2}},x_{1}\frac{m_{B_{c}}}{\sqrt{2}},k_{1\bot}),\quad
k_{2}=(x_{2}\frac{m_{B_{c}}}{\sqrt{2}}r\eta^{+},x_{2}\frac{m_{B_{c}}}{\sqrt{2}}r\eta^{-},k_{2\bot}).$
(3)

For the $D^{*}_{(s)}$ mesons, we define their polarization vector $\epsilon$
as

$\displaystyle\epsilon_{L}=\frac{1}{\sqrt{2}}(\eta^{+},-\eta^{-},0_{\bot}),\quad\epsilon_{T}=(0,0,1),$
(4)

where $\epsilon_{L}$ and $\epsilon_{T}$ denote the longitudinal and transverse
polarization of the $D^{*}_{(s)}$ mesons, respectively.

In this work, we use the same distribution amplitude for the $B_{c}$ meson as
that used in Refs. cpc37-093102 ; epjc45-711 ; prd81-014002 ; prd87-074027 ,

$\displaystyle\Phi_{B_{c}}(x)=\frac{i}{\sqrt{2N_{c}}}[(p/+m_{B_{c}})\gamma_{5}\phi_{B_{c}}(x)]_{\alpha\beta}\;,$
(5)

with

$\displaystyle\phi_{B_{c}}(x)=\frac{f_{B_{c}}}{2\sqrt{2N_{c}}}\delta(x-m_{c}/m_{B_{c}})\exp[-\omega^{2}_{B_{c}}b^{2}/2]\;,$
(6)

where $m_{c}$ is the mass of $c$-quark. Because the $B_{c}$ meson consists of
two heavy quarks $b$ and $c$, just like a heavy quarkonium, the non-
relativistic QCD framework can be applied, which means the leading-order wave
function should be just the zero-point wave function shown in Eq. (6).

For the $D^{(*)}_{(s)}$ mesons, up to twist-3 accuracy, the two-parton light-
cone distribution amplitudes are defined as prd67-054028 ; prd78-014018

$\displaystyle\langle D_{(s)}(p)|q_{\alpha}(z)\bar{c}_{\beta}(0)|0\rangle$
$\displaystyle=$
$\displaystyle\frac{i}{\sqrt{2N_{C}}}\int_{0}^{1}dxe^{ixp\cdot
z}\left[\gamma_{5}\left(p/+m\right)\phi_{D_{(s)}}(x,b)\right]_{\alpha\beta}\;,$
$\displaystyle\langle D_{(s)}^{*}(p)|q_{\alpha}(z)\bar{c}_{\beta}(0)|0\rangle$
$\displaystyle=$
$\displaystyle-\frac{1}{\sqrt{2N_{C}}}\int_{0}^{1}dxe^{ixp\cdot
z}\big{[}\epsilon/_{L}(p/+m)\phi_{D^{*}_{(s)}}(x,b)$ (7)
$\displaystyle+\;\epsilon/_{T}(p/+m)\phi_{D^{*}_{(s)}}(x,b)\big{]}_{\alpha\beta}\;,$

with

$\displaystyle\int_{0}^{1}dx\phi_{D_{(s)}}(x,0)$ $\displaystyle=$
$\displaystyle\frac{f_{D_{(s)}}}{2\sqrt{2N_{c}}}\;,$
$\displaystyle\int_{0}^{1}dx\phi_{D^{*}_{(s)}}(x,0)$ $\displaystyle=$
$\displaystyle\frac{f_{D^{*}_{(s)}}}{2\sqrt{2N_{c}}}\;,$ (8)

as the normalization conditions. We adopt $f_{D}=206.7\pm 8.9$ MeV and
$f_{D_{s}}=260.0\pm 5.6$ MeV in PDG pdg2012  by experimental average for $D$
and $D_{s}$ mesons, respectively. For the $D^{*}$ or $D^{*}_{s}$ meson, we
adopt the same decay constant and distribution amplitude for the longitudinal
and transverse components. Since there is no experimental data, we use
$f_{D^{*}}=270$ MeV and $f_{D^{*}_{s}}=310$ MeV for $D^{*}$ and $D^{*}_{s}$
meson considering of the results in Refs. decay-cons  and assume a $10\%$
uncertainty. The distribution amplitude for the $D_{(s)}$ meson is

$\displaystyle\phi_{D_{(s)}}=\frac{1}{2\sqrt{2N_{C}}}f_{D_{(s)}}6x(1-x)\left[1+C_{D_{(s)}}(1-2x)\right]\exp\left[-\frac{\omega_{D_{(s)}}^{2}b^{2}}{2}\right],$
(9)

which is a $k_{T}$-dependent form with $C_{D}=0.5,\omega_{D}=0.1$ and
$C_{D_{s}}=0.4,\omega_{D_{s}}=0.2$ for $D$ and $D_{s}$ mesons, respectively
prd78-014018 . In this work, we also adopt the same distribution amplitude for
both the vector meson $D^{*}_{(s)}$ and pseudoscalar meson $D_{(s)}$ because
of their small mass difference prd78-014018 .

##  III Form factors of semileptonic decays

The form factors $F_{+}(q^{2}),F_{0}(q^{2})$ for the $B_{c}$ to pseudoscalar
meson $D_{(s)}$ transition induced by the vector current can be defined as
zpc29-637 ; npb592-3

$\displaystyle\langle
D_{(s)}(p_{2})|\bar{q}(0)\gamma_{\mu}b(0)|B_{c}(p_{1})\rangle$
$\displaystyle=$
$\displaystyle\left[(p_{1}+p_{2})_{\mu}-\frac{m_{B_{c}}^{2}-m^{2}}{q^{2}}q_{\mu}\right]F_{+}(q^{2})$
(10) $\displaystyle+$
$\displaystyle\frac{m_{B_{c}}^{2}-m^{2}}{q^{2}}q_{\mu}F_{0}(q^{2}),$

where $q=p_{1}-p_{2}$ is the momentum of the lepton pairs. In order to cancel
the poles at $q^{2}=0$, $F_{+}(0)$ should be equal to $F_{0}(0)$. For the sake
of convenience, we define the auxiliary form factors $f_{1}(q^{2})$ and
$f_{2}(q^{2})$,

$\displaystyle\langle
D_{(s)}(p_{2})|\bar{q}(0)\gamma_{\mu}b(0)|B_{c}(p_{1})\rangle=f_{1}(q^{2})p_{1\mu}+f_{2}(q^{2})p_{2\mu}.$
(11)

In terms of $f_{1}(q^{2})$ and $f_{2}(q^{2})$ the form factors $F_{+}(q^{2})$
and $F_{0}(q^{2})$ are

$\displaystyle F_{+}(q^{2})$ $\displaystyle=$
$\displaystyle\frac{1}{2}\left[f_{1}(q^{2})+f_{2}(q^{2})\right],$
$\displaystyle F_{0}(q^{2})$ $\displaystyle=$
$\displaystyle\frac{1}{2}f_{1}(q^{2})\left[1+\frac{q^{2}}{m_{B_{c}}^{2}-m^{2}}\right]+\frac{1}{2}f_{2}(q^{2})\left[1-\frac{q^{2}}{m_{B_{c}}^{2}-m^{2}}\right].$
(12)

The form factor $F_{T}(q^{2})$ for the $B_{c}\to D_{(s)}$ transition induced
by the tensor current can be defined as npb592-3

$\displaystyle\langle
D_{(s)}(p_{2})|\bar{q}(0)\sigma_{\mu\nu}b(0)|B_{c}(p_{1})\rangle=i\left[p_{2\mu}q_{\nu}-q_{\mu}p_{2\nu}\right]\frac{2F_{T}(q^{2})}{m_{B_{c}}+m}\;.$
(13)

There are seven form factors $V(q^{2}),A_{0,1,2}(q^{2})$ and
$T_{1,2,3}(q^{2})$ needed for the transition of $B_{c}\to D^{*}_{(s)}$ in this
work. The form factors $V(q^{2})$ and $A_{0,1,2}(q^{2})$ are defined by
npb592-3 ; epjc28-515 ; prd65-014007

$\displaystyle\langle
D^{*}_{(s)}(p_{2})|\bar{q}(0)\gamma_{\mu}b(0)|B_{c}(p_{1})\rangle$
$\displaystyle=$
$\displaystyle\epsilon_{\mu\nu\alpha\beta}\epsilon^{\nu*}p_{1}^{\alpha}p_{2}^{\beta}\frac{2V(q^{2})}{m_{B_{c}}+m}\;,$
(14) $\displaystyle\langle
D^{*}_{(s)}(p_{2})|\bar{q}(0)\gamma_{\mu}\gamma_{5}b(0)|B_{c}(p_{1})\rangle$
$\displaystyle=$ $\displaystyle
i\left[\epsilon_{\mu}^{*}-\frac{\epsilon^{*}\cdot
q}{q^{2}}q_{\mu}\right](m_{B_{c}}+m)A_{1}(q^{2})$ (15) $\displaystyle-$
$\displaystyle
i\left[(p_{1}+p_{2})_{\mu}-\frac{m_{B_{c}}^{2}-m^{2}}{q^{2}}q_{\mu}\right](\epsilon^{*}\cdot
q)\frac{A_{2}(q^{2})}{m_{B_{c}}+m}$ $\displaystyle+$ $\displaystyle
i\frac{2m(\epsilon^{*}\cdot q)}{q^{2}}q_{\mu}A_{0}(q^{2}),$

where $\epsilon^{*}$ is the polarization vector of the $D^{*}_{(s)}$ meson.
The form factors $T_{1,2,3}$ are defined by npb592-3 ; prd53-3672

$\displaystyle\langle D^{*}_{(s)}(p_{2})|\bar{q}(0)\sigma_{\mu\nu}q^{\nu}(1$
$\displaystyle+$
$\displaystyle\gamma_{5})b(0)|B_{c}(p_{1})\rangle=i\epsilon_{\mu\nu\alpha\beta}\epsilon^{*\nu}p_{1}^{\alpha}p_{2}^{\beta}2T_{1}(q^{2})$
(16) $\displaystyle+$
$\displaystyle\left[\epsilon^{*}_{\mu}(m^{2}_{B_{c}}-m^{2})-(\epsilon^{*}\cdot
q)(p_{1}+p_{2})_{\mu}\right]T_{2}(q^{2})$ $\displaystyle+$
$\displaystyle(\epsilon^{*}\cdot
q)\left[q_{\mu}-\frac{q^{2}}{m^{2}_{B_{c}}-m^{2}}(p_{1}+p_{2})_{\mu}\right]T_{3}(q^{2})\;,$

with $T_{1}(0)=T_{2}(0)$ implied by the identity

$\displaystyle\sigma_{\mu\nu}\gamma_{5}=-\frac{i}{2}\epsilon_{\mu\nu\alpha\beta}\sigma^{\alpha\beta}\;.$
(17)

In the transverse configuration $b$-space and by including the Sudakov form
factors and the threshold resummation effects, we obtain the $B_{c}\to
D_{(s)}$ form factors $f_{1}(q^{2}),f_{2}(q^{2})$ and $F_{T}(q^{2})$ as
follows

$\displaystyle f_{1}(q^{2})$ $\displaystyle=$ $\displaystyle 16\pi
m_{B_{c}}^{2}rC_{F}\int dx_{1}dx_{2}\int
b_{1}db_{1}b_{2}db_{2}\;\phi_{B_{c}}(x_{1})\phi_{D_{(s)}}(x_{2},b_{2})$ (18)
$\displaystyle\times$ $\displaystyle\Bigl{\\{}\left[1-rx_{2}\right]\cdot
h_{1}(x_{1},x_{2},b_{1},b_{2})\cdot\alpha_{s}(t_{1})\exp\left[-S_{ab}(t_{1})\right]$
$\displaystyle-$ $\displaystyle\left[r+2x_{1}(1-\eta)\right]\cdot
h_{2}(x_{1},x_{2},b_{1},b_{2})\cdot\alpha_{s}(t_{2})\exp\left[-S_{ab}(t_{2})\right]\Bigr{\\}}\;,$
$\displaystyle f_{2}(q^{2})$ $\displaystyle=$ $\displaystyle 16\pi
m_{B_{c}}^{2}C_{F}\int dx_{1}dx_{2}\int
b_{1}db_{1}b_{2}db_{2}\;\phi_{B_{c}}(x_{1})\phi_{D_{(s)}}(x_{2},b_{2})$ (19)
$\displaystyle\times$
$\displaystyle\Bigl{\\{}\left[1-2rx_{2}(1-\eta)\right]\cdot
h_{1}(x_{1},x_{2},b_{1},b_{2})\cdot\alpha_{s}(t_{1})\exp\left[-S_{ab}(t_{1})\right]$
$\displaystyle+$ $\displaystyle\left[2r-x_{1}\right]\cdot
h_{2}(x_{1},x_{2},b_{1},b_{2})\cdot\alpha_{s}(t_{2})\exp\left[-S_{ab}(t_{2})\right]\Bigr{\\}}\;,$
$\displaystyle F_{T}(q^{2})$ $\displaystyle=$ $\displaystyle 8\pi
m_{B_{c}}^{2}C_{F}(1+r)\int dx_{1}dx_{2}\int
b_{1}db_{1}b_{2}db_{2}\;\phi_{B_{c}}(x_{1})\phi_{D_{(s)}}(x_{2},b_{2})$ (20)
$\displaystyle\times$ $\displaystyle\Bigl{\\{}\left[1-rx_{2}\right]\cdot
h_{1}(x_{1},x_{2},b_{1},b_{2})\cdot\alpha_{s}(t_{1})\exp\left[-S_{ab}(t_{1})\right]$
$\displaystyle+$ $\displaystyle\left[2r-x_{1}\right]\cdot
h_{2}(x_{1},x_{2},b_{1},b_{2})\cdot\alpha_{s}(t_{2})\exp\left[-S_{ab}(t_{2})\right]\Bigr{\\}}\;,$

where $C_{F}=4/3$ is the color factor. The functions $h_{1}$ and $h_{2}$, the
scales $t_{1}$, $t_{2}$ and the Sudakov factors $S_{ab}$ are the same as those
given in Refs. prd67-054028 ; cpc37-093102 .

The expressions of form factors $V(q^{2}),A_{0,1,2}(q^{2})$ and
$T_{1,2,3}(q^{2})$ for the $B_{c}\to D^{*}_{(s)}$ transition in the pQCD
approach are:

$\displaystyle V(q^{2})$ $\displaystyle=$ $\displaystyle 8\pi
m_{B_{c}}^{2}C_{F}(1+r)\int dx_{1}dx_{2}\int
b_{1}db_{1}b_{2}db_{2}\;\phi_{B_{c}}(x_{1})\phi^{T}_{D^{*}_{(s)}}(x_{2},b_{2})$
(21) $\displaystyle\times$ $\displaystyle\Bigl{\\{}\left[1-rx_{2}\right]\cdot
h_{1}(x_{1},x_{2},b_{1},b_{2})\cdot\alpha_{s}(t_{1})\exp\left[-S_{ab}(t_{1})\right]$
$\displaystyle+$ $\displaystyle r\cdot
h_{2}(x_{1},x_{2},b_{1},b_{2})\cdot\alpha_{s}(t_{2})\exp\left[-S_{ab}(t_{2})\right]\Bigr{\\}}\;,$
$\displaystyle A_{0}(q^{2})$ $\displaystyle=$ $\displaystyle 8\pi
m_{B_{c}}^{2}C_{F}\int dx_{1}dx_{2}\int
b_{1}db_{1}b_{2}db_{2}\;\phi_{B_{c}}(x_{1})\phi^{L}_{D^{*}_{(s)}}(x_{2},b_{2})$
(22) $\displaystyle\times$
$\displaystyle\Bigl{\\{}\left[1-rx_{2}(r-2\eta)+r\left(1-2x_{2}\right)\right]$
$\displaystyle\times$ $\displaystyle
h_{1}(x_{1},x_{2},b_{1},b_{2})\cdot\alpha_{s}(t_{1})\exp\left[-S_{ab}(t_{1})\right]$
$\displaystyle+$ $\displaystyle\left[r^{2}+x_{1}(1-2r\eta)\right]\cdot
h_{2}(x_{1},x_{2},b_{1},b_{2})\cdot\alpha_{s}(t_{2})\exp\left[-S_{ab}(t_{2})\right]\Bigr{\\}}\;,$
$\displaystyle A_{1}(q^{2})$ $\displaystyle=$ $\displaystyle 16\pi
m_{B_{c}}^{2}C_{F}\frac{r}{1+r}\int dx_{1}dx_{2}\int
b_{1}db_{1}b_{2}db_{2}\;\phi_{B_{c}}(x_{1})\phi^{T}_{D^{*}_{(s)}}(x_{2},b_{2})$
(23) $\displaystyle\times$
$\displaystyle\Bigl{\\{}\left[1+rx_{2}\eta-2rx_{2}+\eta\right]\cdot
h_{1}(x_{1},x_{2},b_{1},b_{2})\cdot\alpha_{s}(t_{1})\exp[-S_{ab}(t_{1})]$
$\displaystyle+$ $\displaystyle\left[r\eta-x_{1}\right]\cdot
h_{2}(x_{1},x_{2},b_{1},b_{2})\cdot\alpha_{s}(t_{2})\exp[-S_{ab}(t_{2})]\Bigr{\\}}\;,$
$\displaystyle A_{2}(q^{2})$ $\displaystyle=$
$\displaystyle\frac{(1+r)^{2}(\eta-r)}{2r(\eta^{2}-1)}\cdot A_{1}(q^{2})-8\pi
m_{B_{c}}^{2}C_{F}\frac{1+r}{\eta^{2}-1}\int dx_{1}dx_{2}\int
b_{1}db_{1}b_{2}db_{2}\;\phi_{B_{c}}(x_{1})$ (24) $\displaystyle\times$
$\displaystyle\phi^{L}_{D^{*}_{(s)}}(x_{2},b_{2})\cdot\Bigl{\\{}\left[\eta(1-r^{2}x_{2})-rx_{2}(1-2\eta^{2}-2r)+(1-r)-r\eta(1+2x_{2})\right]$
$\displaystyle\times$ $\displaystyle
h_{1}(x_{1},x_{2},b_{1},b_{2})\cdot\alpha_{s}(t_{1})\exp\left[-S_{ab}(t_{1})\right]$
$\displaystyle+$
$\displaystyle\left[r(1-x_{1}+2x_{1}\eta^{2})-\eta(r^{2}+x_{1})\right]$
$\displaystyle\times$ $\displaystyle
h_{2}(x_{1},x_{2},b_{1},b_{2})\cdot\alpha_{s}(t_{2})\exp\left[-S_{ab}(t_{2})\right]\Bigr{\\}}\;,$
$\displaystyle T_{1}(q^{2})$ $\displaystyle=$ $\displaystyle 8\pi
m_{B_{c}}^{2}C_{F}\int dx_{1}dx_{2}\int
b_{1}db_{1}b_{2}db_{2}\;\phi_{B_{c}}(x_{1})\phi^{T}_{D^{*}_{(s)}}(x_{2},b_{2})$
(25) $\displaystyle\times$
$\displaystyle\Bigl{\\{}\left[1+r(1-x_{2}(2+r-2\eta))\right]\cdot
h_{1}(x_{1},x_{2},b_{1},b_{2})\cdot\alpha_{s}(t_{1})\exp\left[-S_{ab}(t_{1})\right]$
$\displaystyle+$ $\displaystyle r\left[1-x_{1}\right]\cdot
h_{2}(x_{1},x_{2},b_{1},b_{2})\cdot\alpha_{s}(t_{2})\exp\left[-S_{ab}(t_{2})\right]\Bigr{\\}}\;,$
$\displaystyle T_{2}(q^{2})$ $\displaystyle=$ $\displaystyle 16\pi
m_{B_{c}}^{2}C_{F}\frac{r}{1-r^{2}}\int dx_{1}dx_{2}\int
b_{1}db_{1}b_{2}db_{2}\;\phi_{B_{c}}(x_{1})\phi^{T}_{D^{*}_{(s)}}(x_{2},b_{2})$
(26) $\displaystyle\times$
$\displaystyle\Bigl{\\{}\left[(1-r)(1+\eta)+2rx_{2}(r-\eta)+rx_{2}(2\eta^{2}-r\eta-1)\right]$
$\displaystyle\times$ $\displaystyle
h_{1}(x_{1},x_{2},b_{1},b_{2})\cdot\alpha_{s}(t_{1})\exp\left[-S_{ab}(t_{1})\right]$
$\displaystyle+$ $\displaystyle\left[r(1+x_{1})\eta-r^{2}-x_{1}\right]\cdot
h_{2}(x_{1},x_{2},b_{1},b_{2})\cdot\alpha_{s}(t_{2})\exp\left[-S_{ab}(t_{2})\right]\Bigr{\\}}\;,$
$\displaystyle T_{3}(q^{2})$ $\displaystyle=$
$\displaystyle\frac{r+\eta}{r}\cdot\frac{1-r^{2}}{2(\eta^{2}-1)}\cdot
T_{2}(q^{2})-\frac{1-r^{2}}{(\eta^{2}-1)}$ (27) $\displaystyle\times$
$\displaystyle 8\pi m_{B_{c}}^{2}C_{F}\int dx_{1}dx_{2}\int
b_{1}db_{1}b_{2}db_{2}\;\phi_{B_{c}}(x_{1})\phi^{L}_{D^{*}_{(s)}}(x_{2},b_{2})$
$\displaystyle\times$
$\displaystyle\Bigl{\\{}\left[1+rx_{2}(\eta-2)+\eta\right]\cdot
h_{1}(x_{1},x_{2},b_{1},b_{2})\cdot\alpha_{s}(t_{1})\exp\left[-S_{ab}(t_{1})\right]$
$\displaystyle+$ $\displaystyle\left[x_{1}\eta-r\right]\cdot
h_{2}(x_{1},x_{2},b_{1},b_{2})\cdot\alpha_{s}(t_{2})\exp\left[-S_{ab}(t_{2})\right]\Bigr{\\}}\;.$

One should note that the expressions for the form factors $f_{1,2}(q^{2})$,
$F_{T}(q^{2})$, $V(q^{2})$, $A_{0,1,2}(q^{2})$ and $T_{1,2,3}(q^{2})$ given in
Eqs. (18)-(27) are the results at leading order of the pQCD approach. The
next-to-leading-order contributions to the form factors of
$B\to(\pi,K,\eta^{(\prime)})$ transitions given in Refs. prd85-074004 ;
prd86-114025 ; prd87-097501  are not available here because of the large mass
of $c$-quark and $(D_{(s)},D^{*}_{(s)})$ mesons.

One should note that the pQCD predictions for the considered form factors are
reliable only for the small values of $q^{2}$. For the form factors in the
large-$q^{2}$ region, one has to make an extrapolation for them from the
low-$q^{2}$ region to large-$q^{2}$ region. In this work we make the
extrapolation by using the formula in Refs. cpc37-093102 ; prd79-054012

$\displaystyle F(q^{2})=F(0)\cdot\exp{\left[a\cdot
q^{2}+b\cdot(q^{2})^{2}\right]},$ (28)

where $F$ stands for the form factors $F_{0,+,T},V,A_{0,1,2}$ and $T_{1,2,3}$,
and $a,b$ are the constants to be determined by the fitting procedure.

The $B_{c}^{-}\to\bar{D}^{0}l^{-}\bar{\nu}_{l}$ and
$B_{c}^{-}\to\bar{D}^{*0}l^{-}\bar{\nu}_{l}$ decays are from the quark level
$b\to ul^{-}\bar{\nu}$ charged current transition. The effective Hamiltonian
for such transition is eff-ham

$\displaystyle{\cal H}_{eff}(b\to
ul^{-}\bar{\nu}_{l})=\frac{G_{F}}{\sqrt{2}}V_{ub}\;\bar{u}\gamma_{\mu}(1-\gamma_{5})b\cdot\bar{l}\gamma^{\mu}(1-\gamma_{5})\nu_{l},$
(29)

where $G_{F}=1.16637\times 10^{-5}GeV^{-2}$ is the Fermi-coupling constant and
$V_{ub}$ is one of the CKM matrix elements. With the form factors calculated
in Eqs. (18,19,21-24), one can easily get the differential decay width
expression for $B_{c}^{-}\to\bar{D}^{0}l^{-}\bar{\nu}_{l}$ and
$B_{c}^{-}\to\bar{D}^{*0}l^{-}\bar{\nu}_{l}$.

For those flavor-changing neutral-current one-loop decay modes, such as
$B_{c}\to D^{(*)}l^{+}l^{-}$ and $B_{c}\to D_{s}^{(*)}l^{+}l^{-}$, are
transitions of $b\to dl^{+}l^{-}$ and $b\to sl^{+}l^{-}$ at quark level,
respectively. The effective Hamiltonians and the corresponding differential
decay widths are more complicated, we refer the readers to Refs. prd86-114025
; prd53-3672 ; npb612-25 ; epjc40-565 ; prd61-074024 .

For the decay modes of $B_{c}\to D^{(*)}_{s}\nu\bar{\nu}$, the effective
Hamiltonian is eff-ham

$\displaystyle{\cal H}_{eff}(b\to s\nu\bar{\nu})$ $\displaystyle=$
$\displaystyle\frac{G_{F}}{\sqrt{2}}\frac{\alpha_{em}}{2\pi\sin^{2}(\theta_{W})}V_{tb}V_{ts}^{*}\eta_{X}X(x_{t})\;\left[{\bar{s}}\gamma^{\mu}(1-\gamma_{5})b\right]\left[{\bar{\nu}}\gamma_{\mu}(1-\gamma_{5})\nu\right]$
(30)

where $\theta_{W}$ is the Weinberg angle with $\sin^{2}(\theta_{W})=0.231$
pdg2012 , $V_{tb}$ and $V_{ts}$ are CKM matrix elements and
$\alpha_{em}\approx 1/137$ is the fine structure constant. The function
$X(x_{t})$ can be found in Ref. eff-ham , while $\eta_{X}\approx 1$ is the QCD
radiative correction factor eff-ham . As for the decay modes of $B_{c}\to
D^{(*)}\nu\bar{\nu}$, their effective Hamiltonian can be obtained by a simple
replacement of $s\to d$ in Eq. (30). The corresponding differential decay
widths for $B_{c}\to D_{(s)}\nu\bar{\nu}$ is the same as
$B\to\pi(K)\nu\bar{\nu}$ in Ref. prd86-114025  except the replacements
$m_{B}\to m_{B_{c}}$ and $m_{P}\to m$. While for the decay modes of $B_{c}\to
D^{*}_{s}\nu\bar{\nu}$, the differential decay width is jhep0707-072

$\displaystyle\frac{d\Gamma(B_{c}\to D^{*}_{s}\nu\bar{\nu})}{dq^{2}}$
$\displaystyle=$
$\displaystyle\frac{G^{2}_{F}\alpha^{2}_{em}}{2^{10}\pi^{5}m^{3}_{B_{c}}}\cdot\left|\frac{X(x_{t})}{sin^{2}(\theta_{w})}\right|^{2}\cdot\eta^{2}_{X}\cdot\left|V_{tb}V^{*}_{ts}\right|^{2}\lambda^{\frac{1}{2}}\bigg{\\{}8\lambda
q^{2}\frac{V^{2}}{(m_{B_{c}}+m)^{2}}$ (31)
$\displaystyle+\frac{1}{m^{2}}\bigg{[}\lambda^{2}\frac{A^{2}_{2}}{(m_{B_{c}}+m)^{2}}+(m_{B_{c}}+m)^{2}(\lambda+12m^{2}q^{2})\cdot
A^{2}_{1}$ $\displaystyle-2\lambda(m^{2}_{B_{c}}-m^{2}-q^{2})\cdot
Re[A^{*}_{1}A_{2}]\bigg{]}\bigg{\\}}\;,$

where $V,A_{1}$ and $A_{2}$ are the form factors of $B_{c}\to D^{*}_{s}$
transition, and the phase-space factor

$\displaystyle\lambda=(m^{2}_{B_{c}}+m^{2}-q^{2})^{2}-4m^{2}_{B_{c}}m^{2}\;.$
(32)

##  IV Numerical results and discussions

In the numerical calculations we adopt the following input parameters pdg2012

$\displaystyle m_{B^{-}_{c}}$ $\displaystyle=$ $\displaystyle 6.277~{}{\rm
GeV},\quad m_{\bar{D}^{0}}=1.865~{}{\rm GeV},\quad m_{D^{-}}=1.870~{}{\rm
GeV},$ $\displaystyle m_{\bar{D}^{*0}}$ $\displaystyle=$ $\displaystyle
2.007~{}{\rm GeV},\quad m_{D^{*-}}=2.010~{}{\rm GeV},\quad
m_{D_{s}^{-}}=1.969~{}{\rm GeV},$ $\displaystyle m_{\bar{D}_{s}^{*-}}$
$\displaystyle=$ $\displaystyle 2.112~{}{\rm GeV},\quad m_{\tau}=1.777~{}{\rm
GeV},\quad m_{c}=1.275\pm 0.025~{}{\rm GeV},$ $\displaystyle\tau_{B_{c}}$
$\displaystyle=$ $\displaystyle(0.45\pm 0.04)~{}{\rm ps},$ (33)

For the CKM matrix element $V_{ub}$, we adopt the value in Refs. prd86-114025
; babar-Vub . And we use $|V_{tb}|=0.999,|V_{ts}|=0.0404$ and
$|V_{td}/V_{ts}|=0.211$ pdg2012  in this work. As for the decay constant of
the $B_{c}$ meson, we adopt $0.489~{}{\rm GeV}$ plb651-171  as its central
value, and give it an uncertainty of $0.050~{}{\rm GeV}$.

The numerical values of the $B_{c}\to D$ and $B_{c}\to D_{s}$ transition form
factors $F_{0,+,T}$ at $q^{2}=0$ and their fitted parameters $a,b$ are listed
in Table 1. The numerical values of the form factors $V,A_{0,1,2}$ and
$T_{1,2,3}$ at $q^{2}=0$ for the $B_{c}\to D^{*}$ and $B_{c}\to D^{*}_{s}$
transitions are listed in Table 2. The first error of the pQCD predictions for
the form factors in Table 1 and Table 2 is induced by the $B_{c}$ meson wave
function parameter $\omega_{B_{c}}=1.0\pm 0.1$; the second error comes from
the uncertainty of decay constant $f_{B_{c}}$; the third error comes from the
uncertainty of decay constants of the $D^{(*)}_{(s)}$ mesons; the fourth error
in Tables 1 and 2 comes from the uncertainty of $D_{(s)}^{(*)}$ wave function
$C_{D^{(*)}}=0.5\pm 0.1$ or $C_{D^{(*)}_{s}}=0.4\pm 0.1$; the fifth error
comes from $m_{c}=1.275\pm 0.025$ GeV. The errors from the uncertainty of
$\omega_{D^{(*)}}=0.10\pm 0.02$ or $\omega_{D^{(*)}_{s}}=0.20\pm 0.04$ are
very small that have been neglected.

Unlike the form factors at maximum recoil, the extrapolation parameters $a$,
$b$ of the form factors are less sensitive to the decay constant and wave
function of $D_{(s)}^{(*)}$ meson. In Tables 1 and 2, we only show
uncertainties for the parameter $a$ and $b$ from $B_{c}$ meson wave function
parameter $\omega_{B_{c}}$, and from quark mass uncertainty $m_{c}=1.275\pm
0.025$ GeV. As a comparison, we also present some results obtained by other
authors based on different methods in Table 3. It is easy to see that our
results are consistent with the results in literature.

Table 1: The pQCD predictions for form factors $F_{0},F_{+}$ and $F_{T}$ at $q^{2}=0$ and the parametrization constants $a$ and $b$ for $B_{c}\to D$ and $B_{c}\to D_{s}$ transitions. | $\ \ F(0)$ | $\ \ a$ | $\ \ b$
---|---|---|---
$F_{0}^{B_{c}\to D}$ |  $0.19\pm 0.02\pm 0.02\pm 0.01\pm 0.01\pm 0.01$ |  $0.038\pm 0.001\pm 0.000$ |  $0.0013\pm 0.0001$
$F_{+}^{B_{c}\to D}$ |  $0.19\pm 0.02\pm 0.02\pm 0.01\pm 0.01\pm 0.01$ |  $0.059\pm 0.001\pm 0.001$ |  $0.0020\pm 0.0001$
$F_{T}^{B_{c}\to D}$ |  $0.20\pm 0.02\pm 0.02\pm 0.01\pm 0.01\pm 0.01$ |  $0.070\pm 0.001\pm 0.001$ |  $0.0021^{+0.0000}_{-0.0001}$
$F_{0}^{B_{c}\to D_{s}}$ |  $0.27\pm 0.03\pm 0.03\pm 0.02\pm 0.01\pm 0.01$ |  $0.039\pm 0.002\pm 0.001$ |  $0.0015^{+0.0001}_{-0.0000}$
$F_{+}^{B_{c}\to D_{s}}$ |  $0.27\pm 0.03\pm 0.03\pm 0.02\pm 0.01\pm 0.01$ |  $0.061\pm 0.002\pm 0.001$ |  $0.0023^{+0.0001}_{-0.0000}$
$F_{T}^{B_{c}\to D_{s}}$ |  $0.28\pm 0.03\pm 0.03\pm 0.02\pm 0.01\pm 0.01$ |  $0.073\pm 0.002\pm 0.001$ |  $0.0025^{+0.0000}_{-0.0001}$
Table 2: The pQCD predictions for form factors $A_{0,1,2},V$ and $T_{1,2,3}$ at $q^{2}=0$ and the parametrization constants $a$ and $b$ for $B_{c}\to D^{*}$ and $B_{c}\to D^{*}_{s}$ transitions. | $\ F(0)$ | $\ a$ | $\ b$
---|---|---|---
$A_{0}^{B_{c}\to D^{*}}$ |  $0.17\pm 0.02\pm 0.02\pm 0.02\pm 0.01\pm 0.00$ |  $0.063\pm 0.001\pm 0.001$ |  $0.0024\pm 0.0000\pm 0.0000$
$A_{1}^{B_{c}\to D^{*}}$ |  $0.18\pm 0.02\pm 0.02\pm 0.02\pm 0.01\pm 0.01$ |  $0.043\pm 0.001\pm 0.001$ |  $0.0018\pm 0.0001\pm 0.0001$
$A_{2}^{B_{c}\to D^{*}}$ |  $0.20\pm 0.02\pm 0.02\pm 0.02\pm 0.01\pm 0.01$ |  $0.067\pm 0.001\pm 0.001$ |  $0.0026\pm 0.0001\pm 0.0001$
$V^{B_{c}\to D^{*}}$ |  $0.25\pm 0.03\pm 0.03\pm 0.03\pm 0.01\pm 0.01$ |  $0.073\pm 0.002\pm 0.001$ |  $0.0029\pm 0.0001\pm 0.0001$
$T_{1}^{B_{c}\to D^{*}}$ |  $0.22\pm 0.02\pm 0.02\pm 0.02\pm 0.01\pm 0.01$ |  $0.063\pm 0.001\pm 0.001$ |  $0.0027\pm 0.0001\pm 0.0001$
$T_{2}^{B_{c}\to D^{*}}$ |  $0.22\pm 0.02\pm 0.02\pm 0.02\pm 0.01\pm 0.01$ |  $0.038\pm 0.001\pm 0.001$ |  $0.0017\pm 0.0001\pm 0.0001$
$T_{3}^{B_{c}\to D^{*}}$ |  $0.20\pm 0.02\pm 0.02\pm 0.02\pm 0.01\pm 0.01$ |  $0.077\pm 0.002\pm 0.001$ |  $0.0049\pm 0.0001\pm 0.0001$
$A_{0}^{B_{c}\to D^{*}_{s}}$ |  $0.21\pm 0.02\pm 0.02\pm 0.02\pm 0.01\pm 0.01$ |  $0.064\pm 0.001\pm 0.001$ |  $0.0031\pm 0.0002\pm 0.0001$
$A_{1}^{B_{c}\to D^{*}_{s}}$ |  $0.23\pm 0.02\pm 0.02\pm 0.02\pm 0.01\pm 0.01$ |  $0.044\pm 0.002\pm 0.001$ |  $0.0022\pm 0.0002\pm 0.0001$
$A_{2}^{B_{c}\to D^{*}_{s}}$ |  $0.25\pm 0.03\pm 0.03\pm 0.03\pm 0.01\pm 0.01$ |  $0.069\pm 0.002\pm 0.001$ |  $0.0035\pm 0.0002^{+0.0002}_{-0.0001}$
$V^{B_{c}\to D^{*}_{s}}$ |  $0.33\pm 0.03\pm 0.03\pm 0.03\pm 0.02\pm 0.01$ |  $0.075\pm 0.002\pm 0.001$ |  $0.0039\pm 0.0002\pm 0.0001$
$T_{1}^{B_{c}\to D^{*}_{s}}$ |  $0.28\pm 0.03\pm 0.03\pm 0.03\pm 0.02\pm 0.01$ |  $0.064\pm 0.001\pm 0.001$ |  $0.0035\pm 0.0002\pm 0.0001$
$T_{2}^{B_{c}\to D^{*}_{s}}$ |  $0.28\pm 0.03\pm 0.03\pm 0.03\pm 0.02\pm 0.01$ |  $0.039\pm 0.001\pm 0.001$ |  $0.0023\pm 0.0002\pm 0.0001$
$T_{3}^{B_{c}\to D^{*}_{s}}$ |  $0.27\pm 0.03\pm 0.03\pm 0.03\pm 0.02\pm 0.01$ |  $0.082\pm 0.002\pm 0.001$ |  $0.0068\pm 0.0002\pm 0.0002$
Table 3: Comparison of $B_{c}\to D^{(*)}_{(s)}$ transition form factors at $q^{2}=0$ evaluated in this paper with other methods. $B_{c}\to D^{(*)}$ | $F_{+}(0)=F_{0}(0)$ |  $F_{T}(0)$ |  $A_{0}(0)$ |  $A_{1}(0)$ |  $A_{2}(0)$ |  $V(0)$ | $T_{1}(0)=T_{2}(0)$ |  $T_{3}(0)$
---|---|---|---|---|---|---|---|---
pQCD | $0.19$ | $0.20$ | $0.17$ | $0.18$ | $0.20$ | $0.25$ | $0.22$ | $0.20$
Ref.prd79-054012  | $0.16$ | $-$ | $0.09$ | $0.08$ | $0.07$ | $0.13$ | $-$ | $-$
Ref.prd68-094020  | $0.14$ | $-$ | $0.14$ | $0.17$ | $0.19$ | $0.18$ | $-$ | $-$
Ref.prd63-074010  | $0.189$ | $-$ | $0.284$ | $0.146$ | $0.158$ | $0.296$ | $-$ | $-$
Ref.epjc51-833  | $0.35$ | $-$ | $0.05$ | $0.32$ | $0.57$ | $0.57$ | $-$ | $-$
Ref.hepph-0211021  | $0.32$ | $-$ | $0.35$ | $0.43$ | $0.51$ | $1.66$ | $-$ | $-$
Ref.prd79-304004  | $0.075$ | $-$ | $0.081$ | $0.095$ | $0.11$ | $0.16$ | $-$ | $-$
$B_{c}\to D_{s}^{(*)}$ | $F_{+}(0)=F_{0}(0)$ |  $F_{T}(0)$ |  $A_{0}(0)$ |  $A_{1}(0)$ |  $A_{2}(0)$ |  $V(0)$ | $T_{1}(0)=T_{2}(0)$ |  $T_{3}(0)$
pQCD | $0.27$ | $0.28$ | $0.21$ | $0.23$ | $0.25$ | $0.33$ | $0.28$ | $0.27$
Ref.prd79-054012  | $0.28$ | $-$ | $0.17$ | $0.14$ | $0.12$ | $0.23$ | $-$ | $-$
Ref.hepph-0211021  | $0.45$ | $-$ | $0.47$ | $0.56$ | $0.65$ | $2.02$ | $-$ | $-$
Ref.prd79-304004  | $0.15$ | $-$ | $0.16$ | $0.18$ | $0.20$ | $0.29$ | $-$ | $-$

With the form factors given, it is straightforward to calculate the branching
ratios for all the considered semileptonic decays by performing the numerical
integration over the whole range of $q^{2}$. For the $b\to u$ charged current
process, with $l=(e,\mu)$, the decay rates are the following:

$\displaystyle Br(B^{-}_{c}\to\bar{D}^{0}l^{-}\bar{\nu}_{l})$ $\displaystyle=$
$\displaystyle(3.15^{+0.97}_{-0.72}(\omega_{B_{c}})^{+0.68}_{-0.61}(f_{B_{c}})^{+0.29}_{-0.27}(m_{c})^{+0.31}_{-0.29}(C_{D})^{+0.28}_{-0.27}(f_{D})\pm
0.28(\tau_{B_{c}}))\cdot 10^{-5},$ $\displaystyle
Br(B^{-}_{c}\to\bar{D}^{0}\tau^{-}\bar{\nu}_{\tau})$ $\displaystyle=$
$\displaystyle(2.16^{+0.72}_{-0.52}(\omega_{B_{c}})^{+0.46}_{-0.42}(f_{B_{c}})^{+0.22}_{-0.19}(m_{c})^{+0.20}_{-0.19}(C_{D})^{+0.19}_{-0.18}(f_{D})\pm
0.19(\tau_{B_{c}}))\cdot 10^{-5},$ $\displaystyle
Br(B^{-}_{c}\to\bar{D}^{*0}l^{-}\bar{\nu}_{l})$ $\displaystyle=$
$\displaystyle(1.09^{+0.34}_{-0.26}(\omega_{B_{c}})^{+0.23}_{-0.21}(f_{B_{c}})^{+0.13}_{-0.11}(m_{c})\pm
0.10(C_{D^{*}})^{+0.23}_{-0.21}(f_{D^{*}})\pm 0.10(\tau_{B_{c}}))\cdot
10^{-4},$ $\displaystyle Br(B^{-}_{c}\to\bar{D}^{*0}\tau^{-}\bar{\nu}_{\tau})$
$\displaystyle=$
$\displaystyle(0.64^{+0.20}_{-0.16}(\omega_{B_{c}})^{+0.14}_{-0.12}(f_{B_{c}})^{+0.08}_{-0.07}(m_{c})^{+0.06}_{-0.05}(C_{D^{*}})^{+0.13}_{-0.12}(f_{D^{*}})\pm
0.06(\tau_{B_{c}}))\cdot 10^{-4},$ (34)

where the errors come from the uncertainties of $\omega_{B_{c}}=1.0\pm 0.1$,
$f_{B_{c}}=0.489\pm 0.050$ GeV, $m_{c}=(1.275\pm 0.025)$ GeV,
$C_{D^{(*)}}=0.5\pm 0.1$, $f_{D}=(206.7\pm 8.9)$ MeV or $f_{D^{*}}=(270\pm
27)$ MeV and $\tau_{B_{c}}=(0.45\pm 0.04)$ ps, respectively.

For the flavor-changing neutral-current processes, after making the numerical
integration over the whole range of $4m_{l}^{2}\leq
q^{2}\leq(m_{B_{c}}-m)^{2}$, we get the pQCD predictions for the branching
ratios of considered decay modes which are listed in Table 4. The errors of
the pQCD predictions in Table 4 come from the uncertainties of
$\omega_{B_{c}}$, $m_{c}$, $C_{D^{(*)}}$ or $C_{D^{(*)}_{s}}$, $f_{D^{(*)}}$
or $f_{D_{s}^{(*)}}$ and $\tau_{B_{c}}$, respectively.

Table 4: The pQCD predictions for the branching ratios of the considered decays ($l=e,\mu$). Decay modes | pQCD predictions
---|---
$Br(B^{-}_{c}\to D^{-}l^{+}l^{-})$ | $(3.79^{+1.16}_{-0.86}(\omega_{B_{c}})^{+0.81}_{-0.74}(f_{B_{c}})^{+0.35}_{-0.32}(m_{c})^{+0.37}_{-0.35}(C_{D})^{+0.33}_{-0.32}(f_{D})\pm 0.34(\tau_{B_{c}}))\cdot 10^{-9}$
$Br(B^{-}_{c}\to D^{-}\tau^{+}\tau^{-})$ | $(1.03^{+0.38}_{-0.27}(\omega_{B_{c}})^{+0.22}_{-0.20}(f_{B_{c}})^{+0.12}_{-0.10}(m_{c})^{+0.09}_{-0.08}(C_{D})\pm 0.09(f_{D})\pm 0.09(\tau_{B_{c}}))\cdot 10^{-9}$
$Br(B^{-}_{c}\to D^{-}\bar{\nu}\nu)$ | $(3.13^{+0.96}_{-0.71}(\omega_{B_{c}})^{+0.67}_{-0.61}(f_{B_{c}})^{+0.30}_{-0.26}(m_{c})^{+0.31}_{-0.29}(C_{D})^{+0.28}_{-0.26}(f_{D})\pm 0.28(\tau_{B_{c}}))\cdot 10^{-8}$
$Br(B^{-}_{c}\to D_{s}^{-}l^{+}l^{-})$ | $(1.56^{+0.46}_{-0.36}(\omega_{B_{c}})^{+0.33}_{-0.30}(f_{B_{c}})^{+0.17}_{-0.15}(m_{c})^{+0.13}_{-0.12}(C_{D_{s}})\pm 0.07(f_{D_{s}})\pm 0.14(\tau_{B_{c}}))\cdot 10^{-7}$
$Br(B^{-}_{c}\to D_{s}^{-}\tau^{+}\tau^{-})$ | $(0.38^{+0.13}_{-0.10}(\omega_{B_{c}})^{+0.08}_{-0.07}(f_{B_{c}})^{+0.05}_{-0.04}(m_{c})\pm 0.03(C_{D_{s}})\pm 0.02(f_{D_{s}})\pm 0.03(\tau_{B_{c}}))\cdot 10^{-7}$
$Br(B^{-}_{c}\to D_{s}^{-}\bar{\nu}\nu)$ | $(1.29^{+0.39}_{-0.30}(\omega_{B_{c}})^{+0.28}_{-0.25}(f_{B_{c}})^{+0.14}_{-0.12}(m_{c})^{+0.11}_{-0.10}(C_{D_{s}})\pm 0.06(f_{D_{s}})\pm 0.11(\tau_{B_{c}}))\cdot 10^{-6}$
$Br(B^{-}_{c}\to D^{*-}l^{+}l^{-})$ | $(1.21^{+0.36}_{-0.28}(\omega_{B_{c}})^{+0.26}_{-0.23}(f_{B_{c}})^{+0.14}_{-0.12}(m_{c})\pm 0.11(C_{D^{*}})^{+0.25}_{-0.23}(f_{D^{*}})\pm 0.11(\tau_{B_{c}}))\cdot 10^{-8}$
$Br(B^{-}_{c}\to D^{*-}\tau^{+}\tau^{-})$ | $(0.16^{+0.05}_{-0.04}(\omega_{B_{c}})\pm 0.03(f_{B_{c}})\pm 0.02(m_{c})\pm 0.01(C_{D^{*}})^{\pm}0.03(f_{D^{*}})\pm 0.01(\tau_{B_{c}}))\cdot 10^{-8}$
$Br(B^{-}_{c}\to D^{*-}\bar{\nu}\nu)$ | $(1.10^{+0.34}_{-0.26}(\omega_{B_{c}})^{+0.24}_{-0.21}(f_{B_{c}})^{+0.13}_{-0.11}(m_{c})\pm 0.10(C_{D^{*}})^{+0.23}_{-0.21}(f_{D^{*}})\pm 0.10(\tau_{B_{c}}))\cdot 10^{-7}$
$Br(B^{-}_{c}\to D_{s}^{*-}l^{+}l^{-})$ | $(4.40^{+1.40}_{-1.05}(\omega_{B_{c}})^{+0.95}_{-0.85}(f_{B_{c}})^{+0.72}_{-0.57}(m_{c})^{+0.32}_{-0.31}(C_{D^{*}_{s}})^{+0.92}_{-0.84}(f_{D^{*}_{s}})\pm 0.39(\tau_{B_{c}}))\cdot 10^{-7}$
$Br(B^{-}_{c}\to D_{s}^{*-}\tau^{+}\tau^{-})$ | $(0.52^{+0.18}_{-0.13}(\omega_{B_{c}})^{+0.11}_{-0.10}(f_{B_{c}})^{+0.10}_{-0.08}(m_{c})\pm 0.03(C_{D^{*}_{s}})^{+0.11}_{-0.10}(f_{D^{*}_{s}})\pm 0.05(\tau_{B_{c}}))\cdot 10^{-7}$
$Br(B^{-}_{c}\to D_{s}^{*-}\bar{\nu}\nu)$ | $(4.04^{+1.30}_{-0.97}(\omega_{B_{c}})^{+0.87}_{-0.78}(f_{B_{c}})^{+0.68}_{-0.53}(m_{c})^{+0.29}_{-0.28}(C_{D^{*}_{s}})^{+0.85}_{-0.77}(f_{D_{s}^{*}})\pm 0.36(\tau_{B_{c}}))\cdot 10^{-6}$

From the pQCD predictions for the form factors $F_{0,+,T}$ in Table 1, the
form factors $V,A_{0,1,2}$ and $T_{1,2,3}$ in Table 2 and the pQCD predictions
for the branching ratios as listed in Eq. (IV) and in Table 4, we have the
following points:

  1. (i)

All the form factors for the transitions $B_{c}\to D^{(*)}_{s}$ are larger
than the corresponding values for the transitions $B_{c}\to D^{(*)}$ at
$q^{2}=0$, which characterizes the SU(3) breaking effect.

  2. (ii)

$F_{0}(0)$ equals to $F_{+}(0)$ by definition for the $B_{c}\to D$ or
$B_{c}\to D_{s}$ transition, but they have different $q^{2}$ dependence by the
different parameters $(a,b)$. $T_{1}(0)$ equals to $T_{2}(0)$ for the
$B_{c}\to D^{*}$ or $B_{c}\to D^{*}_{s}$ transition claimed by the Eq. (17) as
they are given in Table 2 although their expressions are different in Eqs.(25,
26).

  3. (iii)

Because of the phase space suppression, the branching ratios of the decay
modes with a $\tau$ in the final product are smaller than those decay modes
with electron or muon in the final product for the the charged current
process. And for the flavor changing neutral current processes, with two
$\tau$’s in the final product, the branching ratios are much smaller than the
corresponding decays with electron or muon pairs in the final product.

  4. (iv)

The branching ratios of the decay modes with $\bar{\nu}\nu$ are almost an
order magnitude larger than the corresponding decays with $l^{+}l^{-}$ after
the summation over the three neutrino generations. Because of the strong
suppression of the CKM factor $|V_{td}/V_{ts}|^{2}=|0.211|^{2}$ pdg2012 , the
branching ratios for the decay modes with $b\to d$ transitions are much
smaller than those decay modes with the $b\to s$ transitions.

In order to reduce the theoretical uncertainty of the form factor
calculations, we define two ratios $R_{D}$ and $R_{D^{*}}$ among the branching
ratios for the the charged-current processes

$\displaystyle R_{D}$ $\displaystyle=$
$\displaystyle\frac{Br(B^{-}_{c}\to\bar{D}^{0}\tau^{-}\bar{\nu}_{\tau})}{Br(B^{-}_{c}\to\bar{D}^{0}l^{-}\nu_{l})}=0.69\pm
0.01(\omega_{B_{c}})^{+0.01}_{-0.00}(m_{c})\;,$ (35) $\displaystyle R_{D^{*}}$
$\displaystyle=$
$\displaystyle\frac{Br(B^{-}_{c}\to\bar{D}^{*0}\tau^{-}\bar{\nu}_{\tau})}{Br(B^{-}_{c}\to\bar{D}^{*0}l^{-}\nu_{l})}=0.59^{+0.00}_{-0.01}(\omega_{B_{c}})^{+0.00}_{-0.01}(m_{c})\;,$
(36)

with $l=(e,\mu)$. These two relations will be tested by experiments.

##  V Summary

In this paper we studed the $B_{c}\to(D_{(s)},D^{*}_{(s)})$ transition form
factors $F_{0,+,T}(q^{2})$ and $V(q^{2}),A_{0,1,2}(q^{2}),T_{1,2,3}(q^{2})$ in
the pQCD factorization approach based on $k_{T}$ factorization. The pQCD
predictions for the values of the $B_{c}\to D_{(s)}$ and $B_{c}\to
D^{*}_{(s)}$ transition form factors agree with those obtained using other
methods. Utilizing these form factors, we calculated the branching ratios for
all the semileptonic decays of $B_{c}^{+}\to
D^{(*)}_{(s)}(l^{+}\nu_{l},l^{+}l^{-},\nu\bar{\nu})$. Because of phase space
suppression, the production ratios of the decay modes with lepton $\tau$ in
the final product are smaller than the corresponding decays with electron or
muon in the final product. The branching ratios of the decay modes with
$\bar{\nu}\nu$ are almost an order magnitude larger than the corresponding
decays with $l^{+}l^{-}$ after the summation over the three neutrino
generations. The branching ratios for the decays with the $b\to d$ transition
are much smaller than those with the $b\to s$ transitions.

In order to reduce the theoretical uncertainty of the pQCD predictions, we
defined two ratios $R_{D}$ and $R_{D^{*}}$ among the branching ratios for the
the charged-current processes. The pQCD predictions are

$\displaystyle R_{D}$ $\displaystyle=$
$\displaystyle\frac{Br(B^{-}_{c}\to\bar{D}^{0}\tau^{-}\bar{\nu}_{\tau})}{Br(B^{-}_{c}\to\bar{D}^{0}l^{-}\nu_{l})}\approx
0.7\;,$ (37) $\displaystyle R_{D^{*}}$ $\displaystyle=$
$\displaystyle\frac{Br(B^{-}_{c}\to\bar{D}^{*0}\tau^{-}\bar{\nu}_{\tau})}{Br(B^{-}_{c}\to\bar{D}^{*0}l^{-}\nu_{l})}\approx
0.6\;,$ (38)

with $l=(e,\mu)$. It would possible to test these predictions by LHCb and the
forthcoming Super-B experiments.

###### Acknowledgements.

This work is supported in part by the National Natural Science Foundation of
China under Grant No. 11375208, 11228512, and 11235005.

## References

  * (1) I.P. Gouz, V.V. Kiselev, A.K. Likhoded, V.I. Romanovsky, and O.P. Yushchenko, Phys. At. Nucl. 67, 1559 (2004).
  * (2) N. Brambilla et al., (Quarkonium Working Group), Report No. CERN-2005-005; M.P. Altarelli and F. Teubert, Int. J. Mod. Phys. A 23, 5117 (2008).
  * (3) R. Dhir and R.C. Verma, Phys. Scr. 82, 065101 (2010).
  * (4) A. Faessler, Th. Gutsche, M.A. Ivanov, J.G. Körner, and V.E. Lyubovitskij, Eur. Phys. J. direct C 4, 1 (2002).
  * (5) C.Q. Geng, C.W. Hwang, and C.C. Liu, Phys. Rev. D 65, 094037 (2002).
  * (6) D. Ebert, R.N. Faustov, and V.O. Galkin, Phys. Rev. D 82, 034032 (2010).
  * (7) I. Ahmed, M. A. Paracha, M. Junaid, A. Ahmed, A. Rehman and M. J. Aslam, arXiv:1107.5694.
  * (8) A. Ahmed, I. Ahmed, M. A. Paracha, M. Junaid, A. Rehman and M. J. Aslam, arXiv:1108.1058.
  * (9) K. Azizi, F. Falahati, V. Bashiry, and S.M. Zebarjad, Phys. Rev. D 77, 114024 (2008).
  * (10) K. Azizi and R. Khosravi, Phys. Rev. D 78, 036005 (2008).
  * (11) T.W. Yeh and H.-n. Li, Phys. Rev. D 56, 1615 (1997); H.-n. Li and H.L. Yu, Phys. Rev. Lett. 74, 4388 (1995); Phys. Lett. B 353, 301 (1995); Phys. Rev. D 53, 2480 (1996).
  * (12) Y.Y. Keum, H.-n. Li, and A.I. Sanda, Phys Lett. B 504, 6 (2001); Phys. Rev. D 63, 054008 (2001); C. -D. Lu, K. Ukai and M. -Z. Yang, Phys. Rev. D 63, 074009 (2001); C. -D. Lu and M. -Z. Yang, Eur. Phys. J. C 23, 275 (2002).
  * (13) G. Sterman, Phys. Lett. B 179, 281 (1986); Nucl. Phys. B281, 310 (1987); S. Catani and L. Trentadue, Nucl. Phys. B327, 323 (1989); H.-n. Li, Phys. Rev. D 66, 094010 (2002).
  * (14) W.F. Wang and Z.J. Xiao, Phys. Rev. D 86, 114025 (2012).
  * (15) W.F. Wang, Y.Y. Fan, M. Liu, and Z.J. Xiao Phys. Rev. D 87, 097501 (2013).
  * (16) R.H. Li, C.D. Lü, and W. Wang Phys. Rev. D 79, 034014 (2009); R.H. Li, C.D. Lü, W. Wang, and X.X. Wang, Phys. Rev. D 79, 014013 (2009);
  * (17) Z.J. Xiao and X. Liu, Phys. Rev. D 84, 074033 (2011); X. Liu, Z.J. Xiao, and C.D. Lü , Phys. Rev. D 81, 014022 (2010); J.F. Cheng, D.S. Du, and C.D. Lü, Eur. Phys. J. C 45,711 (2006); J.F. Sun, D.S. Du, and Y.L. Yang, Eur. Phys. J. C 60, 107 (2009).
  * (18) W.F. Wang, Y.Y. Fan, and Z.J. Xiao, Chin. Phys. C 37, 093102 (2013).
  * (19) T. Kurimoto, H.-n. Li, and A.I. Sanda, Phys. Rev. D 67, 054028 (2003).
  * (20) J.F. Cheng, D.S. Du, and C.D. Lü, Eur. Phys. J. C 45, 711 (2006).
  * (21) X. Liu, Z.J. xiao, and C.D. Lü, Phys. Rev. D 81, 014022 (2010).
  * (22) Z.T. Zou, X. Yu, and C.D. Lü, Phys. Rev. D 87, 074027 (2013).
  * (23) R.H. Li, C.D. Lü, and H. Zou, Phys. Rev. D 78, 014018 (2008).
  * (24) J. Beringer, et al. (Particle Data Group), Phys. Rev. D 86, 010001 (2012).
  * (25) D. Becirevic, Ph. Boucaud, J.P. Leroy, V. Lubicz, G. Martinelli, F. Mescia, and F. Rapuano, Phys. Rev. D 60, 074501 (1999); M.A. Ivanov, J.G. Korner, P. Santorelli, Phys. Rev. D 73, 054024 (2006); G.L. Wang, Phys. Lett. B 633, 492 (2006); D. Becirevic, V. Lubicz, F. Sanfilippo, S. Simula, and C. Tarantino, J. High Energy Phys. 02 (2012) 042; P. Gelhausen, A. Khodjamirian, A.A. Pivovarov, and D. Rosenthal, Phys. Rev. D 88, 014015 (2013);
  * (26) M. Wirbel, B. Stech, and M. Bauer, Z. Phys. C 29, 637 (1985).
  * (27) M. Beneke and Th. Feldmann, Nucl. Phys. B592, 3 (2001).
  * (28) C.D. Lü and M.Z. Yang, Eur. Phys. J. C 28, 515 (2003).
  * (29) T. Kurimoto, H.-n. Li, and A. I. Sanda, Phys. Rev. D 65, 014007 (2001).
  * (30) P. Colangelo, F. De Fazio, P. Santorelli, and E. Scrimieri, Phys. Rev. D 53, 3672 (1996).
  * (31) H.-n. Li, Y.L. Shen, and Y.M. Wang, Phys. Rev. D 85, 074004 (2012).
  * (32) W. Wang, Y.L. Shen, and C.D. Lü, Phys. Rev. D 79, 054012 (2009).
  * (33) A. J. Buras, M. Misiak, M. Munz, and S. Pokorski, Nucl. Phys. B424, 374 (1994);
G. Buchalla, A. J. Buras, and M. E. Lautenbacher, Rev. Mod. Phys. 68, 1125
(1996).

  * (34) M. Beneke, Th. Feldmann, and D. Seidel, Nucl. Phys. B612, 25 (2001).
  * (35) W.J. Li, Y.B. Dai, and C.S. Huang, Eur. Phys. J. C 40, 565 (2005).
  * (36) A. Ali, P. Ball, L. T. Handoko, and G. Hiller, Phys. Rev. D 61, 074024 (2000).
  * (37) T. M. Aliev, A. S. Cornell, and Naveen Gaur, J. High Energy Phys. 07 (2007) 072;
T. Barakat, New J. Phys. 4, 25 (2002).

  * (38) P. del Amo Sanchez et al. (BABAR Collaboration), Phys. Rev. D 83, 032007 (2011); Phys. Rev. D 83, 052011 (2011).
  * (39) T.W. Chiu, T.H. Hsieh, C.H. Huang, and K. Ogawa (TWQCD Collaboration), Phys. Lett. B 651, 171 (2007).
  * (40) D. Ebert, R. N. Faustov, and V. O. Galkin, Phys. Rev. D 68, 094020 (2003).
  * (41) M. A. Ivanov, J. G. Körner, and O, N, Pakhomova, Phys. Lett. B 555, 189 (2003).
  * (42) T. Huang, and F. Zuo, Eur. Phys. J. C 51, 833 (2007).
  * (43) V.V. Kiselev, arXiv:hep-ph/0211021.
  * (44) R. Dhir and R.C. Verma, Phys. Rev. D 79, 034004 (2009).