# Feasibility of searching for the Cabibbo-favored
$D^{\ast}$ ${\to}$ $\overline{K}{\pi}^{+}$, $\overline{K}^{\ast}{\pi}^{+}$,
$\overline{K}{\rho}^{+}$ decays

Yueling Yang  Institute of Particle and Nuclear Physics, Henan Normal
University, Xinxiang 453007, China  Kang Li  Institute of Particle and Nuclear
Physics, Henan Normal University, Xinxiang 453007, China  Zhenglin Li
Institute of Particle and Nuclear Physics, Henan Normal University, Xinxiang
453007, China  Jinshu Huang  School of Physics and Electronic Engineering,
Nanyang Normal University, Nanyang 473061, China  Qin Chang  Institute of
Particle and Nuclear Physics, Henan Normal University, Xinxiang 453007, China
Junfeng Sun  Institute of Particle and Nuclear Physics, Henan Normal
University, Xinxiang 453007, China

###### Abstract

The current knowledge on the $D^{\ast}$ mesons are still inadequate.
Encouraged by the positive development prospects of high-luminosity and high-
precision experiments, the Cabibbo-favored nonleptonic $D^{\ast}$ ${\to}$
$\overline{K}{\pi}^{+}$, $\overline{K}^{\ast}{\pi}^{+}$,
$\overline{K}{\rho}^{+}$ weak decays are studied with the naive factorization
approach. It is found that branching ratios of these processes can reach up to
${\cal O}(10^{-10})$ or more, and can be accessible at STCF, CEPC, FCC-ee and
LHCb@HL-LHC experiments in the future. It might even be possible to search for
the $D^{{\ast}0}$ ${\to}$ $K^{{\ast}-}{\pi}^{+}$ and $K^{-}{\rho}^{+}$ decays
at the running SuperKEKB experiments.

Phys. Rev. D 106, 036029 (2022)

$D^{\ast}$ meson; weak decay; branching ratio.

##  I Introduction

The need for a charm quark as its natural incorporation in a unified
description of the weak and electromagnetic interactions, was first postulated
by Bjorken, Glashow, Iliopoulos, and Maiani Phys.Lett.11.255 ; PhysRevD.2.1285
. The charm flavor is strictly conserved in electromagnetic and strong
interactions. Although the $J/{\psi}$ particle was first discovered by two
vastly different experiments in 1974 Phys.Rev.Lett.33.1404 ;
Phys.Rev.Lett.33.1406 , and interpreted as a bound system of the $c\bar{c}$
pair, the charm flavor is hidden away in the quantity $R$ near the narrow
resonance $J/{\psi}$ peak. The observable quantity $R$ measures the total
hadron production rate relative to the electrodynamic ${\mu}^{+}{\mu}^{-}$
process in $e^{+}e^{-}$ collisions, $R$ ${\equiv}$
${\sigma}(e^{+}e^{-}{\to}{\rm
hadrons})/{\sigma}(e^{+}e^{-}{\to}{\mu}^{+}{\mu}^{-})$, determines the
properties of the species of quark carrying color and fractional electric
charge. The value of $R$ with the lowest order approximation should remain
approximately constant with center-of-mass energy of the $e^{+}e^{-}$ pair
annihilation, as long as no new final state fermion pair appears. From the
view of $R$, the most striking evidence for the existence of a charm quark is
the discovery of charmed mesons in positron-electron annihilation at SPEAR
experiments in 1976 Phys.Rev.Lett.37.255 ; Phys.Rev.Lett.37.569 . The clearly
evident increase in $R$ near the $4$ GeV region in Fig. 1 has been widely
taken by particle physicists to be the characteristic features of charm
threshold.

Figure 1: $R$ near the charm threshold regions versus the center-of-mass
energy $\sqrt{s}$ of $e^{+}e^{-}$ annihilation, corresponding to Fig.53.3 of
Ref. pdg2020 , where the experimental data are available at web page
https://pdg.lbl.gov/2021/hadronic-xsections pdg2020 , the lower and upper
horizontal dashed lines correspond to $R$ $=$ $2$ and $10/3$, respectively;
and the vertical dashed lines denote charm thresholds.

The ordinary charmed mesons consist of the charmed quark plus a light
nonstrange quark forming an isodoublet, $c\bar{u}$ and $c\bar{d}$, symbolized
respectively by $D^{0}$ and $D^{+}$ corresponding to the $1^{1}S_{0}$ state
with spin-parity quantum $J^{P}$ $=$ $0^{-}$ for ground pseudoscalar mesons,
$D^{{\ast}0}$ and $D^{{\ast}+}$ corresponding to the $1^{3}S_{1}$ state with
spin-parity quantum $J^{P}$ $=$ $1^{-}$ for ground vector mesons. In electron-
positron collisions, the charmed quark pair originated from the virtual
timelike photon combines with the light $u$ and $d$ quarks out of the vacuum
with the orbit-spin $L$-$S$ coupling and forms the charmed mesons through
strong interactions. In principle, the charmed meson pairs are equally likely
to be neutral and charged. It is clearly seen from Fig. 1 that there are many
broader resonances above charm threshold, such as ${\psi}(3770)$,
${\psi}(4040)$, ${\psi}(4160)$ and so on. These broader sibling ${\psi}$
resonances presumably decay strongly to pairs of oppositely charmed mesons,
$D\bar{D}$, $D\bar{D}^{\ast}$, and $D^{\ast}\bar{D}^{\ast}$. The relative
production cross sections of charmed mesons close to threshold were once
thought to be, for example,
${\sigma}(D\bar{D})$:${\sigma}(D\bar{D}^{\ast})$:${\sigma}(D^{\ast}\bar{D}^{\ast})$
${\simeq}$ 1:4:7 PhysRevLett.37.477 ; PhysRevLett.37.398 ; Phys.Lett.B.65.55 ;
Phys.Lett.B.66.70  or other ratios Phys.Lett.B.70.106 ; PhysRevD.15.759 ;
PhysRevD.16.1327 ; PhysRevD.17.765  based on different theoretical
calculation. The measured production ratio may differ considerably from those
theoretical estimates. At energies a little far above the
$D^{\ast}\bar{D}^{\ast}$ threshold, it is generally believed that the
production cross sections of $D^{\ast}$ mesons are greater than those of $D$
mesons in $e^{+}e^{-}$ collisions, which is clearly seen from Fig. 2.

Figure 2: Measured production cross sections for charmed meson pairs near the
threshold regions versus the center-of-mass energy $\sqrt{s}$ of $e^{+}e^{-}$
annihilation, where the experimental data are available from the CLEO-c group
in Refs. PhysRevD.80.072001 ; cpc.42.043002 ; the vertical dashed lines denote
charm thresholds.

The center-of-mass energy for $D^{\ast}\bar{D}^{\ast}$ production is larger
than that for $D\bar{D}$ production. Due to the engineering design of
$e^{+}e^{-}$ accelerators and detectors, plus the actual operation energy and
time, more experimental data on the pseudoscalar $D$ mesons rather than the
vector $D^{\ast}$ mesons have been accumulated at CLEC-c and BES experiments.
The properties of $D$ mesons have been extensively and carefully studied by
many experimental groups. To go along with this, the experimental information
on $D^{\ast}$ mesons is still very limited by now pdg2020 . The measurement of
the mass of $D^{\ast}$ mesons which was quoted but not adopted by Particle
Data Group was carried out in 1977 PhysRevLett.39.1301 ; Phys.Lett.B.69.503 ,
not updated for 45 years. In addition, only three decay modes of $D^{{\ast}0}$
($D^{{\ast}+}$) mesons have been observed pdg2020 .

The fitted $D^{\ast}$ masses are $m_{D^{{\ast}0}}$ $=$ $2006.85(5)$ MeV and
$m_{D^{{\ast}+}}$ $=$ $2010.26(5)$ MeV pdg2020 . There is a hierarchical
relationship among hadron mass,

$m_{D^{{\ast}0}}-m_{D^{0}}\,>\,m_{{\pi}^{0}},$ (1)
$m_{D^{{\ast}0}}-m_{D^{+}}\,<\,m_{{\pi}^{-}},$ (2)
$m_{D^{{\ast}+}}-m_{D^{0}}\,>\,m_{{\pi}^{+}},$ (3)
$m_{D^{{\ast}+}}-m_{D^{+}}\,>\,m_{{\pi}^{0}},$ (4)

which results in the $D^{{\ast}0}$ ${\to}$ $D^{+}{\pi}^{-}$ decay being
kinematically forbidden. The $D^{{\ast}}$ mesons can decay mainly through the
strong and electromagnetic interactions, and are dominated by $D^{\ast}$
${\to}$ $D{\pi}$ processes. Kinematically, the $D^{\ast}$ ${\to}$ $D{\pi}$
decays are inhibited by very compact phase space. Dynamically, for the
$D^{\ast}$ ${\to}$ $D{\pi}$ decays, the ${\pi}$ emission processes are
suppressed by the phenomenological Okubo–Zweig–Iizuka rule ozi-o ; ozi-z ;
ozi-i , and the recoiled ${\pi}$ processes are suppressed by a color match
factor owing to the color singlet requirements of quark combinations.
Consequently, the decay width of $D^{\ast}$ mesons is very narrow, for
example, ${\Gamma}_{D^{{\ast}+}}$ $=$ $83.4{\pm}1.8$ keV pdg2020 .
Experimentally, the momentum of the pion in the rest frame of $D^{\ast}$
mesons is very soft, about $40$ MeV. The difficulties in signal reconstruction
result in inefficiency of particle identification.

Besides, the $D^{{\ast}}$ meson decay through weak interactions is also legal
and allowable within the standard model (SM) of elementary particle. As it is
well known, the strong and electromagnetic interactions are related only to
vector currents, while the weak interactions are related to both vector and
axial vector currents. Study of the $D^{{\ast}}$ meson weak decays can not
only improve our knowledge of the properties $D^{{\ast}}$ mesons, but also
test the axial vector current interactions in SM and enrich our understanding
on the decay mechanism of $D^{{\ast}}$ mesons. In addition, the Cabibbo-
Kobayashi-Maskawa (CKM) matrix elements PhysRevLett.10.531 ; PTP.49.652
describing the quark mixing can be determined and overconstrained from
$D^{{\ast}}$ meson weak decays.

It is not hard to visualize that the occurrence probability of $D^{{\ast}}$
meson weak decays is very small, insignificant when compared with that of
strong and electromagnetic decays. Fortunately, the experimental data on
$D^{{\ast}}$ mesons are accumulating increasingly. Now, more than
$5{\times}10^{7}$ $D^{{\ast}{\pm}}$ mesons have been accumulated in energy
regions $\sqrt{s}$ ${\in}$ $[4.1,4.6]$ at BESIII experiments JHEP2022.05.155 .
It is promisingly expected that there will be a total of about
$5{\times}10^{10}$ $c\bar{c}$ pairs at the SuperKEKB PTEP.2019.123C01 . Given
the charm quark fragmentation fractions $f(c{\to}D^{{\ast}+})$ ${\approx}$
$25\%$ and $f(c{\to}D^{{\ast}0})$ ${\approx}$ $23\%$ epjc.76.397 , more than
$2{\times}10^{10}$ $D^{{\ast}0}$ and $D^{{\ast}{\pm}}$ will be available at
SuperKEKB. If it is considered optimistically that the production cross
sections ${\sigma}(e^{+}e^{-}{\to}D^{0}\bar{D}^{{\ast}0})$ ${\approx}$
${\sigma}(e^{+}e^{-}{\to}D^{+}D^{{\ast}-})$ ${\approx}$ $2$ nb and
${\sigma}(e^{+}e^{-}{\to}D^{{\ast}0}\bar{D}^{{\ast}0})$ ${\approx}$
${\sigma}(e^{+}e^{-}{\to}D^{{\ast}+}D^{{\ast}-})$ ${\approx}$ $3$ nb near
${\psi}(4040)$ resonance, as illustrated in Fig. 2 (or Table 1 and 2 of Ref.
JHEP2022.05.155 ), about $8{\times}10^{10}$ $D^{{\ast}0}$ and
$D^{{\ast}{\pm}}$ mesons will be available at the super ${\tau}$-charm factory
(STCF) STCF ; SCTF  with a total integrated luminosity of $10\,{\rm ab}^{-1}$.
It is highly expected that there will be a total of about $10^{12}$ $Z^{0}$
bosons at the Circular Electron Positron Collider (CEPC) cepc  and $10^{13}$
$Z^{0}$ bosons at the Future Circular Collider (FCC-ee) fcc . Considering the
inclusive branching ratios ${\cal B}(Z^{0}{\to}D^{{\ast}{\pm}}X)$ $=$
$(11.4{\pm}1.3)\%$ pdg2020  and the approximation ${\cal
B}(Z^{0}{\to}D^{{\ast}0}X/\bar{D}^{{\ast}0}X)$ ${\approx}$ ${\cal
B}(Z^{0}{\to}D^{{\ast}{\pm}}X)$, more than $10^{11}$ and $10^{12}$
$D^{{\ast}}$ mesons will be available at CEPC and FCC-ee with a total
integrated luminosity of $20\,{\rm ab}^{-1}$, respectively. In addition, the
inclusive production cross sections in hadron-hadron collisions are measured
to be ${\sigma}(pp{\to}D^{{\ast}+}X)$ ${\approx}$ $0.8$ mb at center-of-mass
energy of $\sqrt{s}$ $=$ $13$ TeV by the LHCb group JHEP.2016.03.159  and
${\sigma}_{\rm tot}(D^{{\ast}+})$ ${\approx}$ $2.1$ mb at $\sqrt{s}$ $=$ $7$
TeV by the ALICE group JHEP.2012.07.191 , respectively. A conservative
estimate is that more than $2{\times}10^{14}$ $D^{{\ast}}$ mesons will be
available by the LHCb detector at the future High Luminosity LHC (HL-LHC)
hadron collider with a total integrated luminosity of $300\,{\rm fb}^{-1}$
1808.08865 . A huge amount of experimental data provides a solid foundation
and valuable opportunities for studying and understanding the properties of
$D^{{\ast}}$ mesons, including the search for $D^{{\ast}}$ meson weak decays.

