pion/.style=black,charged scalar, dash pattern=on 3pt off 2pt,very thick
phot/.style=black,photon,very thick omeg/.style=black,fermion,very thick
cutAs/.style=black!80,very thick cutBs/.style=black!80,thin,decoration=pre
length=0.06cm, post length=0.0cm,border,segment
length=0.4mm,amplitude=1mm,angle=270, postaction=decorate,draw

# Khuri-Treiman analysis of $J/\psi\to\pi^{+}\pi^{-}\pi^{0}$

M. Albaladejo  0000-0001-7340-9235 S. Gonzàlez-Solís  0000-0003-1947-5420 Ł.
Bibrzycki  0000-0002-6117-4894 C. Fernández-Ramírez  0000-0001-8979-5660 N.
Hammoud  0000-0002-8395-0647 V. Mathieu  0000-0003-4955-3311 M. Mikhasenko
0000-0002-6969-2063 G. Montaña  0000-0001-8093-6682 R. J. Perry
0000-0002-2954-5050 A. Pilloni  0000-0003-4257-0928 A. Rodas
0000-0003-2702-5286 W. A. Smith  0009-0001-3244-6889 A. Szczepaniak
0000-0002-4156-5492 D. Winney  0000-0002-8076-243X

###### Abstract

We study the decay $J/\psi\to\pi^{+}\pi^{-}\pi^{0}$ within the framework of
the Khuri-Treiman equations. We find that the BESIII experimental di-pion mass
distribution in the $\rho(770)$-region is well reproduced with a once-
subtracted $P$-wave amplitude. Furthermore, we show that $F$-wave
contributions to the amplitude improve the description of the data in the
$\pi\pi$ mass region around 1.5 GeV. We also present predictions for the
$J/\psi\to\pi^{0}\gamma^{*}$ transition form factor.

##  1 Introduction

Decays of the lowest-lying charmonium states provide an excellent environment
to study light hadron spectroscopy, search for exotic mesons, test QCD and
QCD-based models, as well as testing theoretical techniques in a region where
both non-perturbative and perturbative QCD effects play a role.

In this work we analyze the decay $J/\psi\to\pi^{+}\pi^{-}\pi^{0}$, to study
the dynamics of the three-pion system at low and intermediate energies under
rather clean conditions. Here, the final state invariant mass distribution can
contain contributions from the $P$-wave $(J^{PC}=1^{--})$ and $F$-wave
$(J^{PC}=3^{--})$ states of the $\pi\pi$ subsystem. Previous experimental
studies from BESII [1] and BABAR [2] showed that the $P$-wave $\rho(770)\pi$
intermediate state dominates the process, but limited statistics prevented any
detailed study of substructures in the $3\pi$ system. While the dominance of
the $\rho(770)$ resonance can be clearly seen in the Dalitz plot distribution
and projection measurements by the BESIII collaboration obtained with roughly
1.9 million $J/\psi\to\pi^{+}\pi^{-}\pi^{0}$ events [3], there are hints of
contributions other than the $\rho(770)$. For example, the absence of events
in the center of the Dalitz plot indicates the contribution from additional
states and/or partial waves which may interfere destructively with the
$\rho(770)$. Exactly the opposite situation is found for the partner reaction
$\psi(2S)\to\pi^{+}\pi^{-}\pi^{0}$. There, the 7872 events from BESIII [3]
show a completely different shape of the $\pi\pi$ invariant mass distribution
and the Dalitz plot — the $\rho\pi$ contribution is subleading and almost all
events are found in the center of the Dalitz plot, with data indicating that
the main contribution comes from a higher mass resonance, i.e. the
$\rho(2150)$ resonance with $J^{PC}=1^{--}$. The different picture between the
$J/\psi$ and $\psi(2S)$ decays into $\pi^{+}\pi^{-}\pi^{0}$ and the lack of
reasonable explanations within the quark model is known as the $\rho\pi$
puzzle and still remains largely unresolved (see e.g. Refs. [4, 5, 6, 7, 8],
and references therein). New high-statistics BESIII data on $J/\psi$ decays
will soon be available [9, 10], which could be used to greatly improve the
theoretical uncertainties associated to vector charmonium decays. In
particular, they might help clarify the $\rho\pi$ puzzle, as well as provide
access to high-precision $\rho$-$\omega$ mixing effect analyses and motivate
coupled channel studies with the decays $J/\psi\to K^{+}K^{-}\pi^{0}$ and
$J/\psi\to K_{S}K^{\pm}\pi^{\mp}$.

The decay $J/\psi\to\pi^{+}\pi^{-}\pi^{0}$ has previously been studied within
the context of the Veneziano model [11], and using aspects of unitarity and
analyticity constraints [12, 13]. Here, we adapt the Khuri-Treiman (KT)
framework [14], applied extensively in the isospin-violating decay $\eta\to
3\pi$ [15, 16, 17, 18, 19, 20, 21] and in the decay of light vector isoscalar
resonances $\omega,\phi\to 3\pi$ [22, 23, 24], to the analysis of the vector
charmonium decay $J/\psi\to\pi^{+}\pi^{-}\pi^{0}$. We show that one
subtraction in the KT equations satisfactorily describes the BESIII
experimental di-pion mass distribution at the peak of the $\rho(770)$. In
addition, we find that $F$-wave effects are needed to describe the
intermediate energy region around 1.5 GeV. We also apply our analysis
techniques to predict the $J/\psi\to\pi^{0}\gamma^{*}$ transition form factor.
Our study lays the groundwork for a detailed analysis of $J/\psi$ decays using
the large data sample currently being collected at BESIII.

This paper is organized as follows. In Section 2 we review the KT formalism
for the $J/\psi\to 3\pi$ decay. In Section 3 we apply the formalism to the
BESIII data and discuss the results. In Section 4, we present predictions for
the $J/\psi\to\pi^{0}\gamma^{*}$ transition form factor, and we summarize our
findings in Section 5.

##  2 Formalism

###  2.1 Decay amplitude and kinematics

The amplitude for the decay
$J/\psi(p_{V})\to\pi^{0}(p_{0})\;\pi^{+}(p_{+})\;\pi^{-}(p_{-})$ can be
expressed in terms of a kinematic prefactor and a single invariant scalar
function $F(s,t,u)$ containing the dynamical information,

$\mathcal{M}(s,t,u)=i\,\epsilon_{\mu\nu\alpha\beta}\;\epsilon^{\mu}(p_{V})\,p_{+}^{\nu}\,p_{-}^{\alpha}\,p_{0}^{\beta}\,\,F(s,t,u)\,,$
(2.1)

where $\epsilon_{\mu\nu\alpha\beta}$ is the Levi-Civita tensor and
$\epsilon^{\mu}(p_{V})$ is the polarization vector of the $J/\psi$ meson. The
particle momenta are related to the Mandelstam variables through:

$s=(p_{+}+p_{-})^{2}\,,\quad t=(p_{0}+p_{+})^{2}\,,\quad
u=(p_{0}+p_{-})^{2}\,,$ (2.2)

with $s+t+u=m_{J/\psi}^{2}+3m_{\pi}^{2}$. In this paper, we work in the
isospin limit with $m_{\pi}\doteq m_{\pi^{\pm}}=m_{\pi^{0}}$ and
$m_{\pi}=(2m_{\pi^{\pm}}+m_{\pi^{0}})/3$. The scattering angle in the
$s$-channel, defined by the center of mass of the $\pi^{+}\pi^{-}$ pair, is
denoted by $\theta_{s}$ and is given by:

$\cos\theta_{s}(s,t,u)=\frac{t-u}{4\,p(s)\,q(s)}\,,\quad\sin\theta_{s}(s,t,u)=\frac{\sqrt{\phi(s,t,u)}}{2\sqrt{s}\,p(s)\,q(s)}\,,$
(2.3)

where the momenta $p(s)$ and $q(s)$,

$p(s)=\frac{\lambda^{\frac{1}{2}}(s,m_{\pi}^{2},m_{\pi}^{2})}{2\sqrt{s}}\,,\quad
q(s)=\frac{\lambda^{\frac{1}{2}}(s,m_{J/\psi}^{2},m_{\pi}^{2})}{2\sqrt{s}}\,,$
(2.4)

are, respectively, the momenta of the $\pi^{\pm}$ and $\pi^{0}$ in the
$s$-channel. $\lambda(a,b,c)=a^{2}+b^{2}+c^{2}-2ab-2bc-2ca$ is the Källén, or
triangle, function [25]. The zeroes of the well-known Kibble function [26] ,

$\phi(s,t,u)=\left(2\sqrt{s}\;\sin\theta_{s}\;p(s)\,q(s)\right)^{2}=s\,t\,u-m_{\pi}^{2}(m_{J/\psi}^{2}-m_{\pi}^{2})^{2}~{},$
(2.5)

define the boundaries of the physical regions of the process. The Dalitz-plot
boundaries in $t$ for a given value of $s$ for $J/\psi\to 3\pi$ lie within the
interval $[t_{\rm{min}}(s),\;t_{\rm{max}}(s)]$, with

$t_{\rm{max,min}}(s)=\frac{m_{J/\psi}^{2}+3m_{\pi}^{2}-s}{2}\pm
2\,p(s)\,q(s)~{},$ (2.6)

while the allowed range for $s$ is given by $s_{\rm{min}}=4m_{\pi}^{2}$ to
$s_{\rm{max}}=(m_{J/\psi}-m_{\pi})^{2}\,.$

Finally, the measured differential decay width can be written in terms of the
invariant amplitude $F(s,t,u)$ as

$\frac{d^{2}\Gamma}{ds\,dt}=\frac{1}{(2\pi)^{3}}\,\frac{1}{32\,m_{J/\psi}^{3}}\,\frac{1}{3}\frac{\phi(s,t,u)}{4}\;|F(s,t,u)|^{2}\,.$
(2.7)

###  2.2 Khuri–Treiman equations for $J/\psi\to 3\pi$

The KT formalism for the $J/\psi\to 3\pi$ amplitude $F(s,t,u)$ is formally
identical to the well-established one for the $\omega\to 3\pi$ decay amplitude
[22, 23, 24, 27], and has been discussed in Ref. [28] (see also Ref. [29]). As
shown in these references, the $s$-channel partial-wave expansion for
$F(s,t,u)$ is given by

$F(s,t,u)=\sum_{J\,\rm{odd}}^{\infty}(p(s)\,q(s))^{J-1}\;P_{J}^{\prime}(z_{s})\;f_{J}(s)\,,$
(2.8)

where $z_{s}=\cos\theta_{s}$ and $P_{J}^{\prime}(z_{s})$ is the derivative of
the Legendre polynomial. The KT representation of the scalar function
$F(s,t,u)$ in Eq. (2.8) may be obtained by replacing the infinite sum of
partial waves in the $s$-channel with the sum of three so-called isobar
amplitudes, one for each of the $s$-, $t$\- and $u$-channels. By truncating
the partial wave expansion of each isobar amplitude at $J_{\text{max}}=1$ we
obtain the following crossing-symmetric isobar decomposition [22, 23, 30]:

$F(s,t,u)=F_{1}(s)+F_{1}(t)+F_{1}(u)\,,$ (2.9)

where each isobar amplitude, $F_{1}(x)$, has only a right-hand or unitary cut
in its respective Mandelstam variable. The relation between $F_{1}(s)$ and
$f_{1}(s)$ is obtained by projecting Eq. (2.9) onto the $s$-channel partial
wave,

$\displaystyle f_{1}(s)$ $\displaystyle=F_{1}(s)+\hat{F}_{1}(s)\,,$ (2.10)
$\displaystyle\hat{F}_{1}(s)$ $\displaystyle\equiv
3\int_{-1}^{1}\frac{dz_{s}}{2}\;(1-z_{s}^{2})\;F_{1}(t(s,z_{s}))\,,$ (2.11)

