# Assessing and Accelerating Coverage in Deep Reinforcement Learning

Arpan Kusari  University of Michigan Transportation Research Institute
University of Michigan
Ann Arbor, Michigan, USA
<EMAIL_ADDRESS>
###### Abstract

Current deep reinforcement learning (DRL) algorithms utilize randomness in
simulation environments to assume complete coverage in the state space.
However, particularly in high dimensions, relying on randomness may lead to
gaps in coverage of the trained DRL neural network model, which in turn may
lead to drastic and often fatal real-world situations. To the best of the
author’s knowledge, the assessment of coverage for DRL is lacking in current
research literature. Therefore, in this paper, a novel measure, Approximate
Pseudo-Coverage (APC), is proposed for assessing the coverage in DRL
applications. We propose to calculate APC by projecting the high dimensional
state space on to a lower dimensional manifold and quantifying the occupied
space. Furthermore, we utilize an exploration-exploitation strategy for
coverage maximization using Rapidly-Exploring Random Tree (RRT). The efficacy
of the assessment and the acceleration of coverage is demonstrated on standard
tasks such as Cartpole, highway-env.

###### Index Terms:

coverage assessment, assurance, dimensionality reduction, coverage
maximization

##  I Introduction

Sequential decision making is an essential component of robotics, utilizing a
sequence of actions (chosen from a set of possible actions) to navigate an
uncertain environment to best achieve a goal. Reinforcement learning [1] has
been proposed to provide a formal approach to solve the sequential decision
making problem where an agent chooses actions so as to maximize the cumulative
rewards. Over the last decade, deep reinforcement learning (DRL), which
combines RL with deep neural networks, has been applied successfully to to
various applications with continuous state space and discrete/continuous
action space. The wide-ranging applications where DRL has been successfully
applied include board games (e.g., Go [2] and Chess [3]), video games (e.g.,
Atari games [4]), and complex robotics control tasks (e.g. [5]).

Currently DRL tasks utilize reward per episode as the common metric for
measuring the performance of the model and completion of the training process.
However, there is an inherent assumption placed on the simulation, which
generates samples for the training, to be able to sample uniformly from the
entire sample space and provide adequate coverage at the end of the training
process. Coverage has been defined as the measure of how much of the Design
and Verification Environment (DVE) has been covered via testing[6]. There is
extensive literature on coverage driven verification (CDV) for functional
verification of hardware designs [7]. Coverage testing for assurance of
machine learning models has been explored by transforming high-dimensional
data into a lower dimensional manifold in [8]. CDV has also been applied
successfully to Deep Neural Networks (DNN) by developing four different
dataset quality measures and measuring the quality of most common test
datasets on popular DNN models [9]. However, to the best of the author’s
knowledge, there has been no research conducted exploring coverage testing in
the context of DRL. Assessment of coverage appears to be an important factor
which can influence the acceptance of the DRL model in a real-world setting.

Coverage testing in traditional development relies on functional requirements
or model-based techniques, which are, by their very nature, deterministic and
provides a tractable solution. Assessment of coverage is complicated in DRL
for three primary reasons:

  * •

Data for training a DRL is generated by the simulation environment during the
training process which can make coverage testing difficult

  * •

There is no explicit model (or functional requirements) in model-free DRL to
rely on in order to do exhaustive coverage testing

  * •

The DRL tasks are often characterized by a high dimensional continuous state
space which suffers from the ”curse of dimensionality” [10] and makes
calculation of coverage problematic

Since DRL tasks involve high-dimensional state space, reducing the
dimensionality of the state space by projecting the state samples to a lower
dimensional manifold makes the computation of coverage tractable. In this
paper, we reduce the state space to two dimensions using t-distributed
Stochastic Neighbor Embedding (t-SNE) [11]. Thus, we get a pseudo measure of
coverage where more is inherently better without equating it to actual
coverage. Even after reducing the dimensionality to two, the simulation
environment generates a lot of samples during the training process and
computing the exact area of these samples is often a complex problem.
Therefore, we implement an over-approximation of the coverage by utilizing an
occupancy grid-like approach [12]. Using a fine grained 2-D grid, we apply a
voting process of whether a particular point falls within a grid cell. We then
calculate the ratio of number of occupied grid-cells (non-zero grids) to the
total number of grid-cells. We refer this measure to Approximate Pseudo-
Coverage (APC).

For coverage maximization, we note the following two observations:

  * •

Currently, the simulation environments for DRL are blind to the training
process and thus, are not adaptive with the changes in the training process.
Our hypothesis is that providing a feedback to the simulation environment
would help in maximizing the coverage.

  * •

The starting state has a large impact on the evolution of the episode [13].
Searching for unexplored areas of the state space can provide a maximization
in coverage. Also, small perturbations of high scoring start states can help
us in maximizing rewards as well.

With these observations in mind, we introduce an exploration-exploitation
strategy for determining the start state for each episode. Using an epsilon-
greedy policy, we explore using a Rapidly Exploring Random Tree (RRT) [14] to
efficiently search high-dimensional non-convex spaces and exploit using a
small random perturbation of highest scoring start state.

The paper is organized as follows: Section II provides a background on deep RL
and coverage testing, Section III details the coverage assessment and
maximization procedure, Section IV presents the results of the proposed
procedure on standard tasks and Section V gives the discussions and
conclusions.

##  II Background

###  II-A Coverage testing

Coverage is defined as a measure of completeness of set of tests for checking
a model of a design. In coverage based verification, the planning of coverage
is performed before the start of the verification process [15]. Each
measurable task is known as a coverage task and together, they form a coverage
model [16]. Pseudorandom generation techniques are used in test generators in
CDV [7]. CDV uses functional requirements or models to create this test suite.
In model-based test generation, a model is explored or traversed to obtain
abstract tests, i.e. tests at the same level of abstraction as the model [17].

