# The $a_{0}(980)$ in the single Cabibbo-suppressed process
$\Lambda_{c}\to\pi^{0}\eta p$

Hong-Shen Li  School of Physics and Microelectronics, Zhengzhou University,
Zhengzhou, Henan 450001, China  Le-Le Wei  School of Physics and
Microelectronics, Zhengzhou University, Zhengzhou, Henan 450001, China  Man-Yu
Duan  School of Physics and Microelectronics, Zhengzhou University, Zhengzhou,
Henan 450001, China  En Wang<EMAIL_ADDRESS>School of Physics and
Microelectronics, Zhengzhou University, Zhengzhou, Henan 450001, China  De-Min
Li  School of Physics and Microelectronics, Zhengzhou University, Zhengzhou,
Henan 450001, China

###### Abstract

In this work, we have investigated the Cabibbo-suppressed process
$\Lambda_{c}\to\pi^{0}\eta p$, by taking into account the intermediate scalar
state $a_{0}(980)$, which could be dynamically generated from the $S$-wave
pseudoscalar-pseudoscalar interaction within the chiral unitary approach. We
have calculated the $\pi^{0}\eta$ invariant mass distribution, and found that
there is a significant structure associated to the $a_{0}(980)$. We have also
roughly estimated the branching fraction
$\mathcal{B}(\Lambda_{c}\to\pi^{0}\eta p)\approx 10^{-4}$. We encourage our
experimental colleagues to measure the process $\Lambda_{c}\to\pi^{0}\eta p$
for searching for the state $a_{0}(980)$ signal in this reaction, which can be
used to test the molecular nature of $a_{0}(980)$ and our theoretical
predictions.

##  I Introduction

Investigating the dynamic structure of resonances is one of the important
directions in hadron physics, and the decay mechanism of the lightest charm
baryon $\Lambda_{c}$ is of great significance to understand the interplay of
the weak and strong interaction in the charm region Cheng:2015iom ;
Ebert:1983ih ; Lu:2016ogy ; Geng:2018upx ; Oset:2016lyh . In the last decades,
lots of experimental information has been accumulated BESIII:2018mes ;
LHCb:2018mes ; BESIII:2018evi ; Belle:2018obs ; Belle:2017sea ; BESIII:2017obs
; BESIII:2019mes , and there exist many theoretical studies about
$\Lambda_{c}$ decay Chao-Qiang Geng:2019sin ; Hai-Yang Cheng:2018sin ; R.
Pavao:2018pro ; Xie:2017erh ; Xiao-Hai Liu:2019vis ; Jung Keun Ahn:2019hyp ;
Ju-Jun Xie:2016the ; Ju-Jun Xie:2017int .

Recently, the BESIII Collaboration has reported the $K^{+}K^{-}$ invariant
mass distribution of the process $\Lambda_{c}\to pK^{+}K^{-}$, which shows an
enhancement structure near the $K^{+}K^{-}$ threshold BESIII:2018mes . In Ref.
Wang:2020 , we have analyzed the experimental measurement of the
$\Lambda_{c}\to pK^{+}K^{-}$, and concluded that the enhancement structure
near the threshold is mainly due to the resonance $f_{0}(980)$, and
$a_{0}(980)$ provides a small contribution. Since both $a_{0}(980)$ and
$f_{0}(980)$ couple to the $K\bar{K}$, the mechanisms of their productions are
crucial for us to deeply understand about the $K\bar{K}$ enhancement in the
process $\Lambda_{c}\to pK^{+}K^{-}$ and other processes. For instance, the
LHCb Collaboration has argued that the $a_{0}(980)$ plays a more important
role than the $f_{0}(980)$ in decay $\bar{B}^{0}\to J/\psi K^{+}K^{-}$, and
reported the branching fraction $\mathcal{B}(\bar{B}^{0}\to J/\psi
a_{0}(980),a_{0}(980)\to K^{+}K^{-})=(4.70\pm 3.31\pm 0.27)\times 10^{-7}$
Aaij:2013mtm , however, both the $f_{0}(980)$ and $a_{0}(980)$ resonances are
expected to contribute to the $K^{+}K^{-}$ distribution, as discussed in Ref.
Liang:2015qva .

Taking into account the uncertainties of the experimental measurements and the
undefined theoretical parameters, it is still difficult to extract the
relative weight of the $a_{0}(980)$ and $f_{0}(980)$ in the processes
involving $K\bar{K}$ final states, which is necessary for us to understand the
production mechanisms of $a_{0}(980)$ and $f_{0}(980)$ in the $K^{+}K^{-}$
channel. Since the $a_{0}(980)$ mainly decays to $\pi\eta$ channel, we propose
to investigate the process $\Lambda_{c}\to\pi^{0}\eta p$, which should be
useful to understand the production mechanism of the $a_{0}(980)$. The
reaction of $\gamma p\to p\pi^{0}\eta$ has been measured by Crystal Barrel at
the electron stretcher accelerator ELSA, where the $a_{0}(980)$ is clearly
identified Gutz:2014wit . Up to our knowledge, the process
$\Lambda_{c}\to\pi^{0}\eta p$ has not yet been investigated theoretically and
experimentally.

For the nature of the $a_{0}(980)$, there are many theoretical explanations,
such as tetraquark, molecular state (see the review ‘Scalar mesons below 2
GeV’ of Particle Data Group (PDG) PDG2018 ). The internal structure of the
$a_{0}(980)$ is essential to establish the decay mechanism of the process
$\Lambda_{c}\to\pi^{0}\eta p$, and is also crucial for us to understand the
spectrum of the scalar mesons with isospin $I=1$ Wang:2017pxm . Among those
different explanations, the molecular nature of the $a_{0}(980)$, dynamically
generated from the $S$-wave pseudoscalar-pseudoscalar interaction within the
chiral unitary approach Oller:1997chi ; Nieves:1999bet ; Cabrera:2005eva ;
Guo:2006dyn , has been widely studied in literature, such as $B^{0}\rightarrow
D^{0}\pi^{0}\eta$ Liang:2014ama , $\bar{B}^{0}\rightarrow J/\psi\pi^{0}\eta$
Liang:2015qva , $\bar{B}^{0}\rightarrow\eta_{c}\pi^{0}\eta$ Xie:2018rqv ,
$\chi_{c1}\rightarrow\eta\pi^{+}\pi^{-}$ Liang:2016hmr ,
$\eta_{c}\rightarrow\eta\pi^{+}\pi^{-}$ Debastiani:2017 ,
$D^{0}\rightarrow\bar{K}^{0}\pi^{0}\eta$ Xie:2014tma ,
$D^{+}\to\pi^{+}\pi^{0}\eta$ Duan:2020vye ; Ikeno:2021kzf , $D^{0}\to
K^{-}\pi^{+}\eta$ Toledo:2020zxj , $D^{+}_{s}\to\pi^{+}\pi^{0}\eta$
Molina:2019udw ; Hsiao:2019ait ; Ling:2021qzl ,
$\tau^{-}\rightarrow\nu_{\tau}\pi^{-}\pi^{0}\eta$ Dai:2018rra ,
$J/\psi\rightarrow K^{+}K^{-}a_{0}(980)$ Liang:2019jtr ,
$J/\psi\to\gamma\pi^{0}\eta$ Xiao:2019lrj , and $\Lambda_{c}\to K^{+}K^{-}p$
Wang:2020 .

In this work, we will investigate the process $\Lambda_{c}\to\pi^{0}\eta p$,
by considering the picture of $a_{0}(980)$ being dynamically generated from
the $S$-wave pseudoscalar-pseudoscalar interaction in the chiral unitary
approach. By presenting our predictions for the $\pi^{0}\eta$ invariant mass
distribution and the branching fraction of this reaction, we would like to
provide a motivation for experimentalists to measure this process. The
experimental information of this process could be useful to learn about the
nature of $a_{0}(980)$, and also constrain the parameters that is crucial to
understand the mechanism of the $a_{0}(980)$ and $f_{0}(980)$ production in
the $K^{+}K^{-}$ system.

The paper is organized as follows. In Sect. II, we introduce our model for the
process $\Lambda_{c}\to\pi^{0}\eta p$, which can be divided into three steps,
weak process, hadronization, and the final state interaction. The results and
discussions are given in Sect. III, a short summary is given in the last
section.

##  II Formalism

Figure 1: The microscopic quark level diagrams for the process
$\Lambda_{c}\rightarrow\pi^{0}\eta p$, (a) the $W^{+}$ internal emission
($c\to dW^{+}\to d\bar{d}u$) with $\bar{q}q$ produced from the vacuum between
$d$ and $\bar{d}$, (b) the $W^{+}$ internal emission ($c\to sW^{+}\to
s\bar{s}u$) with $\bar{q}q$ produced from the vacuum between $s$ and
$\bar{s}$, (c) the $W^{+}$ internal emission ($c\to dW^{+}\to d\bar{d}u$) with
$\bar{u}u$ produced from the vacuum between $u$ and $ud$, (d) the $W^{+}$
internal emission ($c\to sW^{+}\to s\bar{s}u$) with $\bar{u}u$ produced from
the vacuum between $u$ and $ud$, and (e) the $W^{+}$ external emission ($c\to
sW^{+}\to s\bar{s}u$) with $\bar{u}u$ produced from the vacuum between $s$ and
$ud$.

In analogy to Refs. Wang:2020 ; Wang:2015pcn ; Lu:2016roh ; Dai:2018nmw , the
single Cabibbo-suppressed decay of $\Lambda_{c}\rightarrow\pi^{0}\eta p$ can
proceed via three steps: the weak process, the hadronization, and the final
state interaction. Firstly, the $c$ quark of the initial $\Lambda_{c}$ decays
into a $W^{+}$ boson and a $d$ (or $s$) quark by weak interaction, then
$W^{+}$ boson turns into $u\bar{d}$ (or $u\bar{s}$) pair. In order to give
rise to the final states of the process $\Lambda_{c}\rightarrow\pi^{0}\eta p$,
the $ud$ pair of the initial $\Lambda_{c}$ and the $d\bar{d}$ (or $s\bar{s}$)
quarks from the weak decay must be hadronized with $\bar{q}q$
($\equiv\bar{u}u+\bar{d}d+\bar{s}s$) created from vacuum with $J^{PC}=0^{++}$,
which can be classified as the $W^{+}$ internal emission and $W^{+}$ external
emission, respectively depicted in Figs. 1(a)-(d) and Fig. 1(e).