It is natural to want to know the probability of $D^{{\ast}}$ meson weak
decays, and whether the study on the $D^{{\ast}}$ meson weak decays is
feasible or not in the coming future. Our study Eur.Phys.J.C.81.1110  has
tentatively shown that the purely leptonic $D^{{\ast}+}$ decays might be
measurable at $e^{+}e^{-}$ colliders, such as SuperKEKB, STCF, CEPC and FCC-
ee, and HL-LHC hadron collider 1808.08865 . Purely leptonic $D^{{\ast}+}$
decays are Cabibbo-suppressed, while leptonic $D^{{\ast}0}$ decays are induced
by the flavor-changing neutral weak currents. In this paper, we would like to
investigate the $D^{{\ast}}$ ${\to}$ $\overline{K}{\pi}^{+}$,
$\overline{K}{\rho}^{+}$ and $\overline{K}^{\ast}{\pi}^{+}$ decays within the
SM. The motivation is as follows. Dynamically, these decays are induced by $W$
boson emission, and favored by CKM elements, and thus should have relatively
larger branching ratios and higher priority to do research among $D^{{\ast}}$
meson weak decays. Kinematically, the final states are back-to-back, and have
definitive momentum and discrete energy in the rest frame of $D^{{\ast}}$
mesons. Experimentally, the energetic charged pion and kaon of final states of
nonleptonic $D^{{\ast}}$ weak decays are more easily captured by the sensitive
detectors, when compared with $D^{\ast}$ ${\to}$ $D{\pi}$ decays where the
soft pion identifications are very challenging, and when compared with
leptonic $D^{{\ast}+}$ decays where the final neutrinos will bring the signal
reconstruction into additional complications. Hence, there would be a
relatively higher reconstruction efficiency for nonleptonic $D^{{\ast}}$ weak
decays.

The remainder of this paper is organized as follows. Theoretical framework and
phenomenological approach dealing with nonleptonic $D^{{\ast}}$ weak decays
are presented in Section II. Numerical results of branching ratios and our
comments are given in Section III. Section IV is a brief summary. The analytic
expressions of decay amplitudes are collected in Appendixes A and B.

##  II Theoretical framework

The effective Hamiltonian in charge of nonleptonic $D^{{\ast}}$ ${\to}$
$\overline{K}{\pi}^{+}$, $\overline{K}{\rho}^{+}$ and
$\overline{K}^{\ast}{\pi}^{+}$ decays is written as RevModPhys.68.1125 ,

${\cal H}_{\rm
eff}\,=\,\frac{G_{F}}{\sqrt{2}}\,V_{cs}^{\ast}\,V_{ud}\,\big{\\{}C_{1}\,O_{1}+C_{2}\,O_{2}\big{\\}}+{\rm
h.c.},$ (5)

where both Fermi constant $G_{F}$ ${\approx}$ $1.166{\times}10^{-5}$ ${\rm
GeV}^{-2}$ pdg2020  and Wilson coefficients $C_{1,2}$ are universal effective
couplings. Wilson coefficients are computable with the perturbation theory at
the mass scale of $W$ boson. The scale-dependence of Wilson coefficients can
be obtained with the renormalization group equation RevModPhys.68.1125 . The
nature of the product of $G_{F}\,C_{i}$ is similar to a gauge parameter, for
example, the electric charge $e$ for electromagnetic interactions. The CKM
elements correspond to different effective operators, and have been determined
precisely, ${|}V_{cs}{|}$ $=$ $0.987(11)$, and ${|}V_{ud}{|}$ $=$
$0.97370(14)$ pdg2020 . The explicit expressions of effective four-quark
operators are written as,

$\displaystyle O_{1}$ $\displaystyle=$
$\displaystyle\big{[}\bar{s}_{\alpha}\,{\gamma}^{\mu}\,(1-{\gamma}_{5})\,c_{\alpha}\big{]}\,\big{[}\bar{u}_{\beta}\,{\gamma}_{\mu}\,(1-{\gamma}_{5})\,d_{\beta}\big{]},$
(6) $\displaystyle O_{2}$ $\displaystyle=$
$\displaystyle\big{[}\bar{s}_{\alpha}\,{\gamma}^{\mu}\,(1-{\gamma}_{5})\,c_{\beta}\big{]}\,\big{[}\bar{u}_{\beta}\,{\gamma}_{\mu}\,(1-{\gamma}_{5})\,d_{\alpha}\big{]},$
(7)

where ${\alpha}$ and ${\beta}$ are the color indices. Owing to only one kaon
meson in final states, there is no penguin operators.

Taking the $D^{{\ast}0}$ ${\to}$ $K^{-}{\pi}^{+}$ decay for example, the
amplitude can be written as,

${\cal A}(D^{{\ast}0}{\to}K^{-}{\pi}^{+})\,=\,{\langle}K^{-}{\pi}^{+}{|}{\cal
H}_{\rm
eff}{|}D^{{\ast}0}{\rangle}\,=\,\frac{G_{F}}{\sqrt{2}}\,V_{cs}^{\ast}\,V_{ud}\,\sum\limits_{i=1}^{2}C_{i}\,{\langle}K^{-}{\pi}^{+}{|}O_{i}{|}D^{{\ast}0}{\rangle}.$
(8)

Clearly, the remaining part of decay amplitudes is to reasonably evaluate the
hadronic matrix elements (HMEs)
${\langle}K^{-}{\pi}^{+}{|}O_{i}{|}D^{{\ast}0}{\rangle}$, which sandwich the
initial and final states together with quark operators. HMEs closely relate to
transitions between quarks and hadron states, and involve short- and long-
distance contributions which complicate theoretical calculations.

Phenomenologically, the naive factorization (NF) approach zpc.34.103  was
often used to deal with HMEs, owing to its simple and clear physics picture,
and good performance for nonleptonic $B$ and $D$ weak decays induced by
external $W$ emission processes. Based on the color transparency hypothesis
npbps.11.325  that energetic final hadron states have flown far away from each
other before influence of soft gluons, HMEs of four-quark operators using NF
approach are divided into two HMEs of hadron currents, and the final state
interactions can be dismissed for the time being. HMEs of hadron currents are
conventionally parametrized by decay constants and hadron transition form
factors, which can be obtained from data and nonperturbative methods. In
addition, it has been shown in Ref. PhysRevD.81.074021  that although the
contributions from $W$-exchange and $W$-annihilation diagrams are the same
order as those from tree-emission diagrams due to the large non-perturbative
contributions for charmed meson decays, the tree-emission diagrams always give
the largest contributions in the total amplitudes, thereby give the correct
orders of magnitude of branching fractions. In this paper, we will apply NF
approach to the concerned nonleptonic $D^{\ast}$ weak decays.

The parametrization schemes of hadron current HMEs are generally written as
zpc.29.637 ; ijmpa.30.1550094 ; jhep.1912.102 ,

$\displaystyle{\langle}0{|}\bar{q}_{1}\,{\gamma}_{\mu}\,q_{2}{|}P(k){\rangle}$
$\displaystyle=$ $\displaystyle 0,$ (9)
$\displaystyle{\langle}0{|}\bar{q}_{1}\,{\gamma}_{\mu}\,{\gamma}_{5}\,q_{2}{|}P(k){\rangle}$
$\displaystyle=$ $\displaystyle i\,f_{P}\,k_{\mu},$ (10)
$\displaystyle{\langle}0{|}\bar{q}_{1}\,{\gamma}_{\mu}\,q_{2}{|}V(k,{\epsilon}){\rangle}$
$\displaystyle=$ $\displaystyle f_{V}\,m_{V}\,{\epsilon}_{\mu},$ (11)
$\displaystyle{\langle}0{|}\bar{q}_{1}\,{\gamma}_{\mu}\,{\gamma}_{5}\,q_{2}{|}V(k,{\epsilon}){\rangle}$
$\displaystyle=$ $\displaystyle 0,$ (12)
${\langle}P(p_{2}){|}\,\bar{q}\,{\gamma}_{\mu}\,c\,{|}D^{\ast}({\epsilon}_{1},p_{1}){\rangle}\,=\,{\varepsilon}_{{\mu}{\nu}{\alpha}{\beta}}\,{\epsilon}_{1}^{\nu}\,P^{\alpha}\,q^{\beta}\frac{V^{D^{\ast}P}(q^{2})}{m_{D^{\ast}}+m_{P}},$
(13)
$\displaystyle{\langle}P(p_{2}){|}\,\bar{q}\,{\gamma}_{\mu}\,{\gamma}_{5}\,c\,{|}D^{\ast}({\epsilon}_{1},p_{1}){\rangle}$
(14) $\displaystyle=$
$\displaystyle+i\,2\,m_{D^{\ast}}\,\frac{{\epsilon}_{1}{\cdot}q}{q^{2}}\,q_{\mu}\,A_{0}^{D^{\ast}P}(q^{2})+i\,{\epsilon}_{1,{\mu}}\,(m_{D^{\ast}}+m_{P})\,A_{1}^{D^{\ast}P}(q^{2})$
$\displaystyle+i\,\frac{{\epsilon}_{1}{\cdot}q}{m_{D^{\ast}}+m_{P}}\,P_{\mu}\,A_{2}^{D^{\ast}P}(q^{2})-i\,2\,m_{D^{\ast}}\,\frac{{\epsilon}_{1}{\cdot}q}{q^{2}}\,q_{\mu}\,A_{3}^{D^{\ast}P}(q^{2}),$
$\displaystyle{\langle}V({\epsilon}_{2},p_{2}){|}\,\bar{q}\,{\gamma}_{\mu}\,c\,{|}D^{\ast}({\epsilon}_{1},p_{1}){\rangle}\,=\,-({\epsilon}_{1}{\cdot}{\epsilon}_{2}^{\ast})\,\big{\\{}P_{\mu}\,V_{1}^{D^{\ast}V}(q^{2})-q_{\mu}\,V_{2}^{D^{\ast}V}(q^{2})\big{\\}}$
$\displaystyle+\frac{({\epsilon}_{1}{\cdot}q)\,({\epsilon}_{2}^{\ast}{\cdot}q)}{m_{D^{\ast}}^{2}-m_{V}^{2}}\,\big{\\{}\big{[}P_{\mu}-\frac{m_{D^{\ast}}^{2}-m_{V}^{2}}{q^{2}}\,q_{\mu}\big{]}\,V_{3}^{D^{\ast}V}(q^{2})+\frac{m_{D^{\ast}}^{2}-m_{V}^{2}}{q^{2}}\,q_{\mu}\,V_{4}^{D^{\ast}V}(q^{2})\big{\\}}$
$\displaystyle-({\epsilon}_{1}{\cdot}q)\,{\epsilon}_{2,{\mu}}^{\ast}\,V_{5}^{D^{\ast}V}(q^{2})+({\epsilon}_{2}^{\ast}{\cdot}q)\,{\epsilon}_{1,{\mu}}\,V_{6}^{D^{\ast}V}(q^{2}),$
(15)
$\displaystyle{\langle}V({\epsilon}_{2},p_{2}){|}\,\bar{q}\,{\gamma}_{\mu}\,{\gamma}_{5}\,c\,{|}D^{\ast}({\epsilon}_{1},p_{1}){\rangle}$
(16) $\displaystyle=$
$\displaystyle-i\,{\varepsilon}_{{\mu}{\nu}{\alpha}{\beta}}\,{\epsilon}_{1}^{\alpha}\,{\epsilon}_{2}^{{\ast}{\beta}}\,\big{\\{}\big{[}P^{\nu}-\frac{m_{D^{\ast}}^{2}-m_{V}^{2}}{q^{2}}\,q^{\nu}\big{]}\,A_{1}^{D^{\ast}V}(q^{2})+\frac{m_{D^{\ast}}^{2}-m_{V}^{2}}{q^{2}}\,q^{\nu}\,A_{2}^{D^{\ast}V}(q^{2})\big{\\}}$
$\displaystyle-\frac{i\,{\varepsilon}_{{\mu}{\nu}{\alpha}{\beta}}\,P^{\alpha}\,q^{\beta}}{m_{D^{\ast}}^{2}-m_{V}^{2}}\,\big{\\{}({\epsilon}_{2}^{\ast}{\cdot}q)\,{\epsilon}_{1}^{\nu}\,A_{3}^{D^{\ast}V}(q^{2})-({\epsilon}_{1}{\cdot}q)\,{\epsilon}_{2}^{{\ast},{\nu}}\,A_{4}^{D^{\ast}V}(q^{2})\big{\\}},$

where $m_{P}$ and $f_{P}$ are the mass and decay constant of pseudoscalar $P$
meson, respectively; $m_{V}$, $f_{V}$, and ${\epsilon}_{V}$ are the mass,
decay constant, and polarization vector of vector $V$ meson, respectively;
$A_{i}^{D^{\ast}h}$ and $V_{i}^{D^{\ast}h}$ are mesonic transition form
factors; The momentum $P$ $=$ $p_{1}$ $+$ $p_{2}$ and $q$ $=$ $p_{1}$ $-$
$p_{2}$. There are some relationships among form factors,

$(m_{D^{\ast}}+m_{P})\,A_{1}^{D^{\ast}P}(q^{2})+(m_{D^{\ast}}-m_{P})\,A_{2}^{D^{\ast}P}(q^{2})\,=\,2\,m_{D^{\ast}}\,A_{3}^{D^{\ast}P}(q^{2}),$
(17) $A_{0}^{D^{\ast}P}(0)\,=\,A_{3}^{D^{\ast}P}(0),$ (18)
$A_{1}^{D^{\ast}V}(0)\,=\,A_{2}^{D^{\ast}V}(0),$ (19)
$V_{3}^{D^{\ast}V}(0)\,=\,V_{4}^{D^{\ast}V}(0).$ (20)

The numerical values of form factors have been comprehensively calculated with
the light front quark model in Ref. jhep.1912.102 , and are listed in Table 1.