There are different kinds of coverage which can be broadly classified into two
categories, normal coverage and pseudo-coverage [6]. These constitute:

  * •

Normal coverage - This refers to exact coverage testing where achieving 100%
explains the completeness of the tests for DVE. There are two kinds of
coverage testing which fall in this category:

    * –

Implementation coverage - measuring which lines of code have been executed or
how many state bits have been toggled

    * –

Functional coverage - user defined metric intended to investigate to which
extent the functionality of a given design under test (DUT) has been verified

  * •

Pseudo-coverage - This type of coverage testing relies on measures where more
coverage implies better without having a exact upper limit. Examples include:

    * –

Novelty coverage - measures if a new sample is sufficiently far from the
current samples

    * –

Criticality coverage - measures the closeness of the samples from the terminal
state

Coverage in the context of machine learning models is not a thoroughly studied
problem. With deep learning algorithms being increasingly used in safety
critical applications, assessing and accelerating coverage for such
applications is of paramount importance.

###  II-B Deep Reinforcement Learning

The central idea behind reinforcement learning (RL) is to learn the sequence
of actions that an agent takes to maximize the cumulative reward. In an RL
task, for any given state $s_{t}\in S$, the exact action is not properly
defined a priori and the agent needs to randomly choose actions $a_{t}$ from a
set of actions, $A$, to find the optimal action at that state based on the
q-value of action $Q(s_{t},a_{t})$. The mathematical formulation governing
this process is known as Markov Decision Process (MDP), where the next state
is completely defined by the current state and action (Markov property) and
choice of an action leads to a scalar reward [1].

The RL problem is exactly defined and solved for a discrete state space.
However, for a continuous state-space problem, function approximators have to
be utilized to learn the value function and/or policy function. Deep neural
networks being universal function approximators, can be used to learn such
functions. The problem that appears then is that exact convergence is not
known and thus, practitioners have to resort to using reward per episode as
the common metric to assess if the training is complete.

##  III Methodology

###  III-A Coverage assessment

The goal of coverage assessment is to find the gaps in coverage of test data
for the model under test. The pseudo-code for the coverage assessment
algorithm is given in Algorithm 1. The first issue in getting a true coverage
estimate is that the state space in DRL applications is inherently high-
dimensional. For an efficient computation of coverage, we need to reduce the
dimensionality of the samples from the simulation. The preliminary requirement
in such a dimensionality reduction is that the local neighborhood be preserved
i.e. similar points in the high dimensional space should appear close by in
the low-dimensional embedding. In this paper, we utilize t-SNE, a nonlinear
dimensionality reduction technique proposed by [11] which embeds data from a
high-dimensional setting onto a two or three dimensional space. t-SNE
constructs a probability distribution over pairs of high-dimensional samples
such that neighboring points are assigned higher probability while distant
points are assigned lower probability. Then, a low-dimensional embedding is
chosen such that the Kullback-Leibler divergence (KL divergence) is minimized
with the high-dimensional probability distribution.

The other issue in developing a coverage assessment technique for DRL is that
the data is generated using a simulation concurrently with the training
process. The extents of the coverage are, thus, unknown a priori. We create a
buffer for storing the state samples and after a certain ‘k’ number of
episodes, create a two-dimensional t-SNE embedding from the samples. The
pseudo-coverage for the training process until the current episode can then be
defined as the area covered by these embedded state samples. There are
different approximations for calculation of area of scatter samples. A common
one is estimating convex hulls of the sampled points. However, for points not
in a single cluster, the calculation of convex hulls is complex and not
completely necessary. In this paper, we utilize an over-approximation using a
regularly spaced 2-D grid. The idea of the evenly spaced grid is similar to
occupancy grid mapping [12] where the non-zero grid-cells represent the
coverage of the embedded states. Based on the current extent of the embedded
states, we create a grid with a user-defined cell size. The cell size is an
important parameter, where there is a trade-off between granularity and
computation. The voting process calculates the histogram of the 2-D points and
returns the non-zero grid-cells. The buffer is reset and the process is
repeated for the next ‘k’ episodes with the non-zero grid-cells being added to
the previous count. APC is calculated as the ratio of non-zero grid-cells to
the total number of grid-cells.

Result: Approximate Pseudo-Coverage

$grids\leftarrow Empty$;

while  _training_ do

$statesBuffer\leftarrow Empty$;

Run k episodes and collect statesBuffer;

$statesEmbedded\leftarrow RunTsne(statesbuffer)$;

$limits\leftarrow CalculateLimits(statesembedded)$;

if  _$grids==Empty$_ then

$grids\leftarrow CreateGrid(limits,cellsize)$;

else

$grids\leftarrow AdjustGrid(grids)$;

end if

$voteGrids\leftarrow Histogram(statesEmbedded,grids)$;

end while

$APC=\dfrac{count(NonzeroGrids)}{count(grids)}$;

Algorithm 1 Coverage assessment

###  III-B Coverage maximization

While the process of executing more and more episodes would continue to
increase the coverage asymptotically, it might take a long time to fill the
gaps in coverage, esp. for high-dimensional state space problems. Coverage
maximization automates this process leading to faster coverage attainment. The
idea behind maximizing coverage is to generate a feedback loop to the
simulation, creating an adaptive process. We detail the procedure for coverage
maximization in Algorithm 2. In order to balance exploration of the unknown
initial states with the exploitation of known high-scoring initial states, we
implement an epsilon-greedy policy for determination of initial states.

In the beginning of the training process, we want to find unique initial
states which are in the neighborhood of the known initial states. We implement
a RRT structure to find these unique initial states. RRT is designed to
efficiently search non-convex, high-dimensional spaces such that the tree is
constructed by biasing the exploration towards unexplored portions of the
space [18]. The idea is to cover the entire state space efficiently without
any coverage gaps.