where the inhomogeneity $\hat{F}_{1}(s)$ contains the $s$-channel projection
of the left-hand cut contributions due to the $t$\- and $u$-channels, and its
evaluation in the decay region requires a proper analytical continuation [31].
Assuming elastic unitarity with only two-pion intermediate states, we arrive
at the KT equation for the $J/\psi\to 3\pi$ decay, i.e. the unitarity relation
for the isobar amplitude $F_{1}(s)$:

$\displaystyle{\rm{disc}}\,F_{1}(s)$
$\displaystyle=2i\left(F_{1}(s)+\hat{F}_{1}(s)\right)\;\sin\delta_{1}(s)\;e^{-i\delta_{1}(s)}\;\theta(s-4m_{\pi}^{2})\,,$
(2.12)

where $\delta_{1}(s)$ is the $P$-wave $\pi\pi$ phase shift, which is real.

Given the discontinuity relation in Eq. (2.12), one can write an unsubtracted
dispersion relation for $F_{1}(s)$ as

$\displaystyle F_{1}(s)=\frac{1}{2\pi
i}\int_{4m_{\pi}^{2}}^{\infty}ds^{\prime}\;\frac{{\rm{disc}}\,F_{1}(s^{\prime})}{s^{\prime}-s}\,,$
(2.13)

the solution of which can be written as:

$F_{1}(s)=\Omega_{1}(s)\left(a+\frac{s}{\pi}\int_{4m_{\pi}^{2}}^{\infty}\frac{ds^{\prime}}{s^{\prime}}\frac{\sin\delta_{1}(s^{\prime})\,\hat{F}_{1}(s^{\prime})}{|\Omega_{1}(s^{\prime})|\left(s^{\prime}-s\right)}\right)\,,$
(2.14)

where $\Omega_{1}(s)$ is the usual Omnès function [32],

$\Omega_{1}(s)=\exp\left[\frac{s}{\pi}\int_{4m_{\pi}^{2}}^{\infty}\frac{ds^{\prime}}{s^{\prime}}\frac{\delta_{1}(s^{\prime})}{s^{\prime}-s}\right]\,.$
(2.15)

The subtraction constant $a$ in Eq. (2.14) is the only free parameter in the
model. It is in general complex, $a=|a|\,e^{i\phi_{a}}$. While its modulus
$|a|$ can be fixed from the experimental $J/\psi\to 3\pi$ decay width, no
observable of the decay is sensitive to the overall phase $\phi_{a}$, so we
can set $\phi_{a}=0$. Since it determines the overall normalization of the
amplitude, the constant $a$ can be factored out.

We note that due to the asymptotic behavior of $F_{1}(s)$ in Eq. (2.14), the
amplitude $F(s,t,u)$ satisfies the Froissart-Martin bound [33, 34, 22]. Also
note that, even though $F_{1}(s)/\Omega_{1}(s)$ in Eq. (2.14) looks like a
once-subtracted dispersion relation, $F_{1}(s)$ actually satisfies the
unsubtracted dispersion relation given in Eq. (2.13). Therefore, the energy
dependence of $F_{1}(s)$ is a pure prediction given solely by the phase shift
$\delta_{1}(s)$. Here, we take $\delta_{1}(s)$ from the phase shift
parametrizations of Ref. [35] that are valid roughly up to $\sqrt{s}=2$ GeV.
These phase shifts contain information about inelastic channels, but given
that the inelasticity is found to be rather small until about 1.4 GeV we
refrain to consider them. Therefore, the phase shift that we employ have the
physics of the $\rho(770)$ and also the effects of the higher $\rho(1450)$ and
$\rho(1770)$ resonances. For our analysis, beyond $\sqrt{s}=\Lambda\equiv
1.85$ GeV we smoothly guide the $\delta_{1}(s)$ to $\pi$ through [36, 27]

$\displaystyle\delta_{\infty}(s)\equiv\lim_{s\to\infty}\delta_{1}(s)=\pi-\frac{\alpha}{\beta+\left(s/\Lambda^{2}\right)^{3/2}}\,,$
(2.16)

where $\alpha$ and $\beta$ are parameters introduced so that the phase
$\delta_{1}(s)$ and its first derivative $\delta^{\prime}(s)$ are continuous
at $s=\Lambda^{2}$. Their explicit expressions read

$\displaystyle\alpha=\frac{3\left(\pi-\delta_{1}(\Lambda^{2})\right)^{2}}{2\Lambda^{2}\delta_{1}^{\prime}(\Lambda^{2})}\,,\quad\beta=-1+\frac{3\left(\pi-\delta_{1}(\Lambda^{2})\right)}{2\Lambda^{2}\delta_{1}^{\prime}(\Lambda^{2})}\,.$
(2.17)

This ensures the expected asymptotic $1/s$ behavior of $\Omega_{1}(s)$. The
three phase shifts $\delta_{1}(s)$ from Ref. [35] that we use as an input are
shown in Fig. 1 up to 2.5 GeV. The different solutions come from using
different $\pi\pi$ scattering data sources. As seen, the behavior of the phase
shift solution I suggests a large interference between the $\rho^{\prime}$ and
$\rho^{\prime\prime}$, with a sizable change in the phase in the region
between 1.5 and 1.8 GeV, while solutions II and III looks smoother in this
region. For our analysis, we use solution I as our central input for the phase
and solutions II and III to quantify the systematic uncertainties in our
calculations.

Figure 1: Solutions I, II and III for the $P$-wave phase shift $\delta_{1}(s)$
from Ref. [35] valid roughly up to $\sqrt{s}=2$ GeV. The solution of Ref. [37]
(dotted red line) is valid only up to about $\sqrt{s}=1.3$ GeV, and is shown
for completeness.

We solve Eq. (2.14) following a numerical iterative procedure similar to Refs.
[22, 16, 20, 21, 38]. We use $F_{1}(s)=\Omega_{1}(s)$ as an efficient initial
input to calculate $\hat{F}_{1}(s)$ from Eq. (2.11), which subsequently is
inserted as an input in Eq. (2.14) for the computation of an updated
$F_{1}(s)$. This cyclic calculation is repeated until the solution converges.
In Fig. 2, we show the solutions for $F_{1}(s)$ (normalized to $a=1$) after
each iteration step along with the initial input (dashed blue line). As can be
seen, convergence is achieved after three iterations. The difference between
the final solution (solid black) and the starting point, i.e.
$F_{1}(s)=\Omega_{1}(s)$ (dashed blue), is rather small, hinting at moderate
crossed-channel effects.

Figure 2: Convergence behavior of the iterative procedure for the real (left
plot) and imaginary (right plot) parts of the amplitude $F_{1}(s)$ given in
Eq. (2.14) using solution I of the phase shift $\delta_{1}(s)$ as input. The
vertical line denotes the two-pion threshold.

Note that when the crossed-channel rescattering effects are removed from the
isobar $F_{1}(s)$, i.e. when $\hat{F}_{1}(s)=0$ in Eq. (2.14), $F_{1}(s)$ is
simply the pure Omnès function multiplied by a constant,

$F_{1}(s)=a^{\prime}\Omega_{1}(s)\,,$ (2.18)

which implies the following isobar decomposition of the full amplitude (cf.
Eq. (2.9)):

$F(s,t,u)=a^{\prime}\left(\Omega_{1}(s)+\Omega_{1}(t)+\Omega_{1}(u)\right)\,.$
(2.19)

In this case, a new normalization constant $a^{\prime}$ has to be chosen to
reproduce the $J/\psi\to 3\pi$ decay width. Also note that Eq. (2.14) can be
written in the form

$F_{1}(s)=\Omega_{1}(s)\left(a+b^{\prime}\,s+\frac{s^{2}}{\pi}\int_{4m_{\pi}^{2}}^{\infty}\frac{ds^{\prime}}{(s^{\prime})^{2}}\frac{\sin\delta_{1}(s^{\prime})\,\hat{F}_{1}(s^{\prime})}{|\Omega_{1}(s^{\prime})|\left(s^{\prime}-s\right)}\right)\,,$
(2.20)

where $b^{\prime}$ satisfies the following sum rule [22]:

$b\equiv
b^{\prime}/a=\frac{1}{\pi}\int_{4m_{\pi}^{2}}^{\infty}\frac{ds^{\prime}}{(s^{\prime})^{2}}\frac{\sin\delta_{1}(s^{\prime})\,\hat{F}_{1}(s^{\prime})/a}{|\Omega_{1}(s^{\prime})|}\,.$
(2.21)

The subtraction constant, $b$, is complex due to the presence of the three-
particle cut in the physical region of the decay amplitude. This value is
found to be:

$b_{\text{sum}}\simeq\,0.141\,e^{2.321\,i}\ \text{GeV}^{-2}~{}.$ (2.22)

Had we used solution II or III of the phase shift $\delta_{1}(s)$ (cf. Fig.
1), we would have obtained $b_{\text{sum}}\simeq\,0.129\,e^{2.640\,i}$ GeV-2
and $b_{\text{sum}}\simeq\,0.124\,e^{2.811\,i}$ GeV-2, respectively.

Performing one subtraction on Eq. (2.13) leads to the solution [22, 20, 30]:

$F_{1}(s)=a\left[F_{a}(s)+b\,F_{b}(s)\right]\,,$ (2.23a) where now $b$ is not
constrained to satisfy Eq. (2.21), and the functions $F_{a}(s)$ and $F_{b}(s)$
are given by $\displaystyle F_{a}(s)$ $\displaystyle=$
$\displaystyle\Omega_{1}(s)\left[1+\frac{s^{2}}{\pi}\int_{4m_{\pi}^{2}}^{\infty}\frac{ds^{\prime}}{s^{\prime
2}}\frac{\sin\delta_{1}(s^{\prime})\,\hat{F}_{a}(s^{\prime})}{|\Omega_{1}(s^{\prime})|(s^{\prime}-s)}\right]\,,$
(2.23b) $\displaystyle F_{b}(s)$ $\displaystyle=$
$\displaystyle\Omega_{1}(s)\left[s+\frac{s^{2}}{\pi}\int_{4m_{\pi}^{2}}^{\infty}\frac{ds^{\prime}}{s^{\prime
2}}\frac{\sin\delta_{1}(s^{\prime})\,\hat{F}_{b}(s^{\prime})}{|\Omega_{1}(s^{\prime})|(s^{\prime}-s)}\right]\,.$
(2.23c)

These functions only need to be calculated once since they are independent of
the numerical values of $a$ and $b$ and, as we will discuss in Sec. 3, $a$ and
$b$ will become fit parameters. In Fig. 3, we show the solutions for
$F_{a}(s)$ and $F_{b}(s)$ using a numerical iterative procedure similar to the
one described previously. In this case, nine iterations are needed to obtain
convergent solutions. Strictly speaking, the amplitude $F(s,t,u)$ built from
$F_{1}(s)$ in Eq. (2.23a) does not satisfy the asymptotic Froissart-Martin
bound for an arbitrary value of the parameter $b\neq b_{\text{sum}}$ [cf. Eq.
(2.22)]. The main advantage of introducing one subtraction is that, due to the
additional $1/s^{\prime}$ factor introduced, we reduce the importance of the
high energy region of the dispersion integrals where the phase shift is not
well known. By letting the subtraction constant $b$ be a free parameter, we
can partially absorb our ignorance of the higher energy part of the integral.
This allows us to parametrize some unknown energy dependence of the $J/\psi\to
3\pi$ interaction not directly related to $\pi\pi$ rescattering. As we will
show in the following section, the once-subtracted parametrization provides a
good representation of the data from BESIII in the $\rho(770)$ resonance
region.

Figure 3: Convergence behavior of the iterative procedure for the real (left
plots) and imaginary (right plots) parts of the amplitudes $F_{a}(s)$ (Eq.
(2.23b), upper plots) and $F_{b}(s)$ (Eq. (2.23c), lower plots) using solution
I of the phase shift as $\delta_{1}(s)$ input. The vertical line denotes the
two-pion threshold.

##  3 Results

###  3.1 $P$-wave contribution