Table 1: Numerical values of decay constants pdg2020 ; Eur.Phys.J.C.81.1110  and form factors at the pole $q^{2}$ $=$ $0$ jhep.1912.102 . $f_{\pi}$ $=$ $130.2{\pm}1.2$ MeV |  $f_{K}$ $=$ $155.7{\pm}0.3$ MeV |  $f_{\rho}$ $=$ $207.7{\pm}1.6$ MeV |  $f_{K^{\ast}}$ $=$ $202.5^{+6.5}_{-6.7}$ MeV
---|---|---|---
$V_{1}^{D^{\ast}{\rho}}$ $=$ $0.65$ |  $V_{2}^{D^{\ast}{\rho}}$ $=$ $0.51$ |  $V_{3}^{D^{\ast}{\rho}}$ $=$ $0.29$ |  $V_{5}^{D^{\ast}{\rho}}$ $=$ $1.42$
$V_{6}^{D^{\ast}{\rho}}$ $=$ $0.68$ |  $A_{1}^{D^{\ast}{\rho}}$ $=$ $0.59$ |  $A_{3}^{D^{\ast}{\rho}}$ $=$ $0.22$ |  $A_{4}^{D^{\ast}{\rho}}$ $=$ $0.24$
$V_{1}^{D^{\ast}K^{\ast}}$ $=$ $0.74$ |  $V_{2}^{D^{\ast}K^{\ast}}$ $=$ $0.43$ |  $V_{3}^{D^{\ast}K^{\ast}}$ $=$ $0.26$ |  $V_{5}^{D^{\ast}K^{\ast}}$ $=$ $1.50$
$V_{6}^{D^{\ast}K^{\ast}}$ $=$ $0.78$ |  $A_{1}^{D^{\ast}K^{\ast}}$ $=$ $0.69$ |  $A_{3}^{D^{\ast}K^{\ast}}$ $=$ $0.15$ |  $A_{4}^{D^{\ast}K^{\ast}}$ $=$ $0.22$
$V^{D^{\ast}{\pi}}$ $=$ $0.92$ |  $A_{0}^{D^{\ast}{\pi}}$ $=$ $0.68$ |  $A_{1}^{D^{\ast}{\pi}}$ $=$ $0.74$ |  $A_{2}^{D^{\ast}{\pi}}$ $=$ $0.61$
$V^{D^{\ast}K}$ $=$ $1.04$ |  $A_{0}^{D^{\ast}K}$ $=$ $0.78$ |  $A_{1}^{D^{\ast}K}$ $=$ $0.85$ |  $A_{2}^{D^{\ast}K}$ $=$ $0.68$

##  III Branching ratio

The branching ratio of two-body nonleptonic $D^{\ast}$ decays is defined as

${\cal
B}r\,=\,\frac{G_{F}^{2}}{48\,{\pi}}\,{|}V_{cs}{|}^{2}\,{|}V_{ud}{|}^{2}\,\frac{p_{\rm
cm}}{m_{D^{\ast}}^{2}\,{\Gamma}_{D^{\ast}}}\,\sum\limits_{s_{i}}{\cal M}{\cal
M}^{\dagger}.$ (21)

where $p_{\rm cm}$ is the center-of-mass momentum of final states in the rest
frame of the $D^{\ast}$ meson, and ${\Gamma}_{D^{\ast}}$ is the decay width of
the $D^{\ast}$ meson. ${\cal M}$ denotes the decay amplitude. For $D^{\ast}$
${\to}$ $\overline{K}{\pi}^{+}$ decays, there is only the $p$-wave amplitude.
For $D^{\ast}$ ${\to}$ $\overline{K}^{\ast}{\pi}^{+}$ and
$\overline{K}{\rho}^{+}$ decays, there are $s$-, $p$\- and $d$-wave
amplitudes. The analytic expressions of decay amplitudes with NF approach are
collected in Appendixes A and B.

The decay width of the $D^{{\ast}+}$ meson has been experimentally determined
to be ${\Gamma}_{D^{{\ast}+}}$ $=$ $83.4{\pm}1.8$ keV pdg2020 , while there
only was an upper limit for the decay width of the $D^{{\ast}0}$ meson,
${\Gamma}_{D^{{\ast}0}}$ $<$ $2.1$ MeV pdg2020 . It is widely assumed that
there should be a relation between decay widths for the $p$-wave strong decay
$D^{\ast}$ ${\to}$ $D{\pi}^{0}$ PhysRevD.32.189 ; PhysRevD.32.810 ;
PhysRevD.88.034034 ,

$\frac{{\Gamma}(D^{{\ast}0}{\to}D^{0}{\pi}^{0})}{{\Gamma}(D^{{\ast}+}{\to}D^{+}{\pi}^{0})}\,=\,\frac{p_{{\rm
cm}D^{0}{\pi}^{0}}^{3}}{p_{{\rm
cm}D^{+}{\pi}^{0}}^{3}}\,\frac{m_{D^{{\ast}+}}^{2}}{m_{D^{{\ast}0}}^{2}}.$
(22)

So it is expected that the decay width of $D^{{\ast}0}$,

${\Gamma}_{D^{{\ast}0}}\,=\,{\Gamma}_{D^{{\ast}+}}\frac{{\cal
B}r(D^{{\ast}+}{\to}D^{+}{\pi}^{0})}{{\cal
B}r(D^{{\ast}0}{\to}D^{0}{\pi}^{0})}\,\frac{p_{{\rm
cm}D^{0}{\pi}^{0}}^{3}}{p_{{\rm
cm}D^{+}{\pi}^{0}}^{3}}\,\frac{m_{D^{{\ast}+}}^{2}}{m_{D^{{\ast}0}}^{2}}\,=\,55.9^{+5.9}_{-5.4}\,\text{keV},$
(23)

where the comprehensive uncertainties are conservative estimates, and come
from data error of ${\Gamma}_{D^{{\ast}+}}$, branching ratios and mesonic
mass. The full decay width ${\Gamma}_{D^{{\ast}0}}$ in Eq.(23) is well
consistent with theoretical prediction of Ref. PhysRevD.88.034034 , and will
be used in our calculation to estimate branching ratios for $D^{{\ast}0}$
decays.

Our numerical calculation shows that branching ratios for nonleptonic
$D^{\ast}$ decays are,

$\displaystyle{\cal B}r(D^{{\ast}+}{\to}\overline{K}^{0}{\pi}^{+})$
$\displaystyle{\approx}$ $\displaystyle 1.6{\times}10^{-10},$ (24)
$\displaystyle{\cal B}r(D^{{\ast}+}{\to}\overline{K}^{{\ast}0}{\pi}^{+})$
$\displaystyle{\approx}$ $\displaystyle 4.4{\times}10^{-10},$ (25)
$\displaystyle{\cal B}r(D^{{\ast}+}{\to}\overline{K}^{0}{\rho}^{+})$
$\displaystyle{\approx}$ $\displaystyle 8.3{\times}10^{-10},$ (26)
$\displaystyle{\cal B}r(D^{{\ast}0}{\to}K^{-}{\pi}^{+})$
$\displaystyle{\approx}$ $\displaystyle 7.3{\times}10^{-10},$ (27)
$\displaystyle{\cal B}r(D^{{\ast}0}{\to}K^{{\ast}-}{\pi}^{+})$
$\displaystyle{\approx}$ $\displaystyle 2.0{\times}10^{-9},$ (28)
$\displaystyle{\cal B}r(D^{{\ast}0}{\to}K^{-}{\rho}^{+})$
$\displaystyle{\approx}$ $\displaystyle 2.9{\times}10^{-9}.$ (29)

Our comments on branching ratios are as follows.

(1) The light valence quarks of the $D^{{\ast}0}$ and $D^{{\ast}+}$ meson are
usually call spectator quarks, and they do not take part in the weak
interaction directly, when $W$ annihilation and $W$ exchange contributions are
not considered with the NF approach. The partial width for $D^{{\ast}0}$ and
$D^{{\ast}+}$ meson decays into the similar isospin final states should in
principle be approximately equal. In fact, $D^{{\ast}+}$ meson decays are
dynamically induced by both external and internal $W$ emission interactions,
and the interference between these two contributions is destructive owing to a
combination of Wilson coefficients. This is the main reason for he
hierarchical relationship, ${\cal
B}r(D^{{\ast}+}{\to}\overline{K}^{0}{\pi}^{+})$ $<$ ${\cal
B}r(D^{{\ast}0}{\to}K^{-}{\pi}^{+})$, ${\cal
B}r(D^{{\ast}+}{\to}\overline{K}^{{\ast}0}{\pi}^{+})$ $<$ ${\cal
B}r(D^{{\ast}0}{\to}K^{{\ast}-}{\pi}^{+})$, and ${\cal
B}r(D^{{\ast}+}{\to}\overline{K}^{0}{\rho}^{+})$ $<$ ${\cal
B}r(D^{{\ast}0}{\to}K^{-}{\rho}^{+})$. Similar hierarchical phenomena have
also been observed experimentally in the Cabibbo-favored pseudoscalar $D^{0}$
and $D^{+}$ meson decays, for example, ${\cal
B}r(D^{+}{\to}K^{0}_{S}{\pi}^{+})$ $=$ $1.562(31)\%$ and ${\cal
B}r(D^{0}{\to}K^{-}{\pi}^{+})$ $=$ $3.950(31)\%$ pdg2020 .

(2) There are three partial wave amplitudes for $D^{\ast}$ ${\to}$ $PV$
decays, while only the $p$-wave amplitude contributes to $D^{\ast}$ ${\to}$
$PP$ decays. Hence, there is a hierarchical relationship among branching
ratios, ${\cal B}r(D^{{\ast}+}{\to}PP)$ $<$ ${\cal B}r(D^{{\ast}+}{\to}PV)$
and ${\cal B}r(D^{{\ast}0}{\to}PP)$ $<$ ${\cal B}r(D^{{\ast}0}{\to}PV)$. In
addition, because of decay constants $f_{\rho}$ $>$ $f_{\pi}$, branching ratio
for emission ${\rho}^{+}$ is generally larger than that for emission
${\pi}^{+}$, for either $D^{{\ast}0}$ or $D^{{\ast}+}$ decays. Similar
phenomena have also been seen in pseudoscalar $D$ meson decays, for example,
${\cal B}r(D^{0}{\to}K^{-}{\rho}^{+})$ $=$ $(11.3{\pm}0.77)\%$ and ${\cal
B}r(D^{0}{\to}K^{{\ast}-}{\pi}^{+})$ $=$ $(2.31^{+0.40}_{-0.20})\%$ determined
from $D^{0}$ ${\to}$ $K^{-}{\pi}^{+}{\pi}^{0}$ decays pdg2020 .

(3) As it is well known the vector ${\rho}$ and $K^{\ast}$ mesons are
resonances, and they will decay promptly into pseudoscalar mesons via the
strong interactions, with branching ratios ${\cal B}r({\rho}{\to}{\pi}{\pi})$
${\sim}$ $100\,\%$ and ${\cal B}r(K^{\ast}{\to}K{\pi})$ ${\sim}$ $100\,\%$
pdg2020 . The vector ${\rho}$ ($K^{\ast}$) meson in $D^{{\ast}0}$ ${\to}$
$K^{-}{\rho}^{+}$ ($K^{{\ast}-}{\pi}^{+}$) decay should be reconstructed by
the final pseudoscalar mesons. Besides, both $D^{{\ast}0}$ ${\to}$
$K^{-}{\rho}^{+}$ and $D^{{\ast}0}$ ${\to}$ $K^{{\ast}-}{\pi}^{+}$ decays
contribute to $D^{{\ast}0}$ ${\to}$ $K^{-}{\pi}^{+}{\pi}^{0}$. The branching
ratio of the three-body $D^{{\ast}0}$ ${\to}$ $K^{-}{\pi}^{+}{\pi}^{0}$ decay
can be approximately the sum of the branching ratios of these two two-body
$D^{{\ast}0}$ ${\to}$ $PV$ modes whose interference effect happens in a small
region of the Dalitz plot. A similar case can be seen by the comparison
between the branching ratio ${\cal B}r(D^{0}{\to}K^{-}{\pi}^{+}{\pi}^{0})$ $=$
$(14.4{\pm}0.5)\,\%$ and the sum of partial branching ratios ${\cal
B}r(D^{0}{\to}K^{-}{\rho}^{+})$ $=$ $(11.3{\pm}0.7)\,\%$ and ${\cal
B}r(D^{0}{\to}K^{{\ast}-}{\pi}^{+})$ $=$ $(2.31^{+0.40}_{-0.20})\,\%$ pdg2020
. Therefore, it may be an educated guess that the three-body $D^{{\ast}0}$
${\to}$ $K^{-}{\pi}^{+}{\pi}^{0}$ decay will have a relatively larger
branching ratio, ${\cal B}r(D^{{\ast}0}{\to}K^{-}{\pi}^{+}{\pi}^{0})$ ${\sim}$
$5.0{\times}10^{-9}$, and can be more easily investigated in experiments,
compared to the two-body $D^{{\ast}0}$ ${\to}$ $K^{-}{\rho}^{+}$ and
$K^{{\ast}-}{\pi}^{+}$ decays.

(4) Theoretical predictions on branching ratios are easily influenced by a
number of factors, including final state interactions and other contributions.
For example, studies PhysRevD.67.014001 ; PhysRevD.81.074021  using the flavor
topology diagrammatic approach have shown that for Cabibbo-favored $D$ meson
decays, the $W$ exchange and annihilation contributions should deserve due
attention, although the external $W$ emission contributions always give the
largest contributions in the total amplitudes. We would like to point out that
what we want is whether it is feasible to investigate the nonleptonic
$D^{{\ast}}$ weak decays at the future experiments, so the magnitude order
estimation rather than precise calculation on branching ratio may be enough.
It is generally believed that the NF approach can give a reasonable and
correct magnitude order estimation on branching ratio for nonleptonic heavy
flavored meson decays arising from the external $W$ emission weak
interactions. In this sense, the magnitude order of branching ratios in
Eq.(24-29) seems to be reliable. The potential event numbers of the concerned
$D^{\ast}$ ${\to}$ $PP$, $PV$ decays are listed in Table 2. It is clear that
the Cabibbo-favored $D^{\ast}$ ${\to}$ $\overline{K}{\pi}^{+}$,
$\overline{K}^{\ast}{\pi}^{+}$, $\overline{K}{\rho}^{+}$ decays can be
measurable at the future STCF, CEPC, FCC-ee and LHCb@HL-LHC experiments. The
$D^{{\ast}+}$ ${\to}$ $K^{{\ast}-}{\pi}^{+}$ and $K^{-}{\rho}^{+}$ decays can
also be investigated at SuperKEKB experiments.