For the exploitation phase, we want to explore the vicinity of high-scoring
initial states. We utilize an idea from neuroevolution [19] where the best
performing agent (described as parent) is mutated by applying a small Gaussian
noise to the initial state. Therefore, we pick the highest scoring initial
state (or select randomly from equally highest scoring states) and perturb the
state by a small amount. This provides adequate coverage around the
potentially high-scoring initial states.

Result: initStatesList

initialize $epsilon\leftarrow 1$;

initialize $epsilonDecay\leftarrow 0.998$;

initialize $rrt\leftarrow beginNode$;

initialize $rrtDist$;

initialize $initStateList\leftarrow Empty$;

while  _training_ do

if  _$random\leq epsilon$_ then

$initStateCandidate\leftarrow Random(stateSize)$;

$neighborNode\leftarrow GetNeighborNode(initStateCandidate)$;

$newNode\leftarrow CreateNode(neighborNode,rrtDist)$;

$rrt\leftarrow Append(newNode)$;

else

$maxInitState\leftarrow Max(initStatesList)$;

$initState\leftarrow Perturb(maxInitState)$;

end if

$initStatesList\leftarrow Append(initState)$;

$epsilon\leftarrow epsilon*epsilonDecay$;

end while

Algorithm 2 Coverage Maximization

##  IV Results

The coverage assessment and maximization technique is validated in two example
environments: Cartpole [20] and highway driving [21].

Figure 1: Visualizations of embedded states colored by the q-value of that
state with the first action at every 100 episodes.  Figure 2: Visualizations
of embedded states derived using coverage maximization colored by the q-value
of that state with the first action at every 100 episodes.

Figure 3: Comparison of mean and bounds of (a) pseudo-coverage metrics and (b)
reward per episode for original method and our proposed method with coverage
maximization strategy for 50 runs

###  IV-A Cartpole

The cartpole environment [20], also known as the inverted pendulum problem, is
a well-known problem in RL whereby a pole is attached to a cart. The pole is
unstable and the only way to keep it upright is by moving the cart
horizontally left or right by one unit. The episode ends when the pole is more
than 15∘ from vertical, or the cart moves more than $2.4$ units from the
center. The state vector is composed of the position of the cartpole, the
angle of the cartpole, and their derivatives. The default reward is the amount
of time the pole is upright in a given episode (time limit of 500 timesteps).
Due to the nature of the continuous state space, it is impossible to visit
every state and action pair during training or to know when the training is
complete. OpenAI judges the training to be complete when the agent gets an
average of 475 (out of a possible 500) or above for 100 consecutive episodes.

Deep Q-network (DQN) has been proposed by [22] which provided the first
algorithm which stably combines deep neural networks with RL to solve
continuous state discrete action problems. DQN uses a neural network that
gives the Q-values for every action and uses a buffer to store old states and
actions to sample from which helps to stabilize training. The cartpole
environment is solved using the DQN approach.

In order to provide an intuition of the pseudo-coverage for the state samples,
we plot the embedded state samples colored by the q-value corresponding to the
first action at every 100 episodes in Figure 1. The figure shows that even at
600 episodes, there are gaps in the embedded state space which leads to gaps
in coverage. Also, we find that in the initial episodes, the state samples
form sparse non-convex clusters finding the area of which is a costly and
time-consuming problem. The over-approximation using a grid provides an
inexpensive solution while keeping over-approximation bounded.

Correspondingly, we utilize coverage maximization to determine the initial
state samples keeping the random seed constant for both implementations. We
use the position and orientation of the cartpole as an abstraction of the
model in the RRT structure. This makes the computations efficient while still
allowing for sufficient exploration of the state space. The embedded state
samples for the environment with coverage maximization is shown in Figure 2.
Visually, we can see that the space utilization using coverage maximization
strategy is higher than normal starting from the early episodes and by 400
episode mark, there are almost no gaps present in the embedded state samples.
This provides a visual confirmation that coverage maximization does bias the
system to fill gaps in the coverage. To get some statistics on the efficacy of
the coverage, we run the original implementation and the proposed
implementation 50 times. The mean and standard deviation of the APC and the
reward per episode is given in Figure 3. We observe that the proposed method
is effective in increasing the coverage by 19.5% and the reward per episode by
42.3% at the end of 1000 episodes.

###  IV-B Highway driving

We implement an environment for simulated highway driving and tactical
decision making provided by [21]. We choose this environment as a
representative example of high dimensional environment with a lot of
variability. The environment is built on top of OpenAI gym environment and
features an easy-to-use abstraction layer for training DRL to navigate an
autonomous vehicle through traffic. Figure 4 shows the layout of the
simulation environment. The green rectangle denotes the ego-vehicle while the
blue rectangles denote the other vehicles. The obstacle vehicles are capable
of changing lanes with each vehicle making independent lane-change decisions
at every instant. The ego-vehicle motion is described by kinematic bicycle
model[23]. The longitudinal motion of the other vehicles are governed by the
Intelligent Driver model [24] while the lane change decisions are governed by
the MOBIL model [25].

Figure 4: Highway driving simulation environment

Figure 5: Comparison of (a) pseudo-coverage metrics and (b) reward per episode
for original method and our proposed method with coverage maximization
strategy

The state space is $V\times F$ array comprised of $F$ features for $V$
vehicles. The features included are _S ={presence, x, y, vx, vy}_. The action
space $A$ is discrete with _A ={left-lane, idle, right-lane, faster, slower}_.
The reward function is composed of a velocity term and a collision term:

$R(s,a)=a\times\dfrac{v-v_{min}}{v_{max}-v_{min}}-b\times collision$ (1)

where $v,v_{max}$ and $v_{min}$ are the current speed, maximum speed and
minimum speed respectively and $a$ and $b$ are coefficients of the velocity
and collision penalty terms respectively.

