# Elucidating the internal structure of hadrons through direct photon
production

David F. Rentería-Estrada  Roger J. Hernández-Pinto  German F. R. Sborlini

###### Abstract

The accurate description of the internal structure of hadrons is a very
challenging task. In order to compare the predictions with the highly-accurate
experimental data, it is necessary to control any possible source of
theoretical uncertainties. Thus, we can use the information extracted from
final state measurement to constrain our knowledge about the internal
structure of hadrons. In this work, we describe how direct photon production
can be exploited to unveil details about the partonic distributions inside
protons. Also, we explain how to describe QCD-QED corrections to hadron plus
photon production at colliders, focusing on the accurate reconstruction of the
partonic momentum fractions from experimentally accessible observables.

##  1 Introduction

One of the best strategies for studying the internal structure of the non-
fundamental particles is based on the parton model. It relies on factorization
properties and parton distribution functions (PDFs) that allow the
characterization of the hadronic structure [1]. The PDFs are obtained from the
experiments, for example from the Large Hadron Collider (LHC), through
advanced model fitting methods. However, this methodology is not sufficient to
explain the total spin of the proton. Using the currently available data, only
30% of the total spin can be explained from the quarks contributions [2, 3,
4]. Furthermore, data extracted from deep inelastic scattering (DIS)
experiments are not sufficient to place restrictions on the polarized
distributions of quarks and gluons [5, 6].

There are different methodologies to characterize the internal structure of
hadrons. A method that gives direct access to the dynamics of the partons and
obtain more information about the proton’s spin, consists of detecting the
photons created by the internal interaction of the partons within the hadron:
these are the so-called _hard_ photons. The proton-proton collision creates a
dense medium with very high temperatures: as a consequence, any particle that
couples to QCD partons will interact with this medium. Due to the non-
perturbative nature of this interactions, it is usually preferable to detect
particles that weakly interacts with this medium: this is why we will focus on
the identification of hard photons.

In this work, we analyze the production of one direct photon accompanied by a
hadron, including leading-order (LO) corrections in QED and next-to-leading
order (NLO) corrections in QCD. Specifically, we analyze the effects of QED
corrections on the distribution of the cross-section, for the energies of the
PHENIX (200 GeV) and LHC (13 TeV) experiments [7]. This complements the
phenomenological study presented in Ref. [8] and constitutes a step forward
towards a better understanding of the process and an updated description of
the phenomenology of hadron-photon production at colliders.

##  2 Cross-section calculation

In the parton model, thanks to the factorization theorem [1], the cross-
section can be described by the convolution between parton distribution
functions, fragmentation functions (FFs) and the partonic cross-section. The
non-perturbative effects associated with the low energy interactions are
included in the PDFs and FFs, while the partonic cross-section is calculated
from a perturbative framework, i.e.

$\displaystyle d\sigma_{H_{1}\,H_{2}\rightarrow h\,\gamma}^{\rm DIR}$
$\displaystyle=$ $\displaystyle\sum_{a_{1}a_{2}a_{3}}\int
dx_{1}dx_{2}dz\,f^{(H_{1})}_{a_{1}}(x_{1},\mu_{I})f^{(H_{2})}_{a_{2}}(x_{2},\mu_{I})\,d^{(h)}_{a_{3}}(z,\mu_{F})d\hat{\sigma}^{\rm
DIR}_{a_{1}\,a_{2}\rightarrow a_{3}\,\gamma}\,.$ (1)

The PDF is represented by $f_{a}^{(H)}(x,\mu_{I})$, describing the probability
density of finding a parton $a$ inside a hadron $H$, with momentum fraction
$x$ and initial energy scale $\mu_{I}$. Similarly, the FF,
$d_{a_{3}}^{h}(z,\mu_{F})$, represents the probability density of generating
one hadron $h$ with momentum fraction $z$ at the energy scale $\mu_{F}$.

To estimate the cross-section we considered two possible scenarios: i) the
_direct_ contribution (i.e. the photon is generated directly from the partonic
interaction), and ii) the _resolved_ contribution (the photon comes from the
fragmentation of a hadron). Due to the quantum nature of the process, it is
not possible to identify the true origin of the detected photon. To mitigate
this problem we implement the smooth cone isolation algorithm [9]. This is a
good method that give the efficiently separate the direct from the resolved
component. Due to the high correlation between the production of hard photons
and isolated photons, we implemented the procedure described in Ref.[9] to
identify the hard photons. The smooth cone isolation algorithm allows for soft
gluons throughout space, and ensures an IR-safe cross-section. Finally, under
these assumptions, Eq. (1) can be written as

$\displaystyle d\sigma_{H_{1}\,H_{2}\rightarrow h\,\gamma}$ $\displaystyle=$
$\displaystyle\sum_{a_{1}a_{2}a_{3}}\int
dx_{1}dx_{2}dz\,f^{(H_{1})}_{a_{1}}(x_{1},\mu_{I})f^{(H_{2})}_{a_{2}}(x_{2},\mu_{I})\,d^{(h)}_{a_{3}}(z,\mu_{F})d\hat{\sigma}^{\rm
ISO}_{a_{1}\,a_{2}\rightarrow a_{3}\,\gamma}\,,$ (2)

where $d\hat{\sigma}^{\rm ISO}_{a_{1}\,a_{2}\rightarrow a_{3}\,\gamma}$ is the
partonic cross-section that contains the smooth cone isolation prescription.
There are two partonic channels that contribute at LO,

$q\bar{q}\rightarrow\gamma g\,,\quad qg\rightarrow\gamma q\,,$ (3)

while at NLO we must also consider:

$q\bar{q}\rightarrow\gamma gg\,,\quad qg\rightarrow\gamma gq\,,\quad
gg\rightarrow\gamma q\bar{q}\,,\quad q\bar{q}\rightarrow\gamma
Q\bar{Q}\,,\quad qQ\rightarrow\gamma qQ~{}.$ (4)

To implement NLO QCD corrections we rely on the FKS subtraction algorithm
[10], which divides the real emission phase-space into different regions
containing non-overlapping singularities and cancels them with the local IR
counter-terms. Renormalization was carried out within the $\overline{\rm MS}$
scheme. Also, to estimate LO QED contribution, we add

$\displaystyle d\hat{\sigma}^{\rm ISO,QED}_{a_{1}\,a_{2}\rightarrow
a_{3}\,\gamma}$ $\displaystyle=$
$\displaystyle\frac{\alpha^{2}}{4\pi^{2}}\,\int d{\rm PS}^{2\rightarrow
2}\,\frac{|{\cal
M}^{(0)}_{QED}|^{2}(x_{1}K_{1},x_{2}K_{2},K_{3}/z,K_{4})}{2\hat{s}}\,{\cal
S}_{2}\,,$ (5)

to Eq. (2). The new partonic channels involved in the QED corrections are
$q\gamma\rightarrow\gamma q$ and $q\bar{q}\rightarrow\gamma\gamma$, where the
photon PDF appears.

##  3 Results

We focus on the process $p\,+\,p\,\rightarrow\,h\,+\,\gamma$, specifically
looking for a charged pion in the final state. We implement a Monte-Carlo code
with the following default kinematic cuts:

  1. 1.

The rapidities of the photon and pion are restricted to $|\eta|\leq 0.35$.

  2. 2.

The transverse momentum of the photon is restricted to $5\,{\rm GeV}\leq
p_{T}^{\gamma}\leq 15\,{\rm GeV}$.

  3. 3.

The transverse momenta of the pion must be greater than 2 GeV.

  4. 4.

We do not consider restrictions on the azimuthal angles of the photon and the
pion. We constrain $\Delta\phi=|\phi^{\pi}-\phi^{\gamma}|\geq 2$ (to include
those events where the photon and pion are produced close to a back-to-back
configuration).

These cuts simulate PHENIX experiment at RHIC. We have explored center-of-mass
energies of $E_{CM}=200{\rm GeV}$ (PHENIX) and $E_{CM}=13{\rm TeV}$ (LHC).
Furthermore, we used more modern sets of PDFs and FFs than those in the work
of Ref. [8]. Then, we explored the impact of NLO QCD and LO QED corrections as
a function of the center-of-mass energy. We implemented the LHAPDF interface
[11] to select different PDF sets. For the default configuration, the NLO QCD
correction is fixed using the PDF set NNPDF3.1. To analyse the effects of the
combined LO QED + NLO QCD corrections, we use NNPDF3.1luxQEDNLO [12, 13] .
Finally, we update the FF routines to implement the recent set DSS2014 [14].
In Fig. 1 we compare the process for different center-of-mass energies. As
expected from the behaviour of the electromagnetic coupling, QED corrections
are more relevant with increasing of center-of -mass energy, as well as in the
high $p_{T}$ region.

Figure 1: Comparison of the combined NLO QCD + LO QED (purple line) versus the
pure NLO QCD corrections to the $p_{T}^{\pi}$ distribution of
$p+p\rightarrow\pi^{+}+\gamma$. In the left, we consider PHENIX kinematics,
whilst we switched to LHC energies in the right plot.

##  4 Conclusions

In this work, we have studied the phenomenology of photon-hadron production.
We included LO QED + NLO QCD corrections and performed a careful comparison
with the QCD only contributions. We used the NNPDF3.1luxQEDNLO set, centering
the analysis on the $p_{T}^{\pi}$ spectrum. We found small but non-negligible
corrections: ${\cal O}(2\,\%)$ for PHENIX and ${\cal O}(8\,\%)$ for LHC
center-of-mass energies. As already discussed in other works [8, 7], these
results suggest the hadron plus photon might be an important process to impose
tighter constraints on PDFs. In particular, using novel computational
techniques (such as machine learning and artificial intelligence), we could
expect to have a better understanding of the partonic kinematics by using the
photon as a probe to explore the dynamics of the collisions [15].

## Acknowledgements

We would like to thank P. Zurita for reading this manuscript and fruitful
discussions for future works. This research was partially supported by COST
Action CA16201 (PARTICLEFACE). The work of D. F. R.-E. and R. J. H.-P. is
supported by CONACyT through the Project No. A1- S-33202 (Ciencia Basica).
Besides, R. J. H.-P. is also funded by Ciencia de Frontera 2021-2042 and
Sistema Nacional de Investigadores from CONACyT.

## References

  * [1] J. C. Collins, D. E. Soper and G. F. Sterman, _Factorization of Hard Processes in QCD_ , _Adv. Ser. Direct. High Energy Phys._ 5 (1989) 1–91, [hep-ph/0409313].
  * [2] European Muon collaboration, J. Ashman et al., _An Investigation of the Spin Structure of the Proton in Deep Inelastic Scattering of Polarized Muons on Polarized Protons_ , _Nucl. Phys. B_ 328 (1989) 1.
  * [3] HERMES collaboration, A. Airapetian et al., _Precise determination of the spin structure function g(1) of the proton, deuteron and neutron_ , _Phys. Rev. D_ 75 (2007) 012007, [hep-ex/0609039].
  * [4] COMPASS collaboration, V. Y. Alexakhin et al., _The Deuteron Spin-dependent Structure Function g1(d) and its First Moment_ , _Phys. Lett. B_ 647 (2007) 8–17, [hep-ex/0609038].
  * [5] D. de Florian, R. Sassot, M. Stratmann and W. Vogelsang, _Global Analysis of Helicity Parton Densities and Their Uncertainties_ , _Phys. Rev. Lett._ 101 (2008) 072001, [0804.0422].
  * [6] D. de Florian, R. Sassot, M. Stratmann and W. Vogelsang, _Extraction of Spin-Dependent Parton Densities and Their Uncertainties_ , _Phys. Rev. D_ 80 (2009) 034030, [0904.3821].
  * [7] D. F. Rentería-Estrada, R. J. Hernández-Pinto and G. F. R. Sborlini, _Analysis of the Internal Structure of Hadrons Using Direct Photon Production_ , _Symmetry_ 13 (2021) 942, [2104.14663].
  * [8] D. de Florian and G. F. R. Sborlini, _Hadron plus photon production in polarized hadronic collisions at next-to-leading order accuracy_ , _Phys. Rev. D_ 83 (2011) 074022, [1011.0486].
  * [9] S. Frixione, _Isolated photons in perturbative QCD_ , _Phys. Lett. B_ 429 (1998) 369–374, [hep-ph/9801442].
  * [10] S. Frixione, Z. Kunszt and A. Signer, _Three jet cross-sections to next-to-leading order_ , _Nucl. Phys. B_ 467 (1996) 399–442, [hep-ph/9512328].
  * [11] A. Buckley, J. Ferrando, S. Lloyd, K. Nordström, B. Page, M. Rüfenacht et al., _LHAPDF6: parton density access in the LHC precision era_ , _Eur. Phys. J. C_ 75 (2015) 132, [1412.7420].
  * [12] A. V. Manohar, P. Nason, G. P. Salam and G. Zanderighi, _The Photon Content of the Proton_ , _JHEP_ 12 (2017) 046, [1708.01256].
  * [13] L. Buonocore, P. Nason, F. Tramontano and G. Zanderighi, _Leptons in the proton_ , _JHEP_ 08 (2020) 019, [2005.06477].
  * [14] D. de Florian, R. Sassot, M. Epele, R. J. Hernández-Pinto and M. Stratmann, _Parton-to-Pion Fragmentation Reloaded_ , _Phys. Rev. D_ 91 (2015) 014035, [1410.6027].
  * [15] D. Rentería-Estrada et al., _Reconstructing partonic kinematics at colliders with machine learning_ , to be published.


# Enforcing nonholonomic constraints in Aerobat,
a roosting flapping wing model