We now compare our KT amplitudes defined in the previous section to the
experimental data from the BESIII collaboration [3]. Given that the Dalitz
plot distribution is not publicly available, we are only able to analyze the
di-pion mass projection of Eq. (2.7), computed on the $\sqrt{s}\equiv
m_{\pi\pi}$ invariant mass, shown in Fig. 2 of Ref. [3]. A Poisson
distribution is assumed to assign uncertainty for every bin. High statistics
of the data sample make it challenging to achieve an accurate description of
the data with reasonably simple models. Nevertheless, we will be able to
obtain a qualitative description of the data in the whole energy range. We
start by using the unsubtracted KT amplitude Eq. (2.14). The single free
parameter $a$ only affects the overall normalization of the amplitude and can
be fixed from the $J/\psi\to\pi^{+}\pi^{-}\pi^{0}$ decay width. Using the PDG
values $\Gamma_{J/\psi}=92.6$ keV and
${\rm{BR}}(J/\psi\to\pi^{+}\pi^{-}\pi^{0})=2.10(8)\%$ [39] one finds
$|a|\simeq 0.051$ GeV-3. In Fig. 4, we compare our prediction to the
$m_{\pi\pi}$ distribution by BESIII with proper normalization [cf. Eq. (3.1)].

Figure 4: BESIII (red circles) [3] measurement of the $m_{\pi\pi}$ invariant
mass distribution for the decay $J/\psi\to 3\pi$ as compared to our prediction
without crossed-channel effects (dotted brown line), with the unsubtracted KT
amplitude (dashed green line) and our fit in Eq. (3.2) including one
subtraction (black solid line). The gray band accounts for the systematic
uncertainties attached to our calculations. See main text for details.

In the figure, we also show the result obtained when the crossed-channel
rescattering is neglected [cf. Eq. (2.19)], in which case the global
normalization is found to be $|a^{\prime}|\simeq 0.046$ GeV-3. As can be
observed, the result of the latter solution (dotted brown line) lies below
that of the unsubtracted KT $F_{1}(s)$ solution at the peak of the
$\rho$-meson, and neither reproduce the experimental data in this region. In
addition, both appear to fail at describing the intermediate energy region. In
order to achieve a better description of the data, we next use the more
flexible, once-subtracted amplitude Eqs. (2.23b) and (2.23c), with the
additional subtraction constant $b$ fitted to BESIII data. For our analysis,
we define

$\chi^{2}_{{\rm{data}}}=\sum_{i=1}\left(\frac{N_{\text{ev},i}-\mathcal{N}d\Gamma_{i}^{\rm{th}}/dm_{\pi\pi}}{\sigma_{N_{\text{ev},i}}}\right)^{2}\,,$
(3.1)

where $N_{\text{ev},i}$ and $\sigma_{N_{\text{ev},i}}$ are, respectively, the
experimental number of events distribution and the corresponding error in the
$i$-th bin and $d\Gamma_{i}^{\rm{th}}/dm_{\pi\pi}$ is the theoretical
expression for the decay distribution [cf. Eq. (2.7)]. For
$\sigma_{N_{\text{ev},i}}$ we take $\sqrt{N_{{\rm{ev}},i}}$. The constant
$\mathcal{N}$ is at this stage an arbitrary normalization. Since we are not
determining the branching ratio, we reabsorb the global normalization of the
amplitude $a$ into $\mathcal{N}$ and fix alone this overall constant from the
fit to the BESIII data. The sum in Eq. (3.1) runs over the 80 data points and
we take into account an efficiency of about $0.3$ for the number of events and
the errors in our fits [3].

The $\chi^{2}_{{\rm{data}}}$ minimization yields

$|b|=0.198(35)~{}{\rm{GeV^{-2}}}\,,\quad\phi_{b}=2.675(300)\,,$ (3.2)

which implies $|a|=0.0565(22)$ GeV-3 for the normalization of the amplitude
upon using the ${\rm{BR}}(J/\psi\to\pi^{+}\pi^{-}\pi^{0})$ from the PDG. The
statistical error is negligible and the quoted error is the theoretical
systematic uncertainty attached to our calculations. This is obtained from the
absolute value of the difference between the fits performed with solutions I
(central solution) and III of the phase shift $\delta_{1}^{1}(s)$ (cf. Fig.
1), which gives the largest variation. We observe that the systematic errors
attached are sizable, of about $18\%$ and $11\%$ for $|b|$ and $\phi_{b}$,
respectively. We also note that this value stays close to its sum-rule
prediction given in Eq. (2.22). Therefore, we conclude that the pion-pion
$P$-wave phase shift saturates the sum rule for the $J/\psi\to 3\pi$ partial
wave to about $75\%$. This result is to be compared to similar sum rules for
$\omega\to 3\pi$ in Ref. [38], where the fitted value of $b$ was found to be
quite different than its sum-rule $b_{\rm{sum}}$, and for $\phi\to 3\pi$ in
Ref. [22], where it was observed that the difference between the fitted $b$
and $b_{\rm{sum}}$ was small. The result of the fit is shown in Fig. 4 with
the normalization of the events distribution resulting from the fits,
$\mathcal{N}=7.64(33)\times 10^{8}$ in units of $(2.4\,\rm{MeV})^{-1}$. The
gray error band in the figure accounts for the systematic uncertainties
associated to our fits and is defined as the (symmetrized) difference between
the fit results obtained with solution I of the phase shift with respect to
the ones from solution III, which give the largest difference. It can be seen
that this fit provides a satisfactory description of experimental data up to
$m_{\pi\pi}\sim 1$ GeV (the elastic region). However, we obtain high values of
the $\chi^{2}/$dof of about 200 but this problem is not critical. We shall
come back to discuss this point below. Here we stress that the once-subtracted
KT amplitude is able to reproduce the $\rho(770)$ function shape and note that
contributions of partial waves other than the elastic $P$-wave, which is the
main one, seem to be required to describe the intermediate energy region
around $m_{\pi\pi}\sim 1.5$ GeV. The next allowed partial wave is the
$F$-wave, which we will include in the following subsection. As we will see,
the inclusion of an explicit $F$-wave improves the quality of the fit.

In Fig. 5, we show the Dalitz plot distribution resulting from our fit, which
exhibits unambiguous contributions from $\rho(770)$ resonances which appear as
bands along the Dalitz plot boundaries, with almost no events in the center of
the Dalitz plot. The visual comparison with the corresponding BESIII Dalitz-
plot data shows a good agreement (see Fig. 2 in Ref. [3]).

Figure 5: Dalitz plot distribution $d^{2}\Gamma/ds\,dt$ (in arbitrary units)
resulting from our fit in Eq. (3.2).

###  3.2 Inclusion of the $F$-wave contribution

The isobar decomposition of the amplitude including $F$-waves follows from Eq.
(2.8) and reads [22, 38]:

$\displaystyle\begin{aligned}
F(s,t,u)&=F_{1}(s)+F_{1}(t)+F_{1}(u)\\\\[4.30554pt]
&+(p(s)q(s))^{2}P_{3}^{\prime}(z_{s})F_{3}(s)+(p(t)q(t))^{2}P_{3}^{\prime}(z_{t})F_{3}(t)+(p(u)q(u))^{2}P_{3}^{\prime}(z_{u})F_{3}(u)\,,\end{aligned}$
(3.3)

where $F_{1}(s)$ is the $P$-wave isobar [cf. Eq. (2.23a)], $F_{3}(s)$ is the
$F$-wave isobar amplitude, which as $F_{1}(s)$ only has a right-hand cut, and:

$z_{t}=\frac{s-u}{4p(t)q(t)}~{},\qquad z_{u}=\frac{s-t}{4p(u)q(u)}~{}.$ (3.4)

The discontinuity of the $F$-wave is expressed by:

$\displaystyle{\rm{disc}}\,F_{3}(s)$
$\displaystyle=2i\left(F_{3}(s)+\hat{F}_{3}(s)\right)\;\sin\delta_{3}(s)\;e^{-i\delta_{3}(s)}\;\theta(s-4m_{\pi}^{2})\,,$
(3.5)

where $\delta_{3}(s)$ and $\hat{F}_{3}(s)$ are the $F$-wave phase shift and
inhomogeneity, respectively. Here, we will simplify Eq. (3.5) by neglecting
$\hat{F}_{3}(s)$, as done for instance in Ref. [20]. The solution is then
given by:

$F_{3}(s)=p_{3}(s)\Omega_{3}(s)\,,$ (3.6)

where $\Omega_{3}(s)$ is the $F$-wave Omnès function (cf. Eq. (2.15))

$\Omega_{3}(s)=\exp\left[\frac{s}{\pi}\int_{4m_{\pi}^{2}}^{\infty}\frac{ds^{\prime}}{s^{\prime}}\frac{\delta_{3}(s^{\prime})}{s^{\prime}-s}\right]\,.$
(3.7)

In order to obtain the required input phase $\delta_{3}(s)$, we model the
$F$-wave contribution by a $\rho_{3}(1690)$ resonance ($J^{PC}=3^{--}$). While
the dominant decay mode of the $\rho_{3}(1690)$ is to $4\pi$, we only consider
here its decay to $\pi\pi$ and neglect inelastic channels effects. We use the
following Breit-Wigner representation for $F_{3}(s)$:

$F_{3}(s)|_{\rm{BW}}=\frac{m_{\rho_{3}}^{2}}{m_{\rho_{3}}^{2}-s-im_{\rho_{3}}\Gamma^{\ell=3}_{\rho_{3}}(s)}\,,$
(3.8)

with the energy-dependent width given by

$\displaystyle\Gamma^{\ell}_{R}(s)$ $\displaystyle=$
$\displaystyle\frac{\Gamma_{R}m_{R}}{\sqrt{s}}\left(\frac{p(s)}{p(m_{R}^{2})}\right)^{2\ell+1}\left(F_{R}^{\ell}(s)\right)^{2}\,.$
(3.9)

The $F_{R}^{\ell}(s)$ denotes the Blatt-Weisskopf factor that limits the
growth of the isobar [40]. For $\ell=3$ it is given by:

$\displaystyle F_{R}^{\ell=3}(s)$ $\displaystyle=$
$\displaystyle\sqrt{\frac{z_{0}(z_{0}-15)^{2}+9(2z_{0}-5)^{2}}{z(z-15)^{2}+9(2z-5)^{2}}}\,,\quad
z=r_{R}^{2}p^{2}(s)\,,\quad z_{0}=r_{R}^{2}p^{2}(m_{\rho_{3}}^{2})\,,$ (3.10)

with the hadronic scale $r_{R}=2{\mathrm{\,Ge\kern-1.00006ptV}}^{-1}$. The
phase can then be computed from the relation

$\tan\delta_{3}(s)=\frac{{\rm{Im}}F_{3}(s)|_{\rm{BW}}}{{\rm{Re}}F_{3}(s)|_{\rm{BW}}}\,,$
(3.11)

which completes our representation of the $F$-wave isobar $F_{3}(s)$. Using
$m_{\rho_{3}}=1688$ MeV and $\Gamma_{\rho_{3}}=161$ MeV from the PDG, in Fig.
6 we display the model for the phase $\delta_{3}(s)$ Eq. (3.11) and the output
for the corresponding Omnès function $\Omega_{3}(s)$ Eq. (3.7) that we use for
our analysis.

Figure 6: $F$-wave phase shift $\delta_{3}(s)$ Eq. (3.11) (left plot) and
output for the Omnès function $\Omega_{3}(s)$ Eq. (3.7) (right plot).

Finally, the function $p_{3}(s)$ in Eq. (3.6) is a polynomial that
parametrizes the energy dependence not directly related to the propagation of
the $\rho_{3}(1690)$ resonance and fixes the strength of the $F$-wave
amplitude. In order to achieve a satisfactory description of the data, we take
$p_{3}(s)$ linear in $s$ with parameters relative to the $P$-wave amplitude,
i.e. $p_{3}(s)=a(|c|e^{i\phi_{c}}+|d|e^{i\phi_{d}}\,s)$, such that the overall
normalization of the amplitude $a$ can be factored out in Eq. (3.3) and
absorbed in $\mathcal{N}$ (cf. Eq. (3.1)) as in the previous subsection. By
minimizing Eq. (3.1), we obtain the following values for the fit parameters:

$|b|=0.205(34)~{}{\rm{GeV^{-2}}}\,,\quad\phi_{b}=2.784(298)\,,$ (3.12)

for the $P$-wave subtraction constant, and

$\displaystyle\begin{aligned} |c|\times
10^{2}&=4.38(1.46)~{}{\rm{GeV^{-4}}}\,,\quad&\phi_{c}&=3.80(5)\,,\\\\[4.30554pt]
|d|\times
10^{2}&=1.58(46)~{}{\rm{GeV^{-6}}}\,,\quad&\phi_{d}&=0.65(8)\,,\end{aligned}$
(3.13)

for the parameters of the $F$-wave subtracted polynomial $p_{3}(s)$. Again,
the quoted error in the previous equations is the systematic uncertainty
obtained from using the different $P$-wave phase shifts $\delta_{1}(s)$ as
input. The result of this fit implies $|a|=0.0581(60)$ GeV-3 for the overall
normalization of the amplitude and it is plotted in Fig. 7 as the dash-dotted
blue line using the event distribution normalization from the fits,
$\mathcal{N}=8.09(41)\times 10^{8}$ in units of $(2.4\,\rm{MeV})^{-1}$. In the
figure, the result of the standalone $P$-wave fit [cf. Eq.(3.2)] is also shown
for comparison. As seen, the $\rho_{3}(1690)$-induced $F$-wave contribution
improves the description of the data around 1.5 GeV. Numerically, we find that
the individual $F$-wave contribution is rather small, while the interference
between the $P$\- and $F$-waves gives a correction of a few percent in the
region $m_{\pi\pi}\sim 1.5$ GeV. The $\chi^{2}/$dof remains high (about 100).
However, with the systematic uncertainties associated to our fits (blue error
band in Fig. 7), we conclude that our representation of the amplitude is
capable of describing the two more prominent features shown by the data: the
line shape of the BESIII measurements in the vicinity of the $\rho(770)$
resonance as well as the movement of the function at $m_{\pi\pi}\sim 1.5$ GeV
due to the $F$-wave effects.111We shall wait for the arrival of new Dalitz
distribution experimental data from BESIII to ascribe a strict statistical
meaning to our $\chi^{2}$ fits. As for the Dalitz-plot distribution, the
$F$-wave effects provides no significant change with respect to Fig. 3.2 and
we thus refrain to show them here.

Figure 7: BESIII (red circles) [3] measurement of the $m_{\pi\pi}$ invariant
mass distribution for the decay $J/\psi\to 3\pi$ as compared to our fits in
Eqs. (3.2) (solid black line), (3.12) and (3.13) (dot-dashed blue line). The
blue error band accounts for the systematic uncertainties attached to our
calculations. See main text for details.

##  4 $J/\psi\to\pi^{0}\gamma^{*}$ transition form factor

Figure 8: Diagrammatic representation of the two-pion contribution to the
discontinuity of the $J/\psi\pi^{0}$ transition form factor [cf. Eq. (4.1)].
The blue and red circles represent, respectively, the full $s$-channel
$P$-wave $J/\psi\to 3\pi$ amplitude $f_{1}(s)$ and the pion vector form factor
$F_{\pi}^{V}(s)$.

The $J/\psi\pi^{0}$ transition form factor (TFF), $f_{J/\psi\pi^{0}}(s)$,
governs the $J/\psi\to\pi^{0}\gamma^{\ast}$ amplitude and its energy
dependence is experimentally accessible from the decays
$J/\psi\to\pi^{0}e^{+}e^{-}$ and $J/\psi\to\pi^{0}\mu^{+}\mu^{-}$. At present,
there is no measurement of the shape of the form factor and the only
experimental information on these decays is the measurement of the branching
ratio by the BESIII collaboration, $BR(J/\psi\to\pi^{0}e^{+}e^{-})=(7.56\pm
1.32\pm 0.50)\times 10^{-7}$ [41]. This measurement was obtained subtracting
the $\rho$ resonance contribution and assuming that excited $c\bar{c}$
exchanges, e.g. coming from off-shell $\psi^{\prime}$ contributions, dominate
the energy-dependence of the form factor. Refs. [42, 28] showed that
subtracting this contribution is not well motivated, as the light vector meson
contributions to the form factor actually dominate the decay. Using the
formalism previously employed for the decays of light vector mesons
$\omega/\phi\to\pi^{0}\gamma^{*}$ [43, 24], we present a dispersive
description of $f_{J/\psi\pi^{0}}(s)$ comparable to Ref. [28], but with the
difference that our analysis is driven by the $J/\psi\to 3\pi$ experimental
data analysis presented in Sec. 3.

A dispersive representation of $f_{J/\psi\pi^{0}}(s)$ is fully determined, up
to possible subtractions, by the discontinuity across the right hand cut.
Here, we focus on the light-quark resonance contributions to the
discontinuity, which dominate the form factor at low and intermediate
energies. Additional $c\bar{c}$ contributions can arise close to the upper
limit of the accessible phase space, $\sqrt{s}=m_{J/\psi}-m_{\pi^{0}}$, and in
fact can dominate the transition form factor there [42, 28], but these
contributions appear in a region of the Dalitz decays which are strongly
suppressed by phase space [42, 28], rendering the task of experimentally
observing them nearly impossible. Bearing this in mind, and because of the
absence of experimental data for the form factor, we do not consider them in
our analysis.

In order to be consistent with the elastic approximation in the
$J/\psi\to\pi^{+}\pi^{-}\pi^{0}$ study, we include only the two-pion
intermediate state contribution to the discontinuity (see Fig. 8 for a
diagrammatic interpretation):

${\rm{disc}}f_{J/\psi\pi^{0}}(s)=i\;\frac{p^{3}(s)}{6\pi\sqrt{s}}\;{F_{\pi}^{V}}^{*}(s)\;f_{1}(s)\;\theta(s-4m_{\pi}^{2})\,,$
(4.1)

which requires as input the full $s$-channel $P$-wave $J/\psi\to 3\pi$
amplitude $f_{1}(s)$ given in Eq. (2.10) and the pion vector form factor
complex-conjugate ${F_{\pi}^{V}}^{*}(s)$, which we approximate by the Omnès
function (complex-conjugate) given in Eq. (2.15). Given that we are using a
once-subtracted dispersion relation for the $J/\psi\to 3\pi$ KT equations, an
unsubtracted dispersion relation for the TFF, as used for instance in Ref.
[28], would result in a divergent integral if no cutoff is used. Therefore, we
use a once-subtracted dispersion relation for the TFF itself,

$\displaystyle
f_{J/\psi\pi^{0}}(s)=|f_{J/\psi\pi^{0}}(0)|\,e^{i\phi_{J/\psi\pi^{0}}(0)}+\frac{s}{12\pi^{2}}\int_{4m_{\pi}^{2}}^{\infty}\frac{ds^{\prime}}{(s^{\prime})^{3/2}}\frac{p^{3}(s^{\prime})\;{F_{\pi}^{V}}^{*}(s^{\prime})\;f_{1}(s^{\prime})}{(s^{\prime}-s)}\,,$
(4.2)

where we indicate explicitly the existence of a non-vanishing phase of
$f_{J/\psi\pi^{0}}(s)$ at $s=0$. This is implied by the cross-channel effects,
i.e. the functions ${F_{\pi}^{V}}^{*}(s)$ and $f_{1}(s)$ do not have the same
phase, and the discontinuity of $f_{J/\psi\pi^{0}}(s)$ is in general complex
[43, 24]. The modulus of the subtraction constant $|f_{J/\psi\pi^{0}}(0)|$ can
be fixed from the $J/\psi\to\pi^{0}\gamma$ partial decay width

$\Gamma(J/\psi\to\pi^{0}\gamma)=\frac{e^{2}(m_{J/\psi}^{2}-m_{\pi^{0}}^{2})^{3}}{96\pi
m_{J/\psi}^{3}}\;|f_{J/\psi\pi^{0}}(0)|^{2}\,.$ (4.3)

Using the value of the partial decay width of $J/\psi\to\pi^{0}\gamma$ [39] in
combination with the above equation, one obtains:

$|f_{J/\psi\pi^{0}}(0)|=6.0(3)\times 10^{-4}\quad\rm{GeV}^{-1}\,.$ (4.4)

The phase $\phi_{J/\psi\pi^{0}}(0)$ is a free parameter that can only be
accessed from the transition form factor experimental data (see e.g. Ref.
[24]). Due to the absence of data for $J/\psi\to\pi^{0}\gamma^{\ast}$, we set
$\phi_{J/\psi\pi^{0}}(0)=0$ in our study.

In Fig. 9, we show up to $\sqrt{s}=2$ GeV our prediction for the absolute
value of the transition form factor resulting from Eq. (4.2) and using the
results from Eq. (3.2) (solid black line). This is our central result for the
form factor. In this figure, however, we also show the result of using the
unsubtracted KT solution for $J/\psi\to 3\pi$ (dashed blue line). It is worth
noting that both curves are similar and only a slight difference is observed
at the $\rho$ peak. Additionally, the calculations when an unsubtracted
dispersion relation for the form factor is used are also shown in the figure,
both with an unsubtracted (dotted red line) and once-subtracted (dot-dashed
green line) $J/\psi\to 3\pi$ amplitude. In the latter case, we have cut the
dispersive integral at 4 GeV2 to avoid the dispersion relation to diverge.
Again, both curves are similar. In this case, the value at the real photon
energy can be calculated from the sum rule [43, 28]:

$f_{J/\psi\pi^{0}}(0)=\frac{1}{12\pi^{2}}\int_{4m_{\pi}^{2}}^{\infty}ds^{\prime}\frac{p^{3}(s^{\prime})F_{\pi}^{V*}(s^{\prime})f_{1}(s^{\prime})}{(s^{\prime})^{3/2}}\,.$
(4.5)

This value is found to be $|f_{J/\psi\pi^{0}}(0)|=5.0(2)\times 10^{-4}$ GeV-1
for both versions of the unsubtracted dispersion relation. The quoted
uncertainty is the systematic error from using the different phase shifts as
input. This value is in qualitative agreement with the value extracted from
the measured $J/\psi\to\pi^{0}\gamma$ in Eq. (4.4), indicating that the
normalization is saturated by the two-pion intermediate state contribution by
roughly $85\%$. The difference between the various lines provides an estimate
of the theoretical uncertainty associated to our description. We expect our
study to strengthen the case for new experimental measurements of the shape of
this form factor, which would allow improving the understanding of radiative
$J/\psi$ decays.

Figure 9: Prediction for the absolute value of the transition form factor
$J/\psi\to\pi^{0}\gamma^{*}$ using Eq. (4.2) (solid black line) and variants
of it. See main text for details.

##  5 Summary

We have analyzed the decay $J/\psi\to\pi^{+}\pi^{-}\pi^{0}$ within the
framework of the Khuri-Treiman equations, which satisfy the constraints
imposed by unitarity, analyticity and crossing symmetry. We have included the
$P$-wave effects of the $\pi\pi$ subsystem up to around 2 GeV, which are
controlled by the $\pi\pi$ $P$-wave scattering-phase shift. We have seen that
one subtraction in the $P$-wave amplitude is necessary to achieve a good
description of the experimental data in the $\rho(770)$-region. The
corresponding subtraction constant was fixed from fits to the di-pion
invariant mass distribution from BESIII. We have also seen that the $P$-wave
alone is not capable of reproducing the data in the mass region around
$m_{\pi\pi}\sim 1.5$ GeV, and that the inclusion of an $F$-wave contribution
arising from the $\rho_{3}(1690)$ brings theory closer to data in this region.
In addition, we have provided predictions for the transition form factor
$J/\psi\to\pi^{0}\gamma^{*}$ up to 2 GeV. Our study lays the groundwork for an
event-by-event likelihood fit of high-precision data from $J/\psi$ decays,
which are expected to be available from BESIII in a near future.

## Acknowledgments