We calculate the coverage metric in a similar fashion to the previous example.
Since the initial placement of the vehicles is the most important factor in
randomization of the simulation, for constructing the RRT structure, we
utilize the longitudinal and lateral positions as abstractions of the entire
space for the $V$ vehicles.

In Figure 5(a) and (b), we provide the coverage metrics and reward for a run
of the environment with and without the proposed implementation. The plots
show that for a high-dimensional state space DRL problem, our proposed
approach is able to increase the coverage while simultaneously maximizing the
cumulative reward per episode.

##  V Conclusions

In this work, we studied the necessity of utilizing coverage as a metric for
training of DRL tasks. We proposed a new metric to review the quality of the
sampled states for high-dimensional state space. For the two standard tasks,
we provided the pseudo-measure of coverage by transforming the state space
samples into a 2-D grid and then voting on the occupied grids. We also
proposed using a feedback loop from the training process to the simulation
environment and forcing the simulation environment to explore efficiently in
the bounds of the state space, both of which led to substantial improvements
for the training process. Thus, we utilized RRT for exploration of unexplored
initial states and perturbing high-scoring initial state for exploitation to
search for higher reward states.

For future work, we would like to extend the coverage estimation towards
providing convergence guarantees for DRL tasks based on the proposed coverage
metric.

## Acknowledgment

I would like to thank Dr. Aghapi Mordovanakis for taking the time to edit the
manuscript and for stimulating discussions. Also, I would like to mention Dr.
Md Tawhid Bin Waez for providing invaluable suggestions during the ideation
phase.

## References

  * [1] R. S. Sutton and A. G. Barto, _Reinforcement learning: An introduction_. MIT press, 2018.
  * [2] D. Silver, J. Schrittwieser, K. Simonyan, I. Antonoglou, A. Huang, A. Guez, T. Hubert, L. Baker, M. Lai, A. Bolton, _et al._ , “Mastering the game of go without human knowledge,” _nature_ , vol. 550, no. 7676, pp. 354–359, 2017.
  * [3] D. Silver, T. Hubert, J. Schrittwieser, I. Antonoglou, M. Lai, A. Guez, M. Lanctot, L. Sifre, D. Kumaran, T. Graepel, _et al._ , “A general reinforcement learning algorithm that masters chess, shogi, and go through self-play,” _Science_ , vol. 362, no. 6419, pp. 1140–1144, 2018.
  * [4] V. Mnih, K. Kavukcuoglu, D. Silver, A. Graves, I. Antonoglou, D. Wierstra, and M. Riedmiller, “Playing atari with deep reinforcement learning,” _arXiv preprint arXiv:1312.5602_ , 2013.
  * [5] S. Levine, C. Finn, T. Darrell, and P. Abbeel, “End-to-end training of deep visuomotor policies,” _The Journal of Machine Learning Research_ , vol. 17, no. 1, pp. 1334–1373, 2016.
  * [6] Y. Hollander, “Machine learning for coverage maximization,” Sep 2016. [Online]. Available: https://blog.foretellix.com/2016/09/01/machine-learning-for-coverage-maximization/
  * [7] D. Araiza-Illan, D. Western, A. Pipe, and K. Eder, “Coverage-driven verification—,” in _Haifa Verification Conference_. Springer, 2015, pp. 69–84.
  * [8] T. Byun and S. Rayadurgam, “Manifold for machine learning assurance,” _arXiv preprint arXiv:2002.03147_ , 2020.
  * [9] S. Mani, A. Sankaran, S. Tamilselvam, and A. Sethi, “Coverage testing of deep learning models using dataset characterization,” _arXiv preprint arXiv:1911.07309_ , 2019.
  * [10] R. Bellman, “Dynamic programming,” _Science_ , vol. 153, no. 3731, pp. 34–37, 1966.
  * [11] L. v. d. Maaten and G. Hinton, “Visualizing data using t-sne,” _Journal of machine learning research_ , vol. 9, no. Nov, pp. 2579–2605, 2008.
  * [12] H. Moravec and A. Elfes, “High resolution maps from wide angle sonar,” in _Proceedings. 1985 IEEE international conference on robotics and automation_ , vol. 2. IEEE, 1985, pp. 116–121.
  * [13] A. Tavakoli, V. Levdik, R. Islam, and P. Kormushev, “Prioritizing starting states for reinforcement learning,” _arXiv preprint arXiv:1811.11298_ , 2018\.
  * [14] S. M. LaValle, “Rapidly-exploring random trees: A new tool for path planning,” 1998.
  * [15] A. Piziali, “Coverage-driven verification,” _Functional Verification Coverage Measurement and Analysis_ , pp. 109–137, 2008.
  * [16] M. Benjamin, D. Geist, A. Hartman, G. Mas, R. Smeets, and Y. Wolfsthal, “A study in coverage-driven test generation,” in _Proceedings of the 36th annual ACM/IEEE Design Automation Conference_ , 1999, pp. 970–975.
  * [17] H. Lackner, H. Schlingloff, and A. Berlin, “Modeling for automated test generation-a comparison.” in _MBEES_. Citeseer, 2012, pp. 57–70.
  * [18] S. M. LaValle and J. J. Kuffner, “Rapidly-exploring random trees: Progress and prospects,” _Algorithmic and computational robotics: new directions_ , no. 5, pp. 293–308, 2001.
  * [19] F. P. Such, V. Madhavan, E. Conti, J. Lehman, K. O. Stanley, and J. Clune, “Deep neuroevolution: Genetic algorithms are a competitive alternative for training deep neural networks for reinforcement learning,” _arXiv preprint arXiv:1712.06567_ , 2017.
  * [20] G. Brockman, V. Cheung, L. Pettersson, J. Schneider, J. Schulman, J. Tang, and W. Zaremba, “Openai gym,” _arXiv preprint arXiv:1606.01540_ , 2016.
  * [21] E. Leurent, “An environment for autonomous driving decision-making,” https://github.com/eleurent/highway-env, 2018.
  * [22] V. Mnih, K. Kavukcuoglu, D. Silver, A. A. Rusu, J. Veness, M. G. Bellemare, A. Graves, M. Riedmiller, A. K. Fidjeland, G. Ostrovski, _et al._ , “Human-level control through deep reinforcement learning,” _Nature_ , vol. 518, no. 7540, p. 529, 2015.
  * [23] P. Polack, F. Altché, B. d’Andréa Novel, and A. de La Fortelle, “The kinematic bicycle model: A consistent model for planning feasible trajectories for autonomous vehicles?” in _2017 IEEE Intelligent Vehicles Symposium (IV)_. IEEE, 2017, pp. 812–818.
  * [24] M. Treiber, A. Hennecke, and D. Helbing, “Congested traffic states in empirical observations and microscopic simulations,” _Physical review E_ , vol. 62, no. 2, p. 1805, 2000.
  * [25] A. Kesting, M. Treiber, and D. Helbing, “General lane-changing model mobil for car-following models,” _Transportation Research Record_ , vol. 1999, no. 1, pp. 86–94, 2007.