For the $W^{+}$ internal emission mechanisms of Figs. 1(a) and 1(b), the
$d\bar{d}$ or $s\bar{s}$ quarks from the weak decay are hadronized with
$\bar{q}q$ ($\equiv\bar{u}u+\bar{d}d+\bar{s}s$) created from vacuum. It should
be pointed out that, in Fig. 1(b), the $s\bar{s}$ pair with isospin $I=0$,
together with the created $\bar{q}q$ pair, can not be hadronized into the
$\pi^{0}\eta$ system with isospin $I=1$, which implies that the mechanism of
Fig. 1(b) has no contribution to the process $\Lambda_{c}\to\pi^{0}\eta p$.
For Fig. 1(a), after hadronization we have the flavor wave function of the
final state from the $\Lambda_{c}$ weak decay,

$\displaystyle\Lambda_{c}$ $\displaystyle=$
$\displaystyle\frac{1}{\sqrt{2}}c(ud-du)$ (1) $\displaystyle\Rightarrow$
$\displaystyle
V_{cd}V_{ud}d\left(\bar{u}u+\bar{d}d+\bar{s}s\right)\bar{d}u\frac{1}{\sqrt{2}}\left(ud-
du\right)$ $\displaystyle=$ $\displaystyle
V_{cd}V_{ud}\sum_{i}\left(M_{2i}\times M_{i2}\right)p$ $\displaystyle=$
$\displaystyle V_{cd}V_{ud}\left(M^{2}\right)_{22}p,$

where $V_{cd}$ and $V_{ud}$ are the CKM matrix elements, and the flavor wave
functions of the baryons are $p=u(ud-du)/\sqrt{2}$ and $\Lambda_{c}=c(ud-
du)/\sqrt{2}$. $M$ is the matrix in terms of the pseudoscalar mesons, using
the standard $\eta$-$\eta^{\prime}$ mixing Roca:2003uk ; Bayar:2014qha ;
Liang:2014ama ,

$\displaystyle M$ $\displaystyle=$
$\displaystyle\left(\begin{array}[]{ccc}u\bar{u}&u\bar{d}&u\bar{s}\\\
d\bar{u}&d\bar{d}&d\bar{s}\\\ s\bar{u}&s\bar{d}&s\bar{s}\end{array}\right)$
(5) $\displaystyle\Rightarrow$
$\displaystyle\left(\begin{array}[]{ccc}\frac{\pi^{0}}{\sqrt{2}}+\frac{\eta}{\sqrt{3}}+\frac{\eta^{\prime}}{\sqrt{6}}&\pi^{+}&K^{+}\\\
\pi^{-}&-\frac{\pi^{0}}{\sqrt{2}}+\frac{\eta}{\sqrt{3}}+\frac{\eta^{\prime}}{\sqrt{6}}&K^{0}\\\
K^{-}&\bar{K}^{0}&-\frac{\eta}{\sqrt{3}}+\frac{2\eta^{\prime}}{\sqrt{6}}\end{array}\right).$
(9)

We ignore the $\eta^{\prime}$ component in this work, because the
$\eta^{\prime}$ has a large mass and does not play the role in the generation
of the $a_{0}(980)$. Now, we have the possible components of the final states
after the hadronization,

$\displaystyle H^{(a)}$ $\displaystyle=$ $\displaystyle
V_{cd}V_{ud}\left(M^{2}\right)_{22}p$ (11) $\displaystyle=$ $\displaystyle
V_{cd}V_{ud}\left(\pi^{+}\pi^{-}+\frac{1}{2}\pi^{0}\pi^{0}-\frac{2}{\sqrt{6}}\pi^{0}\eta\right.$
$\displaystyle\left.+\frac{1}{3}\eta\eta+K^{0}\bar{K}^{0}\right)p\,,$

where the coefficients account for the relative strengths of the different
components. It should be noted that $H^{(a)}$ is not the effective
Hamiltonian. Since the channels $\pi^{+}\pi^{-}$, $\pi^{0}\pi^{0}$, and
$\eta\eta$ only couple to the system of isospin $I=0$, and have no
contribution for the $a_{0}(980)$ production in the $\Lambda_{c}$ decay, we
eliminate the components of these three channels. Equation (11) can be
rewritten as,

$\displaystyle H^{(a)}$ $\displaystyle=$ $\displaystyle
V_{cd}V_{ud}\left(M^{2}\right)_{22}p$ (12) $\displaystyle=$ $\displaystyle
V_{cd}V_{ud}\left(-\frac{2}{\sqrt{6}}\pi^{0}\eta+K^{0}\bar{K}^{0}\right)p\,.$

For the mechanisms of Figs. 1(c) and 1(d), where the $\bar{u}u$ pair is
produced from the vacuum between $u$ from the $W^{+}$ decay and $ud$ quarks of
the initial $\Lambda_{c}$, we have,

$\displaystyle\Lambda_{c}$ $\displaystyle\Rightarrow$ $\displaystyle
V_{cd}V_{ud}d\bar{d}u\left(\bar{u}u\right)\frac{1}{\sqrt{2}}\left(ud-
du\right)$ (13) $\displaystyle=$ $\displaystyle V_{cd}V_{ud}\left(M_{22}\times
M_{11}\right)\frac{1}{\sqrt{2}}u\left(ud-du\right)$ $\displaystyle=$
$\displaystyle
V_{cd}V_{ud}\left(-\frac{\pi^{0}}{\sqrt{2}}+\frac{\eta}{\sqrt{3}}\right)\left(\frac{\pi^{0}}{\sqrt{2}}+\frac{\eta}{\sqrt{3}}\right)p$
$\displaystyle=$ $\displaystyle
V_{cd}V_{ud}\left(-\frac{1}{\sqrt{6}}\pi^{0}\eta+\frac{1}{\sqrt{6}}\pi^{0}\eta\right)p=0,$

for Fig. 1(c), and

$\displaystyle\Lambda_{c}$ $\displaystyle\Rightarrow$ $\displaystyle
V_{cs}V_{us}s\bar{s}u\left(\bar{u}u\right)\frac{1}{\sqrt{2}}\left(ud-
du\right)$ (14) $\displaystyle=$ $\displaystyle V_{cs}V_{us}\left(M_{33}\times
M_{11}\right)\frac{1}{\sqrt{2}}u\left(ud-du\right)$ $\displaystyle=$
$\displaystyle
V_{cs}V_{us}\left(-\frac{\eta}{\sqrt{3}}\right)\left(\frac{\pi^{0}}{\sqrt{2}}+\frac{\eta}{\sqrt{3}}\right)p$
$\displaystyle=$ $\displaystyle
V_{cs}V_{us}\left(-\frac{1}{\sqrt{6}}\pi^{0}\eta\right)p,$

for Fig. 1(d), respectively, where we have eliminated the components of
$\pi^{0}\pi^{0}$ and $\eta\eta$ that only couples to the system of isospin
$I=0$. One can find that the mechanism of Fig. 1(c) has no contribution to the
$\pi^{0}\eta p$ production. For Fig. 1(d), we have the components of the final
states,

$\displaystyle H^{(d)}$ $\displaystyle=$ $\displaystyle
V_{cs}V_{us}\left(-\frac{1}{\sqrt{6}}\pi^{0}\eta\right)p.$ (15)

Similar to the Figs. 1(c) and 1(d), the $\bar{q}q$ pair can also be created
between the $ud$ pair of the initial $\Lambda_{c}$, and hadronizes into the
proton ($uud$) and $\pi^{0}$ $\bar{d}d$. However, the $ud$ pair of the initial
$\Lambda_{c}$ is the most attractive “good” diquark, and this mechanism will
destruct the strong diquark correlation, which is not favored, as discussed in
Refs. Miyahara:2015cja ; Aitala:1999uq . Thus, we neglect the contributions
from those mechanisms because of the diquark correlation.

Although the final states $\pi^{0}\eta p$ can not be directly produced via the
$W^{+}$ external emission of the $\Lambda_{c}$ decay, the $\Lambda_{c}\to
pK^{+}K^{-}$ can proceed via the mechanism of the $W^{+}$ external emission as
depicted in Fig. 1(e), then undergoes the final state interaction by the
transition $K^{+}K^{-}\to\pi^{0}\eta$, which finally gives rise to the final
states $\pi^{0}\eta p$. For the process $\Lambda_{c}\to pK^{+}K^{-}$, we have
the possible component of the final states,

$\displaystyle H^{(e)}$ $\displaystyle=$ $\displaystyle
V_{cs}V_{us}K^{+}K^{-}p,$ (16)

where $V_{cs}$ and $V_{us}$ are the CKM matrix elements related to the Cabibbo
angle. We take $V_{cd}=V_{us}=-{\rm sin}\theta_{c}=-0.22534$ and
$V_{ud}=V_{cs}={\rm sin}\theta_{c}=0.97427$ Bayar:2014qha ; Wang:2020 ;
Duan:2020vye .

Figure 2: The mechanisms of the decay $\Lambda_{c}\to\pi^{0}\eta p$, (a) tree
diagram, (b) the $S$-wave final state interactions.

For the process $\Lambda_{c}\to\pi^{0}\eta p$, in addition to the tree diagram
of Fig. 2(a), we also take into account the $S$-wave pseudoscalar-pseudoscalar
interaction, which will dynamically generate the scalar $a_{0}(980)$, as
depicted in Fig. 2(b). Taking into the contributions from Figs. 1(a), 1(d),
and 1(e) The total amplitude for the $\Lambda_{c}\to\pi^{0}\eta p$ can be
expressed as,

$\displaystyle\mathcal{M}$ $\displaystyle=$ $\displaystyle
V_{p}V_{cd}V_{ud}\left(-\frac{3}{\sqrt{6}}-\frac{3}{\sqrt{6}}G_{\pi^{0}\eta}t_{\pi^{0}\eta\rightarrow\pi^{0}\eta}\right.$
$\displaystyle\left.+G_{K^{0}\bar{K}^{0}}t_{K^{0}\bar{K}^{0}\rightarrow\pi^{0}\eta}+C\times
G_{K^{+}K^{-}}t_{K^{+}K^{-}\rightarrow\pi^{0}\eta}\right),$

where $V_{p}$ is the vertex of meson-meson production, which contains the weak
amplitudes and the $q\bar{q}$ pair creation. As discussed in Ref. Geng:2018upx
where three-body charmed baryon decays are explicitly studied within the SU(3)
flavor symmetry, the squared amplitude of the $\Lambda_{c}\to\pi^{0}\eta p$
without resonance effects is almost structureless, and the decay amplitude of
the tree diagram is treated as a constant without energy dependence. Indeed,
this hypothesis is widely used in the weak decays of the heavy hadrons
Miyahara:2015cja ; Wang:2020 ; Xie:2018rqv ; Liang:2014tia ; Sekihara:2015iha
. Thus we take $V_{p}$ as a constant in here.