The authors would like to thank Joshua Jackson and Ryan Mitchell (Indiana
University) for fruitful discussions. MA is supported by Generalitat
Valenciana under Grant No. CIDEGENT/2020/002, and by the Spanish Ministerio de
Ciencia e Innovación (MICINN) under contracts No. PID2020-112777GBI00. The
work of SGS is supported by the Laboratory Directed Research and Development
program of Los Alamos National Laboratory under project number 20210944PRD2,
and by the U.S. Department of Energy through the Los Alamos National
Laboratory. Los Alamos National Laboratory is operated by Triad National
Security, LLC, for the National Nuclear Security Administration of U.S.
Department of Energy (Contract No. 89233218CNA000001). This work was supported
by the U.S. Department of Energy contract DE-AC05-06OR23177, under which
Jefferson Science Associates, LLC operates Jefferson Lab, U.S. Department of
Energy Grants No. DE-FG02-87ER40365 and No. DE-FG02-92ER40735, CONACYT
(Mexico) Grant No. A1-S-21389, and Spanish national Grants
PID2020-118758GB-I00 and PID2019–106080 GB-C21. CFR is supported by Spanish
Ministerio de Educación y Formación Profesional under Grant No. BG20/00133. VM
is a Serra Húnter fellow. The work of MM is funded by the Deutsche
Forschungsgemeinschaft under Germany’s Excellence Strategy-EXC-2094-390783311.
DW is supported by National Natural Science Foundation of China Grant No.
12035007 and the NSFC and the Deutsche Forschungsgemeinschaft (DFG, German
Research Foundation) through the funds provided to the Sino-German
Collaborative Research Center TRR110 “Symmetries and the Emergence of
Structure in QCD” (NSFC Grant No. 12070131001, DFG ProjectID 196253076-TRR
110). This work contributes to the aims of the U.S. Department of Energy
ExoHad Topical Collaboration, contract DE-SC0023598.

## References

  * Bai  _et al._ [2004] J. Z. Bai _et al._ (BES), Phys. Rev. D 70, 012005 (2004), arXiv:hep-ex/0402013.
  * Aubert  _et al._ [2004] B. Aubert _et al._ (BaBar), Phys. Rev. D 70, 072004 (2004), arXiv:hep-ex/0408078.
  * Ablikim  _et al._ [2012] M. Ablikim _et al._ (BESIII), Phys. Lett. B 710, 594 (2012), arXiv:1202.2048 [hep-ex].
  * Chen and Braaten [1998] Y.-Q. Chen and E. Braaten, Phys.Rev.Lett. 80, 5060 (1998), arXiv:hep-ph/9801226.
  * Mo  _et al._ [2007] X.-H. Mo, C.-Z. Yuan, and P. Wang, Chin.Phys. C31, 686 (2007), arXiv:hep-ph/0611214.
  * Wang  _et al._ [2012] Q. Wang, G. Li, and Q. Zhao, Phys.Rev. D85, 074015 (2012), arXiv:1201.1681 [hep-ph].
  * Kivel [2023] N. Kivel, (2023), arXiv:2301.03884 [hep-ph].
  * Yan  _et al._ [2023] L.-W. Yan, Y.-H. Chen, C.-G. Duan, and Z.-H. Guo, (2023), arXiv:2301.03869 [hep-ph].
  * Ablikim  _et al._ [2022] M. Ablikim _et al._ (BESIII), Chin. Phys. C 46, 074001 (2022), arXiv:2111.07571 [hep-ex].
  * Ablikim  _et al._ [2020] M. Ablikim _et al._ (BESIII), Chin. Phys. C 44, 040001 (2020), arXiv:1912.05983 [hep-ex].
  * Szczepaniak and Pennington [2014] A. P. Szczepaniak and M. R. Pennington, Phys.Lett. B737, 283 (2014), arXiv:1403.5782 [hep-ph].
  * Guo  _et al._ [2010] P. Guo, R. Mitchell, and A. P. Szczepaniak, Phys.Rev. D82, 094002 (2010), arXiv:1006.4371 [hep-ph].
  * Guo  _et al._ [2012] P. Guo, R. Mitchell, M. Shepherd, and A. P. Szczepaniak, Phys.Rev. D85, 056003 (2012), arXiv:1112.3284 [hep-ph].
  * Khuri and Treiman [1960] N. N. Khuri and S. B. Treiman, Phys. Rev. 119, 1115 (1960).
  * Anisovich and Leutwyler [1996] A. V. Anisovich and H. Leutwyler, Phys. Lett. B375, 335 (1996), arXiv:hep-ph/9601237 [hep-ph].
  * Guo  _et al._ [2015a] P. Guo, I. Danilkin, and A. P. Szczepaniak, Eur. Phys. J. A 51, 135 (2015a), arXiv:1409.8652 [hep-ph].
  * Guo  _et al._ [2015b] P. Guo, I. V. Danilkin, D. Schott, C. Fernández-Ramírez, V. Mathieu, and A. P. Szczepaniak, Phys. Rev. D92, 054016 (2015b), arXiv:1505.01715 [hep-ph].
  * Colangelo  _et al._ [2017] G. Colangelo, S. Lanz, H. Leutwyler, and E. Passemar, Phys. Rev. Lett. 118, 022001 (2017), arXiv:1610.03494 [hep-ph].
  * Colangelo  _et al._ [2018] G. Colangelo, S. Lanz, H. Leutwyler, and E. Passemar, Eur. Phys. J. C78, 947 (2018), arXiv:1807.11937 [hep-ph].
  * Albaladejo and Moussallam [2017] M. Albaladejo and B. Moussallam, Eur. Phys. J. C77, 508 (2017), arXiv:1702.04931 [hep-ph].
  * Gasser and Rusetsky [2018] J. Gasser and A. Rusetsky, Eur. Phys. J. C78, 906 (2018), arXiv:1809.06399 [hep-ph].
  * Niecknig  _et al._ [2012] F. Niecknig, B. Kubis, and S. P. Schneider, Eur. Phys. J. C72, 2014 (2012), arXiv:1203.2501 [hep-ph].
  * Danilkin  _et al._ [2015] I. V. Danilkin, C. Fernández-Ramírez, P. Guo, V. Mathieu, D. Schott, M. Shi, and A. P. Szczepaniak, Phys. Rev. D91, 094029 (2015), arXiv:1409.7708 [hep-ph].
  * Albaladejo  _et al._ [2020a] M. Albaladejo, I. Danilkin, S. Gonzàlez-Solís, D. Winney, C. Fernández-Ramírez, A. Hiller Blin, V. Mathieu, M. Mikhasenko, A. Pilloni, and A. Szczepaniak, Eur.Phys.J. C80, 1107 (2020a), arXiv:2006.01058 [hep-ph].
  * Källén [1964] G. Källén,  _Elementary particle physics_ (Addison-Wesley, Reading, MA, 1964).
  * Kibble [1960] T. W. B. Kibble, Phys. Rev. 117, 1159 (1960).
  * Albaladejo  _et al._ [2022] M. Albaladejo _et al._ (JPAC), Prog. Part. Nucl. Phys. 127, 103981 (2022), arXiv:2112.13436 [hep-ph].
  * Kubis and Niecknig [2015] B. Kubis and F. Niecknig, Phys. Rev. D 91, 036004 (2015), arXiv:1412.5385 [hep-ph].
  * Stamen  _et al._ [2022] D. Stamen, T. Isken, B. Kubis, M. Mikhasenko, and M. Niehus, (2022), arXiv:2212.11767 [hep-ph].
  * Albaladejo  _et al._ [2020b] M. Albaladejo, D. Winney, I. Danilkin, C. Fernández-Ramírez, V. Mathieu, M. Mikhasenko, A. Pilloni, J. Silva-Castro, and A. Szczepaniak (JPAC), Phys. Rev. D 101, 054018 (2020b), arXiv:1910.03107 [hep-ph].
  * Bronzan and Kacser [1963] J. B. Bronzan and C. Kacser, Phys.Rev. 132, 2703 (1963).
  * Omnes [1958] R. Omnes, Nuovo Cim. 8, 316 (1958).
  * Froissart [1961] M. Froissart, Phys. Rev. 123, 1053 (1961).
  * Martin [1963] A. Martin, Phys. Rev. 129, 1432 (1963).
  * Pelaez  _et al._ [2019] J. R. Pelaez, A. Rodas, and J. Ruiz De Elvira, Eur. Phys. J. C 79, 1008 (2019), arXiv:1907.13162 [hep-ph].
  * Gonzàlez-Solís and Roig [2019] S. Gonzàlez-Solís and P. Roig, Eur. Phys. J. C 79, 436 (2019), arXiv:1902.02273 [hep-ph].
  * Garcia-Martin  _et al._ [2011] R. Garcia-Martin, R. Kaminski, J. R. Pelaez, J. Ruiz de Elvira, and F. J. Yndurain, Phys. Rev. D 83, 074004 (2011), arXiv:1102.2183 [hep-ph].
  * Albaladejo  _et al._ [2020c] M. Albaladejo, I. Danilkin, S. Gonzalez-Solis, D. Winney, C. Fernandez-Ramirez, A. N. H. Blin, V. Mathieu, M. Mikhasenko, A. Pilloni, and A. Szczepaniak (JPAC), Eur. Phys. J. C 80, 1107 (2020c), arXiv:2006.01058 [hep-ph].
  * Workman and Others [2022] R. L. Workman and Others (Particle Data Group), PTEP 2022, 083C01 (2022).
  * Blatt and Weisskopf [1952] J. M. Blatt and V. F. Weisskopf,  _Theoretical nuclear physics_ (Springer, New York, 1952).
  * Ablikim  _et al._ [2014] M. Ablikim _et al._ (BESIII), Phys. Rev. D 89, 092008 (2014), arXiv:1403.7042 [hep-ex].
  * Chen  _et al._ [2015] Y.-H. Chen, Z.-H. Guo, and B.-S. Zou, Phys. Rev. D 91, 014010 (2015), arXiv:1411.1159 [hep-ph].
  * Schneider  _et al._ [2012] S. P. Schneider, B. Kubis, and F. Niecknig, Phys. Rev. D 86, 054013 (2012), arXiv:1206.3098 [hep-ph].


# Electron pair emission from surfaces: Some general experimental
considerations

R. Kamrla, W. Widdra  Institute of Physics, Martin-Luther-Universität Halle-
Wittenberg, Von-Danckelmann-Platz 3, D-06120 Halle (Saale), Germany  C.-T.
Chiang  Institute of Atomic and Molecular Sciences, Academia Sinica, Taipei,
Taiwan  F. O. Schumann  Electronic mail<EMAIL_ADDRESS>Max-Planck-
Institut für Mikrostrukturphysik, Weinberg 2, 06120 Halle, Germany

(August 27, 2024)

###### Abstract

We discuss some experimental facets of electron pair emission from surfaces
using two different experimental approaches. In the first case the instrument
consists of a pair of hemispherical analyzers which are operated with
continuous primary beams of electrons or photons. The second instrument
employs a pair of time-of-flight spectrometers which require a pulsed
excitation source. A key experimental quantity is the ratio of ‘true’ to
‘random’ coincidences which can be determined in different ways. Regardless of
the type of instrument the primary flux has to adopt a much smaller value than
in single electron spectroscopy. We describe different approaches to obtain
the relevant count rates, in particular the concept of operating with a
delayed coincidence circuit. We also address the question on how to compare
the two types of spectrometer in terms of their performance.

###### pacs:

73.20.At, 79.60.-i

##  I Introduction