Table 2: The potential event numbers of $D^{\ast}$ ${\to}$ $\overline{K}{\pi}^{+}$, $\overline{K}^{\ast}{\pi}^{+}$, $\overline{K}{\rho}^{+}$ decays at experiments, where the $D^{\ast}$ meson data available has been estimated in Section I. experiment | SuperKEKB | STCF | CEPC | FCC-ee | LHCb@HL-LHC
---|---|---|---|---|---
$N_{D^{\ast}}$ | $2{\times}10^{10}$ | $8{\times}10^{10}$ | $10^{11}$ | $10^{12}$ | $2{\times}10^{14}$
$N_{D^{{\ast}+}{\to}\overline{K}^{0}{\pi}^{+}}$ | $3$ | $13$ | $16$ | $160$ | $3.2{\times}10^{4}$
$N_{D^{{\ast}+}{\to}\overline{K}^{{\ast}0}{\pi}^{+}}$ | $9$ | $35$ | $44$ | $440$ | $8.8{\times}10^{4}$
$N_{D^{{\ast}+}{\to}\overline{K}^{0}{\rho}^{+}}$ | $17$ | $66$ | $83$ | $830$ | $1.66{\times}10^{5}$
$N_{D^{{\ast}0}{\to}K^{-}{\pi}^{+}}$ | $14$ | $58$ | $72$ | $720$ | $1.44{\times}10^{5}$
$N_{D^{{\ast}0}{\to}K^{{\ast}-}{\pi}^{+}}$ | $40$ | $160$ | $200$ | $2000$ | $4.0{\times}10^{5}$
$N_{D^{{\ast}0}{\to}K^{-}{\rho}^{+}}$ | $57$ | $230$ | $287$ | $2870$ | $5.74{\times}10^{5}$
$N_{D^{{\ast}0}{\to}K^{-}{\pi}^{+}{\pi}^{0}}$ | $100$ | $400$ | $500$ | $5000$ | $1.0{\times}10^{6}$

##  IV Summary

Now, more than 45 years after the discovery of the $D^{\ast}$ mesons, our
knowledge and understanding of the nature of the $D^{\ast}$ mesons is far from
enough, and needs to be substantially improved. One of the major reasons that
excessively hindered experimental investigation on $D^{\ast}$ mesons is that
data are too scarce. We should thank the high-luminosity particle physics
experiments for offering us a huge amount of $D^{\ast}$ meson data and a
tempting opportunity to explore the wanted $D^{\ast}$ meson in the future.
Compared with the dominant $D^{\ast}$ ${\to}$ $D{\pi}$ decays which are
subject to kinematical factors, one advantage of nonleptonic $D^{\ast}$ weak
decays is that the final pion and kaon mesons are energetic and easily
detectable by the sensitive high-resolution detectors. In addition, study of
nonleptonic $D^{\ast}$ weak decays is scientifically significant, and provide
us with a new venue for testing SM. In this paper, the Cabibbo-favored two-
body nonleptonic $D^{\ast}$ ${\to}$ $PP$, $PV$ decays were studied by using
the NF approach within SM. It is found that branching ratios for $D^{\ast}$
${\to}$ $\overline{K}{\pi}^{+}$, $\overline{K}^{\ast}{\pi}^{+}$,
$\overline{K}{\rho}^{+}$ decays can reach up to ${\cal O}(10^{-10})$ or more,
and can be accessible at STCF, CEPC, FCC-ee and LHCb@HL-LHC experiments, which
indicate that study of these weak interaction processes is experimentally
feasible and practicable in the future.

## Acknowledgments

The work is supported by the National Natural Science Foundation of China
(Grant Nos. 11705047, U1632109, 11875122) and Natural Science Foundation of
Henan Province (Grant No. 222300420479), the Excellent Youth Foundation of
Henan Province (Grant No. 212300410010). We would like to acknowledge the
useful help and valuable discussion from Prof. Haibo Li (IHEP@CAS), Prof.
Shuangshi Fang (IHEP@CAS), Prof. Frank Porter (Caltech), Prof. Antimo Palano
(INFN), Prof. Chengping Shen (Fudan University), Dr. Ping Xiao (Fudan
University), Dr. Qingping Ji (Henan Normal University), Dr. Huijing Li (Henan
Normal University) and Ms. Liting Wang (Henan Normal University), and positive
comments and constructive suggestions from referees.

##  Appendix A amplitudes for $V$ ${\to}$ $P_{1}$ $+$ $P_{2}$ decays

With the conventions of Eq.(9), Eq.(10), Eq.(13) and Eq.(14), the general
expression of HMEs with NF approach for $V$ ${\to}$ $P_{1}$ $+$ $P_{2}$
transition can be written as,

${\cal
M}\,=\,{\langle}P_{2}{|}V_{\mu}-A_{\mu}{|}0{\rangle}\,{\langle}P_{1}{|}V^{\mu}-A^{\mu}{|}V{\rangle}\,=\,{\cal
M}_{p}\,({\epsilon}_{V}{\cdot}p_{P_{2}}),$ (30) ${\cal
M}_{p}\,=\,2\,f_{P_{2}}\,m_{V}\,A_{0}^{VP_{1}}(0),$ (31)
$\sum\limits_{s_{V}}{\cal M}\,{\cal M}^{\dagger}\,=\,{|}{\cal
M}_{p}{|}^{2}\,p_{\rm cm}^{2},$ (32) $p_{\rm
cm}\,=\,{\lambda}^{1/2}(m_{V}^{2},m_{P_{1}}^{2},m_{P_{2}}^{2})/2\,m_{V},$ (33)
${\lambda}(a,b,c)\,=\,a^{2}+b^{2}+c^{2}-2\,a\,b-2\,b\,c-2\,c\,a.$ (34)

With the NF approach, there is only one amplitude for the neutral
$D^{{\ast}0}$ meson decay in question, which corresponds to external $W$
emission. There are two amplitudes for the charged $D^{{\ast}{\pm}}$ meson
decay in question, which correspond to external and internal $W$ emissions.
The partial wave amplitudes for $D^{{\ast}0}$ ${\to}$ $K^{-}{\pi}^{+}$ decay
are written as,

${\cal M}_{p}\,=\,2\,m_{D^{{\ast}}}\,a_{1}\,f_{\pi}\,A_{0}^{D^{{\ast}}K},$
(35)

and for $D^{{\ast}+}$ ${\to}$ $\overline{K}^{0}{\pi}^{+}$ decay,

${\cal
M}_{p}\,=\,2\,m_{D^{{\ast}}}\,\big{\\{}a_{1}\,f_{\pi}\,A_{0}^{D^{{\ast}}K}+a_{2}\,f_{K}\,A_{0}^{D^{{\ast}}{\pi}}\big{\\}},$
(36)

where the coefficients $a_{1}$ and $a_{2}$ correspond to external and internal
$W$ emission, respectively; and they are defined as,

$a_{1}\,=\,C_{1}+C_{2}/N_{c},$ (37) $a_{2}\,=\,C_{2}+C_{1}/N_{c}.$ (38)

In practice, it is generally believed that coefficients $a_{1,2}$ are also
influenced by nonfactorizable contributions and final state interactions. In
many phenomenological studies on charmed meson weak decays, such as Refs.
zpc.34.103 ; Stech.1985 ; npb.133.315 ; cpc.27.665 ; epjc.42.391 ;
PhysRevD.84.074019 ; PhysRevD.81.074021 ; PhysRevD.93.114010 ;
PhysRevD.86.036012 , coefficients $a_{1}$ ${\approx}$ $1.2$ and $a_{2}$
${\approx}$ $-0.5$ are often used for charmed meson decays by including
comprehensive contributions.

##  Appendix B amplitudes for $V_{1}$ ${\to}$ $V_{2}$ $+$ $P$ decays

With the conventions of Eqs.(9-16), the general expression of HMEs with NF
approach for $V_{1}$ ${\to}$ $V_{2}$ $+$ $P$ transition can be written as,

$\displaystyle{\cal M}$ $\displaystyle=$
$\displaystyle{\langle}V_{2}{|}V_{\mu}-A_{\mu}{|}0{\rangle}\,{\langle}P{|}V_{\mu}-A_{\mu}{|}V_{1}{\rangle}$
(39) $\displaystyle=$ $\displaystyle{\cal
M}_{s}\,({\epsilon}_{V_{1}}{\cdot}{\epsilon}_{V_{2}}^{\ast})+\frac{{\cal
M}_{d}}{m_{V_{1}}\,m_{V_{2}}}\,({\epsilon}_{V_{1}}{\cdot}p_{V_{2}})\,({\epsilon}_{V_{2}}^{\ast}{\cdot}p_{V_{1}})$
$\displaystyle+$ $\displaystyle\frac{{\cal
M}_{p}}{m_{V_{1}}\,m_{V_{2}}}\,{\varepsilon}_{{\mu}{\nu}{\alpha}{\beta}}\,{\epsilon}_{V_{1}}^{\mu}\,{\epsilon}_{V_{2}}^{{\ast}{\nu}}\,p_{V_{1}}^{\alpha}\,p_{V_{2}}^{\beta},$
$\displaystyle{\cal M}^{\prime}$ $\displaystyle=$
$\displaystyle{\langle}P{|}V_{\mu}-A_{\mu}{|}0{\rangle}\,{\langle}V_{2}{|}V_{\mu}-A_{\mu}{|}V_{1}{\rangle}$
(40) $\displaystyle=$ $\displaystyle{\cal
M}_{s}^{\prime}\,({\epsilon}_{V_{1}}{\cdot}{\epsilon}_{V_{2}}^{\ast})+\frac{{\cal
M}_{d}^{\prime}}{m_{V_{1}}\,m_{V_{2}}}\,({\epsilon}_{V_{1}}{\cdot}p_{V_{2}})\,({\epsilon}_{V_{2}}^{\ast}{\cdot}p_{V_{1}})$
$\displaystyle+$ $\displaystyle\frac{{\cal
M}_{p}^{\prime}}{m_{V_{1}}\,m_{V_{2}}}\,{\varepsilon}_{{\mu}{\nu}{\alpha}{\beta}}\,{\epsilon}_{V_{1}}^{\mu}\,{\epsilon}_{V_{2}}^{{\ast}{\nu}}\,p_{V_{1}}^{\alpha}\,p_{V_{2}}^{\beta},$
$\displaystyle{\cal M}_{s}$ $\displaystyle=$
$\displaystyle-i\,f_{V_{2}}\,m_{V_{2}}\,(m_{V_{1}}+m_{P})\ A_{1}^{V_{1}P}(0),$
(41) $\displaystyle{\cal M}_{d}$ $\displaystyle=$
$\displaystyle-i\,f_{V_{2}}\,m_{V_{1}}\,m_{V_{2}}\,\frac{2\,m_{V_{2}}}{m_{V_{1}}+m_{P}}\,A_{2}^{V_{1}P}(0),$
(42) $\displaystyle{\cal M}_{p}$ $\displaystyle=$ $\displaystyle-
f_{V_{2}}\,m_{V_{1}}\,m_{V_{2}}\,\frac{2\,m_{V_{2}}}{m_{V_{1}}+m_{P}}\,V^{V_{1}P}(0),$
(43) $\displaystyle{\cal M}_{s}^{\prime}$ $\displaystyle=$
$\displaystyle-i\,f_{P}\,(m_{V_{1}}^{2}-m_{V_{2}}^{2})\,V_{1}^{V_{1}V_{2}}(0),$
(44) $\displaystyle{\cal M}_{d}^{\prime}$ $\displaystyle=$
$\displaystyle-i\,f_{P}\,m_{V_{1}}\,m_{V_{2}}\,\big{\\{}V_{4}^{V_{1}V_{2}}(0)-V_{5}^{V_{1}V_{2}}(0)+V_{6}^{V_{1}V_{2}}(0)\big{\\}},$
(45) $\displaystyle{\cal M}_{p}^{\prime}$ $\displaystyle=$
$\displaystyle-2\,f_{P}\,m_{V_{1}}\,m_{V_{2}}\,A_{1}^{V_{1}V_{2}}(0),$ (46)
$\displaystyle\sum\limits_{s_{V_{1}}}\sum\limits_{s_{V_{2}}}{\cal M}\,{\cal
M}^{\dagger}$ $\displaystyle=$ $\displaystyle{|}{\cal
M}_{s}{|}^{2}\,(x^{2}+2)+2\,{\cal R}\big{(}{\cal M}_{s}\,{\cal
M}_{d}^{\ast}\big{)}\,x\,(x^{2}-1)$ (47) $\displaystyle+$
$\displaystyle{|}{\cal M}_{d}{|}^{2}\,(x^{2}-1)^{2}+2\,{|}{\cal
M}_{p}{|}^{2}\,(x^{2}-1),$
$\displaystyle\sum\limits_{s_{V_{1}}}\sum\limits_{s_{V_{2}}}{\cal
M}^{\prime}\,{\cal M}^{{\prime}\dagger}$ $\displaystyle=$
$\displaystyle{|}{\cal M}_{s}^{\prime}{|}^{2}\,(x^{2}+2)+2\,{\cal
R}\big{(}{\cal M}_{s}^{\prime}\,{\cal
M}_{d}^{{\prime}{\ast}}\big{)}\,x\,(x^{2}-1)$ (48) $\displaystyle+$
$\displaystyle{|}{\cal M}_{d}^{\prime}{|}^{2}\,(x^{2}-1)^{2}+2\,{|}{\cal
M}_{p}^{\prime}{|}^{2}\,(x^{2}-1),$
$x\,=\,\frac{p_{V_{1}}{\cdot}p_{V_{2}}}{m_{V_{1}}\,m_{V_{2}}}\,=\,\frac{m_{V_{1}}^{2}+m_{V_{2}}^{2}-m_{P}^{2}}{2\,m_{V_{1}}\,m_{V_{2}}}.$
(49)

For $D^{{\ast}0}$ ${\to}$ $K^{{\ast}-}{\pi}^{+}$ decay, the partial wave
amplitudes are written as,

$\displaystyle{\cal M}_{s}^{\prime}$ $\displaystyle=$
$\displaystyle-i\,f_{\pi}\,a_{1}\,(m_{D^{\ast}}^{2}-m_{K^{\ast}}^{2})\,V_{1}^{D^{\ast}K^{\ast}},$
(50) $\displaystyle{\cal M}_{d}^{\prime}$ $\displaystyle=$
$\displaystyle-i\,f_{\pi}\,a_{1}\,m_{D^{\ast}}\,m_{K^{\ast}}\,\big{\\{}V_{4}^{D^{\ast}K^{\ast}}-V_{5}^{D^{\ast}K^{\ast}}+V_{6}^{D^{\ast}K^{\ast}}\big{\\}},$
(51) $\displaystyle{\cal M}_{p}^{\prime}$ $\displaystyle=$
$\displaystyle-2\,f_{\pi}\,a_{1}\,m_{D^{\ast}}\,m_{K^{\ast}}\,A_{1}^{D^{\ast}K^{\ast}}.$
(52)

For $D^{{\ast}0}$ ${\to}$ $K^{-}{\rho}^{+}$ decay, the partial wave amplitudes
are written as,