# Multiple Charged Meson Production in Exclusive $B_{c}$ Decays:
$K+4\pi,$$KK+3\pi$, $7\pi$ Cases.

A.V. Luchinsky  Institute for High Energy Physics, Protvino, Russia

###### Abstract

Theoretical analysis of three exclusive decays,
$B_{c}\to\psi^{(^{\prime})}+\mathcal{R},$where $R=K+4\pi,$$KK+3\pi$, or
$7\pi$, is given. Using the factorization method and the resonance
approximation (methods proved to be very useful in the analysis of some other
similar decays), we have obtained analytical expressions of the required
amplitudes and created distributions over some interesting mass combinations.
Presented results could be used for comparison with forthcoming experimental
data and better understanding of the nature of heavy quarkonia.

##  I Introduction

Heavy quarkonia mesons, that is particles that in valence approximation are
built from heavy (i.e. $c$ of $b$) pairs can be considered as a unique
laboratory to study strong interaction both in perturbative and non-
perturbative regimes. On the partonic level the processes of their production
and decay can be described in terms of creation or annihilation of the heavy
quark, that hadronizes later into experimentally obesrved particle. Doubly
heavy meson with an opened flavor,$B_{c}=bc$ take an intermediate position
between charmonium $\left(c\bar{c}\right)$ and bottomonium
$\left(b\bar{b}\right)$ particles, that allows one to use it as a test for
models that were developed to describe heacy quarkonium particles with hidden
flavor. Nice theoretical review concerning different properties of this
particle can be found, for example, in paper Gershtein et al. (1995)

The presented paper is devoted to phenomenological analysis of some exclusive
$B_{c}$-meson decay with the vector charmonium production accomplished by
$K+4\pi$, $KK+3\pi$ or $7\pi$. This work is a continuation of series of papers
at were devoted to production of other systems of light mesons in the similar
processes.

General idea used in all these works is simple. It is well known that in
valence approximation $B_{c}$-meson is built from $b$\- and $c$-quarks.
Exclusive decay of this particle into vector charmonium and some set of light
mesons on the partonic level can be described by the weak decay of the
constituent $b$-quark: $b\to Wc$. Produced $W$ boson hadronizes into the final
system of light particles, while $c$-quark with a spectator $\bar{c}$ forms
$J/\psi$ or $\psi(2S)$ . The subprocess $B_{c}\to W\psi^{(^{\prime})}$ can be
described in terms of weak $B_{c}$-meson form-factors. The suitable
theoretical model for the W- boson hadronization, on the other hand, is the
so-called resonance approach.

Described above simple model turned out to be surprisingly powerful. It was
used to describe decays $B_{c}\to J/\psi+3\pi$, etc., and the predictions for
differential distributions and integrated branching fractions are in good
agreement withe experimental results. In the present paper we are extending
the analysis to three new decays.

The rest of the paper is organized as follows. In the next section a short
description of the adopted approach is given and the parametrization of the
$B_{c}$-meson’s form-factors is resented. Section 3 is devoted to theoretical
description of same observed already decays, that mile be used as subprocesses
for new ones. These new decays will be studied in section 4, while the last
section is reserved to some discussion and conclusion.

##  II Exclusive $B_{c}$ decays

In the partonic approximation $B_{c}$ meson is build from $b\bar{c}$ quarks,
so its exclusive decays into vector charmonium and a set of light mesons can
be described as a weak decay of the constituent $b$ quark. Produced $c$ with a
spectator $\bar{c}$ forms the final charmonium, while $W$-boson hadronizes
into the final system of light mesons $\mathcal{R}$. The schematic Feynman
diagram of this process is shown in figure 1.

Figure 1: Feynman diagram for $B_{c}\to\psi\mathcal{R}$ decay

In the factorization approximation this diagram can be written in the form

$\displaystyle\mathcal{M}$ $\displaystyle=$
$\displaystyle\frac{G_{F}V_{bc}}{\sqrt{2}}a_{1}H^{\mu}\epsilon_{\mu}^{(\mathcal{R})},$

where the Wilson coefficient $a_{1}(m_{b})\approx 1.14$ describes the effect
of final state interaction Buchalla et al. (1996), $G_{F}$ is a Fermi coupling
constant, $V_{bc}$ is a corresponding coefficient of the CKM matrix, while
$H^{\mu}$ and $\epsilon_{\mu}^{(\mathcal{R})}$ are the amplitudes of the
$B_{c}\to\psi W$ and $W\to R$ transitions respectively.