In the past decades photoelectron spectroscopy has developed into the most
powerful experimental technique to determine the detailed electronic
properties of matter. On the experimental side, a wide range of energy-,
momentum- and spin-resolving electron spectrometers has been employed. Here
the large majority of spectrometers records single electron spectra, e.g. the
kinetic energy distribution upon photon absorption. However, also the emission
of two electrons excited by absorption of a single photon is possible, which
might be due to a strong electron-electron correlation. This process can be
identified if two time-resolving electron detector are used. Both detectors
are connected via a coincidence logic. Only if both electrons hit the
detectors within an appropriate small relative time window, the two-electron
event is registered. In general, the coincidence rates are often rather small,
in the range of 0.1 to 10 events/s. Experimentally it is challenging to
discriminate these true coincidence events from the large number of single-
electron events originating from different photons.To achieve the
experimentally desired signal-to-noise ratios, adequate event statistics are
required. The resulting measurement time scales with the inverse of the
repetition rate of the pulsed excitation. Therefore high repetition rates are
favored, which are only limited by the requirement that the slowest electrons
have reached the detector within one excitation periode. The latter is
typically in the order of a few MHz and is fulfilled at synchrotron light
sources in special operation modes as well as with pulsed UV laser sources. In
Fig.1 electron pair emission due to single-photon absorption known as double
photoemission (DPE) is illustrated. An unwanted background signal arises, if
two photons are independently absorbed as this will also lead to the emission
of an electron pair. It is customary to refer to the latter process as
‘random’ coincidences whereas the genuine DPE signal is known as ‘true’
coincidence. As we discuss later the process of interest scales linearly with
the primary flux, while the undesirable process scales quadratically with the
flux. This unwanted pathway very quickly overwhelms the genuine signal and the
operation has to commence with much reduced flux. Let us quote a numerical
example for electron pair emission due to primary electron impact, known as
(e,2e). If one utilizes an experiment with energy dispersive elements, as
described later, the sample current is $10^{-13}$ A. This is many orders of
magnitudes lower than for Auger electron spectroscopy using a primary electron
beam for which primary currents are of the order of a few $\mu$A. The
requirement of low primary flux asks for a detection scheme with large
detection probability. This means that the detector solid angle should be as
large as possible while a large electron kinetic energy window should be
covered. There exist different solutions for experiments in the gas phase or
at surfaces.[1, 2, 3, 4, 5, 6, 7, 8]These schemes can be separated into energy
dispersive detection or time-of-flight (ToF) spectroscopy. The latter asks for
a pulsed excitation source which in the case of photon beams is more involved.
For efficient measurement scheme a repetition rate of the order of a few MHz
is needed, because once the slowest electron has been detected the next
excitation can commence. This requirement can be fulfilled at synchrotron
light sources in special operation modes, but the available beamtime is
restricted.More recent developments in laser technology using high-order
harmonic generation (HHG) allow to perform in-house DPE experiments.[9, 10,
11, 12, 13] While it is possible to carry out these type of experiments also
with an energy dispersive element together with a standard vacuum-ultraviolet
VUV light source, the available photon energies are limited to a few lines.
These sources employ the discharge of nobel gases of He, Ar and Ne.[14] We use
in our work a monochromatized He light source.[15, 16] For ToF spectroscopy,
HHG light sources provide more options.

Figure 1: Sketch of the emission of an electron pair via single photon
absorption (top) or absorption of two photons (bottom). The first case is the
process of interest and leads to true coincidences. The other pathway leads to
unwanted background commonly referred to as ‘random’ coincidences.

It is well-established on how to determine the individual contribution of
‘random’ and ‘true’ coincidences. Here we discuss a novel approach for ToF-
spectroscopy similar to a recently reported concept in Auger-photoelectron
coincidences exploring synchrotron radiation on surfaces, for a recent
overview about coincidence spectroscopy we refer to the literature.[17, 18] In
this work we compare this with the common approach of a double logarithmic
analysis. We will demonstrate that both approaches yield the same
characteristic parameters for different DPE experiments in which the intensity
levels are orders of magnitude different.

Figure 2:  We display an example of the sum of ‘true’ and ‘random’ coincidence
count rates as a function of single photoemission count rate in a double
logarithmic manner. The arrow indicates the singles rate for which the number
of ‘true’ and ‘random’ coincidences are equal.

We start by considering how a genuine pair emission can be identified and how
the statistical probabilities for ‘true’ and ‘random’ events have to be
formulated. Second, we introduce the concept of delayed coincidences and how
the relevant quantities can be experimentally determined. Third, we provide
experimental data to prove the concept and how it can be related to a known
procedure. Finally, we provide an outlook on how the ‘random’ events can be
removed from coincidence event spectra.

##  II signature of ‘true’ coincidences

Although our discussion will be valid for a variety of experimental
approaches, we will present data obtained from two different set-ups for the
investigation of solid surfaces. Our ToF instrument consists of a pair of
lenses. The lens axes include an angle of 90∘ and have an angular opening of
$\pm$15∘.[19] As excitation sources a pulsed electron gun or an HHG light
source is available. The light source can operate at a repetition rate of
0.2-1 MHz and provides a photon energy range from 13 to 39 eV.[10, 13, 12] The
energy dispersive experiment utilizes a pair of hemispherical analyzers.[20]
Also in this case the two electron-optical axes include an angle 90∘ with an
angular acceptance of $\pm$15∘ within the scattering plane. This instrument
has proven to be very versatile. We reported on experiments with primary
electrons, positrons, He2e+ ion and photons from laboratory and synchrotron
sources.[21, 22, 23, 24, 25]Let us define the probability to find one photon
in a light pulse by $P_{1}(\lambda)$ while the probability $P_{2}(\lambda)$
describes the probability to find two photons in a pulse. The quantity
$\lambda$ is the average number of photons in a pulse and hence proportional
to the primary flux. If we assume a Poisson distribution and $\lambda\ll 1$ we
obtain the following expressions:

$P_{1}(\lambda)=\lambda\leavevmode\nobreak\ \leavevmode\nobreak\ ,\\\
\leavevmode\nobreak\ \leavevmode\nobreak\
P_{2}(\lambda)=\frac{1}{2}\lambda^{2}$ (1)

The key result is that the probability $P_{2}(\lambda)$ goes quadratically
with $\lambda$ and therefore with the flux, while the term $P_{1}(\lambda)$
displays a linear relation. This explains the flux dependence for ‘random’ and
‘true’ coincidences, respectively.This difference in the flux dependence can
be observed by presenting the coincidence rate as a function of singles rate
in a double logarithmic manner, see Fig.2. The dashed lines indicate the
behavior for low and high flux. In a simplified analysis one can evaluate
these regions separately. The slopes are close to 1 and 2 indicative of the
expected behavior. The intercept of those lines occurs at a primary flux at
which the number of ‘true’ and ‘random’ coincidences are equal. We will
present below an improved procedure to analyze the data. The key point is that
for this study it is proven that‘true’ coincidence exist. Furthermore, one can
adjust the primary flux such that the ratio of true-to-random coincidences
(termed TR ratio) is at an acceptable level. The usual wisdom is to aim for a
value above 1. This type of analysis can be used for energy dispersive
experiments and ToF spectrometer. Implicit is the assumption that the primary
flux can be increased to such a high level that ‘random’ coincidences dominate
the count rate. While this is in many cases not a constraint it can be an
issue. In such a case the log-log plot can not be generated.In the pioneering
experiment of Bothe and Geiger a different way appropriate for energy
dispersive experiments was discussed.[26] It is based on the measurement of
the arrival time difference of the two emitted electrons at their respective
detectors, see Fig.3. The implementation for an instrument which consists of a
pair of hemispherical analyzers has been discussed in detail, here we recall
only the conceptional points. [27, 28, 20] For a valid coincidence event the
arrival time ($t_{left}$ and $t_{right}$) of two electrons at the respective
detector with respect to the coincidence trigger is known. Therefore we can
plot the coincidence intensity as a function of the time difference
$dt$=$t_{left}-t_{right}$. A schematic curve is shown in Fig.3. The emergence
of a peak is evidence of ‘true’ coincidences as discussed in the
literature.[26, 1, 29, 30, 31, 2, 32]. There is no temporal relation between
the ‘random’ coincidences which explains the constant intensity outside the
peak region. Apart from proving the existence of ‘true’ coincidences this
approach allows the straightforward determination of the TR ratio of ‘true’ to
‘random’ events.

Figure 3: Sketch of the arrival time histogram. The emergence of the peak is
proof of the emission of electron pairs. The width $t_{true}$ is a consequence
of the time dispersion of the spectrometer. The histogram has a total width of
$t_{tot}$.

The total width of the $dt$ curve is given by $t_{tot}$. The base width of the
peak is $t_{true}$ as indicated by the pair of vertical dashed lines, see
Fig.3. This result is used to define a region of interest where ‘true’
coincidences are located. The peak width is in accordance with the time
dispersion of the spectrometer.[33, 34, 35] The generation of a log-log plot
with suitable statistics takes a considerable time and will interrupt the
actual measurement, because largely different primary fluxes have to be
used.In this respect the use of a histogram is very powerful. We can
illustrate this point by observing how the histogram evolves. For this we take
a data file from a DPE experiment from a Pb surface. The experiment in
question captured the coincidence of a $5d$ core electron and related Auger
electron upon excitation with 40.8 eV photons, see Fig.8 of Aliaev et al.[36].
The coincidence rate was around 1 cps for this study.

The curve of Fig.4 (a) shows the histogram after data acquisition of 200
coincidence events which for the given coincidence rate translates to 200 s.
Even after this short time it is clear that ‘true’ coincidences exist and the
TR ratio can be computed. The evolution of the TR ratio as a function of the
counted events is displayed in Fig.4 (b) and one can see that the value of the
TR ratio is accurately known after around 2000 counts.

At the beginning of a new study the required intensity level for a
sufficiently good TR ratio is not known, hence the primary flux needs to be
varied. For a log-log plot the picture emerges only after a sufficiently large
flux window of 2-3 orders of magnitude, see Fig.2. This is different for the
histogram measurements where a single measurement gives a result.

Figure 4:  Panel (a) shows the arrival time histograms for 200 coincidence
counts. The evolution of the TR ratio as a function of the accumulated counts
is plotted in panel (b). The data are from a measurement from a Pb surface,
see Fig.8 (a) in Aliaev et al[36].

##  III delayed coincidence for ToF

Our ToF experiment is not set-up to record two electrons per photon pulse on
one of the MCP detectors. This limitation is not as severe as it appears,
because previous work has shown that the emission of two electrons into the
same angular range is strongly reduced. This was related to an important
concept of solid state theory known as exchange-correlation hole.[37, 38, 39,
40, 41, 42] The angular range of the reduced intensity region is of the order
of 60∘ and exceeds the angular acceptance of the entrance lenses which is $\pm
15^{\circ}$.It is highly desirable to be able to perform an equally fast
determination of the TR ratio for a ToF experiment. This becomes possible if
we employ a second coincidence circuit. For simplicity and also reflecting our
experimental situation we assume that the time $t_{rep}$ between two
subsequent primary pulses is longer than the time-of-flight of the slowest
detectable electron. In the first circuit the two multi-channelplate (MCP)
signals of the different spectrometers are used. The second circuit, called B,
receives one MCP signal which is delayed by the time $t_{rep}$ between two
consecutive primary pulses. A schematic view of the two coincidence circuits
is seen in Fig.5. While circuit A records ‘true’ and ‘random’ coincidences
circuit B can only record ‘random’ events. As indicated by Fig.1 we define a
‘random’ event as being caused by two independent photons. This is clearly the
case for circuit B.

Figure 5: Schematic signals from two detectors. The dashed vertical lines
indicate the time when photons from the different pulses arrive at the sample.
The time between subsequent pulses is given by $t_{rep}$. Circuit A triggers
if both detectors supply a signal after photon absorption at the first pulse
and display temporal overlap. Circuit B triggers if the temporal overlap
exists between the delayed signal of detector 1 and a signal from detector 2
from the next pulse.

The probability to find two photons in the same pulse is given by
$P_{2}(\lambda)$. This is in contrast to find one photon in one pulse while
the second photon is in the subsequent pulse. For this scenario the
probability is given by the joint product of probabilities
$P_{\lambda}(1)^{2}$. For $\lambda\ll 1$ we can simplify
$P_{\lambda}(1)^{2}=2\cdot P_{\lambda}(2)$ by using the expansions of Eq.(1).

In the case of a ‘random’ coincidence, circuit A will record the first
electron either by detector 1 or 2. The second electron is then registered in
the opposite detector. In other words there exists two detection combinations
in contrast to circuit B. This means that the ‘random’ rate for both circuits
is the same and we can write:

$c_{A}=r+t\\\ ,\quad c_{B}=r$ (2)

The coincidence rates for circuits A and B are given by $c_{A}$ and $c_{B}$,
while the terms $r$ and $t$ refer to the ‘random’ and ‘true’ coincidence
rates. This terms are low count rate approximations as discussed in the
literature.[43] With this we define $ratio$ as:

$ratio=\frac{c_{A}}{c_{B}}=\frac{r+t}{r}=1+\frac{t}{r}$ (3)

We immediately see that $ratio$ is closely related to the quantity of
interest, namely the TR ratio. We recall that the ‘random’ rate varies
quadratic with flux while the ‘true’ events scale linearly with the flux.
Since the singles rate also depends linearly on the flux, a presentation of
$ratio$ as a function of the inverse singles rate is expected to yield a line
with an intercept of 1. If we label the singles rate with $s$ we can
formulate:

$ratio=a_{1}+a_{2}\cdot\frac{1}{s}$ (4)

Figure 6:  In (a) the measurement of the $ratio$ as a function of inverse
singles rate is shown. These are DPE data come from a SrRuO3 surface with a
photon energy of 30 eV. The solid line is a linear fit to the data and reveals
an intercept $a_{1}$=1.03$\pm 0.35$. The horizontal dashed line indicates
$ratio=2$ which is equivalent to TR=1. In panel (b) we plot the delayed
coincidence rate as a function of the singles rate via a log-log plot. We find
for the slope $b_{2}$=1.98$\pm 0.04$.

We expect the intercept $a_{1}$ to be close to 1 according Eq.(3). In this
case the slope $a_{2}$ can be identified with the singles rate $s$ for which
ratio is close to 2 or equivalently the TR ratio is 1.A single photon pulse
will lead with a probability $p_{s}$ to the detection of an electron in one
spectrometer. The time between subsequent pulses is given by $t_{rep}$.
Therefore the single electron rate is :

$s=\frac{p_{s}}{t_{rep}}$ (5)

Similarly we get for the delayed coincidence rate $c_{b}$:

$c_{b}=\frac{p_{s}^{2}}{t_{rep}}=s^{2}\cdot t_{rep}$ (6)

If we take the logarithm of both sides we obtain:

$log(c_{b})=log(t_{rep})+2\cdot log(s)=b_{1}+b_{2}\cdot log(s)$ (7)

This means that the fitting of the data in a double logarithmic fashion will
yield fit parameters $b_{1}$ and $b_{2}$. For our experiment with a
$t_{rep}=2\cdot 10^{-6}$ sec we expect $b_{1}$=-5.699 and $b_{2}$=2. In order
to explore the validity of the previous statements we have performed
experiments using a ToF instrument with an HHG light source as described
previously.[19, 13, 12] The spectrometer settings were identical to those
reported in our previous work. We selected a SrRuO3 sample as target which we
investigated with a few different photon energies.In Fig.6 (a) we show the
behavior of $ratio$ as a function of the inverse singles rate $s^{-1}$ for a
photon energy of 30 eV. It is apparent that the data are well-described by a
line which has an intercept of $a_{1}$=1.03$\pm$0.35 in agreement with the
expected value of 1. The horizontal dashed line marks the level for $ratio$=2
or TR=1. For the other parameter we find $a_{2}$=2013.5 cps which is singles
rate for which the TR ratio is near 1. In Fig.6 (b) we show the delayed
coincidence rate as a function of the singles rate in double logarithmic
manner. A linear fit reveals a slope of $b_{2}$=1.98$\pm 0.04$ which is within
the error of the expected value for ‘random’ coincidences. The vertical dashed
line marks the singles rate $s=a_{2}$ and the crossing point with the solid
line occurs for a delayed coincidence rate of 7.87 cps which is identical to
the ‘true’ rate. We obtain a value of the fit parameter $b_{1}$=-5.6479 which
is within the expected range.

Figure 7: We summarize the flux dependent measurements for a variety of photon
energies. In the panels (a)-(c) we show the variation of the fit parameters
$a_{1}$, $b_{1}$ and $b_{2}$ with their error bars. The dashed horizontal
lines in each plot are the expected values as discussed in the text.

We summarized our findings for a variety of photon energies in Fig.7. The
three panels show the variation of the fit parameters $a_{1}$, $b_{1}$ and
$b_{2}$ together with the expected values as indicated by the dashed
horizontal lines in each plot. The blue data points for $a_{1}$ of Fig.7(a)
are obtained if only the data with highest primary flux are retained. In this
case $a_{1}$=1 is almost exactly fulfilled. The key point is that the material
independent parameters are essentially in agreement with the prediction.

In the description so far we have employed 4 parameters which can be
determined by fitting the experimental data. It is appropriate to perform
calibration measurements with the aim to keep some parameters fixed. The
intercept $a_{1}$ from Eq.(4) can be easily measured by operating at a high
singles rate. This is best done at lower photon energies, because there we
observe a lower ‘true’ rate for the same singles rate. For the evaluation of
the ‘true’ and singles rate for TR=1 we will set $a_{1}$=1 and $b_{2}$=2. In
other words in our analysis we adjust only two parameters. We compare the
results with the usual double logarithmic analysis of the coincidence rate
$r_{a}$ versus the singles rate $s$ formulated as:

$log(r_{a})=log(a\cdot s+b\cdot s^{2})$ (8)

The parameters $a$ and $b$ refer to the contribution of ‘true’ and ‘random’
coincidences, respectively. Once the parameters $a$ and $b$ are obtained in
the analysis it is straightforward to determine the singles and ‘true’ rate
for TR=1. These are given by $s=a/b$ and $t=a^{2}/b$. The outcome of this
analysis and comparison with the result from the delayed coincidence is
plotted in Fig.8. It is clearly shown that both ways of analysis give
essentially the same result. Furthermore, the variation of the ‘true’ rate
covers more than 2 orders of magnitude and starts at a low rate of 0.02 cps
for a photon energy of 15 eV. The demonstration of consistent results clearly
proves that the delayed coincidence is a valid approach. Therefore, the
evaluation of the TR ratio without the need of a flux variation is possible.
Hence monitoring the TR ratio during a measurement similar to those for the
energy dispersive set-up shown in Fig.4 is possible.

Figure 8: Photon energy dependence of the ‘true’ (a) and singles (b) count
rate for TR=1. The black data points refer to the result from the delayed
coincidences while the red data points stem from a double-logarithmic
analysis.

##  IV Removal of ‘random’ background

Besides the knowledge of the ‘random’ contribution to the count rate it is
also beneficial to be able to remove the ‘random’ background from the energy
spectra without the need to collect an additional spectrum at very high
primary flux. This possibility exists for the energy dispersive instrument and
we have discussed our implementation in the literature.[20, 25] It turns out
that this can also be achieved with a ToF instrument as shown in a recent work
on APECS at a synchrotron light source.[17] Effectively the concept of the
delayed coincidence was used which was facilitated by the possibility to
record the absolute time for each photon bunch. One can make use of a delayed
coincidence circuit to get the equivalent information. The key point is to
change the operation of the coincidence circuit and to add an additional
channel in which the signal level indicates the presence of a delayed
coincidence. Rather than using the AND operation with the two MCP signals one
selects the OR option. Therefore, each time at least one spectrometer records
an electron impact a trigger is generated. At this point the detector signals
are read and stored in a list as for the usual AND operation. In Table 1 we
provide the schematic data structure. The 1st column is the recorded event
number. The 2nd and 3rd columns indicate whether the respective spectrometer,
labeled as $left$ and $right$, has kinetic energy and emission angle
information. Consequently the entry is either $yes$ or $no$ for each of these
two columns.

Table 1: Schematic data structure. event | left spect | right spect | coin | delayed coin
---|---|---|---|---
1 | yes | yes | yes | no
2 | yes | no | no | no
3 | no | yes | no | yes
4 | no | yes | no | no
5 | yes | no | no | no
6 | no | yes | no | yes
.. | .. | .. | .. | ..

If for a given event both spectrometer registered an electron this is a
coincidence event which gives the entry $yes$ in the column $coin$. This is
the case for the 1st event in Table 1. This means that the usual AND condition
of the coincidence circuit can be recovered in a post-experiment analysis and
the energy spectrum containing both ‘true’ and ‘random’ coincidences can be
computed.Alongside the coincidence circuit which was changed to an OR
operation we will operate a delayed coincidence circuit in the AND mode as
sketched in Fig.5 as circuit B. This provides the signal for the last column
of the table. It has an entry if a delayed coincidence circuit has triggered.
This is the case for event 3 where only the $right$ spectrometer records an
electron. However, the $left$ spectrometer has an entry for event 2. This
means that the events 2 and 3 originate from two photon pulses separated
exactly by the time $t_{rep}$. This means that we can in a post-experiment
analysis identify a delayed ‘random’ coincidence and compute the energy
spectrum. If we remove this spectrum from the one containing ‘true’ and
‘random’ we obtain effectively the ‘true’ spectrum. This procedure is allowed,
because we have established above that the number of ‘random’ coincidence in
the regular coincidence circuit is identical to those in the delayed circuit.

##  V revised definition of $\lambda$

In the discussion so far we have introduced the primary intensity via the
average number of photons per pulse $\lambda$. However, we are talking about
electron detection, hence we have to introduce a parameter which describes the
probability for electron emission upon photon absorption and detection within
a given spectrometer. This will be a material and photon energy dependent
entity which we want to abbreviate with $y$ for yield. The focus are ‘random’
events which are mainly determined by single electron emission, hence we
ignore double electron emission at this moment. The value of $y$ will depend
also on the detailed spectrometer settings which determines the solid angle of
detection and the probed energy window. This will lead to a singles count rate
$s$ which is easily determined. Conceptually the photon flux can be measured
in absolute units and with the known time between pulses $t_{rep}$ the
relation between $\lambda$ and $y$ reads as:

$y=\frac{s\cdot t_{rep}}{\lambda}$ (9)

This means that the average number of detected electrons per pulse is simply
$y\cdot\lambda$. In a next step one would like to know the probability
function for emitted electrons. In other words what is the probability that
$m$ electrons are detected. For this evaluation one needs first the
distribution function for $n$ photons in a pulse, which we have abbreviated
with $P_{n}(\lambda)$ in Eq.(1). Further, we do not assume that $\lambda\ll
1$. If $n$ photons are absorbed by the sample, the number of electrons $m$
registered by one spectrometer is $m\leq n$. The probability for this to occur
is determined by the binomial distribution:

$el_{n,m}(y)=\binom{n}{m}\cdot y^{m}\cdot(1-y)^{n-m}$ (10)

This finally leads to the probability for the detection of $m$ electrons:

$P_{m}^{elec}(\lambda,y)=\sum_{i=1}^{n}P_{i}(\lambda)\cdot el_{n,m}(y)$ (11)

If the photon distribution is a Poisson distribution then one can show that
$P_{m}^{elec}$ is also a Poisson distribution with a new parameter
$\lambda_{eff}=\lambda\cdot y$. Hence we need to demand that $\lambda_{eff}\ll
1$ in order to derive the equivalent expressions of Eq.(1). From Fig.6 we can
see that the highest singles rate was 15000 cps. With $t_{rep}$=2000 ns we
obtain for $\lambda_{eff}=$0.03. Clearly the condition $\lambda_{eff}\ll 1$ is
fulfilled in our DPE measurements and the approximations made in Eq.(1) are
valid.In the discussion of the delayed coincidence we have ignored the
possibility that detector 1 has a ‘true’ counter part in detector 2. This will
lead to additional events, however the ‘true’ coincidence rate is more than 2
orders of magnitude smaller than the singles rate, see Fig.8. Hence this
contribution can be safely ignored.This observation is not tied to the current
DPE measurement on SrRuO3, but is a rather general result in our activities on
surfaces. This includes also observations from Auger-Photoelectron
coincidences. This in turn means that the consideration of $y$ to contain only
single electron events is justified. This should not be read as an indication
that electron pair emission from surfaces is a weak effect. A combined (e,2e)
and DPE study revealed that 15-40$\%$ of electrons detected by one
spectrometer have a counterpart emitted somewhere in the halfspace.[44]