$\displaystyle{\cal M}_{s}$ $\displaystyle=$
$\displaystyle-i\,f_{\rho}\,a_{1}\,m_{\rho}\,(m_{D^{\ast}}+m_{K})\
A_{1}^{D^{\ast}K},$ (53) $\displaystyle{\cal M}_{d}$ $\displaystyle=$
$\displaystyle-i\,f_{\rho}\,a_{1}\,m_{D^{\ast}}\,m_{\rho}\,\frac{2\,m_{\rho}\,}{m_{D^{\ast}}+m_{K}}\,A_{2}^{D^{\ast}K},$
(54) $\displaystyle{\cal M}_{p}$ $\displaystyle=$ $\displaystyle-
f_{\rho}\,a_{1}\,m_{D^{\ast}}\,m_{\rho}\,\frac{2\,m_{\rho}\,}{m_{D^{\ast}}+m_{K}}\,V^{D^{\ast}K}.$
(55)

For $D^{{\ast}+}$ ${\to}$ $\overline{K}^{{\ast}0}{\pi}^{+}$ decay, each of
partial wave amplitudes can be divided into two parts,

${\cal M}_{i}\,=\,{\cal M}_{i}^{(1)}+{\cal M}_{i}^{(2)},\quad\text{for}\
i\,=\,s,p,d.$ (56) $\displaystyle{\cal M}_{s}^{(2)}$ $\displaystyle=$
$\displaystyle-i\,f_{K^{\ast}}\,a_{2}\,m_{K^{\ast}}\,(m_{D^{\ast}}+m_{\pi})\
A_{1}^{D^{\ast}{\pi}},$ (57) $\displaystyle{\cal M}_{d}^{(2)}$
$\displaystyle=$
$\displaystyle-i\,f_{K^{\ast}}\,a_{2}\,m_{D^{\ast}}\,m_{K^{\ast}}\,\frac{2\,m_{K^{\ast}}\,}{m_{D^{\ast}}+m_{\pi}}\,A_{2}^{D^{\ast}{\pi}},$
(58) $\displaystyle{\cal M}_{p}^{(2)}$ $\displaystyle=$ $\displaystyle-
f_{K^{\ast}}\,a_{2}\,m_{D^{\ast}}\,m_{K^{\ast}}\,\frac{2\,m_{K^{\ast}}}{m_{D^{\ast}}+m_{\pi}}\,V^{D^{\ast}{\pi}},$
(59)

and expressions of ${\cal M}_{i}^{(1)}$ are the same as those of ${\cal
M}_{s,p,d}^{\prime}$ for $D^{{\ast}0}$ ${\to}$ $K^{{\ast}-}{\pi}^{+}$ decay in
Eq.(50), Eq.(51) and Eq.(52).

For $D^{{\ast}+}$ ${\to}$ $\overline{K}^{0}{\rho}^{+}$ decay, each of partial
wave amplitudes can also be divided into two parts similar to Eq.(56),

$\displaystyle{\cal M}_{s}^{(2)}$ $\displaystyle=$
$\displaystyle-i\,f_{K}\,a_{2}\,(m_{D^{\ast}}^{2}-m_{\rho}^{2})\,V_{1}^{D^{\ast}{\rho}},$
(60) $\displaystyle{\cal M}_{d}^{(2)}$ $\displaystyle=$
$\displaystyle-i\,f_{K}\,a_{2}\,m_{D^{\ast}}\,m_{\rho}\,\big{\\{}V_{4}^{D^{\ast}{\rho}}-V_{5}^{D^{\ast}{\rho}}+V_{6}^{D^{\ast}{\rho}}\big{\\}},$
(61) $\displaystyle{\cal M}_{p}^{(2)}$ $\displaystyle=$
$\displaystyle-2\,f_{K}\,a_{2}\,m_{D^{\ast}}\,m_{\rho}\,A_{1}^{D^{\ast}{\rho}},$
(62)

and expressions of ${\cal M}_{i}^{(1)}$ are the same as those of ${\cal
M}_{s,p,d}$ for $D^{{\ast}0}$ ${\to}$ $K^{-}{\rho}^{+}$ decay in Eq.(53),
Eq.(54) and Eq.(55).

## References

  * (1) J. Bjorken, S. Glashow, Phys. Lett. 11, 255 (1964).
  * (2) S. Glashow, J. Iliopoulos, L. Maiani, Phys. Rev. D 2, 1285 (1970).
  * (3) J. Aubert, U. Becker, P. Biggs et al., Phys. Rev. Lett. 33, 1404 (1974).
  * (4) J. Augustin, A. Boyarski, M. Breidenbach et al., Phys. Rev. Lett. 33, 1406 (1974).
  * (5) G. Goldhaber, F. Pierre, G. Abrams et al., Phys. Rev. Lett. 37, 255 (1976).
  * (6) I. Peruzzi, M. Piccolo, G. Feldman et al., Phys. Rev. Lett. 37, 569 (1976).
  * (7) R. Workman et al. (Particle Data Group), Prog. Theor. Exp. Phys. 2022, 083C01 (2022).
  * (8) K. Lane, E. Eichten, Phys. Rev. Lett. 37, 477 (1976).
  * (9) A. Rújula, H. Georgi, S. Glashow, Phys. Rev. Lett. 37, 398 (1976).
  * (10) F. Close, Phys. Lett. B 65, 55 (1976).
  * (11) S. Matsuda, Phys. Lett. B 66, 70 (1977).
  * (12) J. Korner, M. Kuroda, G. Schierholz, Phys. Lett. B 70, 106 (1977).
  * (13) M. Suzuki, W. Wada, Phys. Rev. D 15, 759 (1977).
  * (14) B. Humpert, R. Clark Phys. Rev. D 16, 1327 (1977).
  * (15) Y. Iwasaki, Phys. Rev. D 17, 765 (1978).
  * (16) D. Cronin-Hennessy et al (CLEO Collaboration), Phys. Rev. D 80, 072001 (2009).
  * (17) X. Dong, L. Wang, C. Yuan, Chin. Phys. C 42, 043002 (2018).
  * (18) I. Peruzzi, M. Piccolo, G. Feldman et al, Phys. Rev. Lett. 39, 1301 (1977).
  * (19) G. Goldhaber, J. Wiss, G. Abrams et al, Phys. Lett. B 69, 503 (1977).
  * (20) S. Okubo, Phys. Lett. 5, 165 (1963).
  * (21) G. Zweig, CERN-TH-401 G. Zweig, CERN-TH-412 (1964).
  * (22) J. Iizuka, Prog. Theor. Phys. Suppl. 37-38, 21 (1966).
  * (23) N. Cabibbo, Phys. Rev. Lett. 10, 531 (1963).
  * (24) M. Kobayashi, T. Maskawa, Prog. Theor. Phys. 49, 652 (1973).
  * (25) M. Ablikim et al. (BESIII collaboration), JHEP 05, 155 (2022).
  * (26) E. Kou et al., Prog. Theor. Exp. Phys. 123C01 (2019); 029201 (2020)(E).
  * (27) M. Lisovyi, A. Verbytskyi, O. Zenaiev, Eur. Phys. J. C 76, 397 (2016).
  * (28) X. Lyu (STCF Working group), PoS(BEAUTY2020), 060 (2021).
  * (29) V. Anashin et al., https://ctd.inp.nsk.su/wiki/images/4/47/CDR2_ScTau_en_vol1.pdf.
  * (30) J. Costa et al., IHEP-CEPC-DR-2018-02, arXiv:1811.10545.
  * (31) A. Abada et al., Eur. Phys. J. C 79, 474 (2019).
  * (32) R. Aaij et al. (LHCb Collaboration), JHEP 03, 159 (2016); 09, 013 (2016)(E); 05, 074 (2017)(E).
  * (33) B. Abelev et al. (ALICE Collaboration), JHEP 07, 191 (2012).
  * (34) I. Bediaga et al. (LHCb Collaboration), arXiv:1808.08865.
  * (35) Y. Yang, Z. Li, K. Li, J. Huang, J. Sun, Eur. Phys. J. C 81, 1110 (2021).
  * (36) G. Buchalla, A. Buras, M. Lautenbacher, Rev. Mod. Phys. 68, 1125, (1996).
  * (37) M. Bauer, B. Stech, M. Wirbel, Z. Phys. C 34, 103 (1987).
  * (38) J. Bjorken, Nucl. Phys. B Proc. Suppl. 11, 325 (1989).
  * (39) H. Cheng, C. Chiang, Phys. Rev. D 81, 074021 (2010).
  * (40) M. Wirbel, B. Stech, M. Bauer, Z. Phys. C 29, 637 (1985).
  * (41) J. Sun, L. Chen, Q. Chang et al., Int. J. Mod. Phys. A 30, 1550094 (2015).
  * (42) Q. Chang, L. Wang, X. Li, JHEP 12, 102 (2019).
  * (43) S. Godfrey, N. Isgur, Phys. Rev. D 32, 189 (1985).
  * (44) R. Thews, A. Kamal, Phys. Rev. D 32, 810 (1985).
  * (45) J. Rosner, Phys. Rev. D 88, 034034 (2013).
  * (46) C. Chiang, Z. Luo, L. Rosner, Phys. Rev. D 67, 014001 (2003).
  * (47) B. Stech, 5th Moriond Workshop: Heavy Quarks, 151, (1985).
  * (48) D. Fakirov, B. Stech, Nucl. Phys. B 133, 315 (1978).
  * (49) J Li, M. Yang, D. Du, Chin. Phys. C 27, 665 (2003).
  * (50) Y. Wu, M. Zhong, Y. Zhou, Eur. Phys. J. C 42, 391 (2005).
  * (51) F. Yu, X. Wang, C. Lu, Phys. Rev. D 84, 074019 (2011).
  * (52) H. Cheng, C. Chiang, A. Kuo, Phys. Rev. D 93, 114010 (2016).
  * (53) H. Li, C. Lu, F. Yu, Phys. Rev. D 86, 036012 (2012).


# Quantum vs. Symplectic Computers

Igor Volovich

(June 2024)

###### Abstract

In this paper, we propose the concept of symplectic computers, which have the
potential to be more powerful than quantum computers. Unlike quantum
computing, which consists of a sequence of unitary transformations (gates) and
projectors (measurements), symplectic computation involves a sequence of
symplectic transformations and measurements.

The proposal to explore symplectic computers is based on the following
quantum-symplectic duality. The Schrödinger equation in its standard complex
form describes the unitary evolution of a quantum system, while its real form
describes the symplectic evolution of a classical mechanical system. This
quantum-symplectic duality can be leveraged to enhance the capabilities of
quantum and symplectic computers. In this symplectic approach, the role of a
quantum bit (qubit) is taken by a symplectic bit (symbit).

##  1 Introduction

Recent years have seen significant progress in the field of quantum computer
technology, see, for example, [1, 2].

In this paper, we propose the concept of symplectic computers, which have the
potential to be more powerful than quantum computers. Unlike quantum
computing, which consists of a sequence of unitary transformations (gates) and
projectors (measurements), symplectic computation involves a sequence of
symplectic transformations and measurements. The main mathematical idea of the
symplectic computers is based on the quantum-symplectic duality. The
Schrödinger equation in its standard complex form describes the unitary
evolution of a quantum system, while its real form describes the symplectic
evolution of a classical mechanical system. This quantum-symplectic duality
can be leveraged to enhance the capabilities of symplectic computers. In this
symplectic approach, the role of a quantum bit (qubit) is taken by a
symplectic bit (symbit).

Above mentioned quantum-symplectic duality can be explained as follows. We
consider the final dimensional Hilbert space ${\mathbb{C}}^{N}$. The
Schrodinger equation in its standard complex form reads

$i\dot{\psi}=H\,\psi.$ (1.1)

Here $\psi$ is a complex valued wave function, the dot means the time
derivative and $H$ is an hermitian matrix. This matrix $H$ can be represented
in the form

$H=K+iL,$ (1.2)

where $K$ is a symmetric real matrix and $L$ is a skew real matrix. Now we
rewrite equation (1.1) in the real form by using the representation of the
wave function

$\psi=q+ip,$ (1.3)

where $p$ and $q$ are the real and imaginary parts of $\psi$. Then the
Schrodinger equation (1.1) takes the form

$\displaystyle\dot{q}$ $\displaystyle=$ $\displaystyle\,\,\,Kp+Lq,$ (1.4)
$\displaystyle\dot{p}$ $\displaystyle=$ $\displaystyle-Kq+Lp.$ (1.5)

The Schrodinger equation in the real form (1.4), (1.5) describes a symplectic
evolution of the corresponding mechanical system [3, 4].

This fact can be called the quantum-symplectic duality. The idea of such
duality comes back to the wave-particle duality suggested by Einstein [5] and
de Broglie [6]. From the modern point of view wave-particle duality is the
concept in quantum mechanics that quantum entities exhibit particle or wave
properties according to the experimental circumstances [7, 8].

The unitary group is a subgroup of a symplectic group and we have the relation

$U(N)=Sp(2N,{\mathbb{R}})\cap{\mathbb{O}}(2N,{\mathbb{R}}).$ (1.6)

Therefore, any quantum computation is a symplectic computation. However, not
every symplectic computation is quantum. For this reason, symplectic computers
are, in principle, generally more powerful than quantum ones.

The paper is organized as follows. In Section 2 we discuss the quantum-
symplectic duality. In Subsection 2.1 the Schrodinger equation is presented as
symplectic mechanics. In Subsection 2.2 unitary and symplectic evolutions are
discussed. In Subsection 2.3 qubits and symbits are compared. In Section 3
unitary and symplectic gates are considered. Possible generalizations of the
proposed approach are mentioned in Section 4.

##  2 Quantum-symplectic duality

###  2.1 Schrodinger equation and symplectic mechanics

The Schrodinger equation actually describes the classical symplectic
Hamiltonian system. We consider the finite dimensional Hilbert space
$\mathcal{H}=\mathbb{C^{N}}$. The Schrodinger equation reads

$i\dot{\psi_{a}}=H_{ab}\psi_{b}$ (2.1)

where $a,b=1...N$, $\psi_{a}\in\mathbb{C}$, $H_{ab}=\overline{H}_{ba}$ is an
Hermitian matrix. We assume the summation over the repeating indices.

We rewrite this equation by using real quantities and set

$\psi_{a}=q_{a}+ip_{a},$ (2.2)

where $q_{a},p_{a}$ are real. We also write

$H_{ab}=K_{ab}+iL_{ab},$ (2.3)

where $K_{ab}=K_{ba}$, $L_{ab}=-L_{ba}$ are matrices with real entries. Then
the Schrodinger equation (2.1) takes the form