Eric Sihite, and Alireza Ramezani, SiliconSynapse Laboratory1 1ECE Department,
Northeastern University, Boston, MA, USA.emails: {e.sihite,
<EMAIL_ADDRESS>
###### Abstract

Flapping wing flight is a challenging dynamical problem and is also a very
fascinating subject to study in the field of biomimetic robotics. A Bat, in
particular, has a very articulated armwing mechanism with high degrees-of-
freedom and flexibility which allows the animal to perform highly dynamic and
complex maneuvers, such as upside-down perching. This paper presents the
derivation of a multi-body dynamical system of a bio-inspired bat robot called
Aerobat which captures multiple biologically meaningful degrees-of-freedom for
flapping flight that is present in biological bats. Then, the work attempts to
manifest closed-loop aerial body reorientation and preparation for landing
through the manipulation of inertial dynamics and aerodynamics by enforcing
nonholonomic constraints onto the system. The proposed design paradigm assumes
for rapidly exponentially stable controllers that enforce holonomic
constraints in the joint space of the model. A model and optimization-based
nonlinear controller is applied to resolve the joint trajectories such that
the desired angular momentum about the roll axis is achieved.

##  I Introduction

Figure 1: The simulated upside-down perching maneuver of Aerobat.

The key hypothesis this work tries to inspect is: ”The feasibility of
performing flip turns, a well known attribute of a bat-like landing maneuver,
through the manipulation of inertial dynamics while including aerodynamic
forces.” Bats (and birds) possess no energy-hungry motors widely used in
flying robots for thrust vectoring yet they are more capable than any of these
systems when agility and energy efficiency of flight are concerned. Flying
vertebrates apply the combination of inertial dynamics and aerodynamics
manipulations to showcase extremely agile maneuvers. Unlike rotary- and fixed-
wing systems wherein aerodynamic surface (i.e., ailerons, rudders, propellers,
etc.) come with the sole purpose of aerodynamic force adjustment, the wings
(also called appendages) in birds and bats possess more sophisticated roles.
It is known that birds perform zero-angular-momentum turns by making
differential adjustments (e.g., collapsing armwings) in the inertial forces
led by one wing versus the other. Bats apply a similar mechanism to perform
sharp banking turns [1].

Among these maneuvers, landing (or perching), which flying vertebrates do it
in one way or another for a variety of reasons (e.g., transition to walking,
resting on a perch, hanging from the ceiling of a cave, etc.), is an
interesting maneuver to take inspiration from for aerial robot designs.
Perching birds rotate their wings so that the aerodynamic drag is increased by
creating a high-pressure region inside of the wings and a low-pressure region
behind the wings. This brings the wings to a stalled condition at which point
the generated lift is equal to zero and the animal falls naturally while
employing the legs as a landing gear. Bats do it in a radically different way.
After the self-created stalled condition, they manifest an acrobatic heels-
above-head maneuver that involves catapulting the lower body in a similar way
that a free style swimmer flip turns. Perching insects and birds have been the
source of inspiration and bio-mimicry of them has led to interesting robot
designs in recent years [2, 3, 4]. Remarkably, the bio-mimicry of bat-like
landing is overlooked mainly because not only the aerodynamics adjustments are
involved but also unique design provisions are required to allow for the
manipulation of inertial dynamics.

A bat-style landing maneuver is extremely rich in dynamics and control and its
characteristics are overlooked. Much of attention has been paid to simpler
dynamics such as hovering and straight flight. While mathematical models of
insect-style, rotary- and fixed-wing robots of varying size and complexity are
relatively well developed, models of airborne, fluidic-based vertebrates
locomotion remain largely open due to the complex body articulation involved
in their flight. The mainstream school of thought inspired by insect flight
has conceptualized wing as a mass-less, rigid structure, which is nearly
planar and translates – as a whole or in two-three rigid parts – through space
[5, 6]. In this view, wings possess no inertial effect, are fast that yield
two-time-scale dynamics, permit quasi-static external force descriptions, and
tractable dynamical system. Unfortunately, these paradigms fail to provide
insight into airborne, vertebrate locomotion and an ingredient of a more
complete and biologically realistic model is missing, that is, the
manipulation of inertial dynamics. The manipulation of inertial dynamics is an
under-appreciated aspect in existing paradigms.

The objective of this work is to manifest closed-loop aerial body
reorientation and preparation for landing through the manipulation of inertial
dynamics and aerodynamics by enforcing nonholonomic constraints onto the
system. The proposed design paradigm assumes for rapidly exponentially stable
controllers that enforce holonomic constraints in the joint space of the
model. Then, a model and optimization-based nonlinear controller is applied to
resolve the joint trajectories such that the desired angular momentum about
the roll axis is achieved. First, a brief overview of our motivation will be
presented followed by a model description of the landing maneuver. Then, a
dynamic model and optimization problem will be derived in detail and the
preliminary simulation results will be reported at the end followed by final
remarks and conclusion.

Figure 2: The bio-inspired bat robot, called Aerobat, uses monolithically
fabricated rigid and soft arm-armwing structure. The white material of the
armwing is rigid while the black material is a flexible living hinge. The
flexible PCB forms the wing membrane which helps reducing the overall weight
of the robot and the kinetic sculpture forms a compliant linkage mechanism
that can conform and articulate flapping motion.

Fig. 2 shows the concept design of our bio-inspired bat robot, called Aerobat
[7], which is a continuation to our previous work in [8, 9, 10, 11, 12, 13,
14]. The control and dynamic modeling of this robot were also investigated in
[15, 16, 17, 18]. This bioinspired robot is designed to be very lightweight
and the flexible armwing mechanism is designed to mimic some of the
biologically meaningful degrees-of-freedom (DoF) in a bat’s flapping gait.
This flexible bat armwing structure, called the kinetic sculpture, uses both
rigid and flexible materials that is monolithically fabricated using PolyJet
3D printing technology. This armwing is articulated through a series of four-
bar linkages and crank mechanism which is driven by a single motor to
articulate the wing expansion and retraction during downstroke and upstroke
respectively. This results in maximum wingspan and lift during the downstroke
while reducing the negative lift during the upstroke, therefore forming an
efficient flapping gait. We have also developed a launching landing apparatus
called Harpoon as shown in [19] which will be used in the actual robot once it
has successfully performed the upside-down maneuver to latch onto the perching
location.

The motivation of developing Aerobat is to mimic the complex and highly
articulated natural bat armwing by deforming the armwing morphology and
achieve a varying range of motion that mimics the flapping motion of a
biological bat. This means of articulation facilitates control through
morphological computation [20], where a simple control action can actuate a
very complex motion or trajectory which is suitable for mimicking the complex
wing articulation of an actual bat. This motivates us to develop a dynamic
model for simulations which can be used to develop the stabilizing controller
for this bat robot. As part of the future work of the current armwing design,
we will extend the biologically meaningful degrees-of-freedom of the armwing
to include the feathering and mediolateral movements. This work aims to find
the optimal gait for this armwing structure to achieve stable flight and other
complex maneuvers, such as the upside-down perching maneuver.

##  II System Dynamic Modeling

The Aerobat system can be modeled as five rotating bodies which are attached
to one another with joints or hinges, as shown in Fig. 3. This armwing
mechanism is designed to follow the biologically meaningful degrees-of-freedom
(DoF) of a bat’s armwing, as illustrated in Fig. 4, where these DoFs are
important for flapping wing flight. The shoulder joint plunge angle
$(\theta_{p})$ controls the wing upstroke and downstroke motion which forms
the core flapping motion. The elbow extension/flexion angle $(\theta_{e})$
expands the wing during the downstroke and retracts it during the upstroke
motion, which improves efficiency by reducing the negative lift during the
upstroke. The mediolateral motion $(\theta_{m})$ extends the wings forward and
the feathering motion $(\theta_{f})$ rotates the wing surface plane with
respect to the arm which has an effect of changing the angle of attack.

Figure 3: Aerobat model using five rotating bodies where joints
$J_{i},i=\\{1,\dots,5\\}$ represent the center of rotations, the linear
position $\bm{p}$ represents the center of mass, and $\bm{l}$ represents the
length vectors relevant to the Aerobat model conformation. Figure 4: Aerobat
biologically meaningful degrees-of-freedom angles: armwing plunging
($\theta_{p}$), mediolateral ($\theta_{m}$), elbow extension ($\theta_{e}$),
and feathering ($\theta_{f}$). $\theta_{p}$ and $\theta_{m}$ are defined with
respect to the body frame while $\theta_{e}$ and $\theta_{f}$ are defined with
respect to the arm frame. The superscript $L$ and $R$ are used to represent
left and right wing angles respectively.

###  II-A Euler-Lagrangian Dynamic Formulation

Let a vector with a superscript notation represents the vector defined in a
non-inertial coordinate frame and the vector without superscript is defined
about the inertial frame, e.g. $\bm{x}^{B}$ is the vector $\bm{x}$ about frame
$B$. The frames of references is illustrated in Fig. 3 and the coordinate
frame rotation of the five bodies can be defined as follows:

$\displaystyle\bm{x}$ $\displaystyle=R_{B}\,\bm{x}^{B},\quad$
$\displaystyle\bm{x}^{B}$
$\displaystyle=R_{A_{L}}\,\bm{x}^{A_{L}}=R_{A_{R}}\,\bm{x}^{A_{R}}$ (1)
$\displaystyle\bm{x}^{A_{L}}$ $\displaystyle=R_{W_{L}}\,\bm{x}^{W_{L}},\quad$
$\displaystyle\bm{x}^{A_{R}}$ $\displaystyle=R_{W_{R}}\,\bm{x}^{W_{R}},$

where $R_{B}$ is the body rotation matrix about the inertial frame,
$R_{A_{L}}$ and $R_{A_{R}}$ are the rotation matrix of the left and right arm
respectively about the body frame, and $R_{W_{L}}$ and $R_{W_{R}}$ are the
rotation matrix of the left and right wing respectively about their respective
arm. The corresponding angular velocities for these rotation matrices must
also be represented in the appropriate coordinate frames. Then we have the
following rotation matrix and angular velocity pairs:
$(R_{B},\bm{\omega}_{B})$, $(R_{A_{L}},\bm{\omega}_{A_{L}}^{B})$,
$(R_{A_{R}},\bm{\omega}_{A_{R}}^{B})$,
$(R_{W_{L}},\bm{\omega}_{W_{L}}^{A_{L}})$, and
$(R_{W_{R}},\bm{\omega}_{W_{R}}^{A_{R}})$. Let $\theta$ angles be the
biologically meaningful flapping angles, as illustrated in Fig. 4, and the
superscript $L$ and $R$ represents the left and right wing joint angles
respectively. Then the left armwing rotation matrices are defined as follows:

$\displaystyle R_{A_{L}}$
$\displaystyle=R_{z}(\theta^{L}_{m})\,R_{x}(\theta^{L}_{p}),\,$ $\displaystyle
R_{W_{L}}$ $\displaystyle=R_{x}(\theta^{L}_{e})\,R_{z}(\theta^{L}_{f}),$ (2)

where $R_{x}(\theta)$ and $R_{z}(\theta)$ are the rotational matrix about $x$
and $z$ axis respectively. The left armwing angular velocities are defined as
follows:

$\displaystyle\bm{\omega}_{A_{L}}^{B}$
$\displaystyle=\begin{bmatrix}0,0,\dot{\theta}_{m}^{L}\end{bmatrix}^{\top}+R_{z}(\theta^{L}_{m})\begin{bmatrix}\dot{\theta}_{p}^{L},0,0\end{bmatrix}^{\top}+\bm{\omega}_{B}^{B}$
(3) $\displaystyle\bm{\omega}_{W_{L}}^{A_{L}}$
$\displaystyle=\begin{bmatrix}\dot{\theta}_{e}^{L},0,0\end{bmatrix}^{\top}+R_{x}(\theta^{L}_{e})\begin{bmatrix}0,0,\dot{\theta}_{f}^{L}\end{bmatrix}^{\top}+\bm{\omega}_{A_{L}}^{A_{L}}.$

The right wing derivations can be derived in a similar fashion. Therefore, for
the rest of the paper, only the left wing components will be derived if the
right side also follow a similar derivation.

As shown in Fig. 3, let $\bm{p}_{B}$ be the linear position of the center of
mass of a body, $\bm{l}_{Lj}$ and $\bm{l}_{Rj}$, $j=\\{1,2,3\\}$, be the
length vectors which represent the Aerobat mechanism morphology that are
constant with respect to their local frame of reference. Then the linear
position of the center of mass of the left armwing can be derived as follows:

$\displaystyle\bm{p}_{A_{L}}$
$\displaystyle=\bm{p}_{B}+R_{B}\,\bm{l}_{L1}^{B}+\tfrac{1}{2}\,R_{B}\,R_{A_{L}}\,\bm{l}_{L2}^{A_{L}}$
(4) $\displaystyle\bm{p}_{W_{L}}$
$\displaystyle=\bm{p}_{A_{L}}+\tfrac{1}{2}\,R_{B}\,R_{A_{L}}\,\bm{l}_{L2}^{A_{L}}+R_{B}\,R_{A_{L}}\,R_{W_{L}}\,\bm{l}_{L3}^{W_{L}},$

The linear velocity of the center of mass can be derived from (4) by
differentiating the linear positions with respect to time. Note that the
linear positions and velocities are defined with respect to the inertial
frame.

The kinetic and potential energy of the system can be derived as follows:

$\displaystyle T$
$\displaystyle=\sum_{F\in\mathcal{F}}\left(m_{F}\,\dot{\bm{p}}_{F}^{\top}\,\dot{\bm{p}}_{F}+(\bm{\omega}_{F}^{F})^{\top}\,\hat{I}_{F}\,\bm{\omega}_{F}^{F}\right)\frac{1}{2}$
(5) $\displaystyle U$
$\displaystyle=\sum_{F\in\mathcal{F}}m_{F}\,[0,0,g]\,\bm{p}_{F},$

where $\mathcal{F}=\\{B,A_{L},A_{R}.W_{L},W_{R}\\}$ is the set containing the
frame of references, $m_{F}$ and $\hat{I}_{F}$ are the mass and inertia matrix
of the corresponding body respectively. $\hat{I}_{F}$ is defined about the
local frame of reference which is diagonal and constant. Then the Lagrangian
of the system, $L=T-U$, can be used to derive the equation of motion.

The body rotation $(R_{B},\omega_{B})$ is derived using the modified Euler-
Lagrangian formulation for a rotation in SO(3). This formulation is not
susceptible to gimbal lock which might happen if we use Tait-Bryan angles
during the upside-down maneuver. The modified Euler-Lagrange equation for
rotation in SO(3) can be derived by using Hamilton’s principle [21], which has
the following form:

$\begin{gathered}\frac{d}{dt}\frac{\partial
L}{\partial\bm{\omega}_{B}^{B}}+\bm{\omega}_{B}^{B}\times\frac{\partial
L}{\partial\bm{\omega}_{B}^{B}}+\sum^{3}_{j=1}\bm{r}_{B,j}\times\frac{\partial
L}{\partial\bm{r}_{Bj}}=\bm{\tau}_{B}^{B}\\\
\dot{R}_{B}=R_{B}\,S(\bm{\omega}_{B}^{B}),\end{gathered}$ (6)

where $S(\cdot)$ is a skew operator,
$R_{B}^{\top}=[\bm{r}_{B1},\bm{r}_{B2},\bm{r}_{B3}]$ and $\bm{\tau}_{B}^{B}$
is the non-conservative torque about the generalized coordinate
$\bm{\omega}_{B}^{B}$. The equation of motion of the remaining states can be
solved by using the Euler-Lagrange equation:

$\begin{gathered}\bm{\theta}_{L}=[\theta_{p}^{L},\theta_{m}^{L},\theta_{e}^{L},\theta_{f}^{L}],\qquad\bm{\theta}_{R}=[\theta_{p}^{R},\theta_{m}^{R},\theta_{e}^{R},\theta_{f}^{R}],\\\
\bm{q}_{e}=[\bm{p}_{B}^{\top},\bm{\theta}_{L}^{\top},\bm{\theta}_{R}^{\top}]^{\top},\qquad\frac{d}{dt}\frac{\partial
L}{\partial\dot{\bm{q}}_{e}}-\frac{\partial
L}{\partial\bm{q}_{e}}=\bm{u}_{e},\end{gathered}$ (7)

where $\bm{u}_{e}$ is the non-conservative force about the generalized
coordinate $\bm{q}_{e}$. Combining (6) and (7), the equation of motion can be
formulated into the following form:

$\begin{gathered}\bm{q}=[\bm{r}_{B}^{\top},\bm{p}_{B}^{\top},\bm{\theta}_{L}^{\top},\bm{\theta}_{R}^{\top}]^{\top},\qquad\bm{q}_{d}=[\bm{\omega}_{B}^{\top},\dot{\bm{p}}_{B}^{\top},\dot{\bm{\theta}}_{L}^{\top},\dot{\bm{\theta}}_{R}^{\top}]^{\top}\\\
M\,\dot{\bm{q}}_{d}+\bm{h}=B_{a}\,\bm{u}_{a}+B_{m}\,\bm{u}_{m},\end{gathered}$
(8)

where $\bm{r}_{B}$ is the rotation matrix $R_{B}$ concatenated into a vector
form, $\bm{u}_{a}$ is the generalized aerodynamic forces and torque.
$\bm{u}_{m}$ is the generalized motor torque acting on the armwing joints
which is selected to directly actuate the joints angles $\bm{\theta}_{L}$ and
$\bm{\theta}_{R}$:

$\displaystyle\bm{u}_{m}$
$\displaystyle=[\tau_{p}^{L},\tau_{m}^{L},\tau_{e}^{L},\tau_{f}^{L},\tau_{p}^{R},\tau_{m}^{R},\tau_{e}^{R},\tau_{f}^{R}]^{\top}$
(9) $\displaystyle B_{m}$ $\displaystyle=[0_{8\times 6},I_{8\times
8}]^{\top},$

where $\tau$ represents the torque acting on each joints.

###  II-B Aerodynamic Modeling

Figure 5: The coordinate system of the wing plate. $w_{r}$ and $w_{c}$ are the
wing span and chord length respectively. $\hat{r}\in[0,1]$ and
$\hat{c}\in[-0.5,0.5]$ are the unitless variables to represent a position on
the wing surface. $\bm{l}_{f_{L}}$ is the length vector to the aerodynamic
force $\bm{f}_{a_{L}}$ acting on the left wing.

The aerodynamic forces can be derived using the virtual displacement defined
at the position of the applied aerodynamics force $\bm{f}_{a_{L}}$. Let
$\bm{p}_{f_{L}}$ be the position of the applied aerodynamics force
$\bm{u}_{L}$ on the left wing, defined as follows:

$\displaystyle\bm{p}_{f_{L}}$
$\displaystyle=\bm{p}_{A_{L}}+\tfrac{1}{2}\,R_{B}\,R_{A_{L}}\,\bm{l}_{L2}^{A_{L}}+R_{B}\,R_{A_{L}}\,R_{W_{L}}\,\bm{l}_{f_{L}}^{W_{L}}$
(10)

where $\bm{l}_{f_{L}}^{W_{L}}$ is the length vector from the left wing elbow
joint to where $\bm{f}_{a_{L}}$ is applied, as illustrated in Fig. 5. Define
the position on the wing surface as follows:

$\bm{l}_{f_{L}}=\hat{r}\,w_{r}\,\bm{e}_{r}+\hat{c}\,w_{c}\,\bm{e}_{c}$ (11)

where $w_{r}$ and $w_{c}$ are the wing span and chord length respectively. The
position is represented using the unitless variables $\hat{r}\in[0,1]$ and
$\hat{c}\in[-0.5,0.5]$, as illustrated in Fig. 5. For the left wing,
$\\{\bm{e}_{c},\bm{e}_{n},\bm{e}_{r}\\}=\\{\bm{e}_{1}^{W_{L}},-\bm{e}_{2}^{W_{L}},\bm{e}_{3}^{W_{L}}\\}$.
Then $\bm{p}_{f_{L}}$ and $\bm{p}_{f_{R}}$ can be represented using the
coordinate $(\hat{c},\hat{r})$, e.g. $\bm{p}_{f_{L}}(\hat{c},\hat{r})$.

The generalized forces can be derived by as follows:

$Q_{j}=\bm{f}_{a_{L}}^{\top}\,\frac{\partial\dot{\bm{p}}_{f_{L}}}{\partial
q_{d,j}}=\left(\frac{\partial\dot{\bm{p}}_{f_{L}}}{\partial
q_{d,j}}\right)^{\top}R_{B}\,R_{A_{L}}\,R_{W_{L}}\,\bm{f}_{a_{L}}^{W_{L}},$
(12)

where $q_{d,j}$ is the $j$’th component of the vector $\bm{q}_{d}$ defined in
(8). Since this is a flapping wing robot, the airfoil speed is variable across
the wingspan which means that we need to integrate the aerodynamic forces
across the wing surface. It is possible to represent the $B_{a}$ matrix
without $\bm{l}_{f_{L}}^{W_{L}}$ by calculating the aerodynamic torque
$\bm{\tau}_{a_{L}}^{A_{L}}$ about the plate joint. Evaluate the component of
(12) that has the $\bm{l}_{f_{L}}^{W_{L}}$ term using the following method:

$\displaystyle Q_{\tau_{j}}$
$\displaystyle=\bm{f}^{\top}_{a_{L}}\frac{\partial}{\partial
q_{d,j}}\left(R_{B}\,R_{A_{L}}\,(\bm{\omega}_{W_{L}}^{A_{L}}\times
R_{W_{L}}\,\bm{l}_{f_{L}}^{W_{L}})\right)$ (13)
$\displaystyle=\frac{\partial\bm{\omega}_{W_{L}}^{A_{L}}}{\partial
q_{d,j}}^{\top}\mathrm{cof}(R_{W_{L}})(\bm{l}_{f_{L}}^{W_{L}}\times\bm{f}_{a_{L}}^{W_{L}}),$

where
$\bm{\tau}_{a_{L}}^{W_{L}}=\bm{l}_{f_{L}}^{W_{L}}\times\bm{f}_{a_{L}}^{W_{L}}$
is the aerodynamic torque acting on the left wing plate about its joint. The
right wing also follows a similar derivation, then combining (12) and (13) for
both sides of the wing forms the matrix $B_{a}\in\mathcal{R}^{14\times 12}$
with the input vector
$\bm{u}_{a}=[\bm{f}_{a_{L}}^{W_{L}\top},\bm{\tau}_{a_{L}}^{W_{L}\top},\bm{f}_{a_{R}}^{W_{R}\top},\bm{\tau}_{a_{R}}^{W_{R}\top}]^{\top}$
for the aerodynamic forces and torques acting on the wings which must be
integrated about the wing surface. The integration can be solved using the
following equations:

$\displaystyle\bm{l}_{f_{L}}(\hat{r})$
$\displaystyle=\hat{r}\,w_{r}\,\bm{e}_{r}+0.25\,w_{c}\,\bm{e}_{c}$ (14)
$\displaystyle\alpha(\hat{r})$
$\displaystyle=\mathrm{atan2}(\bm{e}_{n}^{\top}\,\bm{v}_{w}(\hat{r}),\bm{e}_{c}^{\top}\,\bm{v}_{w}(\hat{r}))$
$\displaystyle f_{l}(\hat{r})$
$\displaystyle=\rho\,w_{c}w_{r}|\bm{v}_{w}|^{2}C_{L}(\alpha(\hat{r}))/2$
$\displaystyle f_{d}(\hat{r})$
$\displaystyle=\rho\,w_{c}w_{r}|\bm{v}_{w}|^{2}C_{D}(\alpha(\hat{r}))\,\mathrm{sgn}(\bm{e}_{c}^{\top}\,\bm{v}_{w})/2$
$\displaystyle\bm{f}_{a_{L}}$
$\displaystyle=\smallint_{0}^{1}(f_{l}(\hat{r})\bm{e}_{n}-f_{d}(\hat{r})\bm{e}_{c})d\hat{r}$
$\displaystyle\bm{\tau}_{a_{L}}$
$\displaystyle=\smallint_{0}^{1}\bm{l}_{f_{L}}(\hat{r})\times(f_{l}(\hat{r})\bm{e}_{n}-f_{d}(\hat{r})\bm{e}_{c})d\hat{r},$

where $\rho$ is the air density, $\bm{v}_{w}$ is the airfoil velocity at
$\bm{l}_{f_{L}}$, $\bm{e}_{n}$ and $\bm{e}_{c}$ are the unit vectors as shown
in Fig. 5. $C_{L}$ and $C_{D}$ are the lift and drag coefficients defined
below:

$\displaystyle C_{L}(\alpha)$
$\displaystyle=0.225+1.58\sin(2.13\alpha-7.2^{\circ})$ (15) $\displaystyle
C_{D}(\alpha)$ $\displaystyle=1.92-1.55\cos(2.04\alpha-9.82^{\circ}),$

where $\alpha$ is the angle of attack in degrees. These coefficients are used
in [22] for a flapping wing MAV based on a fruit fly. The integration about
the wing span (integrate about $\hat{c}$) and the right armwing side follow a
similar derivation. Calculating the sum of all of this forces and torques
forms the aerodynamic actuation vector $\bm{u}_{a}$.

##  III Gait Optimization Framework

TABLE I: List of the parameters and their values Param. | Value | Param. | Value
---|---|---|---
$m_{B}$ | 5 g | $\bm{l}_{L1}$ |  $[0,25,25]^{\top}$ mm
$m_{A}$ | 0.35 g | $\bm{l}_{L2}$ |  $[0,0,50]^{\top}$ mm
$m_{W}$ | 5.6 g | $\bm{l}_{L3}$ |  $[0,0,150]^{\top}$ mm
$\hat{I}_{B,x}$ | 0.625 g.cm2 | $\bm{l}_{R1}$ |  $[0,-25,25]^{\top}$ mm
$\hat{I}_{B,y}$ | 3.65 g.cm2 | $\bm{l}_{R2}$ |  $[0,0,50]^{\top}$ mm
$\hat{I}_{B,z}$ | 3.65 g.cm2 | $\bm{l}_{R3}$ |  $[0,0,150]^{\top}$ mm
$\hat{I}_{A,x}$ | 0.147 g.cm2 | $\hat{I}_{W,x}$ | 1.05 g.cm2
$\hat{I}_{A,y}$ | 0.147 g.cm2 | $\hat{I}_{W,y}$ | 2.11 g.cm2
$\hat{I}_{A,z}$ | 0.040 g.cm2 | $\hat{I}_{W,z}$ | 2.11 g.cm2
$w_{c}$ | 150 mm | $\rho$ | 1 kg/m3
$w_{r}$ | 150 mm | $\Omega$ | 10 Hz
Figure 6: Open loop simulation for zero angular momentum gait using the
optimized parameters in (28). The mean of the angular momentum is
approximately zero and the linear velocity is approximately constant.

This section outlines the optimization framework used in order to find the
gait which stabilizes the robot’s angular momentum and another gait for the
upside-down perching maneuver. In order to define the angular momentum about
the robot’s center of mass, define the system center of mass:

$\bm{p}_{CoM}=\frac{1}{\sum_{F\in\mathcal{F}}m_{F}}\sum_{F\in\mathcal{F}}\bm{p}_{F}\,m_{F}.$
(16)

Then calculate the inertial angular momentum of the bodies about this center
of mass:

$\displaystyle\bm{\Pi}_{B}=$
$\displaystyle\,R_{B}\,\hat{I}_{B}\,\bm{\omega}_{B}^{B}-m_{B}\,(\bm{p}_{B}-\bm{p}_{CoM})\times\dot{\bm{p}}_{B}$
(17) $\displaystyle\bm{\Pi}_{A_{L}}=$
$\displaystyle\,R_{B}\,R_{A_{L}}\,\hat{I}_{A_{L}}\,R_{A_{L}}^{\top}\,\bm{\omega}^{B}_{A_{L}}$
$\displaystyle\qquad-
m_{A}\,(\bm{p}_{A_{L}}-\bm{p}_{CoM})\times\dot{\bm{p}}_{A_{L}}$
$\displaystyle\bm{\Pi}_{W_{L}}=$
$\displaystyle\,R_{B}\,R_{A_{L}}\,R_{W_{L}}\,\hat{I}_{W_{L}}\,R_{W_{L}}^{\top}\,\bm{\omega}^{A_{L}}_{W_{L}}$
$\displaystyle\qquad-
m_{W}\,(\bm{p}_{W_{L}}-\bm{p}_{CoM})\times\dot{\bm{p}}_{W_{L}}.$

Then the total angular momentum of the system about the center of mass is:

$\bm{\Pi}=\bm{\Pi}_{B}+\bm{\Pi}_{A_{L}}+\bm{\Pi}_{A_{R}}+\bm{\Pi}_{W_{L}}+\bm{\Pi}_{W_{R}}.$
(18)

In order to design the optimization problem to find the best gait for a stable
flapping flight, we introduce a nonholonomic constraint at the wing joints
acceleration so we can impose a specific gait for the armwing joints in the
simulation. This can be implemented by using a Lagrangian multiplier to impose
the acceleration constraint on (8), as shown below

$\displaystyle M\,\dot{\bm{q}}_{d}+\bm{h}$
$\displaystyle=B_{a}\,\bm{u}_{a}+B_{m}\,\bm{u}_{m}+J_{c}^{\top}\,\bm{\lambda}$
(19) $\displaystyle J_{c}\,\dot{\bm{q}}_{d}$
$\displaystyle=\ddot{\bm{\theta}}_{c},$

where
$\ddot{\bm{\theta}}_{c}=[\ddot{\bm{\theta}}_{c,L}^{\top},\ddot{\bm{\theta}}_{c,R}^{\top}]^{\top}$
is the joint acceleration constraint, $\bm{\theta}_{c,L}$ and
$\bm{\theta}_{c,R}$ are the constraints for the left and right joint angles
respectively, and $J_{c}=[0_{8\times 6},I_{8\times 8}]$. Then the Lagrangian
multiplier $\bm{\lambda}$ can be solved as follows:

$\bm{\lambda}=(J_{c}\,M^{-1}\,J_{c}^{\top})^{-1}\,\left(J_{c}\,M^{-1}\,(\bm{h}-\bm{u})-\ddot{\bm{\theta}}_{c}\right),$
(20)

where $\bm{u}=B_{a}\bm{u}_{a}+B_{m}\bm{u}_{m}$ is the combined forces acting
on the system. The $\bm{u}_{m}$ value is irrelevant as the consequence of
using the nonholonomic constraint on the joint angles.

###  III-A Zero Angular Momentum Gait Optimization

The following open loop trajectory is used for the the left armwing:

$\displaystyle\bm{\theta}_{c,L}(t)=\begin{bmatrix}A_{p}\,\cos(\Omega\,t)+\accentset{\rule{2.79996pt}{0.8pt}}{\theta}_{p}\\\
A_{m}\,\cos(\Omega\,t+\phi_{m})+\accentset{\rule{2.79996pt}{0.8pt}}{\theta}_{m}\\\
A_{e}\,\cos(\Omega\,t+\phi_{e})+\accentset{\rule{2.79996pt}{0.8pt}}{\theta}_{e}\\\
A_{f}\,\cos(\Omega\,t+\phi_{f})+\accentset{\rule{2.79996pt}{0.8pt}}{\theta}_{f}\end{bmatrix},$
(21)

where $\Omega$ is the flapping frequency,
$\accentset{\rule{2.79996pt}{0.8pt}}{\theta}$ is the mean joint angle, $A$ is
the amplitude, and $\phi$ is the phase shift of the other joints with respect
to the plunging motion. The right armwing trajectory is symmetric to the left
armwing which results in a symmetric aerodynamic force and stable in the roll
and yaw motion. Then, there are 11 parameters to optimize using the following
optimization problem:

$\displaystyle\min_{\bm{k}_{1}}$
$\displaystyle\textstyle\sum_{i}(\mathrm{trace}(\bm{z}_{i}^{\top}\,Q\,\bm{z}_{i}))$
(22) subject to
$\displaystyle\,\,\bm{k}_{1,min}\leq\bm{k}_{1}\leq\bm{k}_{1,max},$
$\displaystyle\bm{k}_{1}$
$\displaystyle=[\accentset{\rule{2.79996pt}{0.8pt}}{\theta}_{p},\accentset{\rule{2.79996pt}{0.8pt}}{\theta}_{m},\accentset{\rule{2.79996pt}{0.8pt}}{\theta}_{e},\accentset{\rule{2.79996pt}{0.8pt}}{\theta}_{f},A_{p},A_{m},A_{e},A_{f},\phi_{m},\phi_{e},\phi_{f}]^{\top}$
(23) $\displaystyle\bm{z}$
$\displaystyle=[\bm{\Pi}^{\top},\,\dot{p}_{B,z}]^{\top},$

where $Q$ is a weighting matrix. The optimization cost function is the
summation of the $\mathrm{trace}(\bm{z}_{i}^{\top}\,Q\,\bm{z}_{i})$ in a
numerical simulation of the wing dynamics where $i$ is the $i$’th time step of
the numerical simulation. This cost function is designed to find the gait that
can steadily keep the robot’s angular momentum constant while also capable of
maintaining zero vertical velocity throughout the simulation.

###  III-B Upside-down Perching Optimization

The trajectory found in (22) will be used to find the upside-down perching
maneuver, which is done by introducing the appropriate force imbalances
through the wing joints articulation. Since the gait found in (22) is momentum
stable, we can simply superimpose a joint trajectory on top of this gait to
create the force imbalance. We choose to superimpose a simple offset
trajectory:

$\displaystyle\tfrac{d}{dt}\accentset{\rule{2.79996pt}{0.8pt}}{\bm{\theta}}(t)=\begin{cases}\bm{d}/T,&\text{if
$t\geq t_{0}$ and $t<t_{0}+T$}\\\ -\bm{d}/T,&\text{if $t\geq t_{0}+T$ and
$t<t_{0}+2T$}\\\ 0,&\text{otherwise}\end{cases},$ (24)

where
$\accentset{\rule{2.79996pt}{0.8pt}}{\bm{\theta}}=[\accentset{\rule{2.79996pt}{0.8pt}}{\theta}_{p},\accentset{\rule{2.79996pt}{0.8pt}}{\theta}_{m},\accentset{\rule{2.79996pt}{0.8pt}}{\theta}_{e},\accentset{\rule{2.79996pt}{0.8pt}}{\theta}_{f}]^{\top}$,
$\bm{d}=[d_{p},d_{m},d_{e},d_{f}]^{\top}$ is the optimization parameters which
generates a triangular offset function with amplitude $\bm{d}$, ramp up and
down period of $T$, and ramp startup time of $t_{0}$. To form a force
imbalance, we set the left and right wing offset trajectory to be equal, i.e.
$\tfrac{d}{dt}\accentset{\rule{2.79996pt}{0.8pt}}{\bm{\theta}}_{L}(t)=\tfrac{d}{dt}\accentset{\rule{2.79996pt}{0.8pt}}{\bm{\theta}}_{R}(t)=\tfrac{d}{dt}\accentset{\rule{2.79996pt}{0.8pt}}{\bm{\theta}}(t)$.
However, since the joint acceleration is the input to our constrained system,
we need to impose a bounded acceleration for a smooth velocity trajectory.
This is done by using a workaround where we close the loop by setting the
acceleration constraint as

$\ddot{\bm{\theta}}_{c}=-120(\bm{\theta}-\bm{\theta}_{r})-120(\dot{\bm{\theta}}-\dot{\bm{\theta}}_{r}),$
(25)

where $\bm{\theta}=[\bm{\theta}_{L}^{\top},\bm{\theta}_{R}^{\top}]^{\top}$,
and $\bm{\theta}_{r}$ is the joint angle trajectories from (21) superimposed
with the offset trajectory from (24). The constraint in (25) forms a simple PD
controller which should give us a bounded joint accelerations. Finally, the
objective function is to follow the trajectory of $\bm{\omega}^{B}_{B}$ that
undertakes a rolling motion. The rolling motion is defined using the following
function:

$\displaystyle\phi_{r}$
$\displaystyle=(\pi/2)\tanh(\eta(t)),\quad\eta(t)=(6/2T)(t-t_{0})-3$ (26)
$\displaystyle\dot{\phi}_{r}$
$\displaystyle=(\pi/2)(3/2T)(1-\tanh(\eta(t))^{2}),$

where $\phi_{r}$ is the target rolling angle. (26) forms a smooth ramp up of
$180^{\circ}$ in roll angle within a time span of $2T$ starting from $t_{0}$.
Then the optimization is setup as follows:

$\displaystyle\min_{\bm{k}_{2}}$
$\displaystyle\,\textstyle\sum_{i}(\omega^{B}_{B,x}-\dot{\phi}_{r})^{2},\qquad\bm{k}_{2}=[t_{0},\bm{d}^{\top}]^{\top}$
(27) subject to
$\displaystyle\,\,\bm{k}_{2,min}\leq\bm{k}_{2}\leq\bm{k}_{2,max},$

where $\omega^{B}_{B,x}$ is the first component of $\bm{\omega}_{B}^{B}$ which
represents the body angular velocity about the roll axis, and $i$ is the
simulation discrete time step between $t_{0}$ and $t_{0}+2T$.

##  IV Simulation and Analysis

Figure 7: Open loop simulation for the upside-down perching using the
optimized parameters from (28) and (29). The robot follows the desired roll
velocity but there is a significant overshoot during the slowing down period.
Figure 8: Closed loop simulation using PID controller to track the
trajectories found in the optimizations. The roll angle reached the -180∘ for
perching which indicates its feasibility. However, the controller failed to
fully track the desired trajectory and the simulation is unstable after the
upside-down maneuver. Figure 9: The Aerobat zero angular momentum flapping
motion. The flapping motion is accompanied by the wing elbow expansion,
mediolateral motion where the shoulder bends forwards, and the feathering
motion which rotates the wing about the elbow joint. During the upstroke, the
elbow joint is retracted and the feathering movement adjusts the angle of
attack which reduces negative lift.

The optimization was run using RK4 algorithm where the robot simply ran in an
open loop using the joint constraints described in (19) and (21) for 2
seconds. The robot is subjected to a constant wind speed
$\bm{v}_{air}=[-2,0,0]^{\top}$m/s about the inertial frame and the system
parameters are listed in Table I. We use the weighting matrix of
$Q=\mathrm{diag}([5,5,5,10^{-5}])$ which gives a cost function ratio
approximately equal between the momentum and the robot’s vertical velocity.
Using the interior-point algorithm, the optimization results in the following
trajectory parameters:

$\displaystyle\bm{k}_{1}=[$
$\displaystyle-75.3^{\circ},-16.2^{\circ},-50.7^{\circ},-9.2^{\circ},45.0^{\circ},17.3^{\circ},$
(28) $\displaystyle
27.2^{\circ},29.1^{\circ},-91^{\circ},-112^{\circ},-92.5^{\circ}]^{\top}.$

The resulting gait and the robot’s states can be seen in Fig. 6 and the
illustration of this gait is shown in Fig. 9. The robot has a limit cycle with
a stable angular and linear velocities with a mean pitch angle angle of
$55^{\circ}$. However, if the force symmetry is lost due to the numerical
error or other means, the robot starts to spin wildly and the open-loop gait
is incapable of stabilizing the robot.

The optimization to find upside-down perching maneuver in (27) utilizes the
gait in (28) to solve the optimization problem. We set $T=0.2$s which is 2
wingbeats at 10 Hz flapping frequency. Then we bound $t_{0}\in[1.0,1.1]$s such
that $t_{0}$ represents the start time of this maneuver within a wingbeat. The
optimizer results in the following parameters:

$\bm{k}_{2}=[1.0724\text{
s},-55.7^{\circ},0^{\circ},0^{\circ},-16.6^{\circ}]^{\top}.$ (29)

The $t_{0}$ of 1.0724 seconds indicates that the maneuver starts at 72.4% of a
wingbeat period. Interestingly, the mediolateral and elbow offsets are both
zero even though we did not constrain the optimizer this way. This indicates
that the plunging and feathering movements are the primary driving force for
the upside-down maneuver. Fig. 1 illustrates the gait found in this
optimization and Fig. 7 shows the simulation states of the robot using the
optimized parameters where it shows that the body’s Euler angles does a
$180^{\circ}$ turn in approximately 0.15 seconds. The body angular velocity
follows the target velocity for a brief period of time but then overshoots
during the ramping down period. The robot spins out of control after 1.3
seconds mark due to its upside-down state, but we can assume that it already
has perched onto something by then using our launching landing gear Harpoon
[19].

Now that we have identified all the trajectories to follow, we can then use a
simple PID controller to track these target trajectory and see if the result
matches what the optimizer predicted. We use PID controller gains of
$k_{p}=k_{d}=0.0012$ and $k_{i}=0.006$. Fig. 8 shows the controller
performance, where the robot achieved the upside-down perching maneuver but
very quickly turned unstable. Fig. 8 also shows that the controller does not
track the desired trajectory well, which is likely caused by the PID
controller’s incapability to track the desired trajectory in the presence of
the changes in aerodynamic forces and gravity. This result indicates that a
constant PID gain is not sufficient to track the trajectory and we need to
develop a better controller. Additionally, a momentum stabilizing controller
also need to be implemented so that the robot can perform the upside-down
maneuver more stably.

##  V Conclusions and Future Work

A dynamic model of a bat robot and the optimization framework using
nonholonomic constraint to find the gait which is momentum stable followed by
an upside-down perching maneuver is presented in this paper. The optimizer has
successfully found the parameters which fulfills our criterion and the
corresponding gaits for the stable flying and upside-down perching has been
found. However, the PID controller we used to track these trajectory is
insufficient and we need to improve the controller’s performance to track
these trajectory. For the future work, we can investigate the implementation
of a robust trajectory tracking controller using a rapidly exponentially
stable controllers, such as [23], and a nonlinear controller which actively
stabilizes the robot’s angular momentum. Additionally, we can design a kinetic
sculpture that incorporate these DoFs and investigate a morphological
computation framework that might be more natural for a bio-inspired flexible
armwing mechanism.

## References

  * [1] D. K. Riskin, A. Bergou, K. S. Breuer, and S. M. Swartz, “Upstroke wing flexion and the inertial cost of bat flight,” _Proceedings of the Royal Society B: Biological Sciences_ , vol. 279, no. 1740, pp. 2945–2950, 2012.
  * [2] A. A. Paranjape, S.-J. Chung, and J. Kim, “Novel dihedral-based control of flapping-wing aircraft with application to perching,” _IEEE Transactions on Robotics_ , vol. 29, no. 5, pp. 1071–1084, 2013.
  * [3] C. E. Doyle, J. J. Bird, T. A. Isom, J. C. Kallman, D. F. Bareiss, D. J. Dunlop, R. J. King, J. J. Abbott, and M. A. Minor, “An avian-inspired passive mechanism for quadrotor perching,” _IEEE/ASME Transactions on Mechatronics_ , vol. 18, no. 2, pp. 506–517, 2012.
  * [4] M. Graule, P. Chirarattananon, S. Fuller, N. Jafferis, K. Ma, M. Spenko, R. Kornbluh, and R. Wood, “Perching and takeoff of a robotic insect on overhangs using switchable electrostatic adhesion,” _Science_ , vol. 352, no. 6288, pp. 978–982, 2016.
  * [5] A. J. Bergou, L. Ristroph, J. Guckenheimer, I. Cohen, and Z. J. Wang, “Fruit flies modulate passive wing pitching to generate in-flight turns,” _Physical review letters_ , vol. 104, no. 14, p. 148101, 2010.
  * [6] T. L. Hedrick, B. Cheng, and X. Deng, “Wingbeat time and the scaling of passive rotational damping in flapping flight,” _Science_ , vol. 324, no. 5924, pp. 252–255, 2009.
  * [7] E. Sihite, P. Kelly, and A. Ramezani, “Computational structure design of a bio-inspired armwing mechanism,” _IEEE Robotics and Automation Letters_ , vol. 5, no. 4, pp. 5929–5936, 2020, in press.
  * [8] A. Ramezani, X. Shi, S. J. Chung, and S. Hutchinson, “Bat bot (b2), a biologically inspired flying machine,” in _Int. Conf. on Robotics and Automation (ICRA)_. IEEE, 2016, pp. 3219–3226.
  * [9] J. Hoff, A. Ramezani, S. J. Chung, and S. Hutchinson, “Synergistic design of a bio-inspired micro aerial vehicle with articulated wings.” in _Robotics: science and systems_ , 2016.
  * [10] J. Hoff, A. Ramezani, S. J. Chung, and S. Hutchinson, “Reducing versatile bat wing conformations to a 1-dof machine,” in _Conf. on Biomimetic and Biohybrid Systems_ , 2017, pp. 181–192.
  * [11] A. Ramezani, S. U. Ahmed, J. Hoff, S. J. Chung, and S. Hutchinson, “Describing robotic bat flight with stable periodic orbits,” in _Conf. on Biomimetic and Biohybrid Systems_ , 2017, pp. 394–405.
  * [12] A. Ramezani, S. J. Chung, and S. Hutchinson, “A biomimetic robotic platform to study flight specializations of bats,” _Science Robotics_ , vol. 2, no. 3, pp. Art–No, 2017.
  * [13] U. A. Syed, A. Ramezani, S. J. Chung, and S. Hutchinson, “From rousettus aegyptiacus (bat) landing to robotic landing: Regulation of cg-cp distance using a nonlinear closed-loop feedback,” in _Int. Conf. on Robotics and Automation (ICRA)_. IEEE, 2017, pp. 3560–3567.
  * [14] J. Hoff, A. Ramezani, S. J. Chung, and S. Hutchinson, “Optimizing the structure and movement of a robotic bat with biological kinematic synergies,” _Int. Journal of Robotics Research_ , vol. 37, no. 10, pp. 1233–1252, 2018.
  * [15] A. Ramezani, X. Shi, S. J. Chung, and S. Hutchinson, “Lagrangian modeling and flight control of articulated-winged bat robot,” in _Int. Conf. on Intelligent Robots and Systems (IROS)_. IEEE, 2015, pp. 2867–2874.
  * [16] A. Ramezani, X. Shi, S. J. Chung, and S. A. Hutchinson, “Modeling and nonlinear flight controller synthesis of a bat-inspired micro aerial vehicle,” in _AIAA Guidance, Navigation, and Control Conf._ , 2016.
  * [17] A. Ramezani, X. Shi, S. J. Chung, and S. Hutchinson, “Nonlinear flight controller synthesis of a bat-inspired micro aerial vehicle,” in _AIAA Guidance, Navigation, and Control Conf._ , 2016.
  * [18] J. Hoff, U. Syed, A. Ramezani, and S. Hutchinson, “Trajectory planning for a bat-like flapping wing robot,” in _Int. Conf. on Intelligent Robots and Systems (IROS)_. IEEE, 2019, pp. 6800–6805.
  * [19] A. Ramezani, “Towards biomimicry of a bat-style perching maneuver on structures: the manipulation of inertial dynamics,” in _Int. Conf. on Robotics and Automation (ICRA)_. IEEE, 2020, in press.
  * [20] H. Hauser, A. J. Ijspeert, R. M. Füchslin, R. Pfeifer, and W. Maass, “Towards a theoretical foundation for morphological computation with compliant bodies,” _Biological cybernetics_ , vol. 105, no. 5-6, pp. 355–370, 2011.
  * [21] T. Lee, M. Leok, and N. H. McClamroch, “Global formulations of lagrangian and hamiltonian dynamics on manifolds,” _Springer_ , vol. 13, p. 31, 2017.
  * [22] S. P. Sane and M. H. Dickinson, “The control of flight force by a flapping wing: lift and drag production,” _Journal of experimental biology_ , vol. 204, no. 15, pp. 2607–2626, 2001.
  * [23] A. D. Ames, K. Galloway, K. Sreenath, and J. W. Grizzle, “Rapidly exponentially stabilizing control lyapunov functions and hybrid zero dynamics,” _IEEE Transactions on Automatic Control_ , vol. 59, no. 4, pp. 876–891, 2014.


# Search for the rare decays $J/\psi\to D^{0}e^{+}e^{-}+c.c.$ and
$\psi(3686)\to D^{0}e^{+}e^{-}+c.c.$

M. Ablikim1, M. N. Achasov9,d, S. Ahmed14, M. Albrecht4, A. Amoroso50A,50C, F.
F. An1, Q. An47,39, J. Z. Bai1, O. Bakina24, R. Baldini Ferroli20A, Y. Ban32,
D. W. Bennett19, J. V. Bennett5, N. Berger23, M. Bertani20A, D. Bettoni21A, J.
M. Bian45, F. Bianchi50A,50C, E. Boger24,b, I. Boyko24, R. A. Briere5, H.
Cai52, X. Cai1,39, O. Cakir42A, A. Calcaterra20A, G. F. Cao1,43, S. A.
Cetin42B, J. Chai50C, J. F. Chang1,39, G. Chelkov24,b,c, G. Chen1, H. S.
Chen1,43, J. C. Chen1, M. L. Chen1,39, S. J. Chen30, X. R. Chen27, Y. B.
Chen1,39, X. K. Chu32, G. Cibinetto21A, H. L. Dai1,39, J. P. Dai35,h, A.
Dbeyssi14, D. Dedovich24, Z. Y. Deng1, A. Denig23, I. Denysenko24, M.
Destefanis50A,50C, F. De Mori50A,50C, Y. Ding28, C. Dong31, J. Dong1,39, L. Y.
Dong1,43, M. Y. Dong1,39,43, O. Dorjkhaidav22, Z. L. Dou30, S. X. Du54, P. F.
Duan1, J. Fang1,39, S. S. Fang1,43, X. Fang47,39, Y. Fang1, R.
Farinelli21A,21B, L. Fava50B,50C, S. Fegan23, F. Feldbauer23, G. Felici20A, C.
Q. Feng47,39, E. Fioravanti21A, M. Fritsch23,14, C. D. Fu1, Q. Gao1, X. L.
Gao47,39, Y. Gao41, Y. G. Gao6, Z. Gao47,39, I. Garzia21A, K. Goetzen10, L.
Gong31, W. X. Gong1,39, W. Gradl23, M. Greco50A,50C, M. H. Gu1,39, S. Gu15, Y.
T. Gu12, A. Q. Guo1, L. B. Guo29, R. P. Guo1, Y. P. Guo23, Z. Haddadi26, A.
Hafner23, S. Han52, X. Q. Hao15, F. A. Harris44, K. L. He1,43, X. Q. He46, F.
H. Heinsius4, T. Held4, Y. K. Heng1,39,43, T. Holtmann4, Z. L. Hou1, C. Hu29,
H. M. Hu1,43, T. Hu1,39,43, Y. Hu1, G. S. Huang47,39, J. S. Huang15, X. T.
Huang34, X. Z. Huang30, Z. L. Huang28, T. Hussain49, W. Ikegami Andersson51,
Q. Ji1, Q. P. Ji15, X. B. Ji1,43, X. L. Ji1,39, X. S. Jiang1,39,43, X. Y.
Jiang31, J. B. Jiao34, Z. Jiao17, D. P. Jin1,39,43, S. Jin1,43, T.
Johansson51, A. Julin45, N. Kalantar-Nayestanaki26, X. L. Kang1, X. S. Kang31,
M. Kavatsyuk26, B. C. Ke5, T. Khan47,39, P. Kiese23, R. Kliemt10, B. Kloss23,
L. Koch25, O. B. Kolcu42B,f, B. Kopf4, M. Kornicer44, M. Kuemmel4, M.
Kuhlmann4, A. Kupsc51, W. Kühn25, J. S. Lange25, M. Lara19, P. Larin14, L.
Lavezzi50C, H. Leithoff23, C. Leng50C, C. Li51, Cheng Li47,39, D. M. Li54, F.
Li1,39, F. Y. Li32, G. Li1, H. B. Li1,43, H. J. Li1, J. C. Li1, Jin Li33, Kang
Li13, Ke Li34, Lei Li3, P. L. Li47,39, P. R. Li43,7, Q. Y. Li34, T. Li34, W.
D. Li1,43, W. G. Li1, X. L. Li34, X. N. Li1,39, X. Q. Li31, Z. B. Li40, H.
Liang47,39, Y. F. Liang37, Y. T. Liang25, G. R. Liao11, D. X. Lin14, B.
Liu35,h, B. J. Liu1, C. X. Liu1, D. Liu47,39, F. H. Liu36, Fang Liu1, Feng
Liu6, H. B. Liu12, H. M. Liu1,43, Huanhuan Liu1, Huihui Liu16, J. B. Liu47,39,
J. P. Liu52, J. Y. Liu1, K. Liu41, K. Y. Liu28, Ke Liu6, L. D. Liu32, P. L.
Liu1,39, Q. Liu43, S. B. Liu47,39, X. Liu27, Y. B. Liu31, Y. Y. Liu31, Z. A.
Liu1,39,43, Zhiqing Liu23, Y. F. Long32, X. C. Lou1,39,43, H. J. Lu17, J. G.
Lu1,39, Y. Lu1, Y. P. Lu1,39, C. L. Luo29, M. X. Luo53, T. Luo44, X. L.
Luo1,39, X. R. Lyu43, F. C. Ma28, H. L. Ma1, L. L. Ma34, M. M. Ma1, Q. M. Ma1,
T. Ma1, X. N. Ma31, X. Y. Ma1,39, Y. M. Ma34, F. E. Maas14, M.
Maggiora50A,50C, Q. A. Malik49, Y. J. Mao32, Z. P. Mao1, S. Marcello50A,50C,
J. G. Messchendorp26, G. Mezzadri21B, J. Min1,39, T. J. Min1, R. E.
Mitchell19, X. H. Mo1,39,43, Y. J. Mo6, C. Morales Morales14, G. Morello20A,
N. Yu. Muchnoi9,d, H. Muramatsu45, P. Musiol4, A. Mustafa4, Y. Nefedov24, F.
Nerling10, I. B. Nikolaev9,d, Z. Ning1,39, S. Nisar8, S. L. Niu1,39, X. Y.
Niu1, S. L. Olsen33, Q. Ouyang1,39,43, S. Pacetti20B, Y. Pan47,39, M.
Papenbrock51, P. Patteri20A, M. Pelizaeus4, J. Pellegrino50A,50C, H. P.
Peng47,39, K. Peters10,g, J. Pettersson51, J. L. Ping29, R. G. Ping1,43, R.
Poling45, V. Prasad47,39, H. R. Qi2, M. Qi30, S. Qian1,39, C. F. Qiao43, J. J.
Qin43, N. Qin52, X. S. Qin1, Z. H. Qin1,39, J. F. Qiu1, K. H. Rashid49,i, C.
F. Redmer23, M. Richter4, M. Ripka23, G. Rong1,43, Ch. Rosner14, X. D. Ruan12,
A. Sarantsev24,e, M. Savrié21B, C. Schnier4, K. Schoenning51, W. Shan32, M.
Shao47,39, C. P. Shen2, P. X. Shen31, X. Y. Shen1,43, H. Y. Sheng1, J. J.
Song34, W. M. Song34, X. Y. Song1, S. Sosio50A,50C, C. Sowa4, S.
Spataro50A,50C, G. X. Sun1, J. F. Sun15, S. S. Sun1,43, X. H. Sun1, Y. J.
Sun47,39, Y. K Sun47,39, Y. Z. Sun1, Z. J. Sun1,39, Z. T. Sun19, C. J. Tang37,
G. Y. Tang1, X. Tang1, I. Tapan42C, M. Tiemens26, B. T. Tsednee22, I. Uman42D,
G. S. Varner44, B. Wang1, B. L. Wang43, D. Wang32, D. Y. Wang32, Dan Wang43,
K. Wang1,39, L. L. Wang1, L. S. Wang1, M. Wang34, P. Wang1, P. L. Wang1, W. P.
Wang47,39, X. F. Wang41, Y. Wang38, Y. D. Wang14, Y. F. Wang1,39,43, Y. Q.
Wang23, Z. Wang1,39, Z. G. Wang1,39, Z. H. Wang47,39, Z. Y. Wang1, Zongyuan
Wang1, T. Weber23, D. H. Wei11, J. H. Wei31, P. Weidenkaff23, S. P. Wen1, U.
Wiedner4, M. Wolke51, L. H. Wu1, L. J. Wu1, Z. Wu1,39, L. Xia47,39, Y. Xia18,
D. Xiao1, H. Xiao48, Y. J. Xiao1, Z. J. Xiao29, Y. G. Xie1,39, Y. H. Xie6, X.
A. Xiong1, Q. L. Xiu1,39, G. F. Xu1, J. J. Xu1, L. Xu1, Q. J. Xu13, Q. N.
Xu43, X. P. Xu38, L. Yan50A,50C, W. B. Yan47,39, W. C. Yan47,39, Y. H. Yan18,
H. J. Yang35,h, H. X. Yang1, L. Yang52, Y. H. Yang30, Y. X. Yang11, M. Ye1,39,
M. H. Ye7, J. H. Yin1, Z. Y. You40, B. X. Yu1,39,43, C. X. Yu31, J. S. Yu27,
C. Z. Yuan1,43, Y. Yuan1, A. Yuncu42B,a, A. A. Zafar49, Y. Zeng18, Z.
Zeng47,39, B. X. Zhang1, B. Y. Zhang1,39, C. C. Zhang1, D. H. Zhang1, H. H.
Zhang40, H. Y. Zhang1,39, J. Zhang1, J. L. Zhang1, J. Q. Zhang1, J. W.
Zhang1,39,43, J. Y. Zhang1, J. Z. Zhang1,43, K. Zhang1, L. Zhang41, S. Q.
Zhang31, X. Y. Zhang34, Y. H. Zhang1,39, Y. T. Zhang47,39, Yang Zhang1, Yao
Zhang1, Yu Zhang43, Z. H. Zhang6, Z. P. Zhang47, Z. Y. Zhang52, G. Zhao1, J.
W. Zhao1,39, J. Y. Zhao1, J. Z. Zhao1,39, Lei Zhao47,39, Ling Zhao1, M. G.
Zhao31, Q. Zhao1, S. J. Zhao54, T. C. Zhao1, Y. B. Zhao1,39, Z. G. Zhao47,39,
A. Zhemchugov24,b, B. Zheng48,14, J. P. Zheng1,39, W. J. Zheng34, Y. H.
Zheng43, B. Zhong29, L. Zhou1,39, X. Zhou52, X. K. Zhou47,39, X. R. Zhou47,39,
X. Y. Zhou1, Y. X. Zhou12, K. Zhu1, K. J. Zhu1,39,43, S. Zhu1, S. H. Zhu46, X.
L. Zhu41, Y. C. Zhu47,39, Y. S. Zhu1,43, Z. A. Zhu1,43, J. Zhuang1,39, L.
Zotti50A,50C, B. S. Zou1, J. H. Zou1
(BESIII Collaboration)
1 Institute of High Energy Physics, Beijing 100049, People’s Republic of China
2 Beihang University, Beijing 100191, People’s Republic of China
3 Beijing Institute of Petrochemical Technology, Beijing 102617, People’s
Republic of China
4 Bochum Ruhr-University, D-44780 Bochum, Germany
5 Carnegie Mellon University, Pittsburgh, Pennsylvania 15213, USA
6 Central China Normal University, Wuhan 430079, People’s Republic of China
7 China Center of Advanced Science and Technology, Beijing 100190, People’s
Republic of China
8 COMSATS Institute of Information Technology, Lahore, Defence Road, Off
Raiwind Road, 54000 Lahore, Pakistan
9 G.I. Budker Institute of Nuclear Physics SB RAS (BINP), Novosibirsk 630090,
Russia
10 GSI Helmholtzcentre for Heavy Ion Research GmbH, D-64291 Darmstadt, Germany
11 Guangxi Normal University, Guilin 541004, People’s Republic of China
12 Guangxi University, Nanning 530004, People’s Republic of China
13 Hangzhou Normal University, Hangzhou 310036, People’s Republic of China
14 Helmholtz Institute Mainz, Johann-Joachim-Becher-Weg 45, D-55099 Mainz,
Germany
15 Henan Normal University, Xinxiang 453007, People’s Republic of China
16 Henan University of Science and Technology, Luoyang 471003, People’s
Republic of China
17 Huangshan College, Huangshan 245000, People’s Republic of China
18 Hunan University, Changsha 410082, People’s Republic of China
19 Indiana University, Bloomington, Indiana 47405, USA
20 (A)INFN Laboratori Nazionali di Frascati, I-00044, Frascati, Italy; (B)INFN
and University of Perugia, I-06100, Perugia, Italy
21 (A)INFN Sezione di Ferrara, I-44122, Ferrara, Italy; (B)University of
Ferrara, I-44122, Ferrara, Italy
22 Institute of Physics and Technology, Peace Ave. 54B, Ulaanbaatar 13330,
Mongolia
23 Johannes Gutenberg University of Mainz, Johann-Joachim-Becher-Weg 45,
D-55099 Mainz, Germany
24 Joint Institute for Nuclear Research, 141980 Dubna, Moscow region, Russia
25 Justus-Liebig-Universitaet Giessen, II. Physikalisches Institut, Heinrich-
Buff-Ring 16, D-35392 Giessen, Germany
26 KVI-CART, University of Groningen, NL-9747 AA Groningen, The Netherlands
27 Lanzhou University, Lanzhou 730000, People’s Republic of China
28 Liaoning University, Shenyang 110036, People’s Republic of China
29 Nanjing Normal University, Nanjing 210023, People’s Republic of China
30 Nanjing University, Nanjing 210093, People’s Republic of China
31 Nankai University, Tianjin 300071, People’s Republic of China
32 Peking University, Beijing 100871, People’s Republic of China
33 Seoul National University, Seoul, 151-747 Korea
34 Shandong University, Jinan 250100, People’s Republic of China
35 Shanghai Jiao Tong University, Shanghai 200240, People’s Republic of China
36 Shanxi University, Taiyuan 030006, People’s Republic of China
37 Sichuan University, Chengdu 610064, People’s Republic of China
38 Soochow University, Suzhou 215006, People’s Republic of China
39 State Key Laboratory of Particle Detection and Electronics, Beijing 100049,
Hefei 230026, People’s Republic of China
40 Sun Yat-Sen University, Guangzhou 510275, People’s Republic of China
41 Tsinghua University, Beijing 100084, People’s Republic of China
42 (A)Ankara University, 06100 Tandogan, Ankara, Turkey; (B)Istanbul Bilgi
University, 34060 Eyup, Istanbul, Turkey; (C)Uludag University, 16059 Bursa,
Turkey; (D)Near East University, Nicosia, North Cyprus, Mersin 10, Turkey
43 University of Chinese Academy of Sciences, Beijing 100049, People’s
Republic of China
44 University of Hawaii, Honolulu, Hawaii 96822, USA
45 University of Minnesota, Minneapolis, Minnesota 55455, USA
46 University of Science and Technology Liaoning, Anshan 114051, People’s
Republic of China
47 University of Science and Technology of China, Hefei 230026, People’s
Republic of China
48 University of South China, Hengyang 421001, People’s Republic of China
49 University of the Punjab, Lahore-54590, Pakistan
50 (A)University of Turin, I-10125, Turin, Italy; (B)University of Eastern
Piedmont, I-15121, Alessandria, Italy; (C)INFN, I-10125, Turin, Italy
51 Uppsala University, Box 516, SE-75120 Uppsala, Sweden
52 Wuhan University, Wuhan 430072, People’s Republic of China
53 Zhejiang University, Hangzhou 310027, People’s Republic of China
54 Zhengzhou University, Zhengzhou 450001, People’s Republic of China
a Also at Bogazici University, 34342 Istanbul, Turkey
b Also at the Moscow Institute of Physics and Technology, Moscow 141700,
Russia
c Also at the Functional Electronics Laboratory, Tomsk State University,
Tomsk, 634050, Russia
d Also at the Novosibirsk State University, Novosibirsk, 630090, Russia
e Also at the NRC ”Kurchatov Institute”, PNPI, 188300, Gatchina, Russia
f Also at Istanbul Arel University, 34295 Istanbul, Turkey
g Also at Goethe University Frankfurt, 60323 Frankfurt am Main, Germany
h Also at Key Laboratory for Particle Physics, Astrophysics and Cosmology,
Ministry of Education; Shanghai Key Laboratory for Particle Physics and
Cosmology; Institute of Nuclear and Particle Physics, Shanghai 200240,
People’s Republic of China
i Government College Women University, Sialkot - 51310. Punjab, Pakistan.

###### Abstract

Using the data samples of $(1310.6\pm 7.2)\times 10^{6}$ $J/\psi$ events and
$(448.1\pm 2.9)\times 10^{6}$ $\psi(3686)$ events collected with the BESIII
detector, we search for the rare decays $J/\psi\to D^{0}e^{+}e^{-}+c.c.$ and
$\psi(3686)\to D^{0}e^{+}e^{-}+c.c.$. No significant signals are observed and
the corresponding upper limits on the branching fractions at the $90\%$
confidence level are determined to be $\mathcal{B}(J/\psi\to
D^{0}e^{+}e^{-}+c.c.)<8.5\times 10^{-8}$ and $\mathcal{B}(\psi(3686)\to
D^{0}e^{+}e^{-}+c.c.)<1.4\times 10^{-7}$, respectively. Our limit on
$\mathcal{B}(J/\psi\to D^{0}e^{+}e^{-}+c.c.)$ is more stringent by two orders
of magnitude than the previous results, and the $\mathcal{B}(\psi(3686)\to
D^{0}e^{+}e^{-}+c.c.)$ is measured for the first time.

###### pacs:

13.20.Gd, 12.38.Qk

In the Standard Model (SM), decays of the charmonium resonances $J/\psi$ and
$\psi(3686)$ (collectively referred to as $\psi$ throughout the text) induced
by Flavor Changing Neutral Current (FCNC) are forbidden at the tree level due
to the Glashow-Iliopoulos-Maiani (GIM) mechanism GIM , but can occur via a
$c\to u$ transition at the loop level, _e.g._ , shown in Fig. 1 for the decay
of $\psi~{}\to D^{0}e^{+}e^{-}$. Such decays can also occur via long-distance
effects on the hadron level, which are, according to Ref. FCNC , expected to
have the same order of magnitude as the FCNC process.

\begin{overpic}[width=138.76157pt,height=65.04034pt]{feynman3.eps}
\end{overpic} Figure 1: Typical Feynman diagram for $\psi\to D^{0}e^{+}e^{-}$.

The decay branching fraction for this kind of rare process is expected to be
of order $10^{-10}$ to $10^{-13}$ in the SM SM1 ; SM2 . However, many new
physics models, such as the Topcolor models TopC , the minimal supersymmetric
standard model with R-parity violation Supersy , and the two Higgs doublet
model Thiggs , predict that the decay branching fractions can be enhanced by
two or three orders of magnitude. Searching for experimental evidence for
these FCNC processes offers an ideal opportunity to study non-perturbative QCD
effects and their underlying dynamics, and serves as a probe to search for new
physics beyond the SM newphy1 ; newphy2 .

The semileptonic decay $J/\psi$ $\to D^{0}e^{+}e^{-}$ (the charged conjugate
channel is always implied throughout the text if not mentioned explicitly) is
an interesting decay to study FCNC-induced processes. The BESII experiment has
searched for the decay $J/\psi\to D^{0}e^{+}e^{-}$ by using $58\times 10^{6}$
$J/\psi$ events, and an upper bound on the branching fraction was placed at
the order of $10^{-5}$ besII , which is far away from the theoretical
prediction. A more precise measurement with larger statistics is desirable to
test the theoretical predictions more stringently.

The Beijing Spectrometer (BESIII) detector detector1 , located at the Beijing
Electron-Positron Collider (BEPCII), has collected ($1310.6\pm 7.2$) $\times$
$10^{6}$ $J/\psi$ events jpsinum1 ; jpsinum2 and ($448.1\pm 2.9$) $\times$
$10^{6}$ $\psi(3686)$ events psipnum1 ; psipnum2 , which are the world largest
samples collected with electron-positron collision at threshold. The high
quality and large statistics data samples provide a unique opportunity to
search for physics beyond the SM. In this paper, we present a search for the
rare processes $\psi\to D^{0}e^{+}e^{-}$.

BEPCII is a double-ring $e^{+}e^{-}$ collider running at center-of-mass (c.m.)
energies between 2.0 and 4.6 GeV, reaching a peak luminosity of $1.0\times
10^{33}$ cm-2s-1 at a c.m. energy of 3770 MeV. The cylindrical BESIII detector
has an effective geometrical acceptance of $93\%$ of 4$\pi$ and is divided
into a barrel section and two endcaps. It contains a small-cell, helium-based
(40$\%$ He, 60$\%$ C3H8) main drift chamber (MDC) which provides momentum
measurements for charged particles with a resolution of $0.5\%$ at a momentum
of 1 GeV/$c$ in a magnetic field of 1 Tesla. The energy loss measurement
($dE/dx$) provided by the MDC has a resolution better than $6\%$. A time-of-
flight system (TOF) consisting of 5-cm-thick plastic scintillators can measure
the flight time of charged particles with a time resolution of 80 ps in the
barrel and 110 ps in the endcaps. An electromagnetic calorimeter (EMC)
consisting of 6240 CsI(Tl) crystals in a cylindrical structure and two endcaps
is used to measure the energies of photons and electrons. The energy
resolution of the EMC is $2.5\%$ in the barrel and $5.0\%$ in the endcaps for
photons and electrons with an energy of 1 GeV. The position resolution of the
EMC is 6 mm in the barrel and 9 mm in the endcaps. A muon counter (MUC) system
consisting of more than 1200 m2 of Resistive Plate Chambers (RPC) is used to
identify muons and provides a spatial resolution better than 2 cm. A detailed
description of the BESIII detector can be found in Ref. detector1 .

Large samples of Monte Carlo (MC) simulated events are used to optimize the
event selection criteria, estimate the background contaminations and determine
the selection efficiencies. The MC samples are generated using a geant4-based
geant4 simulation software package BESIII object oriented simulation tool
(BOOST) Deng , which includes the description of geometry and material, the
detector response and the digitization model, as well as tracking for the
detector running conditions and performances. ’Inclusive’ MC samples of
$1225\times 10^{6}$ $J/\psi$ and $506\times 10^{6}$ $\psi(3686)$ generic decay
events are generated to study the backgrounds. In the simulation, production
of the $\psi$ resonances is simulated with the kkmc kkmc1 generator, while the
decays are generated by evtgen evtgen for the known decay modes, setting the
branching fractions according to the Particle Data Group (PDG) PDG , or
lundcharm lundcharm for the remaining unknown decays. The QED Final State
Radiation (FSR) effect is simulated with photos photos . Signal samples
$\psi\to D^{0}e^{+}e^{-}$ of $2.0\times 10^{5}$ events are generated according
to a theoretical calculation based on the vector meson dominance (VMD) model
DIYg , where $e^{+}e^{-}$ is generated from a virtual photon decay.

To study the decay $\psi\to D^{0}e^{+}e^{-}$, we reconstruct the $D^{0}$
signal through its three prominent exclusive hadronic decay modes,
$K^{-}\pi^{+}$ (mode I), $K^{-}\pi^{+}\pi^{0}$ (mode II), and
$K^{-}\pi^{+}\pi^{+}\pi^{-}$ (mode III), which have relatively large branching
fractions, and suffer from relatively low background.

Charged tracks are reconstructed from hits in the MDC. A good charged track is
required to have a polar angle $\theta$ that satisfies $|\cos\theta|<0.93$ and
a distance of closest approach point to the interaction point (IP) within 10
cm along the beam direction and 1 cm in the plane perpendicular to the beam.
The measured ionizing energy loss $dE/dx$ in the MDC and flight time in the
TOF are combined to form a particle identification (PID) confidence level
(C.L.) for each particle hypothesis $i$ ($i$ = $e$, $\pi$, $K$, $p$). Each
track is assigned to the particle type with the highest C.L.

Photons are reconstructed from clusters of energy depositions in the EMC
crystals. The energy deposited in the nearby TOF counters is included to
improve the reconstruction efficiency and energy resolution. Good photons are
required to have energy greater than 25 MeV in the barrel region
($|\cos\theta|<0.80$) or 50 MeV in the endcap region
($0.86<|\cos\theta|<0.92$). The showers in the transition region between the
barrel and endcap are poorly reconstructed and are excluded from the analysis.
To exclude showers from charged particles, a photon must be separated from any
charged track by more than $10^{\circ}$. A requirement on the EMC time $t$
with respect to the event start time of $0\leq\rm{t}\leq 700$ ns is used to
suppress electronic noise and energy deposits in the EMC unrelated to the
events.

The candidate events are selected by requiring an electron-positron pair, a
kaon, one or three pions depending on the $D^{0}$ decay mode as well as two
photons for decay mode II. To suppress backgrounds, the electron or positron
are required to satisfy $E/p>0.8$, where $E$ and $p$ are the corresponding
energy deposited in the EMC and momentum measured in the MDC, respectively. A
vertex fit is performed on the selected charged tracks to ensure that the
events originate at the IP. To improve resolution and reduce backgrounds, a
four-constraint (4C) kinematic fit imposing energy-momentum conservation is
carried out under the hypothesis of $\psi\to K^{-}\pi^{+}e^{+}e^{-}$ (mode I)
or $\psi\to K^{-}\pi^{+}\pi^{+}\pi^{-}e^{+}e^{-}$ (mode III), and
$\chi^{2}_{4C}<60$ is required. For decay mode II, a 5C kinematic fit is
performed under the hypothesis of $\psi\to K^{-}\pi^{+}\gamma\gamma
e^{+}e^{-}$ with an additional constraint on the mass of the $\gamma\gamma$
pair to the $\pi^{0}$ nominal mass (5C), and $\chi^{2}_{5C}<70$ is required.
For events with more than two photon candidates, the $\gamma\gamma$
combination with the least $\chi^{2}_{5C}$ is retained for further analysis.

With the above selection criteria, the dominant backgrounds are from processes
with similar hadronic final states as the signal but with an $e^{+}e^{-}$ pair
which comes from a $\gamma$ conversion from interactions with the detector
material. This effect occurs primarily at the beam pipe, which has an inner
diameter of 63 mm, or the inner MDC wall, which has a diameter of 108 mm. To
suppress these backgrounds, a $\gamma$ conversion finder algorithm xuzr was
developed to reconstruct the $\gamma$ conversion vertex. A set of parameters
is defined to study the properties of the $\gamma$ conversion vertex. These
include the distance $R_{xy}$ from the IP to the reconstructed vertex point of
the $e^{+}e^{-}$ pair in the $\it{x}$-$\it{y}$ plane, the distance
$\Delta_{\it{xy}}$ between the two intersection points of the two circles
describing the trajectories of the electron and positron in the
$\it{x}$-$\it{y}$ plane and the line connecting their centers, the invariant
mass $M(e^{+}e^{-})$ of the $e^{+}e^{-}$ pair, the angle $\theta_{eg}$ between
the photon momentum vector (in the $e^{+}e^{-}$ system) and the direction from
the IP (run averaged) to the reconstructed vertex in the $\it{x}$-$\it{y}$
plane, and the angle $\theta_{e^{+}e^{-}}$ between the momenta of $e^{+}$ and
$e^{-}$ in the $\it{x}$-$\it{y}$ plane. The $\gamma$ conversion vertex is
identified with the criteria $R_{\it{xy}}>2$ cm and
$\Delta_{\it{xy}}<6\sigma_{\Delta_{\it{xy}}}$ as well as
$\cos\theta_{eg}>0.99$, and the events satisfying these criteria are removed.
This selection removes 95.0% of all $\gamma$ conversion events while losing
less than 5.0% efficiency for the signal decays. The criteria are determined
by studying a control sample of $J/\psi\to\gamma\pi^{+}\pi^{-}$ with the
$\gamma$ conversion into a $e^{+}e^{-}$ pair. The resolution
$\sigma_{\Delta_{\it{xy}}}$ of $\Delta_{\it{xy}}$ depends on the angle
$\theta_{e^{+}e^{-}}$ and is also determined from the study of the control
sample.

After the $\gamma$ conversion suppression criteria are applied, the inclusive
$J/\psi$ and $\psi(3686)$ MC samples are used to study the remaining
background contamination. In decay mode I, the dominant backgrounds are
$\psi\to(\gamma)\pi^{+}\pi^{-}e^{+}e^{-}$ $\it{etc}$., due to $K$/$\pi$ mis-
identification for high momentum tracks. The 4C kinematic fits with hypotheses
$\psi\to(\gamma)\pi^{+}\pi^{-}e^{+}e^{-}$ are performed, and the corresponding
$\chi^{2}_{4C}(K^{-}\pi^{+}e^{+}e^{-})<\chi^{2}_{4C}((\gamma)\pi^{+}\pi^{-}e^{+}e^{-})$
is required. Another potential background in the $\psi(3686)$ data sample is
$\psi(3686)\to\pi^{+}\pi^{-}J/\psi$ with subsequent decay $J/\psi\to
e^{+}e^{-}$, where a low momentum $\pi^{\pm}$ is mis-identified to be an
electron while a high momentum electron (positron) is mis-identified as a
kaon. The requirement $|M^{\rm
Recoil}(\pi^{+}\pi^{-})-m(J/\psi)|<0.02~{}\mbox{GeV/$c^{2}$}$ is implemented
to reduce this background, where $M^{\rm Recoil}(\pi^{+}\pi^{-})$ is the
recoil mass of the two low momentum tracks with opposite charges using the
$\pi^{\pm}$ hypothesis, and $m(J/\psi)$ is the nominal $J/\psi$ mass. In decay
mode II, the dominant background is $\psi\to K_{S}K^{-}\pi^{+}$ with
subsequent decays $K_{S}\to\pi^{0}\pi^{0}$ and with a $\pi^{0}$ Dalitz decay.
The requirement $|M^{\rm
Recoil}(K^{-}\pi^{+})-m(K_{S})|>0.06~{}\mbox{GeV/$c^{2}$}$ is applied to
suppress this background, where $M^{\rm Recoil}(K^{-}\pi^{+})$ is the recoil
mass of the $K^{-}\pi^{+}$ system, and $m(K_{S})$ is the $K_{S}$ nominal mass.
Another background is $\psi\to\omega\pi^{+}\pi^{-}$
($\omega\to\pi^{0}e^{+}e^{-}$) with a pion mis-identified as a kaon. This
background is suppressed by a requirement on the kinematic fit quality with
$\chi^{2}_{5C}(K^{-}\pi^{+}\pi^{0}e^{+}e^{-})<\chi^{2}_{5C}(\pi^{+}\pi^{-}\pi^{0}e^{+}e^{-})$.
In decay mode III, the dominant background is $\psi\to e^{+}e^{-}K_{S}K\pi$
with the subsequent $K_{S}$ decay into $\pi^{+}\pi^{-}$. This background is
rejected if any of the $\pi^{+}\pi^{-}$ invariant mass $M(\pi^{+}\pi^{-})$
satisfies $|M(\pi^{+}\pi^{-})-m(K_{S})|<0.025~{}\mbox{GeV/$c^{2}$}$.

\begin{overpic}[width=99.73074pt,height=86.72267pt]{combinefitkpi_jpsi4.eps}
\end{overpic}\begin{overpic}[width=99.73074pt,height=86.72267pt]{combinefitkpi_psip5.eps}
\end{overpic}\begin{overpic}[width=99.73074pt,height=86.72267pt]{combinefitkpipi0_jpsi3.eps}
\end{overpic}\begin{overpic}[width=99.73074pt,height=86.72267pt]{combinefitkpipi0_psip4.eps}
\end{overpic}\begin{overpic}[width=99.73074pt,height=86.72267pt]{combinefitk3pi_jpsi3.eps}
\end{overpic}\begin{overpic}[width=99.73074pt,height=86.72267pt]{combinefitk3pi_psip4.eps}
\end{overpic}

Figure 2:  Distributions of $K^{-}\pi^{+}$ (upper row), $K^{-}\pi^{+}\pi^{0}$
(middle row) and $K^{-}\pi^{+}\pi^{+}\pi^{-}$ (bottom row) invariant masses.
The left and right columns are for the $J/\psi$ and $\psi(3686)$ samples,
respectively. Dots with error bars are data, the solid and dashed curves are
for the signal shape and the total best fit to data, respectively.

After applying the above selection criteria, the distributions of the
$K^{-}\pi^{+}$, $K^{-}\pi^{+}\pi^{0}$, and $K^{-}\pi^{+}\pi^{+}\pi^{-}$
invariant masses for the surviving events in the three $D^{0}$ meson decay
modes are shown in Fig. 2. No $D^{0}$ signals are observed, and therefore
upper limits on the branching fractions at the 90% C.L. are determined.

In the measurements of the branching fractions, the sources of systematic
uncertainty include the detection efficiencies of charged tracks and photons,
the PID efficiency, the kinematic fit, $\gamma$ conversion veto, mass window
requirements, the fit procedure, the decay branching fractions of intermediate
states, as well as the total numbers of $\psi$ events. The individual
systematic uncertainties are estimated and described in detail as follows:

(a) Tracking efficiency: The tracking efficiencies for $\pi^{\pm}$ and
$K^{\pm}$ are studied using control samples of
$J/\psi\to\rho\pi\to\pi^{+}\pi^{-}\pi^{0}$, $J/\psi\to
p\bar{p}\pi^{+}\pi^{-}$, and $J/\psi\to K^{0}_{S}K^{-}\pi^{+}$ trackerror1 ;
trackerror2 . The tracking efficiency for electrons (positrons) is studied
with a control sample of radiative Bhabha events. The differences in tracking
efficiencies between data and MC simulation are 1% per track for $K$, $\pi$
and $e$, respectively, and are taken as the systematic uncertainties.

(b) PID: The PID efficiencies of $\pi^{\pm}$ and $K^{\pm}$ are studied with
the same control samples as in the study of the tracking efficiency
trackerror1 ; trackerror2 . The PID efficiency from the data sample agrees
with that of the MC simulation within 1% for each track. The uncertainty of
the PID efficiency for electrons (positrons) is studied with the control
sample of radiative Bhabha events, and 1.0% is assigned for each electron
(positron). The uncertainty of the $E/p$ requirement for electrons (positrons)
is studied with the control sample $J/\psi\to\pi^{+}\pi^{-}\pi^{0}$
($\pi^{0}\to\gamma e^{+}e^{-}$), and an uncertainty of 2% is assigned.

(c) Photon detection efficiency: The photon detection efficiency is studied
with the control samples $J/\psi\to\pi^{+}\pi^{-}\pi^{0}$, and a weighted
average uncertainty, according to the energy distribution, is determined to be
0.6% per photon.

(d) Kinematic fit: The uncertainty associated with the kinematic fit arises
from the inconsistency of the track helix parameters between data and MC
simulation. Therefore, the three track parameters $\phi_{0}$, $\kappa$ and
$\tan\lambda$ are corrected for the signal MC samples, where the correction
factors are obtained by comparing the pull distributions of the control
samples described in detail in Ref. kinematic . The resulting difference in
the detection efficiencies between the samples with and without the helix
correction is taken as the systematic uncertainty.

(e) $\gamma$ conversion veto: The effect of the $\gamma$ conversion veto is
studied using a control sample of $J/\psi\to\pi^{+}\pi^{-}\pi^{0}$ with the
subsequent Dalitz decay $\pi^{0}\to\gamma e^{+}e^{-}$. A clean control sample
is selected, and the corresponding MC sample is generated with the RhoPi
generator based on a formalism of helicity coupling amplitudes for the process
$J/\psi\to\pi^{+}\pi^{-}\pi^{0}$ rhopi , while a generator for the decay
$\pi^{0}\to\gamma e^{+}e^{-}$ adopts a simple pole approximation in the form
factor $|F(q^{2})|=1+\alpha q^{2}/m^{2}_{\pi^{0}}$ with $\alpha=0.032$ PDG .
The efficiency of the $\gamma$ conversion veto is the ratio of signal yields
with and without the $\gamma$ conversion veto, where the signal yields are
extracted by fitting the $e^{+}e^{-}$ invariant mass. The resulting difference
between data and MC, 1.7%, is taken as the systematic uncertainty.

(f) Mass window requirements: Various requirements of mass window by widening
5 MeV/$c^{2}$ are applied to veto the different backgrounds, the corresponding
uncertainties are studied by changing the appropriate values. The resulting
changes in the final results are taken as the systematic uncertainties.

(g) Branching fractions of intermediate states:  The uncertainties of the
decay branching fractions of intermediate states in the cascade decays are
quoted from the PDG PDG .

(h) Total number of $\psi$ events: The uncertainties on the total numbers of
$J/\psi$ and $\psi(3686)$ events are 0.55% and 0.62%, respectively, which are
determined by studying the inclusive hadron events jpsinum1 ; jpsinum2 ;
psipnum1 ; psipnum2 .

All the individual systematic uncertainties are summarized in Table 1, where
the sourcees of the uncertainties tagged with ${}^{\prime}{\ast}^{\prime}$ are
assumed to be 100% correlated among the three different $D^{0}$ decay modes.
The efficiencies for other selection criteria, the trigger simulation, the
event start time determination and the FSR simulation are quite high
($>99\%$), their systematic uncertainties are estimated to be less than 1%
otherserrors . The total systematic uncertainties are given by the quadratic
sum of the individual uncertainties, assuming all sources to be independent.
The uncertainty due to the fit procedure is considered during the upper limit
determination described in the following.

Table 1: Summary of systematic uncertainties (in %) for $J/\psi\to D^{0}e^{+}e^{-}$ and $\psi\to D^{0}e^{+}e^{-}$, the sources tagged with ${}^{\prime}{\ast}^{\prime}$ are correlated among the different $D^{0}$ decay modes. The hyphen ($-$) indicates the source does not contribute to the channel. |  $D^{0}\to K^{-}\pi^{+}$ |  $D^{0}\to K^{-}\pi^{+}\pi^{0}$ |  $D^{0}\to K^{-}\pi^{+}\pi^{+}\pi^{-}$
---|---|---|---
$J/\psi$ |  $\psi(3686)$ |  $J/\psi$ |  $\psi(3686)$ |  $J/\psi$ |  $\psi(3686)$
Tracking∗ | 4.0 | 4.0 | 4.0 | 4.0 | 6.0 | 6.0
PID∗ | 6.0 | 6.0 | 6.0 | 6.0 | 8.0 | 8.0
$\gamma$ detection | – | – | 1.2 | 1.2 | – | –
Kinematic fit | 1.7 | 1.6 | 1.1 | 1.8 | 2.2 | 2.0
Veto $\gamma$ conversion∗ | 1.7 | 1.7 | 1.7 | 1.7 | 1.7 | 1.7
Veto $K_{S}\to\pi^{0}\pi^{0}$ | – | – | 0.6 | – | – |
Veto $K_{S}\to\pi^{+}\pi^{-}$ | – | – | – | – | 2.1 | 2.2
Veto $J/\psi\to e^{+}e^{-}$ | – | 0.1 | – | – |  | –
Branching fraction | 1.3 | 1.3 | 3.6 | 3.6 | 2.6 | 2.6
$\psi$ total number∗ | 0.55 | 0.62 | 0.55 | 0.62 | 0.55 | 0.62
Others | 1.0 | 1.0 | 1.0 | 1.0 | 1.0 | 1.0
Total | 7.8 | 7.8 | 8.5 | 8.7 | 11.0 | 10.9

Since no significant signal for $\psi\to D^{0}e^{+}e^{-}$ is observed, upper
limits at the 90% C.L. on the branching fractions are determined.
Simultaneous, unbinned maximum likelihood fits on the distributions of
invariant masses $M(K^{-}\pi^{+})$, $M(K^{-}\pi^{+}\pi^{0})$, and
$M(K^{-}\pi^{+}\pi^{+}\pi^{-})$, are carried out for the $J/\psi$ and
$\psi(3686)$ samples. In the fit, the signal shapes are described by the
corresponding signal MC samples and the background shapes are described by 2nd
order polynomial functions. The expected number of signal events in the
$i^{\rm th}$ decay mode is calculated with
$N_{i}=N_{\psi}~{}\cdot~{}\mathcal{B}~{}\cdot~{}\mathcal{B}^{\rm
inter}_{i}~{}\cdot~{}\epsilon_{i}$, where $N_{\psi}$ is the total number of
$\psi$ events, $\mathcal{B}_{i}^{\rm inter}$ is the product of the decay
branching fractions of $D^{0}$ mesons and subsequent intermediate states,
taken from the PDG PDG , and $\epsilon_{i}$ is the detection efficiency from
the signal MC samples. The decay branching fraction $\mathcal{B}$ of $\psi\to
D^{0}e^{+}e^{-}$ is a common parameter among the three $D^{0}$ decay modes.
The overall likelihood values ($\mathcal{L}$) are the products of those of the
three $D^{0}$ decay modes, incorporating systematic uncertainties, which are
separated as correlated and uncorrelated upper1 ; upper2 . The likelihood fits
are carried out with the MINUIT package minuit .

We compute the upper limits on the branching fraction at the 90% C.L. using a
Bayesian method PDG with a flat prior. The optimized likelihoods,
$\mathcal{L}$, are presented as a function of branching fraction
$\mathcal{B}(\psi\to D^{0}e^{+}e^{-})$. The upper limits on the branching
fractions $\mathcal{B}^{\rm UP}$ at the 90% C.L. are the values that yield 90%
of the likelihood integral over $\mathcal{B}$ from zero to infinity:
$\int_{0}^{\mathcal{B}^{\rm
UP}}{\mathcal{L}d\mathcal{B}}/\int_{0}^{\infty}{\mathcal{L}d\mathcal{B}}=0.9$.
To take into account the systematic uncertainties related to the fit process,
two alternative fit scenarios are considered: (1) changing the fit range on
the invariant masses by 10 MeV/$c^{2}$; or (2) replacing the 2nd order
polynomial function with a 3rd order polynomial function for the background.
We try all combinations of the different scenarios. The one with the maximum
upper limits on the branching fractions is taken as the conservative result.
The upper limits at the 90% C.L. on the branching fractions are
$\mathcal{B}(J/\psi\to D^{0}e^{+}e^{-})<8.5\times 10^{-8}$ and
$\mathcal{B}(\psi(3686)\to D^{0}e^{+}e^{-})<1.4\times 10^{-7}$, respectively.
The corresponding normalized likelihood distributions are shown in Fig. 3 and
the best fit curves are shown in Fig. 2.

\begin{overpic}[width=99.73074pt,height=86.72267pt]{likelihoodjpsidrafhl.eps}
\end{overpic}\begin{overpic}[width=99.73074pt,height=86.72267pt]{likelihoodpsidraftmhl.eps}
\end{overpic}

Figure 3: Normalized combined likelihood as a function of the decay branching
fraction $\mathcal{B}(\psi\to D^{0}e^{+}e^{-})$ for $J/\psi$ (left) and
$\psi(3686)$ (right) samples, where the correlated and un-correlated
systematic uncertainties are incorporated. The likelihood function is
normalized with the maximum to be 1. The blue arrow denotes the 90% C.L.

In summary, we perform a search for the rare decays of $J/\psi\to
D^{0}e^{+}e^{-}$ and $\psi(3686)\to D^{0}e^{+}e^{-}$ using samples of
($1310.6\pm 7.2$) $\times$ $10^{6}$ $J/\psi$ events and ($448.1\pm 2.9$)
$\times$ $10^{6}$ $\psi(3686)$ events collected with the BESIII detector. No
significant signal is observed and upper limits at the 90% C.L. for the
branching fractions are determined to be $\mathcal{B}(J/\psi\to
D^{0}e^{+}e^{-})<8.5\times 10^{-8}$ and $\mathcal{B}(\psi(3686)\to
D^{0}e^{+}e^{-})<1.4\times 10^{-7}$, respectively. The limit on
$\mathcal{B}(J/\psi\to D^{0}e^{+}e^{-})$ is more stringent by two orders in
magnitude compared to the previous results, and the $\mathcal{B}(\psi(3686)\to
D^{0}e^{+}e^{-})$ is set for the first time. Though the upper limits are
larger than the SM predictions, they may help to discriminate the different
new physics models or to constrain the parameters in the different physics
models. Additionally, higher statistics $J/\psi$ and $\psi(3686)$ samples may
help to improve the sensitivity of the measurements.

The BESIII collaboration thanks the staff of BEPCII, the IHEP computing center
and the Supercomputing center of USTC for their strong support. This work is
supported in part by National Key Basic Research Program of China under
Contract No. 2015CB856700; National Natural Science Foundation of China (NSFC)
under Contracts Nos. 11235011, 11322544, 11335008, 11425524, 11322544,
11375170, 11275189, 11475169, 11475164, 11575133, 11625523, 11635010; the
Chinese Academy of Sciences (CAS) Large-Scale Scientific Facility Program; the
CAS Center for Excellence in Particle Physics (CCEPP); the Collaborative
Innovation Center for Particles and Interactions (CICPI); Joint Large-Scale
Scientific Facility Funds of the NSFC and CAS under Contracts Nos. 11179007,
U1332201, U1532257, U1532258; CAS under Contracts Nos. KJCX2-YW-N29,
KJCX2-YW-N45, QYZDJ-SSW-SLH003; 100 Talents Program of CAS; National 1000
Talents Program of China; INPAC and Shanghai Key Laboratory for Particle
Physics and Cosmology; German Research Foundation DFG under Contracts Nos.
Collaborative Research Center CRC 1044, FOR 2359; Istituto Nazionale di Fisica
Nucleare, Italy; Koninklijke Nederlandse Akademie van Wetenschappen (KNAW)
under Contract No. 530-4CDP03; Ministry of Development of Turkey under
Contract No. DPT2006K-120470; National Natural Science Foundation of China
(NSFC) under Contracts Nos. 11505034, 11575077; National Science and
Technology fund; The Swedish Research Council; U. S. Department of Energy
under Contracts Nos. DE-FG02-05ER41374, DE-SC-0010118, DE-SC-0010504, DE-
SC-0012069; University of Groningen (RuG) and the Helmholtzzentrum fuer
Schwerionenforschung GmbH (GSI), Darmstadt; WCU Program of National Research
Foundation of Korea under Contract No. R32-2008-000-10155-0

## References

  * (1) S. L. Glashow, J. Iliopoulos and L. Maiani, Phys. Rev. D 2, 1285 (1970).
  * (2) Y. M. Wang et al., J. Phys. G 36, 105002 (2009).
  * (3) M. A. Sanchis-Lonzano, Z. Phys. C 62, 271 (1994).
  * (4) Y. M. Wang et al., Eur. Phys. J. C 54, 107 (2008) .
  * (5) C. Hill, Phys. Lett. B 345, 483 (1995).
  * (6) C. S. Aulakh and R. N. Mohapatra, Phys. Lett. B 119, 136 (1982).
  * (7) S. Glashow and S. Weinberg, Phys. Rev. D 15, 1858 (1977).
  * (8) X. Zhang, HEP&NP, 25, 461 (2001).
  * (9) A. Datta, P. J. Odonnel, S. Pakvasa and X. Zhang, Phys. Rev. D 60, 014011 (1999).
  * (10) M. Ablikim et al. [BES Collaboration], Phys. Lett. B 639, 418 (2006).
  * (11) M. Ablikim et al. [BESIII Collaboration], Nucl. Instrum. Meth. A 614, 3 (2010).
  * (12) M. Ablikim et al. [BESIII Collaboration], Chin. Phys. C 36, 915 (2012).
  * (13) M. Ablikim et al. [BESIII Collaboration], Chin. Phys. C 41, 013001 (2017).
  * (14) M. Ablikim et al. [BESIII Collaboration], Chin. Phys. C 37, 063001 (2013).
  * (15) M. Ablikim et al. [BESIII Collaboration], arXiv:1709.03653 [hep-ex], submitted to Chin. Phys. C.
  * (16) S. Agostinelli et al. [geant4 Collaboration], Nucl Instrum. Meth. A 506, 250 (2003).
  * (17) Z. Y. Deng et al., Chin. Phys. C 30, 371 (2006).
  * (18) S. Jadach, B. F. L. Ward and Z. Was, Comput. Phys. Commun. 130, 260 (2000); Phys. Rev. D 63, 113009 (2001).
  * (19) D. J. Lange, Nucl. Instrum. Meth. A 462, 152 (2001); R. G. Ping et al., Chin. Phys. C 32, 599 (2008).
  * (20) C. Patrignani et al. [Particle Data Group], Chin. Phys. C, 40, 100001 (2016).
  * (21) J. C. Chen et al., Phys. Rev. D 62, 034003 (2000).
  * (22) E. Barberio and Z. Was, Comput. Phys. Commun. 79, 219 (1994).
  * (23) V. M. Budnev and V. A. Karnakov, Pisma Zh. Eksp. Teor. Fiz. 29, 439 (1979).
  * (24) Z. R. Xu and K. L. He, Chin. Phys. C 36, 742 (2012).
  * (25) M. Ablikim et al. [BESIII Collaboration], Phys. Rev. D 83, 112005 (2011).
  * (26) M. Ablikim et al. [BESIII Collaboration], Phys. Rev. D 85, 092012 (2012).
  * (27) M. Ablikim et al. [BESIII Collaboration], Phys. Rev. D 87, 012002 (2013).
  * (28) R. G. Ping, Gang Li and Z. Wang, Formalism of Helicity Coupling Amplitudes for $J/\psi\to\pi^{+}\pi^{-}\pi^{0}$. Commun. Theor. Phys. (Beijing, China) 47 (2007).
  * (29) M. Ablikim et al. [BESIII Collaboration], Phys. Rev. Lett. 110, 252002 (2013).
  * (30) M. R. Convery, Incorporating multiplicative systematic errors in branching ratio limits. SLAC-TN-03-001, 2003.
  * (31) K. Stenson, A more exact solution for incorporating multiplicative systematic uncertainties in branching ratio limits. arXiv:0605236 (2006).
  * (32) F. James and M. Roos, Comput. Phys. Commun. 10, 343 (1975).