We also include the color factor $C$ to account for the relative weight of the
$W^{+}$ external emission mechanism with respect to the one of the $W^{+}$
internal emission mechanism, and the value of $C$ should be around 3 because
we take the number of colors $N_{c}=3$ Wang:2020 ; Dai:2018nmw ; Zhang:2020rqr
. 111 For the $W^{+}$ external emission of Fig. 1(e), the $u\bar{s}$ quarks
from the $W^{+}$ decay will hadronize into the color-singlet $K^{+}$, and the
$u$ and $\bar{s}$ quarks could have three choices of the colors ($N_{C}=3$)
because $W^{+}$ has no color. However, for the $W^{+}$ internal emission of
Fig. 1(a), the $d$ quark from the $c$ decay has the same color as the $c$
quark, and the $\bar{d}$ and $u$ quarks from the $W^{+}$ decay should have the
opposite and fixed colors, because the $u$ quark, together with the $ud$ of
the initial $\Lambda_{c}$, has to match to the color singlet proton. It
implies that the colors of all the final quarks are determined for the $W^{+}$
internal emission mechanism. Thus, the mechanism of the $W^{+}$ external
emission (Fig. 1(e)) is color-favored with respect to the ones of the $W^{+}$
internal emission (Figs. 1(a)-(d)), and we introduce the color factor $C$ for
the mechanism of the $W^{+}$ external emission. The loop function $G_{l}$ is a
diagonal matrix given by,

$\displaystyle G_{l}$ $\displaystyle=$ $\displaystyle
i\int\frac{d^{4}\vec{q}}{(2\pi)^{4}}\frac{1}{(P-q)^{2}-m_{1}^{2}+i\epsilon}\frac{1}{q^{2}-m^{2}_{2}+i\epsilon}$
(18) $\displaystyle=$ $\displaystyle\int_{0}^{q_{\rm
max}}\frac{|\vec{q}\,|^{2}d|\vec{q}\,|}{(2\pi)^{2}}\frac{\omega_{1}+\omega_{2}}{\omega_{1}\omega_{2}}\frac{1}{s-(\omega_{1}+\omega_{2})^{2}+i\epsilon},\,$

where $s$ is the invariant mass squared of the meson-meson system, and
$m_{1,2}$ are the masses of the mesons. The meson energies
$\omega_{1}=\sqrt{(\vec{q}\,)^{2}+m^{2}_{1}}$,
$\omega_{2}=\sqrt{(\vec{q}\,)^{2}+m^{2}_{2}}$. Since the $G_{l}$ is
logarithmically divergent, we use the cut-off method to solve this singular
integral, and choose the cut-off parameter ${q}_{\rm max}=600$ MeV, as Refs.
Wang:2020 ; Dias:2016gou ; Liang:2014tia ; Xie:2014tma ; Ahmed:2020qkv . The
transition amplitude $t_{ij}$ in Eq. (LABEL:eq:amp) is obtained by solving the
Bethe-Salpeter equation in coupled channels,

$T=[1-VG]^{-1}V,$ (19)

where three coupling channels $\pi^{0}\eta$, $K^{+}K^{-}$, and
$K^{0}\bar{K}^{0}$ are included. $V_{ij}$ represents the transition potential
from $i$-channel to $j$-channel, as Refs. Xie:2014tma ; Wang:2020 ,

$\displaystyle
V_{K^{+}K^{-}\to\pi^{0}\eta}=\frac{-\sqrt{3}}{12f^{2}}(3s-\frac{8}{3}m^{2}_{K}-\frac{1}{3}m^{2}_{\pi}-m^{2}_{\eta}),$
(20) $\displaystyle
V_{K^{0}\bar{K}^{0}\to\pi^{0}\eta}=-V_{K^{+}K^{-}\to\pi^{0}\eta},$ (21)
$\displaystyle V_{\pi^{0}\eta\to\pi^{0}\eta}=-\frac{1}{3f^{2}}m^{2}_{\pi},$
(22) $\displaystyle V_{K^{+}K^{-}\to K^{+}K^{-}}=-\frac{1}{2f^{2}}s,$ (23)
$\displaystyle V_{K^{+}K^{-}\to K^{0}\bar{K}^{0}}=-\frac{1}{4f^{2}}s,$ (24)
$\displaystyle V_{K^{0}\bar{K}^{0}\to K^{0}\bar{K}^{0}}=-\frac{1}{2f^{2}}s,$
(25)

with the decay constant $f=f_{\pi}=93$ MeV. $m_{\pi}$ and $m_{K}$ are the
isospin averaged masses of the pion and kaon, respectively.

With the full amplitude of Eq. (LABEL:eq:amp), we can write the $\pi^{0}\eta$
invariant mass distribution for the process $\Lambda_{c}\to\pi^{0}\eta p$,

$\frac{d\Gamma}{dM_{\rm\pi^{0}\eta}}=\frac{1}{(2\pi)^{3}}\frac{p_{p}\tilde{p}_{\pi^{0}}}{4M^{2}_{\Lambda_{c}}}|\mathcal{M}|^{2},$
(26)

where $p_{p}$ is the proton momentum in the $\Lambda_{c}$ rest frame,

$p_{p}=\frac{\lambda^{1/2}(M^{2}_{\Lambda_{c}},M^{2}_{p},M^{2}_{\pi^{0}\eta})}{2M_{\Lambda_{c}}},$
(27)

$\tilde{p}_{\pi^{0}}$ is the $\pi^{0}$ momentum in the $\pi^{0}\eta$ rest
frame.

$\displaystyle\tilde{p}_{\pi^{0}}=\frac{\lambda^{1/2}(M^{2}_{\pi^{0}\eta},m^{2}_{\pi^{0}},m^{2}_{\eta})}{2M_{\pi^{0}\eta}},$
(28)

with Källen function $\lambda(x,y,z)=x^{2}+y^{2}+z^{2}-2xy-2yz-2zx$. All the
masses of the mesons and baryons involved in our calculations are taken from
PDG PDG2018 .

For the process $\Lambda_{c}\to\pi^{0}\eta p$, the rescattering of $\pi^{0}p$
and $\eta p$ systems can give some contributions. Since there are many $N^{*}$
and $\Delta^{*}$ resonances, which could couple to the $\pi^{0}p$ and $\eta p$
PDG2018 , and most of their couplings are unknown, we expect that those
intermediate resonances will give a background contribution around 1 GeV in
the $\pi^{0}\eta$ invariant mass distribution, and not significantly change
the position of $a_{0}(980)$. Thus, we neglect the rescattering of the
$\pi^{0}p$ and $\eta p$ systems in this work, and will consider them when the
experimental data are available in future.

##  III Results and Discussion

Figure 3: The modulus squared of the amplitudes: (a)
$|t_{K^{+}K^{-}\rightarrow\pi^{0}\eta}|^{2}$, (b)
$|t_{K^{0}\bar{K}^{0}\rightarrow\pi^{0}\eta}|^{2}$, (c)
$|t_{\pi^{0}\eta\rightarrow\pi^{0}\eta}|^{2}$. Figure 4: The $M_{\pi^{0}\eta}$
invariant mass distribution of the process $\Lambda_{c}\to\pi^{0}\eta p$. The
purple dash-dotted curve shows the contribution from the $S$-wave
pseudoscalar-pseudoscalar interaction as depicted in Fig. 2(b), the blue
dotted curve stands for the the contribution from the tree level diagram as
depicted in Fig. 2(a), and the red solid curve corresponds to the total
results given by Eq. (LABEL:eq:amp). Figure 5: The $M_{\pi^{0}\eta}$ invariant
mass distribution for the process $\Lambda_{c}\to\pi^{0}\eta p$ with different
values of $C=$ 3, 2, -2, -3.

In this section, we will present our results. In order to show the production
of the $a_{0}(980)$ in the $S$-wave $K^{+}K^{-}$, $K^{0}\bar{K}^{0}$, and
$\pi^{0}\eta$ interactions, we respectively plot the modulus squared
$|t_{ij}|^{2}$ for the transitions $K^{+}K^{-}\rightarrow\pi^{0}\eta$,
$K^{0}\bar{K}^{0}\rightarrow\pi^{0}\eta$, and
$\pi^{0}\eta\rightarrow\pi^{0}\eta$, as shown in Fig. 3. From Fig. 3, one can
see a clear cusp structure around 980 MeV, which is associated to the
$a_{0}(980)$. The cusp structure manifested by the $a_{0}(980)$ is also
supported by the experimental measurements Debastiani:2017 .

In our model, we have two free parameters, $V_{p}$ and $C$. $V_{p}$ is a
global factor and its value does not affect the shape of the $\pi^{0}\eta$
invariant mass distribution. As mentioned above, the value of $C$ should be
around 3. In the first step, we take $C=3$, and show the $\pi^{0}\eta$
invariant mass distribution up to an arbitrary normalization in Fig. 4.
Instead of the cusp structure as shown in Fig. 3, one can find a dip structure
around 980 MeV in Fig. 4, which could be associated to the $a_{0}(980)$.
Although the $a_{0}(980)$, as a dynamically generated state from the $S$-wave
pseudoscalar-pseudoscalar, manifests as a cusp structure in many processes,
such as $\eta_{c}\to\eta\pi^{+}\pi^{-}$ Debastiani:2017 ; Ablikim:2017 , it
shows a dip structure in the process $\Lambda_{c}\to\pi^{0}\eta p$. As we
known, hadron resonances are observed as narrow or broad peaks in the
invariant mass distribution in many cases, a resonance may even show up as a
dip, depending on the interference between the different contributions, as
discussed in Ref. Guo:2019twa . Indeed, this behaviour is relatively common in
hadron physics. For example, the $f_{0}(980)$ manifests itself as a clear peak
in the $\pi^{+}\pi^{-}$ invariant mass distribution of the $J/\psi\to\phi
K^{+}K^{-}$ Wu:2001vz ; Augustin:1988ja , and $B_{s}\to J/\psi\pi^{+}\pi^{-}$
Aaij:2011fx  reactions, but shows up as a dip in the $S$-wave $\pi\pi$
scattering amplitudes Pelaez:2015qba . The dip structures have also been found
in experiments Frabetti:2001ah ; Aubert:2006jq  and theoretical studies
Lu:2016roh ; Wang:2015pcn ; Wang:2015qta ; liu:2017efp .

On the other hand, the $N_{c}$ scaling only indicates the relative strength of
the absolute values, and the relative sign of $C$ is not fixed Wang:2020 ;
Zhang:2020rqr , thus we show the $\pi^{0}\eta$ invariant mass distribution
with different values of $C=3$, $2$, $-2$, $-3$ in Fig. 5, where one can see a
significant cusp structure around 980 MeV for the negative $C$, and a small
dip structure for the positive $C$, which is due to the interference of the
contribution from the final state interaction and the contribution from the
tree diagram. It should be stressed that the positive $C$ is supported by many
studies, such as Refs. Wang:2020 ; Dai:2018nmw . Thus, one dip or cusp
structure in the $\pi^{0}\eta$ invariant mass distribution is expected to be
observed in the more precision measurements for the process
$\Lambda_{c}\to\pi^{0}\eta p$.