##  VI Comparison energy dispersive versus ToF performance

A variety of coincidence spectrometer have been developed addressing different
scientific questions. There is no rule which determines the best approach.
Here we only want to provide a simple estimate on the performance relation
between an energy dispersive and ToF approach. For simplicity we assume that
the emission of pairs is isotropic and constant within an energy window. We
assume at the moment that the energy dispersive and ToF instrument cover the
same solid angle $\Omega$ and energy window. Later we will discuss the
deviations from this description.An example resembling this situation is a
pair of hemispherical analyzers compared to a pair of ToF spectrometers which
are based on the entrance lens of a hemispherical analyzer. In our activities
these approaches have been realized.[45, 20, 19, 13, 12] The key number
characterizing the performance of the hemispheres is the time dispersion. This
is essentially the quantity $t_{true}$, see Figs.3 and 4. A ToF instrument
requires a pulsed excitation source and the time between subsequent pulses is
$t_{rep}$. Let us assume that the average number of primary particles per
pulse is given by $\alpha$ which yields the desired ratio of ‘true’ to
‘random’ . This means that the primary flux for the ToF experiment is:

$f_{ToF}=\frac{\alpha}{t_{rep}}$ (12)

For a proper comparison with an energy dispersive instrument we have to allow
the same number of primary particles within the coincidence window. The
temporal width of this is characterized by the value of $t_{true}$, see Fig.3.
This will ensure the same ratio of ‘true’ to ‘random’ and we obtain for the
primary flux $f_{disp}$:

$f_{disp}=\frac{\alpha}{t_{true}}$ (13)

This allows now to compare the primary flux of the two types of experiments
and we obtain:

$\frac{f_{disp}}{f_{ToF}}=\frac{t_{trep}}{t_{true}}$ (14)

Let us use some typical values for our instruments which are $t_{true}$=20 ns
and $t_{rep}$=1000 ns. This tells us that the energy dispersive system is a
factor of 50 more efficient or equivalently a higher coincidence rate by this
factor is possible.However, the assumption of the same solid angle for both
types of spectrometer is not valid. While our ToF instrument captures 3$\%$ of
the half sphere this value is reduced to 0.7$\%$ for the hemisphere set-up.
Since the solid angle enters quadratic into the detection efficiency it
reduces the performance advantage to about 3.The energy window captured by the
ToF and energy dispersive instrument can be selected. The requirement of the
operation with a low primary flux makes asks for efficient detection hence a
sizable window is required at the expense of energy resolution. In the context
of this work we operated the hemispherical analyzers with a pass energy of 150
eV which provides an energy window of 13.5 eV. The typical settings of the ToF
instrument for coincidence measurements have been described.[46] It is
possible to cover an energy window of 30 eV, but only over a range of 10 eV
the spectrometer captures electrons with the maximum solid angle of 3$\%$.
Outside this range the solid angle drops significantly. This means that the
detection efficiency for a ToF is about a factor of 1.25 larger than for the
hemisphere. For a coincidence experiment one has to consider the square of
this value and this gives a factor 1.56. From this one concludes that the
performance advantage of the hemisphere reduces to a factor 2. In this
estimate other experimental considerations have been ignored, but the main
conclusion is that both approaches provide performance on equal terms.

##  VII summary

For coincidence experiment it is vital to know the ratio of ‘true’ to ‘random’
events which is sensitive to the primary flux. We reviewed the established
procedures on how to identify ‘true’ electron pair emission. A powerful tool
known for energy dispersive coincidence experiments was found to have its
counter part in ToF instrumentation. This approach we termed delayed
coincidence circuit. We discussed the relevant parameters of this procedure
and how they can be determined in an experiment. In order to support this
notion we performed photon energy dependent DPE studies. These resulted in
intensity levels varying by more than 2 order of magnitude. The characteristic
values of the delayed coincidence are in agreement with the predicted behavior
and with the concurrent obtained double-logarithmic presentation.This
experimental proof of the concept allows to monitor the TR ratio during the
actual measurement. Therefore there is no need to resort to the time consuming
flux dependency measurement for the conventional double-logarithmic analysis
which interrupts the actual measurement. Furthermore, it is not required to
have such a high primary flux in which the coincidence rate is dominated by
‘random’ coincidences.We have also outlined how it becomes possible to remove
the ‘random’ contribution from the spectra by using the novel approach.

## References

  * [1] E. Jensen, R. A. Bartynski, S. L. Hulbert, and E. Johnson,  Rev. Sci. Instrum. 63, 3013 (1992).
  * [2] S. Thurgate,  urf. Interface Anal. 20, 627 (1993).
  * [3] R. Herrmann, S. Samarin, H. Schwabe, and J. Kirschner,  Phys. Rev. Lett. 81, 2148 (1998).
  * [4] R. Gotter, A. Ruocco, A. Morgante, D. Cvetko, L. Floreano, F. Tommasini, and G. Stefani,  Nucl. Instrum. Methods A 467-468, 467 (2001).
  * [5] J. Ullrich, R. Moshammer, A. Dorn, R. Dörner, L. P. H. Schmidt, and H. Schmidt-Böcking,  Rep. Prog. Phys. 66, 1463 (2003).
  * [6] M. Hattass, T. Jalowy, A. Czasch, T. Weber, T. Jahnke, S. Schössler, L. P. Schmidt, O. Jagutzki, R. Dörner, and H. Schmidt-Böcking,  Rev. Sci. Instrum. 75, 2373 (2004).
  * [7] J. H. D. Eland, O. Vieuxmaire, T. Kinugawa, P. Lablanquie, R. I. Hall, and F. Penent,  Phys. Rev. Lett. 90, 053003 (2003).
  * [8] F. Penent, P. Lablanquie, R. I. Hall, J. Palaudoux, K. Ito, Y. Hikosaka, T. Aoto, and J. H. D. Eland,  J. Electron. Spectrosc. Relat. Phenom. 144-147, 7 (2005).
  * [9] C.-T. Chiang, A. Blättermann, M. Huth, J. Kirschner, and W. Widdra,  Appl. Phys. Lett. 101, (2012).
  * [10] C.-T. Chiang, M. Huth, A. Trützschler, M. Kiel, F. O. Schumann, J. Kirschner, and W. Widdra,  New J. Phys. 17, 013035 (2015).
  * [11] C.-T. Chiang, M. Huth, A. Trützschler, F. O. Schumann, J. Kirschner, and W. Widdra,  J. Electron. Spectrosc. Relat. Phenom. 200, 15 (2015).
  * [12] C.-T. Chiang, A. Trützschler, M. Huth, R. Kamrla, F. O. Schumann, and W. Widdra,  Prog. Surf. Sci. 95, 100572 (2020).
  * [13] A. Trützschler, M. Huth, C.-T. Chiang, R. Kamrla, F. O. Schumann, J. Kirschner, and W. Widdra,  Phys. Rev. Lett. 118, 136401 (2017).
  * [14] G. Schönhense and U. Heinzmann,  J. Phys. E 16, 74 (1983).
  * [15] Light source MBS L-1 from MB Scientific.
  * [16] Monochromator VUV5040 from VG Scienta.
  * [17] T. Leitner, A. Born, I. Bidermane, R. Ovsyannikov, F.O.L. Johansson, Y. Sassa, A. Föhlisch, A. Lindblad, F.O. Schumann, S. Svensson and N. Mårtensson,  J. Electron. Spectrosc. Relat. Phenom. 250, 147075 (2021).
  * [18] T. Arion and U. Hergenhahn,  J. Electron. Spectrosc. Relat. Phenom. 200, 222 (2015).
  * [19] M. Huth, C.-T. Chiang, A. Trützschler, F. O. Schumann, J. Kirschner, and W. Widdra,  Appl. Phys. Lett. 104, 061602 (2014).
  * [20] F. O. Schumann, R. S. Dhaka, G. A. van Riessen, Z. Wei, and J. Kirschner,  Phys. Rev. B 84, 125106 (2011).
  * [21] I. S. Brandt, Z. Wei, F. O. Schumann, and J. Kirschner,  Phys. Rev. Lett. 113, 107601 (2014).
  * [22] Z. Wei, F. O. Schumann, C. H. Li, L. Behnke, G. Di Filippo, G. Stefani, and J. Kirschner,  Phys. Rev. Lett. 113, 267603 (2014).
  * [23] G. Di Filippo, F. O. Schumann, S. Patil, Z. Wei, G. Stefani, G. Fratesi, M. I. Trioni, and J. Kirschner,  J. Electron. Spectrosc. Relat. Phenom. 211, 32 (2016).
  * [24] C. H. Li, C. Tusche, F. O. Schumann, and J. Kirschner,  Phys. Rev. Lett. 118, 136402 (2017).
  * [25] F. O. Schumann, Y. Aliaev, I. Kostanovskiy, and J. Kirschner,  Phys. Rev. B 101, 115104 (2020).
  * [26] W. Bothe and H. Geiger,  Zeitschrift für Physik A Hadrons and Nuclei 32, 639 (1925).
  * [27] G. A. van Riessen, F. O. Schumann, M. Birke, C. Winkler, and J. Kirschner,  Journal of Physics: Conference Series 185, 012051 (2009).
  * [28] G. A. van Riessen, Z. Wei, R. S. Dhaka, C. Winkler, F. O. Schumann, and J. Kirschner,  J. Phys.: Condens. Matter 22, 092201 (2010).
  * [29] U. Amaldi Jr., A. Egidi, R. Marconero, and G. Pizzella,  Rev. Sci. Instrum. 40, 1001 (1969).
  * [30] I. E. McCarthy and E. Weigold,  Physics Reports 27, 275 (1976).
  * [31] P. Hayes, M. Bennett, J. Flexman, and J. Williams,  Rev. Sci. Instrum. 59, 2445 (1988).
  * [32] J. Kirschner, O. M. Artamonov, and A. Terekhov,  Phys. Rev. Lett. 69, 1711 (1992).
  * [33] M. Völkel and W. Sandner,  J. Phys. E 16, 456 (1983).
  * [34] R. E. Imhof, A. Adams, and G. C. King,  J. Phys. E 9, 138 (1976).
  * [35] O. Kugeler, S. Marburger, and U. Hergenhahn,  Rev. Sci. Instrum. 74, 3955 (2003).
  * [36] Y. Aliaev, I. Kostanovskiy, J. Kirschner, and F. O. Schumann,  Surf. Sci. 677, 167 (2018).
  * [37] F. O. Schumann, J. Kirschner, and J. Berakdar,  Phys. Rev. Lett. 95, 117601 (2005).
  * [38] F. O. Schumann, C. Winkler, G. Kerherve, and J. Kirschner,  Phys. Rev. B 73, 041404(R) (2006).
  * [39] M. Hattass et al.,  Phys. Rev. B 77, 165432 (2008).
  * [40] F. O. Schumann, C. Winkler, and J. Kirschner,  Phys. Rev. Lett. 98, 257604 (2007).
  * [41] F. O. Schumann, C. Winkler, J. Kirschner, F. Giebels, H. Gollisch, and R. Feder,  Phys. Rev. Lett. 104, 087602 (2010).
  * [42] F. O. Schumann and J. Kirschner,  Phys. Rev. B 103, 085435 (2021).
  * [43] H. Kossmann,  Measurement Science and Technology 4, 16 (1993).
  * [44] F. O. Schumann, Y. Aliaev, I. Kostanovskiy, G. Di Filippo, Z. Wei, and J. Kirschner,  Phys. Rev. B 93, 235128 (2016).
  * [45] G. A. van Riessen, F. O. Schumann, M. Birke, C. Winkler, and J. Kirschner,  J. Phys.: Condens. Matter 20, 442001 (2008).
  * [46] M. Huth, A. Trützschler, C.-T. Chiang, R. Kamrla, F. O. Schumann, and W. Widdra,  J. Appl. Phys. 124, 164504 (2018).