$\dot{q_{a}}=K_{ab}p_{b}+L_{ab}q_{b}$ (2.4)
$\dot{p_{a}}=-K_{ab}q_{b}+L_{ab}p_{b}$ (2.5)

These equations are classical symplectic Hamiltonian equations

$\displaystyle\dot{q}_{a}=\frac{\partial H_{sym}}{\partial
p_{a}},\qquad\dot{p}_{a}=-\frac{\partial H_{sym}}{\partial q_{a}},$ (2.6)

with the Hamiltonian

$H_{sym}=\frac{1}{2}(p_{a}K_{ab}\,p_{b}+q_{a}K_{ab}\,q_{b})+p_{a}L_{ab}\,q_{b}$
(2.7)

Remark 1. If the Hamiltonian equation in the Schrodinger equations depends on
time

$i\dot{\psi_{a}}=H_{ab}(t)\psi_{b}$ (2.8)

with

$H_{ab}(t)=K_{ab}(t)+iL_{ab}(t),$ (2.9)

then we get the following equation

$\dot{q_{a}}=K_{ab}(t)p_{b}+L_{ab}(t)q_{b}$ (2.10)
$\dot{p_{a}}=-K_{ab}(t)q_{b}+L_{ab}(t)p_{b}$ (2.11)
$H_{sym}(t)=\frac{1}{2}(p_{a}K_{ab}(t)\,p_{b}+q_{a}K_{ab}(t)\,q_{b})+p_{a}L_{ab}\,q_{b}.$
(2.12)

Note in particular that it seems the problem of quantum control can be reduced
the problem of classical control theory by using quantum-symplectic duality,
see [10] for a discussion of the quantum control theory.

###  2.2 Unitary and symplectic evolution

In this subsection we discuss the relation between representations of
solutions of the Schrodinger/Hamiltonian equations in the complex unitary form
and the real symplectic form.

The solution of the N-component Schrodinger equation

$\psi(t)=U_{t}\,\psi(0),\qquad\psi(t)=\left(\psi_{1}(t),...\psi_{N}(t)\right),$
(2.13)

can be rewritten in term of real N-component vectors

$\psi(t)=q(t)+ip(t).$ (2.14)

Representing the unitary matrix $U_{t}$ in the form

$U_{t}=X_{t}+iY_{t},$ (2.15)

where $X_{t}$ and $Y_{t}$ are real matrices, we have

$\displaystyle\psi(t)=X_{t}q(0)-Y_{t}p(0)+i(Y_{t}q(0)+X_{t}p(0)),$ (2.16)

i.e.

$\displaystyle q(t)$ $\displaystyle=$ $\displaystyle X_{t}q(0)-Y_{t}p(0)$
(2.17) $\displaystyle p(t)$ $\displaystyle=$ $\displaystyle
Y_{t}q(0)+X_{t}p(0),$ (2.18)

These relations can be represented in the form

$\begin{pmatrix}q(t)\\\ p(t)\end{pmatrix}=S_{t}\begin{pmatrix}q(0)\\\
p(0)\end{pmatrix},$ (2.19)

where $S_{t}$ given by

$S_{t}=\left(\begin{array}[]{ll}\,\,\,X_{t}&Y_{t}\\\
-Y_{t}&X_{t}\end{array}\right).$ (2.20)

The unitarity of matrix $U_{t}$,

$U_{t}^{*}\,U_{t}=(X_{t}^{T}-iY_{t}^{T})(X_{t}+iY_{t})=X_{t}^{T}X_{t}+Y_{t}^{T}Y_{t}-i(Y_{t}^{T}X_{t}-X_{t}^{T}Y_{t})=I$
(2.21)

means that

$\displaystyle X_{t}^{T}X_{t}+Y_{t}^{T}Y_{t}$ $\displaystyle=$ $\displaystyle
I,$ (2.22) $\displaystyle Y_{t}^{T}X_{t}-X_{t}^{T}Y_{t}$ $\displaystyle=$
$\displaystyle 0$ (2.23)

These relation imply that $S_{t}$ is a symplectic matrix, i.e. satisfies the
relation

$S^{T}_{t}JS_{t}=J,$ (2.24)

here $S^{T}$ means the transposed matrix and

$J=\left(\begin{array}[]{ll}\,\,\,0&I\\\ -I&0\end{array}\right)$ (2.25)

Indeed,

$\displaystyle S^{T}_{t}JS_{t}$ $\displaystyle=$
$\displaystyle\left(\begin{array}[]{ll}\,\,\,X_{t}^{T}Y_{t}-Y_{t}^{T}X_{t}&\,\,Y_{t}^{T}Y_{t}+X_{t}^{T}X_{t}\\\
-Y_{t}^{T}Y_{t}-X_{t}^{T}X_{t}&\,-Y_{t}^{T}X_{t}+X_{t}^{T}Y_{t}\end{array}\right)$
(2.28)

and taking into account (2.22) and (2.23) we get

$\displaystyle S^{T}_{t}JS_{t}$ $\displaystyle=$
$\displaystyle\left(\begin{array}[]{ll}\,\,\,0&I\\\ -I&0\end{array}\right),$
(2.31)

i.e. equation (2.24).

###  2.3 Qubits and Symbits

In quantum computing, the fundamental unit of information is the qubit. The
qubit is a two-level quantum system, mathematically represented by the complex
Hilbert space $\mathbb{C}^{2}$. A general state of a qubit can be written as:

$|\psi\rangle=\alpha|0\rangle+\beta|1\rangle,\quad\alpha,\beta\in\mathbb{C},\quad|\alpha|^{2}+|\beta|^{2}=1,$
(2.32)

where $|0\rangle$ and $|1\rangle$ form an orthonormal basis for
$\mathbb{C}^{2}$. The coefficients $\alpha$ and $\beta$ are complex numbers
such that the total probability (norm) is equal to one.

In the symplectic computing framework, the analogous concept to a qubit is a
symplectic bit (symbit). A symbit is represented by a two-dimensional real
vector space $\mathbb{R}^{2}$. A general state of a symbit can be written as:

$|\phi)=\alpha|0)+\beta|1),\quad\alpha,\beta\in\mathbb{R},\quad\alpha^{2}+\beta^{2}=1,$
(2.33)

where $|0)$ and $|1)$ form an orthonormal basis for $\mathbb{R}^{2}$. The
coefficients $\alpha$ and $\beta$ here are real numbers such that the total
probability is equal to one. We can interpret $\mathbb{R}^{2}$ as the phase
plane for a dynamical system with coordinates $p$ and $q$.

One can consider tensor products of $N$ sbits:

$\mathbb{R}^{2}\otimes\mathbb{R}^{2}\otimes\cdots\otimes\mathbb{R}^{2}=\mathbb{R}^{2^{N}}.$
(2.34)

This is equivalent to $\mathbb{C}^{N}$. On the other hand, the tensor product
of $N$ qubits is $\mathbb{C}^{2^{N}}$. Therefore, by using tensor products of
sbits, one can obtain more general spaces than those obtained from tensor
products of qubits.

Operations on symbits are performed using symplectic transformations, which
are elements of the symplectic group $Sp(2,\mathbb{R})$. For instance, the
simplest symplectic transformation in $\mathbb{R}^{2}$ can be represented by a
$2\times 2$ matrix $S$ that preserves the symplectic form:

$S^{T}JS=J,\quad\text{where}\quad J=\begin{pmatrix}0&1\\\ -1&0\end{pmatrix}.$

An example of a symplectic transformation is the rotation matrix:

$R(\theta)=\begin{pmatrix}\cos\theta&\sin\theta\\\
-\sin\theta&\cos\theta\end{pmatrix},$

which preserves the symplectic form and hence is a member of
$Sp(2,\mathbb{R})$.

In summary, while qubits are the basic units of quantum information in
$\mathbb{C}^{2}$ and are manipulated using unitary transformations, symbits
are the basic units of symplectic information in $\mathbb{R}^{2}$ and are
manipulated using symplectic transformations. This duality provides a bridge
between quantum and symplectic computing, offering a new perspective on
computational processes.

##  3 Unitary and symplectic gates

###  3.1 Unitary group $U(N)$ and symplectic group $Sp(2N,{\mathbb{R}})$

Let us show that unitary group $U(N)$ is a subgroup of symplectic group
$Sp(2N,{\mathbb{R}})$. Let us remind that a symplectic matrix $S$ of
symplectic group $Sp(2N,\mathcal{R})$ satisfies the relation

$S^{T}JS=J$ (3.1)

where

$J=\left(\begin{array}[]{ccc}0\leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ I\\\ -I\leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ 0\end{array}\right).$ (3.2)

Let $S$ is a $2N\times 2N$ matrix in the form

$S=\left(\begin{array}[]{ll}A&B\\\ C&D\end{array}\right)$

where $A,B,C,D$ are $N\times N$ matrices. The conditions on $A,B,C$ and $D$
are

  * •

$A^{\mathrm{T}}C,B^{\mathrm{T}}D$ symmetric, and
$A^{\mathrm{T}}D-C^{\mathrm{T}}B=I$

  * •

$AB^{\mathrm{T}}$, $CD^{\mathrm{T}}$ symmetric, and
$AD^{\mathrm{T}}-BC^{\mathrm{T}}=I$

Let us show that there is a canonical mapping $\gamma$ of the unitary group
$U(N)$ to $Sp(2N,{\mathbb{R}})$. Let $V\in U(N)$, we present $V$ in the form

$V=X+iY,\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\text{where }X,Y\text{ are }N\times N\text{matrices with real entries}$ (3.3)

We define $\gamma$ by the following formula

$\gamma(V)=\gamma(X+iY)=\left(\begin{array}[]{ll}X&Y\\\
-Y&X\end{array}\right)$ (3.4)

One can check that the unitarity conditions $VV^{*}=V^{*}V=1$ lead to the
conditions for the symplectic matrices.

Remark 2. One can also check that there is a relation
$\gamma(V_{1})\gamma(V_{2})=\gamma(V_{1}V_{2})$.

Remark 3. The Schrodinger equation in the form

$i\dot{\psi}=H\,\psi,$ (3.5)

where

$H=-\Delta+V(x)$ (3.6)

$\Delta$ is the Laplace operator in $\mathbb{R}$, in the real formulation
looks as follows

$\dot{q}=Hp,\qquad\dot{p}=-Hq.$ (3.7)

Spectrum and scattering theory for in this approach are considered in [9]

###  3.2 NOT and $\gamma$(NOT) gates

Quantum computation is a sequence of unitary matrices (gates) of the simple
form. We put into correspondence into such unitary gate a corresponding
symplectic gate using the $\gamma$ defined by (3.4).

See the simplest gate is NOT. It is a unitary matrix
$NOT=\left(\begin{array}[]{ll}0&1\\\ -1&0\end{array}\right)$ defined on the
qubit $\mathbb{C}^{2}$. By using the mapping $\gamma$ we define the
corresponding symplectic gate

$\gamma(NOT)=\left(\begin{array}[]{llll}0&1&0&0\\\ -1&0&0&0\\\ 0&0&0&1\\\
0&0&-1&0\end{array}\right)$ (3.8)

One can construct the symplectic analog of the controlling CNOT gate.

##  4 Conclusion and Discussion

In this paper, we have introduced the concept of symplectic computers and
proposed that they have the potential to offer computational advantages over
quantum computers. By utilizing the quantum-symplectic duality, we
demonstrated that the Schrödinger equation, when expressed in real terms,
corresponds to a classical symplectic mechanical system. This allows for the
description of quantum systems in the framework of symplectic mechanics, with
symplectic bits (symbits) serving as analogs to quantum bits (qubits).

We have remind how the classical symplectic Hamiltonian equations can be
derived from the Schrodinger equation, enabling the use of symplectic
transformations and projectors in symplectic computation, similar to unitary
transformations and measurements in quantum computation. Our examination of
the relationship between unitary and symplectic evolutions revealed that
unitary matrices can be mapped to symplectic matrices, ensuring that quantum
gates have corresponding symplectic counterparts. This suggests that
symplectic computing can replicate the capabilities of quantum computing while
potentially providing additional benefits. Additionally, we discussed the
implications of symplectic mechanics for quantum computation, suggesting that
quantum control problems might be addressed using classical control theory
through quantum-symplectic duality. This could lead to new insights and
methods in both quantum and symplectic computing.

Let us note that our discussion has focused on linear symplectic computation.
Future research could explore non-linear symplectic transformations and the
application of dynamics on symplectic manifolds. Note that the complexity of
some manifolds, including black holes, has already been discussed [12]. It
would also be interesting to explore computational devices based on non-
quadratic Hamiltonian systems, such as the nonlinear Schrodinger equation.
This could provide a connection with approaches to quantum theory hidden in
classical probability [13]. Here, we discussed the finite-dimensional Hilbert
space. We expect that our considerations can be extended to the infinite-
dimensional Hilbert space in line with the considerations in the papers [3, 4,
11].

In this paper, we have discussed quantum and symplectic computers, which
consist of products of unitary or symplectic matrices. It seems natural to
extend the concept of computations to include sequences of elements from
arbitrary Lie groups. From this perspective, unitary representations of Lie
groups would be natural extensions of quantum computers. As a further
generalization, the $q$-deformations of Lie groups could be considered.
Finally, one can study computations on various number fields, including
$p$-adic, see [14, 15, 16].

In summary, the introduction of symplectic computers presents an interesting
development in the field of computation. By further exploring the theoretical
foundations and practical implementations of symplectic computing, we may
discover new computational techniques and applications that enhance both
quantum and classical computing paradigms.

##  5 Acknowledgements

I am very grateful to I. Ya. Aref’eva, V. Zh. Sakbaev, and D. Stepanenko for
fruitful discussions. This work is supported by the Russian Science Foundation
(project 24-11-00039, V.A. Steklov Mathematical Institute)