Let us define the the first amplitude first. It is clear, that it can depend
only on the momenta of the initial and final heavy mesons $p_{1,2}$. There are
several ways how it can be written in the Lorentz-invariant form and in our
paper we will use the parametrization adopted, for example, in paper Kiselev
et al. (2000a):

$\displaystyle H_{\mu}$ $\displaystyle=$
$\displaystyle\left[2M_{\psi}A_{0}\left(q^{2}\right)\frac{q^{\mu}q^{\nu}}{q^{2}}+\left(M_{B_{c}}-M_{\psi}\right)A_{1}\left(q^{2}\right)\left(g^{\mu\nu}-\frac{q^{\mu}q^{\nu}}{q^{2}}\right)-\right.$
$\displaystyle\left.A_{2}\left(q^{2}\right)\left(p_{1}^{\mu}+p_{2}^{\mu}-\frac{M_{B_{c}}^{2}-M_{\psi}^{2}}{q^{2}}q^{\mu}\right)-\frac{2iV\left(q^{2}\right)}{M_{B_{c}+M_{\psi}}}e^{\mu\nu\alpha}p_{1\alpha}p_{2}\right]\epsilon_{\nu}^{(\psi)}.$

Here $\epsilon^{(\psi)}$ is the polarization vector of final charmonium meson
and $A_{0,1,2}\left(q^{2}\right)$, $V\left(q^{2}\right)$ are axial and vector
form factors of $B_{c}\to\psi W$ transition. It is clear that these functions
cannot be determined from perturbative theory, so some other approach should
be used, such as for example QCD sum rules Huang and Zuo (2007); Kiselev et
al. (2000b, c, a); Kiselev (2002), different potential quark models Kiselev et
al. (1993); Gershtein et al. (1995, 1997); Colangelo and De Fazio (2000);
Ivanov et al. (2005), light-front models Anisimov et al. (1999); Choi and Ji
(2009a, b), etc. In our paper we will use form-factors sets presented in paper
Kiselev et al. (2000a).

These form factors can be parametrized as

$\displaystyle F\left(q^{2}\right)$ $\displaystyle=$
$\displaystyle\frac{F\left(0\right)}{1-q^{2}/m_{fit}^{2}},$

where $q^{2}=\left(p_{1}-p_{2}\right)^{2}$ is the transferred momentum squared
and the values of function at zero and maximal agruments are listed in table
1.

| $F(0)$ | $F\left(q_{max}^{2}=M_{B_{c}}^{2}-M_{J/\psi}^{2}\right)$ |  | $F(0)$ | $F\left(q_{max}^{2}=M_{B_{c}}^{2}-M_{J/\psi}^{2}\right)$
---|---|---|---|---|---
$A_{0}$ | 0.60 | 1.6 | $A_{2}$ | 0.69 | 1.4
$A_{1}$ | 0.63 | 1.3 | $V$ | 1.0 | 2.1
Table 1: Parameters of the $B_{c}\to\psi^{(^{\prime})}W$ transition form
factors

All information about the final system $\mathcal{R}$ is hidden in the
effective polarization vector $\epsilon^{(\mathcal{R})}$. Its explicit form
and numerical values depend on the number of final particles, their types and
momenta. The only general issue is that in the limit of isospin conservation
the relation $q^{\mu}\epsilon_{\mu}^{(\mathcal{R})}$ should hold, to the
contribution of $A_{0}$ form factor vanishes. More detailed investigation
requires some model assumptions on the physics of the underlying processes.
Our previous works it was shown, that the resonance model gives the results
that are in a good agreement with the experiment. In the framework of this
model the amplitude of the process under consideration is written in the terms
of hadronic resonances with suitable quantum numbers and the form of the
amplitudes is chosen accordingly. In the present paper we will continue to use
such an approach to describe the production of $\mathcal{R}=K+4\pi,$
$KK+3\pi$, and $7\pi$ states. It turns out that the corresponding processes
can be explained in terms of previously considered and tested on experiment
reactions with $KK\pi$, $3\pi$, $5\pi$ production. For this reason we will
first consider these processes.

##  III Known Decays

Figure 2: Feynman diagrams for old decay processes

Let us first consider $B_{c}\to\psi+3\pi$ decay. In resonance approximation it
can be described by the diagram shown on figure 2(a). The corresponding
amplitude can be written in the form

$\displaystyle\epsilon_{(3\pi)}^{\mu}\left(Q\to k_{1}k_{2}k_{3}\right)$
$\displaystyle\sim$ $\displaystyle
BW_{a1}\left(Q^{2}\right)\left(g_{\mu\alpha}-\frac{Q_{\mu}Q_{\alpha}}{Q^{2}}\right)q_{[13]}^{\alpha}BW_{\rho}\left(q_{\\{13\\}}^{2}\right)+\left\\{q_{1}\Leftrightarrow
q_{2}\right\\},$

where $q_{1,2,3}$ are the momenta of two final $\pi^{+}$ mesons and $\pi^{-}$
meson, $Q=q_{1}+q_{2}+q_{3}$ is the total momentum of virtual $a_{1}$, and
notations

$\displaystyle q_{[ij]}^{\alpha}$ $\displaystyle=$ $\displaystyle
q_{i}^{\alpha}-q_{j}^{\alpha},\qquad
q_{\\{ij\\}}^{\alpha}=q_{i}^{\alpha}+q_{j}^{\alpha}$

were introduced. As you can see, we are using usual Feynman rules to write
this relation and general from of the vertices is determined by quantum
numbers of the participating particles. For example, the interaction
$(q_{1}-q_{3})_{\alpha}$ was selected for $\rho\to\pi\pi$ decay since we are
dealing with vector particle decay into two scalars, so P-wave is in place.
The only interesting thing in the above equation is the propagators of the
virtual $a_{1}$ and $\rho$ mesons. Following Flatte (1976) they are written
using Flatte parametrization, where the energy dependence of the resonance
width is taken into account:

$\displaystyle BW_{a}\left(q^{2}\right)$ $\displaystyle=$
$\displaystyle\frac{M_{a}^{2}}{M_{a}^{2}-q^{2}-iM_{a}\Gamma_{a}(q)},$
$\displaystyle BW_{\rho}\left(q^{2}\right)$ $\displaystyle=$
$\displaystyle\frac{1}{1+\beta}\left[\frac{M_{\rho}^{2}}{M_{\rho}^{2}-q^{2}-iM_{\rho}\Gamma_{\rho}\left(q^{2}\right)}+\beta\frac{M_{\rho^{\prime}}^{2}}{M_{\rho^{\prime}}^{2}-q^{2}-iM_{\rho^{\prime}}\Gamma_{\rho^{\prime}}\left(q^{2}\right)}\right]$

This process was studied theoretically in works Berezhnoy et al. (2011);
Luchinsky (2012); Kuhn and Was (2008)and later experimental confirmation by
the LHCb collaboration followed Aaij et al. (2012) As you can see from figure
3(a), the agreement between theory and experiment is pretty good.

Figure 3: Comparison of theoretical predictions (blue line) with experiment
(red dots).LHCb results shown in the left, center and right figures were
published in papers Aaij et al. (2022), Aaij et al. (2014) and Aaij et al.
(2022) respectively.

This process can be used as a building block to describe the production of the
higher number of charged $\pi$ mesons, the corresponding Feynman diagram is
shown in figure 2(b). The amplitude of this diagram is equal to

$\displaystyle\epsilon_{(5\pi)}^{\mu}$ $\displaystyle=$ $\displaystyle
BW_{a}\left(Q^{2}\right)BW_{f}\left(\left(q_{4}+q_{5}\right)^{2}\right)\epsilon_{(3\pi)}^{\mu}\left(q_{1},q_{2},q_{3}\right)+\text{symmetrization},$

and the comparison of theoretical Luchinsky (2012) and experimental Aaij et
al. (2017) distributions in shown in figure 3(b). As you can see, the
agreement is again quite reasonable.

Let is consider now production of $KK\pi$ state. Two Feynman diagrams, that
were used to describe it, are shown in figure 2(c), (d). The corresponding
amplitude can be written in the following form

$\displaystyle\epsilon_{KK\pi}^{\mu}$ $\displaystyle=$ $\displaystyle
BW_{a}\left(Q^{2}\right)BW_{K^{*}}\left(q_{\\{12\\}}^{2}\right)q_{[12]}^{\mu}$

Figure 3(c) shows that there is a good agreement between theoretical Luchinsky
(2013) and experimental Aaij et al. (2022) results.

##  IV New Decays

In this section three new decays will be described, we will present amplitudes
of these processes and theoretical predictions for some typical distributions.

Figure 4: Diagram for $W\to\mathcal{R}$ transition Figure 5: Distributions in
$B_{c}\to J/\psi+\mathcal{R}$ decays:

Let us consider $\mathcal{R}=K^{+}\pi^{+}\pi^{+}\pi^{-}\pi^{-}$ final state
first. In the resonance approximation $W$-boson hadronization into this system
of light particles can be described by the diagram shown in figure 4(a). As
you can see, this diagram includes $a_{1}\to 3\pi$ transition as a subprocess,
so we can use presented above parametrization to describe it. The form of the
other vertices’ amplitudes is easily determined by the quantum numbers of the
interacting particles. As a result, the following form of the transition
amplitude was used later:

$\displaystyle\epsilon_{K4\pi}^{\mu}$ $\displaystyle=$ $\displaystyle
BW_{a}\left(Q^{2}\right)BW_{K^{*}}\left(q_{\\{14\\}}^{2}\right)e^{\mu\nu\alpha\beta}q_{[14]}^{\nu}\left(q_{\\{14\\}}-q_{\\{234\\}}\right)^{\alpha}\epsilon_{(3\pi)}^{\beta}\left(q_{2,}q_{3,}q_{5}\right).$

In figure 5(a), (b) some distributions are shown. On the last figure in the
$m_{K\pi}$ system you can clearly see a peak caused by $K^{*0}$ resonance.

The subprocess $a_{1}\to KK\pi$, on the other hand, can be used in production
of $KK+3\pi$ final state. The diagram is shown in figure 4(b) and the
corresponding can be written in the form

$\displaystyle\epsilon_{KK3\pi}^{\mu}$ $\displaystyle=$ $\displaystyle
BW_{a}\left(Q^{2}\right)BW_{f}\left(q_{\\{45\\}}^{2}\right)\epsilon_{KK\pi}^{\mu}\left(q_{1}q_{2}q_{3}\right)+\text{symmetrization}.$

In figure 5(c), (d) some distributions of the corresponding hadronic process
are shown. Note that in this case there are no peaks in the $m_{K\pi}$
spectrum since there are no resonances in $K^{+}\pi^{+}$ channel.

The final process to be considered in our article is
$B_{c}\to\psi^{(^{\prime})}+7\pi$ decay. The corresponding $W$-boson
hadronization diagram is shown in figure 4(c).. As you can see, this reaction
includes $a_{1}\to 5\pi$ transition as a subprocess, to its amplitude is
written as

$\displaystyle\epsilon_{7\pi}^{\mu}$ $\displaystyle=$ $\displaystyle
BW_{a}\left(Q^{2}\right)BW_{f}\left(q_{\\{47\\}}^{2}\right)\epsilon_{5\pi}^{\mu}\left(q_{1,}q_{2},q_{3},q_{5},q_{6}\right)+\text{symmetrization}.$

In figure 5(e), (f) distributions over the transferred momentum and mass of
the $\pi^{+}\pi^{-}$ pair are shown. In the latter distribution you can see a
small $\rho$-meson peak, which is almost hidden by the broad $f_{0}$
resonance.