As we knows, in the energy region around 1 GeV, there is only the state
$a_{0}(980)$ that could decay into $\pi\eta$ PDG2018 . Thus if there is a
peak/dip structure observed in the $\pi\eta$ mass distribution around 1 GeV,
it should result from the $a_{0}(980)$.

Since the process $\Lambda_{c}\to\eta\pi^{0}p$ has not yet been measured
experimentally up to our knowledge, we will predict the branching fraction of
the process $\Lambda_{c}\to\pi^{0}\eta p$, which is important for the
measurements of this process in the future. The first step of the
$\Lambda_{c}\to\pi^{0}\eta p$, weak process, is the same as the ones of the
process $\Lambda_{c}\to p\pi^{+}\pi^{-}$ and $\Lambda_{c}\to pK^{+}K^{-}$
Wang:2020 , the normalization factor $V_{p}$, containing all dynamical factors
of the weak process, is expected to be identical for these three processes.
The branching fraction of the process $\Lambda_{c}\to\pi^{0}\eta p$ can be
expressed as,

$\displaystyle\mathcal{B}(\Lambda_{c}\to\pi^{0}\eta p)$ $\displaystyle=$
$\displaystyle\frac{\int\left(\frac{d\Gamma}{dM_{\pi^{0}\eta}}\right)dM_{\pi^{0}\eta}}{\Gamma_{\Lambda_{c}}},$
(29)

where the integrating range of the $M_{\pi^{0}\eta}$ is taken to be
$(m_{\pi}+m_{\eta})<M_{\pi^{0}\eta}<1200$ MeV, because the chiral unitary
approach can only make reliable predictions in the region up to $1200$ MeV.
With the $V^{2}_{p}/\Gamma_{\Lambda_{c}}=0.2$ MeV-1 from Ref. Wang:2020 , we
can roughly estimate the branching fraction
$\mathcal{B}(\Lambda_{c}\to\pi^{0}\eta p)=(1.13\sim 1.26)\times 10^{-4}$ for
$2<C<3$, and $(5.96\sim 8.19)\times 10^{-4}$ for $-3<C<-2$. The order of the
magnitude $10^{-4}$ is expected to be accessible in the BESIII and Belle II
experiments.

##  IV Conclusions

In this paper, the process $\Lambda_{c}\rightarrow\pi^{0}\eta p$ is
investigated within the chiral unitary approach, where the contributions from
the tree level diagram and the final state interactions of $K^{+}K^{-}$,
$K^{0}\bar{K}^{0}$, and $\pi^{0}\eta$ are taken into account.

In our model, the interactions of the meson-meson pairs $K^{+}K^{-}$,
$K^{0}\bar{K}^{0}$, and $\pi^{0}\eta$, preliminarily produced via the
$\Lambda_{c}$ weak decay, result in the intermediate $a_{0}(980)$. We
calculate the modulus square of the transition amplitudes for the
$K^{+}K^{-}\rightarrow\pi^{0}\eta$, $K^{0}\bar{K}^{0}\rightarrow\pi^{0}\eta$,
and $\pi^{0}\eta\rightarrow\pi^{0}\eta$ transitions, and find an obvious cusp
near $980$ MeV, which corresponds to the $a_{0}(980)$ resonance. Up to an
arbitrary normalization, we predict the $\pi^{0}\eta$ invariant mass
distribution, and find a dip structure for the positive values of $C$ and a
cusp structure for the negative values of $C$.

As so far, no measurements of the process $\Lambda_{c}\rightarrow\pi^{0}\eta
p$ were done, we roughly estimate the branching fraction
$\mathcal{B}(\Lambda_{c}\to\pi^{0}\eta p)=(1.13\sim 1.26)\times 10^{-4}$ for
the positive $C$, and $(5.96\sim 8.19)\times 10^{-4}$ for the negative $C$,
which are expected to be accessible in the BESIII and Belle II experiments.
Thus, we would like to call the attention of the experimentalists to measure
this reaction, which should be useful to understand the nature of the
$a_{0}(980)$, and the enhancement structure near the $K^{+}K^{-}$ threshold in
the process $\Lambda_{c}\to pK^{+}K^{-}$ observed by the BESIII Collaboration
BESIII:2018mes .

###### Acknowledgements.

This work is partly supported by the National Natural Science Foundation of
China under Grants No. 11505158. It is also supported by the Key Research
Projects of Henan Higher Education Institutions under No. 20A140027, the
Project of Youth Backbone Teachers of Colleges and Universities of Henan
Province (2020GGJS017), the Fundamental Research Cultivation Fund for Young
Teachers of Zhengzhou University (JC202041042), and the Academic Improvement
Project of Zhengzhou University.