## References

  * [1] 40 years of quantum computing. Nat. Rev. Phys. 4, 1 (2022).
  * [2] Ohya, M., & Volovich, I. (2011), Mathematical foundations of quantum information and computation and its applications to nano-and bio-systems, Springer Science & Business Media, 2011.
  * [3] I.V. Volovich, “Complete Integrability of Quantum and Classical Dynamical Systems“, P-Adic Num Ultrametr Anal Appl 11, 328–334 (2019).
  * [4] I.V. Volovich, “Remarks on the complete integrability of quantum and classical dynamical systems“, arXiv:1911.01335 [math-ph].
  * [5] A. Einstein, Albert (1905). ”Über einen die Erzeugung und Verwandlung des Lichtes betreffenden heuristischen Gesichtspunkt”, Annalen der Physik. (1905) 17, 132–148.
  * [6] L. V. P. R. de Broglie, “Recherches sur la théorie des quanta,” Annals Phys. 2 (1925), 22-128
  * [7] R. P. Feynman, R. B. Leighton, and M. Sands, “The Feynman Lectures on Physics, Vol. 3,” Addison-Wesley, 1965.
  * [8] P. A. M. Dirac, “The Principles of Quantum Mechanics,” Oxford University Press, 4th edition, 1958.
  * [9] V. Zh. Sakbaev, I. V. Volovich, In preparation
  * [10] A. N. Pechen, D. J. Tannor, “Are there traps in quantum control landscapes?”, Phys. Rev. Lett., (2011), v. 106, 120402, 3 pp.
  * [11] V. Zh. Sakbaev, I. V. Volovich, “Analogues of Jacobi and Weyl Theorems for Infinite-Dimensional Tori“, Journal of Stochastic Analysis, (2024) v. 5, No. 1, Article 4, arXiv:2305.05929
  * [12] L. Susskind, “Computational Complexity and Black Hole Horizons,” Fortsch. Phys. 64, 24-43 (2016)
  * [13] Accardi, L., Dhahri, A. Orthogonal Polynomial Decomposition for Random Fields with All Moments. Milan J. Math. 87, 21–56 (2019)
  * [14] I.V. Volovich, “Number theory as the ultimate physical theory“, P-Adic Numbers, Ultrametric Analysis, and Applications, (2010). v. 2, 77-87; CERN-TH.4781/87
  * [15] B. Dragovich, A. Yu. Khrennikov, S. V. Kozyrev, I. V. Volovich, E. I. Zelenov, “p-Adic Mathematical Physics: The First 30 Years”, P-Adic Numbers Ultrametric Anal. Appl., (2017) v. 9, 87–121
  * [16] V. Anashin, “Toward the (non-cellular) automata interpretation of quantum mechanics: Volovich postulates as a roadmap“, Int. J. Mod. Phys. A (2022) v. 37, p. 2243003


capbtabboxtable[][]

11institutetext: Institute of Medical Informatics, University of Lübeck,
Germany
<EMAIL_ADDRESS>
# Voxelmorph++

Going beyond the cranial vault with keypoint supervision and multi-channel
instance optimisation

Mattias P. Heinrich 11 0000-0002-7489-1972
Lasse Hansen 11 0000-0003-3963-7052

###### Abstract

The majority of current research in deep learning based image registration
addresses inter-patient brain registration with moderate deformation
magnitudes. The recent Learn2Reg medical registration benchmark has
demonstrated that single-scale U-Net architectures, such as VoxelMorph that
directly employ a spatial transformer loss, often do not generalise well
beyond the cranial vault and fall short of state-of-the-art performance for
abdominal or intra-patient lung registration. Here, we propose two
straightforward steps that greatly reduce this gap in accuracy. First, we
employ keypoint self-supervision with a novel network head that predicts a
discretised heatmap and robustly reduces large deformations for better
robustness. Second, we replace multiple learned fine-tuning steps by a single
instance optimisation with hand-crafted features and the Adam optimiser.
Different to other related work, including FlowNet or PDD-Net, our approach
does not require a fully discretised architecture with correlation layer. Our
ablation study demonstrates the importance of keypoints in both self-
supervised and unsupervised (using only a MIND metric) settings. On a multi-
centric inspiration-exhale lung CT dataset, including very challenging COPD
scans, our method outperforms VoxelMorph by improving nonlinear alignment by
77% compared to 19% - reaching target registration errors of 2 mm that
outperform all but one learning methods published to date. Extending the
method to semantic features sets new stat-of-the-art performance on inter-
subject abdominal CT registration.

###### Keywords:

registration heatmaps deep learning.

##  1 Introduction

Medical image registration aims at finding anatomical and semantic
correspondences between multiple scans of the same patient (intra-subject) or
across a population (inter-subject). The difficulty of this task with
plentiful clinical applications lies in discriminating between changes in
intensities due to image appearance changes (acquisition protocol, density
difference, contrast, etc.) and nonlinear deformations. Advanced similarity
metrics may help in finding a good contrast-invariant description of local
neighbourhoods, e.g. normalised gradient fields [5] or MIND [12]. Due to the
ill-posedness of the problem some form of regularisation is often employed to
resolve the disambiguity between several potential local minima in the cost
function. Powerful optimisation frameworks that may comprise iterative
gradient descent, discrete graphical models or both (see [24] for an overview)
aim at solving for a global optimum that best aligns the overall scans (within
the respective regions of interest). Many deep learning (DL) registration
frameworks (e.g. DLIR [26] and VoxelMorph [1]) rely on a spatial transformer
loss that may be susceptible to ambiguous optimisation landscapes - hence
multiple resolution or scales levels need to be considered. The focus of this
work is to reflect such local minima more robustly in the loss function of DL-
registration. We propose to predict probabilistic displacement likelihoods as
heatmaps, which can better capture multiple scales of deformation within a
single feed-forward network.

Related work: Addressing large deformations with learning based registration
is generally approached by either multi-scale, label-supervised networks [13,
20, 19] or by employing explicitly discretised displacements [7, 9]. Many
variants of U-Net like architectures have been proposed that include among
others, dual-stream [16], cascades [28] and embeddings [3]. Different to those
works, we do neither explicitly model a discretised displacement space,
multiple scales or warps, nor modify the straightforward feed-forward U-Net of
DLIR or VoxelMorph.

Point-cloud registration is another research field of interest (FlowNet3d
[18]) that has however so far been restricted to lung registration in the
medical domain [6]. Stacked hourglass networks that predict discretised
heatmaps of well-defined anatomical landmarks are commonly used in human pose
estimation [22]. They are, however, restricted to datasets and registration
applications where not only pairwise one-to-one correspondences can be
obtained as training objective but generic landmarks have to be found across
all potential subjects. Due to anatomical variations this restriction often
prevents their use in medical registration. Combining an initial robust
deformation prediction with instance optimisation or further learning fine-
tuning steps (Learn-to-optimise, cf. [25]) has become a new trend in learning-
based registration, e.g. ([23, 14]). This paradigm, which expects coarse but
large displacements from a feed-forward prediction, can shift the focus away
from sub-pixel accuracy and towards avoiding failure cases in the global
transformation due to local minima. It is based on the observation that local
iterative optimisers are very precise when correctly initialised and have
become extremely fast due to GPU acceleration.

Contributions: We demonstrate that a single-scale U-Net without any bells and
whistles in conjunction with a fast MIND-based instance optimisation can
achieve or outperform state-of-the-art in large-deformation registration. This
is achieved by focussing on coarse scale global transformation by introducing
a novel heatmap prediction network head. In our first scenario we employ weak
self-supervision, through automatic keypoint correspondences [10]. Here, the
heatmap enables a discretised integral regression of displacements to directly
and explicitly match the keypoint supervision. Second, we incorporate the
heatmap prediction into a non-local unsupervised metric loss. This enables a
direct comparison within the same network architecture to the commonly used
spatial transformer (warping) loss in unsupervised DL registration and
highlights the importance of providing better guidance to avoid local minima.
Our extensive ablation experiments with and without instance optimisation on a
moderately large and challenging inspiration-exhale lung dataset demonstrate
state-of-the-art performance.

Our code is publicly available at:

https://www.github.com/mattiaspaul/VoxelMorphPlusPlus

Figure 1: Overview of method and qualitative result for held-out case #1 from
[21]. The key new element of our approach is the heatmap prediction head that
is appended to a standard VoxelMorph. It helps overcome local minima through a
probabilistic loss with either automatic keypoint correspondences or non-
locally weighted MIND features.

##  2 Method

Keypoint correspondences are an excellent starting point to explore the
benefits of incorporating a heatmap prediction within DL-registration. Our
method can be either trained when those automatically computed displacements
at around $|K|\approx 2000$ locations per scan are available for training
data, or we can use a non-local unsupervised metric loss (see details below).
In both scenarios, we use Förstner keypoints in the fixed scan to focus on
distinct locations within a region of interest. We will first describe the
baseline single-scale U-Net backbone, followed by our novel heatmap prediction
head and the non-local MIND loss (which is an extension of [11] to 3D).

Baseline backbone: Given two input CT scans, $F\rightarrow\mathbb{R}^{3}$
fixed image and $M\rightarrow\mathbb{R}^{3}$ moving image and a region of
interest $\Omega\in\mathbb{R}^{3}$, we firstly define a feed-forward U-Net [4]
$\Theta(F,M,\Omega,\theta)\rightarrow\mathbb{R}^{C}$ with trainable parameters
$\theta$ that maps the concatenated input towards a shared $C$-dimensional
feature representation $\mathbf{z}$ (we found C$\approx$64 is expressive
enough to represent displacements). This representation may have a lower
spatial resolution than $F$ or $M$ and is the basis for predicting a (sparse)
displacement field $\varphi$ that spatially aligns $F$ and $M$ within
$\Omega$. $\Theta$ comprises in our implementation a total of eleven 3D
convolution blocks, each consisting of a $3\times 3\times 3$ convolution,
instance normalisation (IN), ReLU, a $1\times 1\times 1$ convolution, and
another IN+ReLU. Akin to VoxelMorph, we use $2\times 2\times 2$ max-pooling
after each of the four blocks in the encoder and nearest neighbour upsampling
to restore the resolution in the decoder, but use a half-resolution output.
The network has up to $C=64$ hidden feature channels and 901’888 trainable
parameters.

Due to the fact that this backbone already contains several convolution blocks
on the final resolution at the end of the decoder, it is directly capable of
predicting a continuous displacement field $\varphi$ by simply appending three
more $1\times 1\times 1$ convolutions (and IN+ReLU) with a number of output
channels equal to 3.

Discretised heatmap prediction head:  The aim of the heatmap prediction head
is to map a $C$-dimensional feature vector (interpreted as a $1\times 1\times
1$ spatial tensor with $|K|$ being the batch dimension) into a discretised
displacement tensor $\mathbf{y}\in\mathcal{Q}$ with predefined size and
spatial range $R$ (see Fig. 1). Here we chose $R=0.3$, in a coordinate system
that ranges from $-1$ to $+1$, which captures even large lung motion between
respiratory states. We set $\mathcal{Q}\rightarrow\mathbb{R}^{11\times
11\times 11}$ to balance computational complexity and limit quantisation
effects. This means we need to design another nested decoder that increases
the spatial resolution from 1 to 11. Our heatmap network comprises a transpose
3D convolution with kernel size $7\times 7\times 7$, six further 3D
convolution blocks (kernel size $7\times 7\times 7$ and IN+ReLU) once
interleaved with a single trilinear upsampling to $11\times 11\times 11$. It
has 462’417 trainable parameters and its number of output channels is equal to
1.

Next, we can define a probabilistic displacements tensor
$\mathcal{P}\rightarrow\mathbb{Z}^{6}$ using a softmax operation along the
combined displacement dimensions as:

$\mathcal{P}(\mathbf{x},\Delta\mathbf{x})=\frac{\exp(y(\mathbf{x},\Delta\mathbf{x}))}{\sum_{\Delta\mathbf{x}}\exp(y(\mathbf{x},\Delta\mathbf{x}))},$
(1)

where $\mathbf{x}\rightarrow\mathbb{Z}^{3}$ are global spatial 3D coordinates
and $\Delta\mathbf{x}\rightarrow\mathbb{Z}^{3}$ local 3D displacements. In
order to define a continuous valued displacement field, we apply a weighted
sum:

$\varphi(\mathbf{x})=\sum_{\Delta\mathbf{x}}\mathcal{P}(\mathbf{x},\Delta\mathbf{x})\cdot\mathcal{Q}(\Delta\mathbf{x})$
(2)

This output is used during training to compute a mean-squared error between
predicted and pre-computed keypoint displacements. Since, the training
correspondences are regularised using a graphical model, we require no further
penalty.

Non-local MIND loss: To avoid the previously described pitfalls of directly
employing a spatial transformer (warping) loss, we can better employ the
probabilistic heatmap prediction and compute the discretely warped MIND
vectors of the moving scan implicitly by a weighted average of the underlying
features within pre-defined capture region (where $c$ describes one of the 12
MIND channels) as:

$\mathrm{MIND}_{\mathrm{warped}}(c,\mathbf{x})=\sum_{\Delta\mathbf{x}}\mathcal{P}(\mathbf{x},\Delta\mathbf{x})\cdot\mathrm{MIND}(c,\mathbf{x}+\Delta\mathbf{x})$
(3)

.

Implementation details: Note that the input to both the small regression
network (baseline) and our proposed are feature vectors sampled at the
keypoint locations, which already improves the baseline architecture slightly.
We use trilinear interpolation in all cases where the input and output grids
differ in size to obtain off-grid values. All predicted sparse displacements
$\varphi$ are extrapolated to a dense field using thin-plate-splines with
$\lambda=0.1$ that yields $\varphi^{*}$.

For the baseline regression setup (VoxelMorph) we employ a common MIND warping
loss and a diffusion regularisation penalty that is computed based on the
Laplacian of a kNN-graph ($k=7$) between the fixed keypoints. The weighting of
the regularisation was empirically set to $\alpha=0.25$. We found that using
spatially aggregated CT and MIND tensors the former using average pooling with
kernel size 2, the latter two of those pooling steps, leads to stabler
training in particular for the regression baseline.

Multi-channel instance optimisation: We directly follow the implementation
described in [23]111https://github.com/multimodallearning/convexAdam. It is
initialised with $\varphi^{*}$, runs for 50 iterations, employs a combined
B-spline and diffusion regularisation coupled with a MIND metric loss and a
grid spacing of 2. This step is extremely fast, but relies on a robust
initialisation as we will demonstrate in our experiments. The method can also
be employed when semantic features, e.g. segmentation predictions from an
nnUNet [17], are available in the form of one-hot tensors.

##  3 Experiments and Results