##  V Conclusion

The presented paper is devoted to theoretical analysis of some exclusive
$B_{c}$-meson decays. This work is a continuation of a series of papers, that
are using the same approach for description of some other decays.

The theoretical model used in all these works is pretty simple. The reactions
under consideration are represented as two-step process. The first step is
weak $B_{c}$-meson decay into final victor charmonium and virtual $W$-boson,
which then hadronizes into a system of light mesons. The first stage was
described in terms of $B_{c}$-meson form-factors, while for the second stage
the resonance approximation was used.

This simple method turned out to be surprisingly useful and able to produce
theoretical predictions for a number of Bc-meson decays, that are in good
agreement with experimental results. In the presented paper it was used to
obtain analytical expressions for the amplitudes of three more decays:
$B_{c}\to J/\psi+K+4\pi$, $B_{c}\to J/\psi+KK+3\pi$, and $B_{c}\to
J/\psi+7\pi$. In addition, the numerical analysis of these reactions was
performed and same interesting differential distributions can be found in
thins work.

There is, of course, lot of work to be done in this field. For example, there
are same unknown normalization constants in the used model, that cannot be
determined from experimental data. For this reason we do not make predictions
for the branching fractions of the considered decays, only normalized
distributions are presented. Calculation of the branching fractions will be
the topic of our future work.

The author would like to thank Dr. A. Likhoded, Dr. I. Belyaev, and Dr. D.
Pereima for useful, stimulating, and motivating discussions. This research was
done with support of RFBR grant 20-02-00154 A.

## References

  * Gershtein et al. (1995) S. S. Gershtein, V. V. Kiselev, A. K. Likhoded, and A. V. Tkabladze, Phys. Usp. 38, 1 (1995), eprint hep-ph/9504319.
  * Buchalla et al. (1996) G. Buchalla, A. J. Buras, and M. E. Lautenbacher, Rev. Mod. Phys. 68, 1125 (1996), eprint hep-ph/9512380.
  * Kiselev et al. (2000a) V. V. Kiselev, A. E. Kovalsky, and A. K. Likhoded, Nucl. Phys. B 585, 353 (2000a), eprint hep-ph/0002127.
  * Huang and Zuo (2007) T. Huang and F. Zuo, Eur. Phys. J. C 51, 833 (2007), eprint hep-ph/0702147.
  * Kiselev et al. (2000b) V. Kiselev, A. Likhoded, and A. Onishchenko, Nucl.Phys. B569, 473 (2000b), eprint hep-ph/9905359.
  * Kiselev et al. (2000c) V. V. Kiselev, A. E. Kovalsky, and A. K. Likhoded, in _5th International Workshop on Heavy Quark Physics_ (2000c), eprint hep-ph/0006104.
  * Kiselev (2002) V. V. Kiselev (2002), eprint hep-ph/0211021.
  * Kiselev et al. (1993) V. V. Kiselev, A. K. Likhoded, and A. V. Tkabladze, Phys. Atom. Nucl. 56, 643 (1993).
  * Gershtein et al. (1997) S. S. Gershtein, V. V. Kiselev, A. K. Likhoded, A. V. Tkabladze, A. V. Berezhnoy, and A. I. Onishchenko, in _4th International Workshop on Progress in Heavy Quark Physics_ (1997), eprint hep-ph/9803433.
  * Colangelo and De Fazio (2000) P. Colangelo and F. De Fazio, Phys. Rev. D 61, 034012 (2000), eprint hep-ph/9909423.
  * Ivanov et al. (2005) M. A. Ivanov, J. G. Korner, and P. Santorelli, Phys. Rev. D 71, 094006 (2005), [Erratum: Phys.Rev.D 75, 019901 (2007)], eprint hep-ph/0501051.
  * Anisimov et al. (1999) A. Y. Anisimov, P. Y. Kulikov, I. M. Narodetsky, and K. A. Ter-Martirosian, Phys. Atom. Nucl. 62, 1739 (1999), eprint hep-ph/9809249.
  * Choi and Ji (2009a) H.-M. Choi and C.-R. Ji, Phys. Rev. D 80, 114003 (2009a), eprint arXiv:0909.5028 [hep-ph].
  * Choi and Ji (2009b) H.-M. Choi and C.-R. Ji, Phys. Rev. D 80, 054016 (2009b), eprint arXiv:0903.0455 [hep-ph].
  * Flatte (1976) S. M. Flatte, in _11th Rencontres de Moriond: new fields in hadronic physics_ (1976), pp. 83–92.
  * Berezhnoy et al. (2011) A. V. Berezhnoy, A. K. Likhoded, and A. V. Luchinsky (2011), eprint arXiv:1104.0808 [hep-ph].
  * Luchinsky (2012) A. V. Luchinsky, Phys. Rev. D 86, 074024 (2012), eprint arXiv:1208.1398 [hep-ph].
  * Kuhn and Was (2008) J. H. Kuhn and Z. Was, Acta Phys. Polon. B 39, 147 (2008), eprint hep-ph/0602162.
  * Aaij et al. (2012) R. Aaij et al. (LHCb), Phys. Rev. Lett. 108, 251802 (2012), eprint arXiv:1204.0079 [hep-ex].
  * Aaij et al. (2022) R. Aaij et al. (LHCb), JHEP 01, 065 (2022), eprint arXiv:2111.03001 [hep-ex].
  * Aaij et al. (2014) R. Aaij et al. (LHCb), JHEP 05, 148 (2014), eprint arXiv:1404.0287 [hep-ex].
  * Aaij et al. (2017) R. Aaij et al. (LHCb), Eur. Phys. J. C 77, 72 (2017), eprint arXiv:1610.01383 [hep-ex].
  * Luchinsky (2013) A. V. Luchinsky (2013), eprint arXiv:1307.0953 [hep-ph].