## References

  * (1) H. Y. Cheng, Charmed baryons circa 2015, Front. Phys. (Beijing) 10, 101406 (2015).
  * (2) D. Ebert and W. Kallies, Nonleptonic Decays of Charmed Baryons in the MIT Bag Model, Phys. Lett. 131B, 183 (1983). Erratum: [Phys. Lett. 148B, 502 (1984)].
  * (3) C. D. Lü, W. Wang and F. S. Yu, Test flavor $SU(3)$ symmetry in exclusive $\Lambda_{c}$ decays, Phys. Rev. D 93, 056008 (2016).
  * (4) C. Q. Geng, Y. K. Hsiao, C. W. Liu and T. H. Tsai, Three-body charmed baryon Decays with SU(3) flavor symmetry, Phys. Rev. D 99, 073003 (2019).
  * (5) E. Oset et al. Weak decays of heavy hadrons into dynamically generated resonances, Int. J. Mod. Phys. E 25 (2016), 1630001.
  * (6) M. Ablikim et al. [BESIII Collaboration], Measurement of Singly Cabibbo Suppressed Decays $\Lambda_{c}\rightarrow p\pi^{+}\pi^{-}$ and $\Lambda_{c}\rightarrow pK^{+}K^{-}$, Phys. Rev. Lett. 117, 232002(2016). Addendum: [Phys. Rev. Lett. 120, 029903 (2018)].
  * (7) R. Aaij et al. [LHCb], Measurements of the branching fractions of $\Lambda_{c}^{+}\rightarrow p\pi^{-}\pi^{+}$, $\Lambda_{c}^{+}\rightarrow pK^{-}K^{+}$, and $\Lambda_{c}^{+}\rightarrow p\pi^{-}K^{+}$, JHEP 03 (2018), 043.
  * (8) M. Ablikim et al. [BESIII], Evidence for the singly-Cabibbo-suppressed decay $\Lambda_{c}^{+}\to p\eta$ and search for $\Lambda_{c}^{+}\to p\pi^{0}$, Phys. Rev. D 95 (2017) 111102.
  * (9) K. Abe et al. [Belle], Observation of Cabibbo suppressed and $W$ exchange $\Lambda^{+}_{c}$ baryon decays, Phys. Lett. B 524 (2002), 33-43.
  * (10) B. Pal et al. [Belle], Search for $\Lambda_{c}^{+}\to\phi p\pi^{0}$ and branching fraction measurement of $\Lambda_{c}^{+}\to K^{-}\pi^{+}p\pi^{0}$, Phys. Rev. D 96 (2017) 051102.
  * (11) M. Ablikim et al. [BESIII], Observation of $\Lambda^{+}_{c}\to nK^{0}_{S}\pi^{+}$, Phys. Rev. Lett. 118 (2017) 112001.
  * (12) M. Ablikim et al. [BESIII], Measurement of the absolute branching fractions of $\Lambda_{c}^{+}\to\Lambda\eta\pi^{+}$ and $\Sigma(1385)^{+}\eta$, Phys. Rev. D 99 (2019) 032010.
  * (13) C. Q. Geng, C. W. Liu and T. H. Tsai, Singly Cabibbo suppressed decays of $\Lambda_{c}^{+}$ with SU(3) flavor symmetry, Phys. Lett. B 790 (2019), 225-228.
  * (14) H. Y. Cheng, X. W. Kang and F. Xu, Singly Cabibbo-suppressed hadronic decays of $\Lambda_{c}^{+}$, Phys. Rev. D 97 (2018) 074028.
  * (15) R. Pavao, S. Sakai and E. Oset, Production of $N^{*}(1535)$ and $N^{*}(1650)$ in $\Lambda_{c}\rightarrow\bar{K}^{0}\eta p$ $(\pi N)$ decay,hys. Rev. C 98 (2018) 015201.
  * (16) J. J. Xie and L. S. Geng, Role of the $N^{*}(1535)$ in the $\Lambda^{+}_{c}\to\bar{K}^{0}\eta p$ decay, Phys. Rev. D 96, 054009 (2017).
  * (17) X. H. Liu, G. Li, J. J. Xie and Q. Zhao, Visible narrow cusp structure in $\Lambda_{c}^{+}\to pK^{-}\pi^{+}$ enhanced by triangle singularity, Phys. Rev. D 100 (2019) 054006.
  * (18) J. K. Ahn, S. Yang and S. I. Nam, Hyperon production in $\Lambda_{c}^{+}\to K^{-}p\pi^{+}$ and $\Lambda_{c}^{+}\to K^{0}_{S}p\pi^{0}$, Phys. Rev. D 100 (2019) 034027.
  * (19) J. J. Xie and L. S. Geng, The $a_{0}(980)$ and $\Lambda(1670)$ in the $\Lambda^{+}_{c}\to\pi^{+}\eta\Lambda$ decay, Eur. Phys. J. C 76 (2016) 496.
  * (20) J. J. Xie and L. S. Geng, $\Sigma^{*}_{1/2^{-}}(1380)$ in the $\Lambda^{+}_{c}\to\eta\pi^{+}\Lambda$ decay, Phys. Rev. D 95 (2017) 074024.
  * (21) Z. Wang, Y. Y. Wang, E. Wang, D. M. Li and J. J. Xie, The scalar $f_{0}(500)$ and $f_{0}(980)$ resonances and vector mesons in the single Cabibbo-suppressed decays $\Lambda_{c}\to pK^{+}K^{-}$ and $p\pi^{+}\pi^{-}$, Eur. Phys. J. C 80, 842 (2020).
  * (22) R. Aaij et al. [LHCb], First observation of $\bar{B}^{0}\to J/\psi K^{+}K^{-}$ and search for $\bar{B}^{0}\to J/\psi\phi$ decays, Phys. Rev. D 88 (2013) , 072005.
  * (23) W. H. Liang, J. J. Xie and E. Oset, ${\bar{B}}^{0}$ , $B^{-}$ and ${\bar{B}}^{0}_{s}$ decays into $J/\psi$ and $K{\bar{K}}$ or $\pi\eta$, Eur. Phys. J. C 75, 609 (2015).
  * (24) E. Gutz et al. [CBELSA/TAPS], High statistics study of the reaction $\gamma p\to p\pi^{0}\eta$, Eur. Phys. J. A 50 (2014), 74.
  * (25) P. A. Zyla et al. [Particle Data Group], Review of Particle Physics, PTEP 2020 (2020) 083C01.
  * (26) G. Y. Wang, S. C. Xue, G. N. Li, E. Wang and D. M. Li, Strong decays of the higher isovector scalar mesons, Phys. Rev. D 97, 034030 (2018).
  * (27) J. A. Oller and E. Oset, Chiral symmetry amplitudes in the S wave isoscalar and isovector channels and the $\sigma$, $f_{0}$(980), $a_{0}$(980) scalar mesons, Nucl. Phys. A 620 (1997), 438-456.
  * (28) J. Nieves and E. Ruiz Arriola, Bethe-Salpeter approach for meson meson scattering in chiral perturbation theory, Phys. Lett. B 455 (1999), 30-38.
  * (29) D. Cabrera, E. Oset and M. J. Vicente Vacas, Evaluation of the $\pi\pi$ scattering amplitude in the sigma-channel at finite density, Phys. Rev. C 72 (2005), 025207.
  * (30) F. K. Guo, P. N. Shen, H. C. Chiang, R. G. Ping and B. S. Zou, Dynamically generated $0^{+}$ heavy mesons in a heavy chiral unitary approach, Phys. Lett. B 641 (2006), 278-285.
  * (31) W. H. Liang, J. J. Xie and E. Oset, $\bar{B}^{0}$ decay into $D^{0}$ and $f_{0}(500)$, $f_{0}(980)$, $a_{0}(980)$, $\rho$ and $\bar{B}^{0}_{s}$ decay into $D^{0}$ and $\kappa(800)$, $K^{*0}$, Phys. Rev. D 92, 034008 (2015).
  * (32) J. J. Xie and G. Li, The decays of $\bar{B}^{0}$ , $\bar{B}^{0}_{s}$ and $B^{-}$ into $\eta_{c}$ plus a scalar or vector meson, Eur. Phys. J. C 78, 861(2018).
  * (33) W. H. Liang, J. J. Xie and E. Oset, $f_{0}(500)$ , $f_{0}(980)$ , and $a_{0}(980)$ production in the $\chi_{c1}\rightarrow\eta\pi^{+}\pi^{-}$ reaction, Eur. Phys. J. C 76, 700(2016).
  * (34) V. R. Debastiani, W. H. Liang, J. J. Xie and E. Oset, Predictions for $\eta_{c}\to\eta\pi^{+}\pi^{-}$ producing $f_{0}(500)$, $f_{0}(980)$ and $a_{0}(980)$, Phys. Lett. B 766 (2017), 59-64.
  * (35) J. J. Xie, L. R. Dai and E. Oset, The low lying scalar resonances in the $D^{0}$ decays into $K^{0}_{s}$ and $f_{0}(500)$, $f_{0}(980)$, $a_{0}(980)$, Phys. Lett. B 742, 363(2015).
  * (36) M. Y. Duan, J. Y. Wang, G. Y. Wang, E. Wang and D. M. Li, Role of scalar $a_{0}(980)$ in the single Cabibbo suppressed process $D^{+}\rightarrow\pi^{+}\pi^{0}\eta$, Eur. Phys. J. C 80 (2020), 1041.
  * (37) N. Ikeno, M. Bayar and E. Oset, Combined theoretical study of the $D^{+}\to\pi^{+}\eta\eta$ and $D^{+}\to\pi^{+}\pi^{0}\eta$ reactions, [arXiv:2102.01650 [hep-ph]].
  * (38) G. Toledo, N. Ikeno and E. Oset, Theoretical study of the $D^{0}\to K^{-}\pi^{+}\eta$ reaction, [arXiv:2008.11312 [hep-ph]].
  * (39) R. Molina, J. J. Xie, W. H. Liang, L. S. Geng and E. Oset, Theoretical interpretation of the $D^{+}_{s}\to\pi^{+}\pi^{0}\eta$ decay and the nature of $a_{0}(980)$, Phys. Lett. B 803 (2020), 135279.
  * (40) Y. K. Hsiao, Y. Yu and B. C. Ke, Resonant $a_{0}(980)$ state in triangle rescattering $D_{s}^{+}\rightarrow\pi^{+}\pi^{0}\eta$ decays, Eur. Phys. J. C 80, 895 (2020).
  * (41) X. Z. Ling, M. Z. Liu, J. X. Lu, L. S. Geng and J. J. Xie, Can the nature of $a_{0}(980)$ be tested in the $D_{s}^{+}\to\pi^{+}\pi^{0}\eta$ decay?, [arXiv:2102.05349 [hep-ph]].
  * (42) L. R. Dai, Q. X. Yu and E. Oset, Triangle singularity in $\tau^{-}\to\nu_{\tau}\pi^{-}f_{0}(980)$ ($a_{0}(980)$) decays, Phys. Rev. D 99, 016021(2019).
  * (43) W. H. Liang, H. X. Chen, E. Oset and E. Wang, Triangle singularity in the $J/\psi\rightarrow K^{+}K^{-}f_{0}(980)a_{0}(980)$ decays, Eur. Phys. J. C 79, 411(2019).
  * (44) C. W. Xiao, U. G. Meißner and J. A. Oller, Investigation of $J/\psi\to\gamma\,\pi^{0}\eta(\pi^{+}\pi^{-},\pi^{0}\pi^{0})$ radiative decays including final-state interactions, Eur. Phys. J. A 56 (2020), 23.
  * (45) E. Wang, H. X. Chen, L. S. Geng, D. M. Li and E. Oset, Hidden-charm pentaquark state in $\Lambda^{0}_{b}\to J/\psi p\pi^{-}$ decay, Phys. Rev. D 93, 094001 (2016).
  * (46) J. X. Lu, E. Wang, J. J. Xie, L. S. Geng and E. Oset, The $\Lambda_{b}\rightarrow J/\psi K^{0}\Lambda$ reaction and a hidden-charm pentaquark state with strangeness, Phys. Rev. D 93, 094009 (2016).
  * (47) L. R. Dai, G. Y. Wang, X. Chen, E. Wang, E. Oset and D. M. Li, The $B^{+}\rightarrow J/\psi\omega K^{+}$ reaction and $D^{\ast}\bar{D}^{\ast}$ molecular states, Eur. Phys. J. A 55, 36 (2019).
  * (48) M. Bayar, W. H. Liang and E. Oset, $B^{0}$ and $B^{0}_{s}$ decays into $J/\psi$ plus a scalar or vector meson, Phys. Rev. D 90 (2014) 114004.
  * (49) L. Roca, J. E. Palomar and E. Oset, Decay of axial vector mesons into $VP$ and $P\gamma$, Phys. Rev. D 70 (2004), 094006.
  * (50) K. Miyahara, T. Hyodo and E. Oset, Weak decay of $\Lambda_{c}^{+}$ for the study of $\Lambda$(1405) and $\Lambda$(1670), Phys. Rev. C 92 (2015) no.5, 055204.
  * (51) E. M. Aitala et al. [E791], Multidimensional resonance analysis of $\Lambda^{+}_{c}\to pK^{-}\pi^{+}$, Phys. Lett. B 471 (2000), 449-459.
  * (52) W. H. Liang and E. Oset, $B^{0}$ and $B^{0}_{s}$ decays into $J/\psi$ $f_{0}(980)$ and $J/\psi$ $f_{0}(500)$ and the nature of the scalar resonances, Phys. Lett. B 737 (2014), 70-74.
  * (53) T. Sekihara and E. Oset, Investigating the nature of light scalar mesons with semileptonic decays of $D$ mesons, Phys. Rev. D 92 (2015) no.5, 054038.
  * (54) Y. Zhang, E. Wang, D. M. Li and Y. X. Li, Search for the $D^{*}\bar{D}^{*}$ molecular state $Z_{c}(4000)$ in the reaction $B^{-}\rightarrow J/\psi\rho^{0}K^{-}$, Chin. Phys. C 44, 093107 (2020).
  * (55) J. M. Dias, F. S. Navarra, M. Nielsen and E. Oset, $f_{0}(980)$ production in $D_{s}^{+}\rightarrow\pi^{+}\pi^{+}\pi^{-}$ and $D_{s}^{+}\rightarrow\pi^{+}K^{+}K^{-}$ decays, Phys. Rev. D 94, 096002(2016).
  * (56) H. A. Ahmed, Z. Y. Wang, Z. F. Sun and C. W. Xiao, Study $B^{0}_{(s)}$ decays into $\phi$ and a scalar or vector meson, [arXiv:2011.08758 [hep-ph]].
  * (57) M. Ablikim et al. [BESIII], Amplitude analysis of the $\chi_{c1}\to\eta\pi^{+}\pi^{-}$ decays, Phys. Rev. D 95 (2017) 032002.
  * (58) F. K. Guo, X. H. Liu and S. Sakai, Threshold cusps and triangle singularities in hadronic reactions, Prog. Part. Nucl. Phys. 112, 103757 (2020).
  * (59) N. Wu, BES R measurements and $J/\psi$ decays, [arXiv:hep-ex/0104050 [hep-ex]].
  * (60) J. E. Augustin et al. [DM2], Study of the $J/\psi$ Decay Into Five Pions, Nucl. Phys. B 320, 1-19 (1989)
  * (61) R. Aaij et al. [LHCb], First observation of $B^{0}_{s}\to J/\psi f_{0}(980)$ decays, Phys. Lett. B 698, 115-122 (2011)
  * (62) J. R. Pelaez, From controversy to precision on the sigma meson: a review on the status of the non-ordinary $f_{0}(500)$ resonance, Phys. Rept. 658, 1 (2016).
  * (63) P. L. Frabetti et al. [E687], Evidence for a Narrow Dip Structure at 1.9-GeV/$c^{2}$ in 3$\pi^{+}3\pi^{-}$ Diffractive photoproduction, Phys. Lett. B 514 (2001), 240-246.
  * (64) B. Aubert et al. [BaBar], The $e^{+}e^{-}\to 3(\pi^{+}\pi^{-})$, $2(\pi^{+}\pi^{-}\pi^{0})$ and $K^{+}K^{-}2(\pi^{+}\pi^{-})$ cross sections at center-of-mass energies from production threshold to 4.5-GeV measured with initial-state radiation, Phys. Rev. D 73 (2006), 052003.
  * (65) E. Wang, J. J. Xie and E. Oset, $\chi_{c0}(1P)$ decay into $\bar{\Sigma}~{}\Sigma\pi$ in search of an $I=1$, $1/2^{-}$ baryon state around $\bar{K}N$ threshold, Phys. Lett. B 753, 526-532 (2016).
  * (66) L. J. Liu, E. Wang, J. J. Xie, K. L. Song and J. Y. Zhu, $\Lambda(1405)$ production in the process $\chi_{c0}(1P)\to\bar{\Lambda}\Sigma\pi$, Phys. Rev. D 98, 114017 (2018).