We perform extensive experiments on inspiration-exhale registration of lung CT
scans - arguably one of the most challenging tasks in particular for learning-
based registration [14]. A dataset of 30 scan pairs with large respiratory
differences is collected from EMPIRE10 (8 scan pairs #1, #7, #8, #14, #18,
#20, #21 and #28) [21], Learn2Reg Task 2 [15] and DIR-Lab COPD [2] (10 pairs).
The exhale and scans are resampled to $1.75\times 1.25\times 1.75~{}$mm and
$1.75\times 1.00\times 1.25~{}$mm respectively to account for average overall
volume scaling and a fixed region with dimensions $192\times 192\times 208$
voxels was cropped that centres the mass of automatic lung masks. Note that
this pre-processing approximately halves the initial target registration error
(TRE) of the COPD dataset. Lung masks are also used to define a region-of-
interest for the loss evaluation and to mask input features for the instance
optimisation. We split the data into five folds for cross-validation that
reflect the multi-centric data origin (i.e. approx. two scans per centre are
held out for validation each).

[6.3cm] [5.3cm]  UNet Heatmap Keypoints Error w/o Adam Error w/ Adam Initial /
Adam 10.04 vx 7.41 vx VoxelMorph ✓ 8.17 vx 4.80 vx VM + Heatmap ✓ ✓ 6.49 vx
3.18 vx VM + Keypoints ✓ ✓ 6.30 vx 2.79 vx VoxelMorph++ ✓ ✓ ✓ 5.31 vx 2.34 vx

Figure 2: Cumulative keypoint error of proposed model compared to a VoxelMorph
baseline and using only Adam instance optimisation with MIND. Figure 3:
Results of ablation study on lung CT: VoxelMorph++ improves error reduction of
nonlinear alignment from 18% to 77%.

Keypoint Self-supervision: To create correspondence as self-supervision for
our proposed VoxelMorph++ method, we employ the
corrField[10]222http://www.mpheinrich.de/code/corrFieldWeb.zip, which is
designed for lung registration and based on a discretised displacement search
and a Markov random field optimisation with multiple task specific
improvements. It runs within a minute per scan pair and creates $\approx
2000=|K|$ highly accurate ($\approx 1.68$ mm) correspondences at Förstner
keypoints.

As additional experiment we extend our method to the inter-subject alignment
of 30 abdominal CTs [27] that was also part of the Learn2Reg 2020 challenge
(Task 3) and provides 13 difficult anatomical organ labels for training and
validation.

Ablation study: We consider a five-fold cross-validation for all ablation
experiments with an initial error of 10.04 vx (after translation and scaling)
across 30 scan pairs computed based on keypoint correspondences. Employing
only the Adam instance optimisation with MIND features results in an error of
8.17 vx, with default settings of grid spacing $=2$ voxels, 50 iterations and
$\lambda_{Adam}=0.65$. Note that a dense displacement is estimated with a
parametric B-spline model. We start from the slightly improved VoxelMorph
baseline with MIND loss, diffusion regularisation and increased number of
convolution operations described above. This yields a keypoint error of 8.17
vx that represents an error reduction of 19% and can be further improved to
4.80 vx when adding instance optimisation. A weighting parameter
$\lambda=0.75$ for diffusion regularisation was empirically found with $k=7$
for the sparse neighbourhood graph of keypoints. Replacing the traditional
spatial transformer loss with our proposed heatmap prediction head that uses
the nonlocal MIND loss much improves the performance to 6.49 vx and 3.18 vx
(with and without Adam respectively). But the key improvement can be gained
when including the self-supervised keypoint loss. Using our baseline
VoxelMorph architecture that regresses continuous 3D vectors, we reach 6.30
and 2.79 vx. See Table 3 and Fig. 3 for numerical and cumulative errors. Our
heatmap-based network and the instance optimisation require around 0.43 and
0.41 seconds inference time, respectively. The complexity of transformations
measured as standard deviation of log-Jacobian determinants is on average
0.0554. The number of negative values is zero (no folding) in 8 out of 10 COPD
cases and negligible ($<10^{-4}$) in the others.

Comparison to state-of-the-art: When evaluation the target registration error
(TRE) in mm for the 10 pairs of DIR-Lab COPD [2] one of the most challenging
benchmarks in medical registration with an initial misalignment of 23.36 mm
(and 12.0 mm after pre-alignment), we reach 2.16 mm. This compares very
favourable to VoxelMorph+ with 7.98mm and LapIRN with 4.76mm. Of all published
DL-methods only GraphRegNet [7] is superior with 1.34 mm. The high visual
quality of our registration is shown in Fig. 4.

Figure 4: Exemplary results of our proposed method before (top row) and after
(bottom row) registration. Fixed exhale scans are shown in blue and
inspiration in orange shades respectively (adding up to grayscale when
aligned). For abdominal alignment transformed segmentation labels are shown,
here: right kidney $\blacksquare$, left kidney $\blacksquare$, gallbladder
$\blacksquare$, liver $\blacksquare$, stomach $\blacksquare$, aorta
$\blacksquare$ and pancreas $\blacksquare$ are visible.

Limitations and further potential: We have not yet considered more advanced
network architectures as backbone, e.g. two-stream or multi-level, which are
likely to yield further improvements. However, based on our experiments we
expect that it could merely reduce the reliance on instance optimisation.

Inter-subject abdominal CT registration: We apply our proposed VoxelMorph++
model with nonlocal loss and no architectural modifications to another
challenging task of inter-subject abdominal CT registration with initial Dice
overlap for 13 organs of only 25.9% and weakly supervised learning (45
registration pairs). Following [8], we decouple the semantic feature
extraction and directly train an nnUNet model [17]. The best published
VoxelMorph model that was trained with label-supervision and extended to a
two-stream architecture reached 43.9% [23], the two top-ranked methods of the
Learn2Reg challenge yield 65.7% (ConvexAdam [23]) and 67% (LapIRN [20])
respectively. Directly employing instance optimisation with 25 iterations on
the nnUNet features achieves 62.9%. We use 2048 keypoints that are sampled
inversely proportional to the predicted label maps and employ two warps (and
inverse consistency for the first of them). Our model substantially
outperforms VoxelMorph with 52.3% and sets a new state-of-the-art performance
after instance optimisation reaching 69.6% with a total run time of less than
a second.

##  4 Discussion and conclusions

: Our results demonstrate that contrary to previous belief, a simple single-
scale U-Net architecture can provide large deformation estimation that is
robust enough to reach high accuracy with a subsequent instance optimisation.
The key insight of our work is the importance to predict a discretised heatmap
to alleviate the problematic direct regression and use strong self-supervision
either using automatic keypoint correspondences or a nonlocal multichannel
loss together with a straightforward instance optimisation. Our work is
related to mlVIRN [13], which also incorporates a keypoint loss for lung
registration in addition to lobe segmentations, but has to be trained with
several hundreds of paired CTs and did not report TRE values for DIRlab-COPD.
Our network can be trained within 17 minutes on a single RTX A4000 requiring
less than 2 GByte of VRAM, indicating the improved training efficiency with
fewer scans when using heatmaps. GraphRegNet [7] is similar in that it also
employs heatmaps (integral regression) but more explicitly by defining the
exact same discretised displacement grid beforehand and computing an SSD cost
tensor based on hand-crafted features as input. While it outperforms our
method with a TRE of 1.34mm it appears to be more tailored towards the
specific task and might not be easily extendable to end-to-end feature
learning or abdominal registration.

## References

  * [1] Balakrishnan, G., Zhao, A., Sabuncu, M.R., Guttag, J., Dalca, A.V.: Voxelmorph: a learning framework for deformable medical image registration. IEEE transactions on medical imaging 38(8), 1788–1800 (2019)
  * [2] Castillo, R., Castillo, E., Fuentes, D., Ahmad, M., Wood, A.M., Ludwig, M.S., Guerrero, T.: A reference dataset for deformable image registration spatial accuracy evaluation using the copdgene study archive. Physics in Medicine & Biology 58(9), 2861 (2013)
  * [3] Estienne, T., Vakalopoulou, M., Battistella, E., Henry, T., Lerousseau, M., Leroy, A., Paragios, N., Deutsch, E.: Mics : Multi-steps, inverse consistency and symmetric deep learning registration network (2021)
  * [4] Falk, T., Mai, D., Bensch, R., Çiçek, Ö., Abdulkadir, A., Marrakchi, Y., Böhm, A., Deubner, J., Jäckel, Z., Seiwald, K., et al.: U-net: deep learning for cell counting, detection, and morphometry. Nature methods 16(1), 67–70 (2019)
  * [5] Haber, E., Modersitzki, J.: Intensity gradient based registration and fusion of multi-modal images. Methods of information in medicine 46(03), 292–299 (2007)
  * [6] Hansen, L., Dittmer, D., Heinrich, M.P.: Learning deformable point set registration with regularized dynamic graph cnns for large lung motion in copd patients. In: International Workshop on Graph Learning in Medical Imaging. pp. 53–61. Springer (2019)
  * [7] Hansen, L., Heinrich, M.P.: Graphregnet: Deep graph regularisation networks on sparse keypoints for dense registration of 3d lung cts. IEEE Transactions on Medical Imaging 40(9), 2246–2257 (2021)
  * [8] Hansen, L., Heinrich, M.P.: Revisiting iterative highly efficient optimisation schemes in medical image registration. In: International Conference on Medical Image Computing and Computer-Assisted Intervention. pp. 203–212. Springer (2021)
  * [9] Heinrich, M.P.: Closing the gap between deep and conventional image registration using probabilistic dense displacement networks. In: International Conference on Medical Image Computing and Computer-Assisted Intervention. pp. 50–58. Springer (2019)
  * [10] Heinrich, M.P., Handels, H., Simpson, I.J.: Estimating large lung motion in copd patients by symmetric regularised correspondence fields. In: International conference on medical image computing and computer-assisted intervention. pp. 338–345. Springer (2015)
  * [11] Heinrich, M.P., Hansen, L.: Highly accurate and memory efficient unsupervised learning-based discrete ct registration using 2.5 d displacement search. In: International Conference on Medical Image Computing and Computer-Assisted Intervention. pp. 190–200. Springer (2020)
  * [12] Heinrich, M.P., Jenkinson, M., Bhushan, M., Matin, T., Gleeson, F.V., Brady, M., Schnabel, J.A.: Mind: Modality independent neighbourhood descriptor for multi-modal deformable registration. Medical image analysis 16(7), 1423–1435 (2012)
  * [13] Hering, A., Häger, S., Moltz, J., Lessmann, N., Heldmann, S., van Ginneken, B.: Cnn-based lung ct registration with multiple anatomical constraints. Medical Image Analysis p. 102139 (2021)
  * [14] Hering, A., Hansen, L., Mok, T.C.W., Chung, A.C.S., Siebert, H., Häger, S., Lange, A., Kuckertz, S., Heldmann, S., Shao, W., Vesal, S., Rusu, M., Sonn, G., Estienne, T., Vakalopoulou, M., Han, L., Huang, Y., Brudfors, M., Balbastre, Y., Joutard, S., Modat, M., Lifshitz, G., Raviv, D., Lv, J., Li, Q., Jaouen, V., Visvikis, D., Fourcade, C., Rubeaux, M., Pan, W., Xu, Z., Jian, B., Benetti, F.D., Wodzinski, M., Gunnarsson, N., Sjölund, J., Qiu, H., Li, Z., Großbröhmer, C., Hoopes, A., Reinertsen, I., Xiao, Y., Landman, B., Huo, Y., Murphy, K., Lessmann, N., van Ginneken, B., Dalca, A.V., Heinrich, M.P.: Learn2reg: comprehensive multi-task medical image registration challenge, dataset and evaluation in the era of deep learning (2021)
  * [15] Hering, A., Murphy, K., van Ginneken, B.: Learn2Reg Challenge: CT Lung Registration - Training Data (May 2020). https://doi.org/10.5281/zenodo.3835682, https://doi.org/10.5281/zenodo.3835682
  * [16] Hu, X., Kang, M., Huang, W., Scott, M.R., Wiest, R., Reyes, M.: Dual-stream pyramid registration network. In: International Conference on Medical Image Computing and Computer-Assisted Intervention. pp. 382–390. Springer (2019)
  * [17] Isensee, F., Jaeger, P.F., Kohl, S.A., Petersen, J., Maier-Hein, K.H.: nnu-net: a self-configuring method for deep learning-based biomedical image segmentation. Nature methods 18(2), 203–211 (2021)
  * [18] Liu, X., Qi, C.R., Guibas, L.J.: Flownet3d: Learning scene flow in 3d point clouds. In: Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition. pp. 529–537 (2019)
  * [19] Mok, T.C., Chung, A.: Conditional deformable image registration with convolutional neural network. In: International Conference on Medical Image Computing and Computer-Assisted Intervention. pp. 35–45. Springer (2021)
  * [20] Mok, T.C., Chung, A.C.: Large deformation diffeomorphic image registration with laplacian pyramid networks. In: International Conference on Medical Image Computing and Computer-Assisted Intervention. pp. 211–221. Springer (2020)
  * [21] Murphy, K., Van Ginneken, B., Reinhardt, J.M., Kabus, S., Ding, K., Deng, X., Cao, K., Du, K., Christensen, G.E., Garcia, V., et al.: Evaluation of registration methods on thoracic ct: the empire10 challenge. IEEE transactions on medical imaging 30(11), 1901–1920 (2011)
  * [22] Newell, A., Yang, K., Deng, J.: Stacked hourglass networks for human pose estimation. In: European conference on computer vision. pp. 483–499. Springer (2016)
  * [23] Siebert, H., Hansen, L., Heinrich, M.P.: Fast 3d registration with accurate optimisation and little learning for learn2reg 2021 (2021)
  * [24] Sotiras, A., Davatzikos, C., Paragios, N.: Deformable medical image registration: A survey. IEEE transactions on medical imaging 32(7), 1153–1190 (2013)
  * [25] Teed, Z., Deng, J.: Raft: Recurrent all-pairs field transforms for optical flow. In: European conference on computer vision. pp. 402–419. Springer (2020)
  * [26] de Vos, B.D., Berendsen, F.F., Viergever, M.A., Sokooti, H., Staring, M., Išgum, I.: A deep learning framework for unsupervised affine and deformable image registration. Medical image analysis 52, 128–143 (2019)
  * [27] Xu, Z., Lee, C.P., Heinrich, M.P., Modat, M., Rueckert, D., Ourselin, S., Abramson, R.G., Landman, B.A.: Evaluation of six registration methods for the human abdomen on clinically acquired ct. IEEE Transactions on Biomedical Engineering 63(8), 1563–1572 (2016)
  * [28] Zhao, S., Lau, T., Luo, J., Eric, I., Chang, C., Xu, Y.: Unsupervised 3d end-to-end medical image registration with volume tweening network. IEEE journal of biomedical and health informatics 24(5), 1394–1404 (2019)