# Fast Tracking of 6D Particle Phase Space Using the FC2CT Algorithm for rf
Cavities

M. Southerby<EMAIL_ADDRESS>R. Apsimon
<EMAIL_ADDRESS>Engineering Department, Lancaster University,
Lancaster, LA1 4YW, UK  Cockcroft Institute, Daresbury Laboratory, Warrington,
WA4 4AD, UK

###### Abstract

In this paper, a fast cell to cell tracking algorithm (FC2CT) is developed,
that determines the change in 6D phase space of a particle beam through rf
accelerating cavities operated in both Standing and Traveling Wave modes. The
performance is compared to well trusted tracking codes - ASTRA and RF-Track -
for proton beams with relativistic beta $\approx$ 0.5. The FC2CT algorithm
produces a complete analytical method to determine the change in particle
momentum, with the constant particle velocity approximation over the
integration region. Additional developments for FC2CT are discussed, such as
altering the integration length. FC2CT can be used as a fast method class in
common tracking codes to perform fast simulations through rf cavities.

##  I Introduction

The use of particle accelerators are paramount in multiple industries,
examples including the health care, cargo scanning and research industries. As
a result, there is constant requirement for additional developments within the
particle accelerator industry, ensuring optimal efficiency of machines or
increased beam energies. A large subset of accelerators are linear
accelerators (linacs) that accelerate particles using radio-frequency (rf)
cavities that contain electromagnetic (em) fields for particle acceleration.
Beam dynamics simulations are a vital aspect of accelerator design, however
can be computationally expensive, for example beam dynamics studies for
circular machines requiring many thousands of turns.

In this paper, we outline a fast tracking algorithm for beam dynamics through
coupled cavity linacs (CCL), and is called the Fast Cell to Cell Tracking
algorithm (FC2CT), calculating the change in 6D phase space of a particle over
one rf cell per iteration. The model approximates the particle velocity is
constant over an rf cell, and space charge effects are ignored. The first
approximation is accurate for protons of all energy that are to be accelerated
by a CCL. The proton transition energy for more efficient acceleration by CCL
is often taken to be $\approx$ 80-100 MeV, below which drift tube linacs may
have increased shunt impedance shunt_impedance_diff_strucs . Nevertheless,
some CCL accelerate protons from as low as 37.5 MeV light_linac . FC2CT is
only accurate for electrons over $\approx$ 3 MeV. The tracking algorithm
requires the on axis longitudinal electric field, from which the remaining
non-zero field components ($E_{r},B_{\theta}$) are approximated by determining
the Taylor expansion coefficients in the small radial limit. This method has
been utilised before, in order to quickly provide a method for 3D particle
tracking in well-known tracking codes, such as ASTRA astra . The method is
accurate for both Standing Wave (SW) and Traveling wave (TW) cavity types,
however require slightly different algorithms. In addition, the algorithm
described is purely for the TM010 cavity mode. The accuracy of the FC2CT model
is compared with well trusted tracking codes, ASTRA (SW) and RF-Track (TW),
rftrack .

###  I.1 Lorentz Force law and Maxwell Equations for a TM010 Mode

A particle in a region of non-zero em field will experience a force given by
the Lorentz force griffiths . The force is the rate of change of momentum;

$\Delta\vec{P}=\int
d\vec{P}=q\int\left(\vec{E}(t)+\vec{v}(t)\times\vec{B}(t)\right)dt.$ (1)

Where $\vec{v}(t)$ is the particle velocity. In general, the em field and
particle velocity are a function of time. The em fields in a given region will
be described by Maxwell equations griffiths . For an rf single cell, the TM010
mode is often used for acceleration wangler . The longitudinal electric field
($E_{z}$) and azimuthal magnetic field ($B_{\theta}$) in a TM010 mode for a
simple pillbox cavity are described by the zeroth and first order Bessel
functions of the first kind bessel . rf single cells are often designed with
nose-cones, to concentrate $E_{z}$ on-axis and improve shunt impedance and
require a beam pipe, to allow power coupling and particle transmission. As a
result, the initially separable em field components in the case of a pillbox
cavity become functions of the longitudinal displacement, $z$.

As a particle traversing an rf single cell will be at a small radial
displacement relative to the cell radius, it is possible to describe the non-
zero em field components as a function of the on-axis $E_{z}$ component,
denoted $E_{z}(r=0,t)$, using a Taylor expansion about $r$ = 0. The TM010 mode
in a given rf cell will have negligible $E_{\theta}$, $B_{r}$, $B_{z}$
components, as in a pillbox cavity.

The non-zero em field components observed by a particle at position, $z$,
time, $t$ can be described as follows,

$E_{z}(r,z,t)=\\\
\left[E_{z}(0,z)-\frac{r^{2}}{4}\left(\frac{\partial^{2}E_{z}}{\partial
z^{2}}+\frac{\omega^{2}}{c^{2}}E_{z}(0,z)\right)\right]\textrm{cos}\left(\phi_{0}+\omega
t\right),$ (2)

$E_{r}(r,z,t)=\\\ \left[-\frac{r}{2}\frac{\partial E_{z}}{\partial
z}+\frac{r^{3}}{16}\left(\frac{\partial^{3}E_{z}}{\partial
z^{3}}+\frac{\omega^{2}}{c^{2}}\frac{\partial E_{z}}{\partial
z}\right)\right]\textrm{cos}\left(\phi_{0}+\omega t\right),$ (3)

$B_{\theta}(r,z,t)=\\\
\frac{\omega}{c^{2}}\left[\frac{r}{2}E_{z}(0,z)-\frac{r^{3}}{16}\left(\frac{\partial^{2}E_{z}}{\partial
z^{2}}+\frac{\omega^{2}}{c^{2}}E_{z}(0,z)\right)\right]\\\
\times\textrm{cos}\left(\phi_{0}-\frac{\pi}{2}+\omega t\right).$ (4)

Where $\phi_{0}$ is the rf phase when the particle enters the cell (at $z=0$).
The magnetic field lags the electric field by $\pi/2$. $E_{z}(0,z)$ is the
longitudinal on-axis electric field component, and can be exported from em
solver codes, such as CST microwave solver CST  for a given rf single cell
geometry. The terms $\frac{\partial^{n}E_{z}}{\partial z^{n}}$ describe the
$n_{th}$ derivative of $E_{z}(0,z)$.

It is required to describe the on-axis $E_{z}$ component as a Fourier series
fourier , as $E_{z}$ is periodic over some period, $P$, defined by the phase
advance per cell of the cavity. As $E_{z}(0,z)$ can be described by a Fourier
series, the partial derivatives of $E_{z}$ become straight-forward to compute.
The Fourier terms required to construct $E_{z}$ for a SW cell is shown in Fig.
1.

(a) Figure 1: Normalised on-axis $E_{z}$ component, G(s) and the individual
Fourier components.

For a SW cavity, the longitudinal electric field, and the derivatives, are
given below at time $t$;

$E_{z}(r=0,z,t)=\sum_{n=1}^{N}b_{n}\hskip
2.84526pt\textrm{sin}\left(\frac{2n\pi z}{P}\right)\cos\left(\omega
t+\phi_{0}\right),,$ (5) $\frac{\partial E_{z}}{\partial
z}=\left(\frac{2\pi}{P}\right)\sum_{n=1}^{N}nb_{n}\hskip
2.84526pt\textrm{cos}\left(\frac{2n\pi z}{P}\right)\cos\left(\omega
t+\phi_{0}\right)$ (6) $\frac{\partial^{2}E_{z}}{\partial
z^{2}}=-\left(\frac{2\pi}{P}\right)^{2}\sum_{n=1}^{N}n^{2}b_{n}\hskip
2.84526pt\textrm{sin}\left(\frac{2n\pi z}{P}\right)\cos\left(\omega
t+\phi_{0}\right).$ (7)

Where $b_{n}$ are the Fourier series coefficients. The required number of
coefficients varies for different field profiles, and should be the minimum
number such that the field is accurately represented.

When the operating mode is a $\frac{\pi}{2}$-mode, the field period is two
cell lengths, $P=2L_{cell}$. Thus, $E_{z}$ can be described using only $b_{n}$
terms, as $a_{n}=0$ for all $n$. For TW cavities, the longitudinal electric
field has both real ($A(z)$) and imaginary ($B(z)$) components, and thus both
components (over one period) must be imported from an em field solver. Both
the real and imaginary fields can also be described as a Fourier series,
however will require both $a_{n}$ and $b_{n}$ coefficients. The field seen by
a particle in a TW cavity is the real component of the complex field, and is
given;

$E_{z}(r=0,z,t)=A(z)\cos\left(\omega t+\phi_{0}\right)-B(z)\sin\left(\omega
t+\phi_{0}\right).$ (8)

Where A(z) and B(z) are described by Fourier series;

$A(z)=\sum_{n=1}^{N}b_{n}^{\text{real}}\sin\left(\frac{2n\pi
z}{P}\right)+a_{n}^{\text{real}}\cos\left(\frac{2n\pi z}{P}\right),$ (9)

and

$B(z)=\sum_{n=1}^{N}b_{n}^{\text{imag}}\sin\left(\frac{2n\pi
z}{P}\right)+a_{n}^{\text{imag}}\cos\left(\frac{2n\pi z}{P}\right).$ (10)

From Eqn. 1, the momentum change in the $z$ direction is given;

$\Delta P_{z}=q\int(E_{z}(r,z,t)+v_{x}B_{y}(r,z,t)-v_{y}B_{x}(r,z,t))dt.$ (11)

As the particle velocity over the rf cell ($v_{z}$) is approximately constant,
$dt=\frac{dz}{v_{z}}$. It is also assumed that terms including
$\frac{v_{x}}{v_{z}}$ are negligible, thus;

$\Delta P_{z}=\frac{q}{\beta_{z}c}\int E_{z}(r,z,t)dz.$ (12)

The momentum change in the transverse planes are also given by the Lorentz
force, and are a function of the particles azimuthal angle, $\theta$, depicted
in Fig. 2.

(a) Figure 2: Schematic showing azimuthal angle, $\theta$.

$\Delta P_{x}=qc\cos\left(\theta\right)\left(\frac{1}{\beta_{z}c}\int
E_{r}dz+\int B_{\theta}dz\right)$ (13) $\Delta
P_{y}=qc\sin\left(\theta\right)\left(\frac{1}{\beta_{z}c}\int E_{r}dz+\int
B_{\theta}dz\right).$ (14)

###  I.2 Derivation of single cell SW

Combining Eqns. 2 and 12 the change in $P_{z}$ is given;

$\Delta P_{z}=\\\
\frac{q}{\beta_{z}c}\int_{0}^{L_{cell}}\left[E_{z}(0,z)-\frac{r^{2}}{4}\left(\frac{\partial^{2}E_{z}}{\partial
z^{2}}+\frac{\omega^{2}}{c^{2}}E_{z}(0,z)\right)\right]\\\
\times\textrm{cos}\left(\phi_{0}+\frac{\omega z}{\beta_{z}c}\right)dz,$ (15)

where the integration is defined over a single rf cell (the integration range
can be over multiple cells, or sections of a single rf cell). Equation 15 can
be simplified by noting all terms are of the form;

$F(b_{n},i,\beta_{s}/\beta_{z},\phi_{0})=\\\
\int_{0}^{L_{cell}}\frac{\partial^{i}}{\partial
z^{i}}[E_{z}]\cos\left(\phi_{0}+\frac{\omega z}{\beta_{z}c}\right)dz.$ (16)

Therefore, using the standard result of Eqn. 16, Eqn. 15 is simplified to

$\Delta P_{z}=\\\
\frac{q}{\beta_{z}c}\left(1-\frac{\omega^{2}}{c^{2}}\frac{r^{2}}{4}\right)F(b_{n},0,\beta_{s}/\beta_{z},\phi_{0})-\\\
\frac{q}{\beta_{z}c}\frac{r^{2}}{4}F(b_{n},2,\beta_{s}/\beta_{z},\phi_{0}).$
(17)

Where

$\begin{split}F\biggl{(}b_{n},{\stackanchor{\emph{i}even}{\emph{i}odd}},\beta_{s}/\beta_{z},\phi_{0}\biggr{)}=\frac{u(i)}{2}\biggl{(}\frac{\pi}{L_{cell}}\biggr{)}^{(i-1)}\times\\\
\sum_{n=1}^{\infty}n^{i}b_{n}\Biggl{[}\frac{{\stackanchor{cos}{sin}}\biggl{(}\frac{\beta_{s}}{\beta_{z}}\pi-n\pi+\phi_{0}\biggr{)}-{\stackanchor{cos}{sin}}\bigl{(}\phi_{0}\bigr{)}}{(\frac{\beta_{s}}{\beta_{z}}-n)}{\mp}\\\
\frac{{\stackanchor{cos}{sin}}\biggl{(}\frac{\beta_{s}}{\beta_{z}}\pi+n\pi+\phi_{0}\biggr{)}-{\stackanchor{cos}{sin}}\bigl{(}\phi_{0}\bigr{)}}{(\frac{\beta_{s}}{\beta_{z}}+n)}\Biggr{]}\end{split}$
(18)

and

$u(i)=\begin{cases}+1,&\text{if}\ i=0,1\\\ -1,&\text{if}\ i=2,3.\end{cases}$
(19)

Calculating similar forms for Eqns. 13 and 14 can be completed via
substitution of $E_{r}$ and $B_{\theta}$ as functions of purely $E_{z}$ and
its derivatives, using Eqns. 3 and 4. Finally, the standard result of Eqn. 16
can be used. For the integral of $B_{\theta}$, the input rf phase, $\phi_{0}$
becomes $\phi_{0}-\frac{\pi}{2}$.

###  I.3 Derivation of single cell TW

The derivation of the momentum gain for TW cavities is more involved, as both
the real and imaginary field maps of $E_{z}(r=0)$ are required. This is
because the field behaves as described in Eqn. 8. The longitudinal momentum
change, $\Delta P_{z}$, is therefore found by computing the following
integral;

$\Delta P_{z}=\\\ \frac{q}{\beta_{z}c}\int\Bigl{[}A(z)\cos\left(\frac{\omega
z}{\beta_{z}c}+\phi_{0}\right)-B(z)\sin\left(\frac{\omega
z}{\beta_{z}c}+\phi_{0}\right)\Bigr{]}-\\\
\frac{r^{2}}{4}\Biggl{(}\Bigl{[}\frac{\partial^{2}A(z)}{\partial
z^{2}}\cos\left(\frac{\omega
z}{\beta_{z}c}+\phi_{0}\right)-\frac{\partial^{2}B(z)}{\partial
z^{2}}\sin\left(\frac{\omega z}{\beta_{z}c}+\phi_{0}\right)\Bigr{]}+\\\
\frac{\omega^{2}}{c^{2}}\Bigl{[}A(z)\cos\left(\frac{\omega
z}{\beta_{z}c}+\phi_{0}\right)-B(z)\sin\left(\frac{\omega
z}{\beta_{z}c}+\phi_{0}\right)\Bigr{]}\Biggr{)}dz.$ (20)

Using Eqn. 16, $\Delta P_{z}$ can be calculated for the TW case, noting that
$\sin(x)=\cos(x-\frac{\pi}{2})$. In general, an iteration over a TW single
cell requires four times the number of calculations than a SW cell, due to
requiring both $a_{n}$ and $b_{n}$ coefficients for describing both the real
and imaginary TW field components, $A(z)$ and $B(z)$. As a result, the FC2CT
method is less effective for TW cavities, as the speed benefit is reduced.

As before, functional forms for $\Delta P_{x}$ and $\Delta P_{y}$ can be
determined using Eqns. 3, 4, 16 and 18.

###  I.4 Calculation of Position in SW and TW cavities

The 3D particle position is also updated every rf cell. Firstly, the time of
flight, $t_{cell}$, is calculated as follows;

$t_{cell}=\overline{\beta_{z}}L_{cell},$ (21)

where the average particle velocity along the $i_{th}$ axis is defined

$\overline{\beta_{i}}=\frac{\overline{P_{i}}}{\sqrt{m^{2}+\overline{P_{x}^{2}}+\overline{P_{y}^{2}}+\overline{P_{z}^{2}}}}.$
(22)

The updated positions are therefore;

$z=z+L_{cell},\hskip 5.69054ptx=x+{t_{cell}}c\overline{\beta_{x}},\hskip
5.69054pty=y+{t_{cell}}c\overline{\beta_{y}}.$ (23)

For both the SW and TW methods, the rf phase, $\phi_{0}$, is updated every
cell due to phase slippage. The phase slippage is a function of the difference
between particle and design velocity;

$\Delta\phi=\omega\Delta t\approx\frac{\omega
L_{cell}}{\beta_{s}c}\left(1-\frac{\beta_{s}}{\overline{\beta_{z}}}\right).$
(24)

A script was written in Python python , that takes the on-axis $E_{z}$ field
map as an input, and extracts the Fourier coefficients. The script tracks a
particle beam through an rf cavity with a specified number of the cells, with
a specified gradient. The 6D phase space of each particle is updated after
every cell, as is the rf phase. The position of the particle is increased by
$L_{cell}$ in the longitudinal direction and determined transversely using a
calculated time of flight and transverse velocity. This allows a single
iteration to be performed each rf cell, as opposed to multiple iterations, as
used in conventional tracking methods. Section II describes the accuracy of
the FC2CT method for both SW for protons with kinetic energy 37.5 and 150 MeV,
and TW structures for protons at 150 MeV.

##  II Performance

For the comparison of SW FC2CT, ASTRA tracking code was used. The input
distribution was created using the ASTRA generator function, producing protons
at 37.5 and 150 MeV. Protons at 37.5 MeV are simulated as below this energy
region, proton acceleration is often accomplished with RFQ and drift tube
linac wangler , as the effective shunt impedance is higher than for CCL. Thus,
exploring both 37.5 and 150 MeV proton energies assesses FC2CT for different
values of relativistic Lorentz factor, $\gamma_{r}$. In addition, showing
FC2CT is accurate for proton energies of 150 MeV, demonstrates the method will
be accurate for higher energy protons (as the only approximation is constant
relativistic $\beta$ over a cell).

(a) (b)

Figure 3: Transverse and longitudinal phase space of proton beam (initial
energy 37.5 MeV) after traversing 20 SW cells, as calculated by ASTRA and
FC2CT.

Side-coupled single cells were designed in CST Microwave Solver CST . The cell
lengths were defined such as to keep the field synchronous with protons of
energy of 37.5 and 150 MeV. Figure 3 shows the transverse and longitudinal
phase space of a proton beam after 20 SW cells. The tracked phase space is
shown as calculated by both ASTRA and FC2CT. The two distributions agree
strongly, suggesting the constant $\beta$ approximation is valid for protons
at 37.5 MeV.

(a) (b)

Figure 4: Transverse and longitudinal phase space of proton beam (initial
energy 150 MeV) after traversing 50 SW cells, as calculated by ASTRA and
FC2CT.

Figure 4 displays the output distributions for an initial proton energy of 150
MeV, after 50 SW cells. The two distributions agree well. The largest
discrepancy between distributions occurs on the ‘neck’ of the phase space
plot. Potential reasons for this are discussed later.

(a) Figure 5: Displays the difference in calculated $P_{z}$ as a function of
difference in calculated $z$ (defined as $X_{ASTRA}-X_{FC2CT}$), after 40 SW
cells. The initial proton energy is 150 MeV.

The difference between phase space components after a 40 cell SW structure, as
calculated by ASTRA and FC2CT, are plotted in Fig. 5. The plot displays the
maximum difference in $P_{z}$ occurs at maximum (positive) difference in $z$.
In this region, FC2CT underestimates both longitudinal particle momentum and
position. This result is intuitive, as lower momentum particles will travel
shorter distances. The maximum longitudinal momentum error is 1 MeV/c,
resulting in a maximum percentage error of approximately 0.016 $\%$ in
$P_{z}$. The ASTRA time step was set to 0.001 ns, resulting in approximately
150 iterations per cell.

(a) (b)

Figure 6: Transverse and longitudinal phase space of proton beam (initial
energy 150 MeV) after traversing 45 TW cells, as calculated by RF-Track and
FC2CT.

Particles in this region reside on the ‘neck’ of the phase space plot, see
Fig. 4 (b). These particles begin with gentle acceleration, due to a slightly
accelerating initial phase. As the particles begin to phase slip, the phase
tends to a zero-crossing phase. In the following cells, particles will
experience very small changes in $P_{z}$, and particles become physically
spread out, as bunching occurs. The input distributions of Figures 3 and 4
were uniform in rf phase, however particles in the ‘neck’ are spread out,
implying thus small changes in initial conditions (rf phase) result in large
changes in output conditions. As a result, FC2CT will perform less accurately
on particles in this region.

Very small discrepancies in calculated $P_{z}$ occurs at maximum (negative)
error in $z$. In this region, particle position is overestimated by FC2CT.
However, the discrepancy in $z$ is very small, with a maximum percentage error
of approximately 0.04$\%$. Particles with maximum and minimum $P_{z}$ are most
accurately tracked using FC2CT. This is due to the particles residing at an rf
phase that remains far from the zero-crossing phase, where FC2CT performs most
poorly.

Figure 6 displays the longitudinal and transverse phase space as calculated by
RF-Track and FC2CT for a TW cavity comprised of 45 cells. The longitudinal
electric field component at some time, $t$, is shown in Fig. 9 (b). The FC2CT
method closely agrees with the well trusted tracking code, RF-Track. The
output distribution has traversed a 45 cell TW cavity with $4\pi/5$ phase
advance per cell. For protons at 150 MeV, the length of the cavity is
approximately 1 m. The RF-Track time step was set to 0.003 ns, resulting in
approximately 44 iterations per cell. It is shown later that 44 iterations is
beyond the limit of convergence.

Figure 7 displays the difference in calculated phase space components between
FC2CT and RF-Track. Fig. 7 (a) shows a similar error relationship to the SW
case shown in Fig. 5. The majority of particles have small discrepancies in
$P_{z}$, with the phase space of ‘neck’ particles being less well determined.
The maximum difference in $P_{z}$ is approximately 0.02 $\%$, similar to the
error in the SW case. The main difference is that FC2CT overestimates both
$P_{z}$ and $z$ for all particles, as $\Delta P_{z}$ and $\Delta z$ are always
negative. Similar to the SW case, maximum $P_{z}$ particles are well
approximated, with small errors in both position and momenta. The maximum
position error is $\sim$ 0.0012 m, approximately 6$\%$ of a cell length.

Figure 7 (b) displays the discrepancy in transverse momentum as a function of
transverse position discrepancy. The relationship between errors appears to be
a function of $P_{z}$. For low $P_{z}$ particles (deep purple), $\Delta P_{x}$
is a straight line with respect to $\Delta X$. For high $P_{z}$ particles, the
relationship is also linear, albeit with a shallower gradient. For all other
particles, $\Delta P_{x}$ is generally larger, with a maximum discrepancy of
approximately 6$\%$. The maximum position error is $\sim$ 1.2 $\%$ of the beam
pipe radius.

(a) (b)

Figure 7: Sub-figure a/b displays the difference in calculated $P_{z}/P_{x}$
as a function of difference in calculated $z/x$ (defined as
$X_{ASTRA}-X_{FC2CT}$), after 45 TW cells. The initial proton energy is 150
MeV.

From results displayed, the FC2CT is a highly efficient and accurate tracking
method. The FC2CT method has the ability to vastly reduce the computational
expense of large tracking jobs. Whilst FC2CT has been limited above to
iterations every rf cell, the method allows iterations over any distance, by
changing the bounds of integration in Eqn. 15.

###  II.1 Changes to FC2CT

FC2CT is based upon the approximation of constant particle $\beta$ over the
cell. This approximation is highly effective for multiple energy ranges for
both electrons and protons. The accuracy of the approximation is directly
proportional to the change in particle $\beta$ over an rf cell, $\Delta\beta$.
However, for very low energy protons, $\Delta\beta$ may be small over an rf
cell, however $\beta$ may be changing by a factor of two, thus the accuracy of
the approximation must be normalised to the particle $\beta$, giving
$\frac{\Delta\beta}{\beta}$.

Figure 8 shows the ratio $\frac{\Delta\beta}{\beta}$ for electrons and protons
at different energy. The value of $\frac{\Delta\beta}{\beta}$ is calculated by
estimating the energy gain over an rf cell of appropriate cell length with a
fixed gradient. The quantitative value of $\frac{\Delta\beta}{\beta}$ is thus
irrelevant, but the relative value is important for inferring particle energy
where the constant $\beta$ approximation is valid. The eye guides shown at
37.5 and 150 MeV display the region of the phase space plots shown in Figures
3, 4 and 6, where the FC2CT model is valid and accurate. The value of
$\frac{\Delta\beta}{\beta}$ is similar for electrons at 2.5 MeV and protons at
37.5 MeV. As previously discussed, CCL are often appropriate for protons with
energy beyond $\sim$ 30 MeV, as below RFQ and DTL are superior.

(a) Figure 8: $\frac{\Delta\beta}{\beta}$ as a function of particle energy.
Eye guides are shown at 37.5 and 150 MeV. The starting energy is 2.5 MeV.

(a) $\frac{2\pi}{3}$ disk-loaded TW cavity. (b) $\frac{4\pi}{5}$ small
aperture TW cavity.

Figure 9: $E_{z}$ field for two TW cavity types.

The FC2CT method takes approximately four times longer for tracking through TW
structures, as previously described. For certain TW structures, the beam
aperture is small, in order to maximise the shunt impedance. For these
structures, the cut-off frequency in the beam pipe is very low, and the field
is highly evanescent in this region between cells, as is the case for the TW
field in Fig. 9 (b). In these cases, the TW cell can be approximated as a SW
cell, utilizing the same method outlined in Section I.2. The field has a fixed
phase advance per cell and a symmetric on-axis $E_{z}$ field about the centre
of the rf cell, thus is not a true TW. However, the results can be highly
accurate, whilst faster by a factor of $\sim$ 4.

Figure 10 displays the phase space of a proton distribution after traversing a
45 cell TW structure. The true TW FC2CT method (slow) is shown alongside the
SW FC2CT method (quick). Whilst the true TW FC2CT agrees with RF-Track more
closely, there is not a huge drop in accuracy for the quick method.

(a) (b)

Figure 10: Transverse and longitudinal phase space of proton beam (initial
energy 150 MeV) after traversing 45 TW cells, as calculated by RF-Track and
FC2CT as TW (slow) and SW (quick).

Not all TW structures can be treated as SW structures, however. The more
asymmetric the $E_{z}$ field is inside a TW cell, the less accurate the
approach becomes. Figure 9 (a) shows $E_{z}$ for a disk-loaded electron
cavity. For this type of cavity, approximating the field as a symmetric SW is
not appropriate, and the true TW FC2CT is required. However, the number of
Fourier terms required in the Fourier series description of A(z) and B(z) is
far fewer than for the case in Fig. 9 (b), saving computation time.

Accuracy improvements for the FC2CT model have been explored by 1) running an
iteration over an rf cell twice, updating the particle $\beta$ with a more
accurate (mean) value for the second iteration or, 2) reducing the integration
length of the iteration.

Figure 11 displays the change in initial $P_{z}$ calculation (in a 1D scheme)
due to changes in FC2CT parameters, such as number of iterations and
integration length, as described in points 1) and 2). The integration length
is given as $\frac{L_{iteration}}{L_{cell}}$, so an integration length of 0.5
computes two iterations per cell.

(a) Figure 11: The change in initial 1D calculation of $P_{z}$ for protons at
150 MeV due to 1) updating particle $\beta$ with a mean value, or, 2)
increasing the number of iterations per rf cell by decreasing the integration
length.

Figure 11 shows that updating the particle $\beta$ with a mean value does not
change the value of $P_{z}$, relative to changing the integration region. The
green line shows the change in $P_{z}$ from running over half an rf cell
twice, requiring the same number of calculations as the orange line (two
iterations, but updating $\beta$ for the second iteration). A much larger
change occurs when reducing the integration step, relative to updating
$\beta$.

The change in $P_{z}$ converges for an integration length of 0.01. Therefore,
separating the cell into ten separate calculations is as accurate as the model
can be made, beyond which the change is negligible.

(a) (b)

Figure 12: Error in longitudinal and transverse phase space between RF-Track
and FC2CT. Two versions of FC2CT are shown, with one (blue) and two (red)
iterations per cell.

Figure 12 displays the error in phase space between RF-Track and FC2CT, with
one and two iterations per cell. As expected, computing two iterations per rf
cell produced an output distribution closer to RF-Track.

It is possible to reduce the computation time for FC2CT considerably for both
SW and TW structures by using only the principal harmonic.

As FC2CT was designed as fast tracking code, decreasing the integration length
is not advised, as the increase in accuracy is minor, with a maximum change in
$P_{z}$ of $\approx$ 0.002 $\%$, requiring the computation time to increase by
a factor of ten. In addition to reducing the integration length, computing
time could be further reduced by iterating over multiple cells in one
calculation. In general, the integration length must be chosen to reflect the
particle rest mass and energy, in addition to the required accuracy. Another
time-saving option for FC2CT is to use only the first harmonic in the Fourier
Series expansion. This is possible as higher order spatial harmonics cancel
exactly for relativistic particles, and are very small for non-relativistic
particles.

##  III Conclusion and Next Steps

In this paper, the FC2CT method for fast particle tracking was displayed, with
its performance compared relative to two well trusted tracking codes, ASTRA
and RF-Track. FC2CT demonstrated strong agreement with the codes for both SW
and TW cavities. Additional developments for FC2CT were discussed, such as
altering the integration length, and treating certain TW cavities as SW
cavities.

The SW FC2CT method works very effectively, as the number of computations is
low, relative to TW FC2CT. Whilst TW FC2CT was highly accurate as a particle
tracker, the computation time is approximately four times longer than the SW
case. Currently, the effect of fringe fields at the entrance/exit of an rf
cavity is not accounted for. Further work can incorporate the radial impulse
due to these fringe fields to better approximate a real life situation.

In order to quantitatively measure the reduction in computation time, an FC2CT
mode must be constructed in a compiled coding language. In this study, FC2CT
was written in an interpretive language, the script is not written with
minimising computation time as a major objective.

###### Acknowledgements.

The studies presented have been funded through the Cockcroft Core Grant by
STFC Grants No. ST/P002056/1 and ST/V001612/1.

## References

  * [1] Frank Bowman.  Introduction to Bessel functions.  Courier Corporation, 2012.
  * [2] Dassault Systemes.  CST STUDIO SUITE, R2021x.  https://www.3ds.com/products-services/simulia/products/cst-studio-suite/, 2021\.  Computer code.
  * [3] Veliko Dimov, Michele Caldara, Alberto Degiovanni, L Esposito, DA Fink, Marina Giunta, Adam Jeff, A Valloni, AM Lombardi, SJ Mathot, et al.  Beam commissioning of the 750 mhz proton rfq for the light prototype.  In 9th International Particle Accelerator Conference (IPAC2018), Vancouver, BC, Canada, JACoW, Geneva, Switzerland, pages 658–660. sn, 2018.
  * [4] Klaus Flottmann.  Astra-a space charge tracking algorithm.  Technical report, DESY, Hamburg, 2000.
  * [5] David J Griffiths.  Introduction to Electrodynamics Fourth Edition.  2021\.
  * [6] Andrea Latina et al.  Rf-track: Beam tracking in field maps including space-charge effects. features and benchmarks.  In Proceedings of the 28th Linear Accelerator Conference-LINAC16, 2016.
  * [7] N Pichoff.  Beam dynamics basics in RF linacs.  2006\.
  * [8] Georgi P Tolstov.  Fourier series.  Courier Corporation, 2012.
  * [9] Guido Van Rossum and Fred L Drake Jr.  Python tutorial.  Centrum voor Wiskunde en Informatica Amsterdam, The Netherlands, 1995\.
  * [10] Thomas P Wangler.  RF Linear accelerators.  John Wiley & Sons, 2008.