#  $0\nu\beta\beta$ to the first $2^{+}$ state with two-nucleon mechanism for
L-R symmetric model

Dong-Liang Fanga,b and Amand Faesslerc aInstitute of Modern Physics, Chinese
Academy of Science, Lanzhou, 730000, China  bUniversity of Chinese Academy of
Sciences, Beijing, 100049,China  cInstitute for theoretical physics, Tübingen
University, D-72076, Germany

###### Abstract

We develop the formalism for calculating the decay rate of neutrinoless double
beta decay to the $2^{+}$ excited states within L-R symmetric model. We
consider the effects from induced hadronic currents up to NLO. The QRPA method
in a spherical basis is adopted for the nuclear many-body calculation and the
corresponding nuclear matrix elements are given. Also, the phase space factors
are obtained with numerical electron wave functions. Our results suggest that
the nuclear matrix elements are nucleus dependent and they are generally
smaller than that of the decay to the ground states. And finally, we give a
naive analysis of how current experiment data constrains the L-R symmetric
model.

###### pacs:

14.60.Lm,21.60.-n, 23.40.Bw, 23.40.-s

##  I Introduction

New physics beyond the standard model is always the hottest topic in particle
physics. And the origin of the neutrino mass could be the one among the most
important questions in this area. Perhaps, the most promising explanation for
a small neutrino mass is the see-saw mechanism. The see-saw mechanism can
usually be divided into different categories. To realize such mechanisms,
different new physics models are proposed. One of the appealing proposals is
the L-R symmetric model MS81  with the underlying gauge symmetry
SU(2)${}_{L}\times$SU(2)${}_{R}\times$U(1)B-L. With two more Higgs particles
introduced, one can naturally incorporate the seesaw mechanism into this
model. There are different phenomenologies related to this model. One of the
important consequences is the existence of the so-called neutrinoless double
beta decay ($0\nu\beta\beta$-decay).

However, the underlying mechanism of $0\nu\beta\beta$-decay may not be
unique.Thus, it is important to find a way to identify the decay mechanisms by
future measurements. There are several proposals for this purpose, such as
comparing the ratios of the decay rates for different candidates LRS15  or
measuring the spectra and the angular correlations of the emitted electrons
GLS22 . As an alternative, one also suggests to compare the decay rates of
decays to ground and excited states DLZ11 , especially between the ground
states and the $2^{+}$ excited states GLS22 .

In minimal L-R symmetric model, for the decay to ground state, due to neutrino
propagator’s helicity suppression by the mass mechanism, the non-helicity-
suppressed $q$-term may play a dominant role DKT85 . Under such scenario, the
emission of P-wave electrons will surely lead to visible effect on the angular
correlation of the double $\beta$ spectrum, this is investigated in SDS15 .
Also, for another decay mode – neutrinoless double beta decay to the $2^{+}$
excited states (Hereafter $0\nu\beta\beta(2^{+})$), contributions from such
P-wave electrons will become dominant. While different to the decay to ground
states, the helicity suppressed $m_{e}$ terms are negligible since these terms
come from the NLO parts of the hadronic current. So with such new physics
model, we may probe the underlying mechanism by comparing these different
decay modes. Unfortunately, $0\nu\beta\beta\beta(2^{+})$ is actually rarely
investigated, and the only calculation available is done with the Projected
HFB approach Tom88 . PHFB calculation suggests that the decay is highly
suppressed since the NME is several orders of magnitude smaller than that of
the decay to ground states (Hereafter $0\nu\beta\beta(0^{+})$). This smallness
is caused partly by the suppression from each part of the NME and partly by
the cancellation among them. On the other hand, in Tom99 , one finds that
although the neutrino mass mechanism can contribute to this decay mode, their
NME is about two order of magnitude smaller than the $q$ mechanism. These
together make it impossible to observe $0\nu\beta\beta(2^{+})$.

Nevertheless, our recent calculations FF21  suggest that the NME in Tom88  is
underestimated, our results are orders of magnitude larger than that of Tom88
and cancellations among different components are not observed especially for
$M_{\eta^{\prime}}$. These new results may suggest that decay to the first
$2^{+}$ states is not that heavily suppressed as previously expected. In all
previous calculations, only the vector and axial-vector parts of the hadronic
current are considered. Another important component, namely the pseudo-scalar
part from the pion pole GT58  is not taken into account. As suggested in CDV17
, this pseudo-scalar piece is accounted as a LO contribution like the vector
and axial-vector parts. In most $0\nu\beta\beta$ calculations, the NLO weak-
magnetism contribution is also taken into account. In this work, we
incorporate all these parts into the calculation and study their effects on
the NME’s.

The calculation of NME relies on various nuclear many-body approaches, and we
limit our discussions to traditional ones with effective nuclear forces,
leaving out those of ab initio methods using a nuclear force starting from
bare nucleon-nucleon interaction. For the decay to the ground states with the
standard neutrino mass mechanism, several calculations have been done over the
decades. First we mention the time-consuming large scale shell model (LSSM)
calculations CNP99 ; Men17 , which takes advantage of the existence of the
shell gap, and separates the particles of the core and the nucleons of the
valence part.

Apart from LSSM calculations, which are applicable to limited cases due to the
large computation requirement, several other methods can be applied to more
occassions, e.g. the IBM-2 method BKI13 , the DFT methods of non-relativistic
RM10  and relativistic versions SYR14 , the project HFB method RCC13 , as well
as the QRPA method SPV99 ; SRF13 ; HS15 ; FFS18 , which takes intermediate
states into account. For recent reviews of these calculations, we refer to
YMN21 ; ABD22 .

Meanwhile, the NME calculations for the non-standard LR-symmetric model with
the inclusion of the non-helicity-suppressed $q$ terms are less frequently
considered for $0\nu\beta\beta(0^{+})$. Recent calculations have been done by
QRPA MBK89 ; SDS15 ; SSD17 , LSSM HN15 ; HN18 ; AH19 ; SIR20 ; IS21  as well
as by PHFB RCC19 . Especially for CNP96 ; HN15 ; HN18 ; AH19 , contributions
from SM effective field theory has been thoroughly analyzed besides the LR-
symmetric model. There are also other calculations starting from SM effective
field theories which have different expressions as traditional LR symmetric
modelsCDV17 ; CDV18 .

Compared to above calculations of decay to the ground states, the decay to the
excited $2^{+}$ state is rarely discussed, and the most recent calculations
are about two decades ago Tom99  as mentioned above. To investigate this
special decay mode, we adopt the QRPA method with realistic nuclear forces
FF21 . Our previous results suggest that the final NMEs are larger than
expected with only the vector and axial vector parts of the hadronic currents
included. In this work, we include more components up to NLO to make a more
thorough investigation.

This article is arranged as follows: At first we present the formalism we use.
It is followed by the results for the phase space factors and the nuclear
matrix elements. Then we discuss constraints on the L-R model parameters from
current results. Finally we give a conclusion as well as an outlook.

##  II Formalism

In the L-R symmetric models, such as
SU(2)${}_{L}\times$SU(2)${}_{R}\times$U(1)B-L, after two successive
spontaneous symmetry breakings, the left- and right- handed gauge bosons
acquire masses through the Higgs mechanism, and in general the left- and
right- handed gauge bosons are mixed MS81 ; DKT85 :

$\displaystyle\left(\begin{array}[]{c}W_{L}\\\
W_{R}\end{array}\right)=\left(\begin{array}[]{cc}\cos\xi&-\sin\xi\\\
\sin\xi&\cos\xi\end{array}\right)\left(\begin{array}[]{c}W_{1}\\\
W_{2}\end{array}\right)$ (7)

Here $\xi$ is the mixing angle and $W_{1}$, $W_{2}$ the mass eigenstates of W
bosons.

The neutrinos acquire masses through their Yukawa coupling with Higgs bosons
MS81 :

$\displaystyle\left(\begin{array}[]{c}\nu_{L}\\\
N_{R}\end{array}\right)=\left(\begin{array}[]{cc}U&U^{\prime}\\\
V^{\prime}&V\end{array}\right)\left(\begin{array}[]{c}\nu_{M}\\\
N_{M}\end{array}\right)$ (14)

Where $\nu_{L}^{T}=(\nu_{e},\nu_{\mu},\nu_{\tau})$ and
$N_{R}^{T}=(N_{e},N_{\mu},N_{\tau})$ are the three flavor left- and right-
handed neutrinos. $\nu_{M}$ and $N_{M}$ are their light and heavy mass
eigenstates. The see-saw mechanism can be naturally fulfilled in this model
MS81 .

Starting from the left- and right- gauge-fermion interactions, the effective
weak Hamiltonian can be written following the definition in DKT85 :

$\displaystyle H_{\rm
eff}=\sqrt{\frac{1}{2}}G_{F}\cos\theta_{C}(j_{L\mu}\tilde{J}^{\mu}_{L}+j_{R\mu}\tilde{J}^{\mu}_{R})+h.c.$
(15)

Where the lepton currents are:

$\displaystyle
j_{L(R)}^{\mu}(\vec{x})=\bar{\psi}_{e}(\vec{x})\gamma^{\mu}P_{L(R)}\psi_{\nu}(\vec{x})$
(16)

With $P_{L}=(1-\gamma_{5})/2$ and $P_{R}=(1+\gamma_{5})/2$ respectively.

In current model, the hadronic currents have the form:

$\displaystyle\tilde{J}_{L\mu}$ $\displaystyle\approx$ $\displaystyle
J_{L\mu}$ $\displaystyle\tilde{J}_{R\mu}$ $\displaystyle\approx$
$\displaystyle\eta J_{L\mu}+\lambda J_{R\mu}$ (17)

Here

$\displaystyle\eta\equiv-(g_{R}/g_{L})\tan\xi[1-(M_{1}/M_{2})^{2}]/$ (18)
$\displaystyle[1+\tan^{2}\xi(M_{1}/M_{2})^{2}]\hskip 28.45274pt$
$\displaystyle\lambda\equiv(g_{R}/g_{L})^{2}[(M_{1}/M_{2})^{2}+\tan^{2}\xi]/$
(19) $\displaystyle[1+\tan^{2}\xi(M_{1}/M_{2})^{2}]\hskip 28.45274pt$

Here $M_{1}$ and $M_{2}$ are the mass eigenvalues of $W_{1}$ and $W_{2}$ gauge
bosons respectively.

Within the non-relativistic impulse approximation, under the Breit frame, the
left- or right- handed hadronic currents have the form:

$\displaystyle J_{L\mu}(\vec{x})$ $\displaystyle=$
$\displaystyle(J_{0}(\vec{x}),\vec{J}_{L}(\vec{x}))$ $\displaystyle
J_{R\mu}(\vec{x})$ $\displaystyle=$
$\displaystyle(J_{0}(\vec{x}),\vec{J}_{R}(\vec{x}))$ (20)

Where

$\displaystyle J_{0}(\vec{x})$ $\displaystyle=$
$\displaystyle\sum_{n=1}^{A}g_{V}(q^{2})\delta(\vec{x}-\vec{r}_{n})$
$\displaystyle\vec{J}_{L}$ $\displaystyle=$
$\displaystyle\sum_{n=1}^{A}-[g_{A}(q^{2})\vec{\sigma}_{n}-g_{P}(q^{2})(\vec{\sigma}_{n}\cdot\vec{q})\vec{q}$
$\displaystyle+$ $\displaystyle
i\frac{g_{M}(q^{2})}{2m_{p}}(\vec{\sigma}_{n}\times\vec{q})]\delta(\vec{x}-\vec{r}_{n})$
$\displaystyle\vec{J}_{R}$ $\displaystyle=$
$\displaystyle\sum_{n=1}^{A}[g_{A}(q^{2})\vec{\sigma}_{n}-g_{P}(q^{2})(\vec{\sigma}_{n}\cdot\vec{q})\vec{q}$
(21) $\displaystyle-$ $\displaystyle
i\frac{g_{M}(q^{2})}{2m_{p}}(\vec{\sigma}_{n}\times\vec{q})]\delta(\vec{x}-\vec{r}_{n})$

According to angular momentum conservation, for $0\nu\beta\beta(2^{+})$, the
emitted electrons must be coupled to total angular momentum $J=2$, this
suggests that the dominant contribution comes from the combination of
decomposed partial waves $s_{1/2}$-$p_{3/2}$ DKT85 . Substituting the hadronic
currents into the S-matrix, the decay width can then be written as DKT85 ;
Tom88 :

$\displaystyle\Gamma=G_{1}|M_{\lambda}\langle\lambda\rangle-
M_{\eta}\langle\eta\rangle|^{2}+G_{2}|M^{\prime}_{\eta}\langle\eta\rangle|^{2}$
(22)

Here G’s are phase space factors and M’s are the NMEs. $\langle\lambda\rangle$
and $\langle\eta\rangle$ are new physics parameters connected to $\lambda$ and
$\eta$ defined above asSDS15 :
$\langle\eta\rangle=\eta|\sum_{j}U_{ej}V^{\prime*}_{ej}|$ and
$\langle\lambda\rangle=\lambda|\sum_{j}U_{ej}V^{\prime*}_{ej}(g^{\prime}_{V}/g_{V})|$,
here $g_{V}$ and $g^{\prime}_{V}$ are the vector coupling constant for left-
and right-handed currents respectively. U and V’ are neutrino mass mixing
matrix elements. And j sums over the light neutrino mass eigenstates.

The phase space factor(PSF) can be expressed as Tom88 :

$\displaystyle G_{i}$ $\displaystyle=$ $\displaystyle\frac{4\pi}{\ln
2~{}R_{\rm n}^{4}}\int
dE_{1}dE_{2}\frac{(G\cos\theta_{C})^{4}}{32\pi^{6}}f_{i}p_{1}p_{2}E_{1}E_{2}$
(23) $\displaystyle\times$
$\displaystyle\delta(E_{1}+E_{2}-2m_{e}-Q_{\beta\beta}(2^{+}))$

Here

$\displaystyle
f_{1}=3[|f^{-2-1}|^{2}+|f_{21}|^{2}+|f^{-1-2}|^{2}+|f_{12}|^{2}]$
$\displaystyle
f_{2}=3[|f^{-2}{}_{1}|^{2}+|f^{-1}{}_{2}|^{2}+|f^{-2}{}_{1}|^{2}+|f^{-1}{}_{2}|^{2}]$
(24)

For these $f$ functions such as $f_{12}$, we follow the convention in KI12 .
In our definition, the phase space factor are with the unit of $y^{-1}$, this
is obtained by dividing $R_{\rm n}^{2}$ to eq.(11) of Tom88 . Here $R_{\rm
n}=1.2A^{1/3}$fm is the conventionally defined nuclear radius. By deriving
this, the no Finite de Broglie wave length correction (no FBWC) approximation
is used DKT85 .

The expressions for the NME are much more complicated, and we follow the
conventions in Tom88 , divide the NME into seven parts:

$\displaystyle M_{\lambda}=\sum_{i=1}^{5}C_{\lambda i}M_{i}$ $\displaystyle
M_{\eta}=\sum_{i=1}^{5}C_{\eta i}M_{i}$ $\displaystyle
M^{\prime}_{\eta}=\sum_{i=6}^{7}C^{\prime}_{\eta i}M_{i}$

For the coefficients $C$’s, we follow the definition of Tom88  too and they
are tabulated in Table.1, where we have absorbed $g_{A}$ into the NME unlike
the conventional treatment where coupling constants $g$’s are included in
PSFs.

The NME can be further expressed as $M_{i}=\langle
2_{f}^{+}||\mathcal{M}_{i}||0_{i}^{+}\rangle$, and these operators
$\mathcal{M}_{i}$ can be expressed in a general form:

$\displaystyle\mathcal{M}_{i}=\frac{2R_{\rm
n}}{\pi}\int\frac{qdq}{q+E_{N}}h_{i}(q,r)\mathcal{O}_{i}$ (25)

Here $E_{N}=E_{mx}+M_{m}-(M_{i}+M_{f}+E_{2^{+}})/2$ is the intermediate state
excitation energy relative to the initial and final states. And $h$ is usually
called the neutrino potential and $\mathcal{O}$ is the angular transition
operator.

Where for different NME components, the detailed forms for neutrino potential
$h_{i}$ is as follows:

$\displaystyle h_{1}(q,r)$ $\displaystyle=$
$\displaystyle\frac{1}{g_{A}^{2}(0)}j_{1}(qr)[g_{A}^{2}(q^{2})+2\frac{g_{A}(q^{2})g_{P}(q^{2})q^{2}}{2m_{p}}-\frac{g_{P}^{2}(q^{2})q^{4}}{(2m_{p})^{2}}-2\frac{g^{2}_{M}(q^{2})q^{2}}{(2m_{p})^{2}})]$
$\displaystyle h_{2}(q,r)$ $\displaystyle=$
$\displaystyle\frac{1}{g_{A}^{2}(0)}j_{1}(qr)[g_{A}^{2}(q^{2})-\frac{g_{A}(q^{2})g_{P}(q^{2})q^{2}}{2m_{p}}+\frac{1}{5}\frac{g_{P}^{2}(q^{2})q^{4}}{(2m_{p})^{2}}+\frac{1}{5}\frac{g^{2}_{M}(q^{2})q^{2}}{(2m_{p})^{2}}]$
$\displaystyle h_{3}(q,r)$ $\displaystyle=$
$\displaystyle\frac{1}{g_{A}^{2}(0)}\\{j_{1}(qr)[g_{A}^{2}(q^{2})-\frac{g_{P}(q^{2})g_{A}(q^{2})q^{2}}{2m_{p}}+\frac{1}{5}\frac{g_{P}^{2}(q^{2})q^{4}}{(2m_{p})^{2}}+\frac{1}{5}\frac{g^{2}_{M}(q^{2})q^{2}}{(2m_{p})^{2}}]-j_{3}(qr)[\frac{3}{35}\frac{g_{P}^{2}(q^{2})q^{4}}{(2m_{p})^{2}}+\frac{3}{35}\frac{g^{2}_{M}(q^{2})q^{2}}{(2m_{p})^{2}}]\\}$
$\displaystyle h^{\prime}_{3}(q,r)$ $\displaystyle=$
$\displaystyle\frac{1}{g_{A}^{2}(0)}j_{3}(qr)[\sqrt{\frac{3}{5}}\frac{g_{P}^{2}(q^{2})q^{4}}{(2m_{p})^{2}}+\sqrt{\frac{3}{5}}\frac{g^{2}_{M}(q^{2})q^{2}}{(2m_{p})^{2}}]$
$\displaystyle h_{4}(q,r)$ $\displaystyle=$
$\displaystyle\frac{g_{V}^{2}(q^{2})}{g_{V}^{2}(0)}j_{1}(qr),\quad
h_{5}(q,r)=\frac{g_{V}(q^{2})g_{A}(q^{2})}{g_{V}(0)g_{A}(0)}j_{1}(qr),\quad
h_{6}(q,r)=h_{7}(q,r)=\frac{g_{V}(q^{2})g_{A}(q^{2})}{g_{A}^{2}(0)}\frac{r_{+}j_{1}(qr)}{r}$
(26)

Here $q=|\vec{q}|$ is the exchange momentum carried by the neutrino
propagator. The $g_{\alpha}(q^{2})$’s are the form factors. It can be written
in an empiric dipole form in general:
$g_{V}(q^{2})=g_{V}(0)/(1+q^{2}/\Lambda_{V})^{2}$ and
$g_{A}(q^{2})=g_{A}(0)/(1+q^{2}/\Lambda_{A})^{2}$. Here $g_{V}(0)=1$ and
$g_{A}(0)=1.27$, and we take $\Lambda_{V}=0.85$GeV and $\Lambda_{A}=1.1$GeV
for the energy cutoff. In nuclear environment, $g_{A}(0)$ is usually quenched
with a not definitely known origin. Therefore in the current work we adopt two
values for this coupling constant: the bare one and a quenched one with a
quenching factor $g_{A}=0.75g_{A0}$. And
$g_{M}(q^{2})=(1+\kappa_{1})g_{V}(q^{2})$SPV99  with
$\kappa_{1}=\mu_{p}-\mu_{n}$. Also
$g_{P}(q^{2})=2m_{p}g_{A}(q^{2})/(q^{2}+m_{\pi}^{2})$ is given by the PCAC
hypothesis.

And the angular operators $\mathcal{O}$s have the forms:

$\displaystyle\mathcal{O}_{1}$ $\displaystyle=$
$\displaystyle\vec{\sigma}_{1}\cdot\vec{\sigma}_{2}[\hat{r}\otimes\hat{r}]^{(2)}$
$\displaystyle\mathcal{O}_{2}$ $\displaystyle=$
$\displaystyle[\vec{\sigma}_{1}\otimes\vec{\sigma}_{2}]^{(2)}$
$\displaystyle\mathcal{O}_{3}$ $\displaystyle=$
$\displaystyle[[\vec{\sigma}_{1}\otimes\vec{\sigma}_{2}]^{(2)}\otimes[\hat{r}\otimes\hat{r}]^{(2)}]^{(2)}$
$\displaystyle\mathcal{O}^{\prime}_{3}$ $\displaystyle=$
$\displaystyle[[\vec{\sigma}_{1}\otimes\vec{\sigma}_{2}]^{(2)}\otimes[[\hat{r}\otimes\hat{r}]^{(2)}\otimes[\hat{r}\otimes\hat{r}]^{(2)}]^{(4)}]^{(2)}$
$\displaystyle\mathcal{O}_{4}$ $\displaystyle=$
$\displaystyle[\hat{r}\otimes\hat{r}]^{(2)}$ $\displaystyle\mathcal{O}_{5}$
$\displaystyle=$
$\displaystyle[(\vec{\sigma}_{1}+\vec{\sigma}_{2})\otimes[\hat{r}\otimes\hat{r}]^{(2)}]^{(2)}$
(27)

for the non-primed NMEs $M_{\lambda}$ and $M_{\eta}$.

And

$\displaystyle\mathcal{O}_{6}$ $\displaystyle=$
$\displaystyle[(\vec{\sigma}_{1}-\vec{\sigma}_{2})\otimes[\hat{r}\otimes\hat{r}_{+}]^{(1)}]^{(2)}$
$\displaystyle\mathcal{O}_{7}$ $\displaystyle=$
$\displaystyle[(\vec{\sigma}_{1}-\vec{\sigma}_{2})\otimes[\hat{r}\otimes\hat{r}_{+}]^{(2)}]^{(2)}$
(28)

For the primed NME $M_{\eta}^{\prime}$.

Compared to the expression in Tom88 , $M_{3}$ has one extra term induced by
the hadronic current, we denote it by $M_{3}^{\prime}$. For $M_{4}\sim M_{7}$,
no corrections from the induced hadronic current are presented.

For the sake of comparison with the decay to the ground state cases, we find
that we have also AA, AP, PP and MM components for the space-space current-
current interactions ($M_{1}$, $M_{2}$ and $M_{3}$) coming from different
components of induced hadronic current. Except for $M_{1}$, one finds that PP
and MM components are suppressed by a factor smaller than $1/5$ compared to AA
and AP components.

The NME are calculated in our case with the nuclear many-body approach, the
so-called pn-QRPA as well as charge conserving QRPA methods both with
realistic nuclear forces FF21 . The detailed expression can be found in FF21 :

$\displaystyle M_{i}$ $\displaystyle=$
$\displaystyle\sum_{pnp^{\prime}n^{\prime}}^{J^{\pi}m}\langle
2^{+}_{f}||\widetilde{[c_{p}^{\dagger}\tilde{c}_{n}]}_{J^{\prime}}||J^{\pi}m_{f}\rangle\langle
J^{\pi}m_{f}||J^{\pi}m_{f}\rangle$ (29) $\displaystyle\times$
$\displaystyle\langle
J^{\pi}m_{i}||[c_{p^{\prime}}^{\dagger}\tilde{c}_{n^{\prime}}]_{J}||0^{+}_{i}\rangle$

The expressions for one body densities as well as the overlap of the initial
and final intermediate states can also be found in FF20 , where the final
$2^{+}$ states are obtained by charge conserving QRPA SC93 ; FF20 .

Table 1: The decomposition coefficients of NMEs for $M_{\lambda}$, $M_{\eta}$ and $M_{\eta}^{\prime}$. Here $g_{A}=g_{A}(0)$. i | 1 | 2 | 3 | 4 | 5 | 6 | 7
---|---|---|---|---|---|---|---
$C_{\Lambda i}$ | $\frac{1}{3}g_{A}^{2}$ | $-\frac{2}{3}g_{A}^{2}$ | $\sqrt{\frac{7}{3}}g_{A}^{2}$ | $1$ | $-\sqrt{\frac{3}{2}}g_{A}$ |  |
$C_{\eta i}$ | $\frac{1}{3}g_{A}^{2}$ | $-\frac{2}{3}g_{A}^{2}$ | $\sqrt{\frac{7}{3}}g_{A}^{2}$ | $-1$ | $0$ |  |
$C_{\eta i}^{\prime}$ |  |  |  |  |  | $\sqrt{\frac{1}{2}}g_{A}$ | $-\sqrt{\frac{3}{2}}g_{A}$

##  III Results and discussions

###  III.1 Phase space factors

Table 2: PSFs for various nuclei from current work and that of DKT85  as well as current experimental limits. For results from DKT85 , we have converted to our convention. We also present the constraints on new physics parameter on two special cases. | Q(MeV) |  $G_{1}$ ($10^{-15}y^{-1}$) |  $G_{2}$ ($10^{-15}y^{-1}$) |  $t_{1/2}^{\rm limit}$(yr) | $|\lambda|\gg|\eta|,m_{\beta\beta}$ | $|\eta|\gg|\lambda|,m_{\beta\beta}$
---|---|---|---|---|---|---
|  | this work | DKT85  | this work | DKT85  |  | $2^{+}$ |  $0^{+}$SDS15  | $2^{+}$ |  $0^{+}$SDS15
76Ge | 1.480 | 6.86 | 7.37 | 4.77 | 5.12 |  $>$2.1$\times$1024MAJ20  |  $<$2.13$\times 10^{-5}$ |  $<$5.07$\times 10^{-7}$ |  $<$4.56$\times 10^{-6}$ |  $<$2.81$\times 10^{-9}$
82Se | 2.219 | 50.12 | 55.28 | 40.11 | 44.30 |  $>$1.0$\times$ 1022Bar17  |  $<$8.35$\times 10^{-5}$ |  |  $<$5.63$\times 10^{-5}$ |
96Zr | 2.572 | 140.4 |  | 117.2 |  |  $>$9.1$\times$1020Bar17  |  $<$2.76$\times 10^{-3}$ |  |  $<$2.38$\times 10^{-4}$ |
100Mo | 2.495 | 134.2 | 149.8 | 111.3 | 124.3 |  $>$1.6$\times$1023Bar17  |  $<$2.52$\times 10^{-4}$ |  |  $<$1.73$\times 10^{-5}$ |
116Cd | 1.520 | 18.18 |  | 12.98 |  |  $>$6.2$\times$1022Bar17  |  $<$1.75$\times 10^{-4}$ |  |  $<$2.32$\times 10^{-5}$ |
128Te | 0.423 | 0.225 | 0.267 | 0.0825 | 0.0985 |  |  |  |  |
130Te | 1.991 | 76.41 | 91.37 | 60.04 | 71.84 |  $>$1.4$\times$1023Bar17  |  $<$2.35$\times 10^{-5}$ |  |  $<$5.11$\times 10^{-6}$ |
136Xe | 1.639 | 35.61 | 44.99 | 26.33 | 33.42 |  $>$2.6$\times$1025Bar17  |  $<$2.43$\times 10^{-6}$ |  $<$4.35$\times 10^{-7}$ |  $<$4.17$\times 10^{-6}$ |  $<$2.12$\times 10^{-9}$

For the calculation of PSFs, we use the numerical package Radial SFW95  for
the electron wave-functions and we follow the convention in KI12 . We use a
uniform charge distribution for the calculations of the nuclear static charge
potential and we choose the charge radius to be the same as the nuclear
radius. We neglect the screening effect from the orbital electrons since it
gives minor corrections to PSF KI12  in the case of the decay to the ground
state.

The results of PSFs of various nuclei are presented in Table 2, we also list
these nuclei’s Q values. For most nuclei except 128Te, PSFs span a range for
about two orders of magnitude. Three nuclei have Q values larger than $2$MeV,
and of which 96Zr has the largest PSF for the decay into the excited state.
100Mo has almost the same values for the PSF as 96Zr, since their Q values are
close. While 82Se has a larger Q value, its PSFs are somehow smaller than that
of 130Te mostly due to its smaller atomic number $Z$.

We also present results from an earlier calculation DKT85 , where a Taylor
expansion of the electron wave function is used. For 76Ge, we have also
results from Tom88 , $G_{1}=7.34\times 10^{-15}y^{-1}$ and $G_{2}=5.10\times
10^{-15}y^{-1}$, which are basically the same to results in DKT85  with a
deviation less than 1% since they use similar treatment for electron wave
functions. Our current numerical results generally agree with their results
quantitatively. We find the deviations for ours and theirs are about
10%$\sim$30%, the current numerical results are smaller than their
predictions. There is a strong trend of an increase for the deviations as $Z$
increases. This is reasonable, with the growth of Z, the nuclear radius R also
increases. Since the conventional Taylor expansion method uses $\alpha Z$ and
$WR$ as variables, the errors will grow with the increase of these variables.
One finds the largest phase space factor for 100Mo from their calculations. If
this nucleus has also the largest NME, then it can be one of the most
promising candidate for a future experimental search. To explore such
possibility, we need high precision nuclear many-body calculations. We will
proceed into this direction in the next part.

###  III.2 Nuclear matrix elements

Table 3: NMEs for decay to $2^{+}_{1}$ with $g_{A}=g_{A0}=1.27$, here $a$ and $b$ refer to two different cases with AV-18 and CD-Bonn short range correlations adopted.  |  |  76Ge | 82Se |  96Zr |  100Mo |  116Cd |  128Te |  130Te |  136Xe
---|---|---|---|---|---|---|---|---|---
|  | a | b | a | b | a | b | a | b | a | b | a | b | a | b | a | b
$M_{1}$ | AA | 0.641 | 0.640 | 0.790 | 0.789 | 0.027 | 0.027 | 0.205 | 0.205 | 0.198 | 0.199 | 0.700 | 0.700 | 0.643 | 0.643 | 0.334 | 0.333
AP | 0.481 | 0.480 | 0.649 | 0.648 | 0.046 | 0.046 | 0.141 | 0.141 | -0.058 | -0.057 | 0.425 | 0.425 | 0.449 | 0.448 | 0.563 | 0.562
PP | -0.087 | -0.087 | -0.149 | -0.149 | -0.022 | -0.022 | -0.022 | -0.022 | 0.067 | 0.067 | -0.044 | -0.043 | -0.070 | -0.070 | -0.206 | -0.206
MM | -0.033 | -0.033 | -0.098 | -0.098 | -0.027 | -0.028 | -0.007 | -0.007 | 0.071 | 0.072 | 0.018 | 0.018 | -0.016 | -0.016 | -0.182 | -0.183
tot. | 1.001 | 1.000 | 1.193 | 1.191 | 0.024 | 0.023 | 0.317 | 0.317 | 0.279 | 0.280 | 1.099 | 1.100 | 1.006 | 1.006 | 0.510 | 0.506
$M_{2}$ | AA | -0.249 | -0.252 | -0.005 | -0.006 | 0.238 | 0.239 | -0.052 | -0.054 | -0.992 | -0.996 | -0.268 | -0.269 | -0.175 | -0.175 | 0.405 | 0.408
AP | 0.207 | 0.210 | -0.028 | -0.027 | -0.229 | -0.230 | 0.265 | 0.267 | 0.564 | 0.568 | 0.142 | 0.143 | 0.059 | 0.060 | -0.416 | -0.419
PP | -0.034 | -0.034 | -0.002 | -0.003 | 0.029 | 0.030 | -0.043 | -0.043 | -0.066 | -0.067 | -0.022 | -0.022 | -0.012 | -0.012 | 0.044 | 0.045
MM | -0.021 | -0.023 | -0.005 | -0.006 | 0.014 | 0.014 | -0.022 | -0.023 | -0.035 | -0.037 | -0.014 | -0.015 | -0.008 | -0.009 | 0.022 | 0.023
tot. | -0.097 | -0.098 | -0.041 | -0.042 | 0.052 | 0.052 | 0.148 | 0.147 | -0.530 | -0.532 | -0.162 | -0.163 | -0.136 | -0.136 | 0.056 | 0.058
$M_{3}$ | AA | -0.049 | -0.049 | 0.118 | 0.118 | -0.008 | -0.008 | 0.158 | 0.158 | 0.015 | 0.015 | -0.098 | -0.098 | -0.009 | -0.009 | 0.407 | 0.407
AP | -0.031 | -0.031 | -0.059 | -0.059 | -0.005 | -0.005 | -0.034 | -0.034 | -0.053 | -0.053 | -0.039 | -0.039 | -0.044 | -0.044 | -0.056 | -0.056
PP | 0.004 | 0.004 | -0.000 | -0.000 | -0.001 | -0.001 | -0.001 | -0.001 | 0.001 | 0.001 | 0.006 | 0.006 | 0.003 | 0.003 | -0.009 | -0.009
MM | 0.001 | 0.001 | -0.003 | -0.003 | -0.001 | -0.001 | -0.002 | -0.002 | -0.002 | -0.002 | -0.001 | -0.001 | -0.001 | -0.001 | -0.005 | -0.005
tot. | -0.074 | -0.074 | 0.056 | 0.056 | -0.015 | -0.015 | 0.121 | 0.121 | -0.040 | -0.040 | -0.132 | -0.132 | -0.052 | -0.052 | 0.336 | 0.336
$M^{\prime}_{3}$ | PP | -0.003 | -0.003 | -0.023 | -0.023 | 0.004 | 0.004 | -0.021 | -0.021 | -0.012 | -0.012 | -0.019 | -0.019 | -0.019 | -0.019 | -0.019 | -0.019
MM | -0.001 | -0.001 | -0.012 | -0.012 | 0.001 | 0.001 | -0.010 | -0.010 | -0.007 | -0.007 | -0.010 | -0.010 | -0.010 | -0.010 | -0.009 | -0.009
tot. | -0.004 | -0.004 | -0.035 | -0.035 | 0.005 | 0.005 | -0.031 | -0.031 | -0.019 | -0.019 | -0.029 | -0.029 | -0.030 | -0.030 | -0.028 | -0.028
$M_{4}$ | -0.147 | -0.147 | -0.123 | -0.123 | -0.072 | -0.071 | 0.121 | 0.120 | 0.142 | 0.142 | -0.188 | -0.188 | -0.139 | -0.139 | 0.119 | 0.119
$M_{5}$ | -0.061 | -0.061 | -0.216 | -0.216 | -0.257 | -0.257 | 0.045 | 0.045 | -0.124 | -0.124 | -0.130 | -0.130 | -0.161 | -0.161 | -0.280 | -0.280
$M_{\lambda}$ | 0.398 | 0.399 | 0.952 | 0.952 | 0.262 | 0.261 | 0.284 | 0.285 | 0.909 | 0.911 | 0.384 | 0.385 | 0.597 | 0.597 | 1.527 | 1.523
$M_{\eta}$ | 0.597 | 0.598 | 0.862 | 0.862 | 0.005 | 0.004 | 0.112 | 0.114 | 0.431 | 0.435 | 0.558 | 0.559 | 0.626 | 0.626 | 0.854 | 0.849
$M_{6}$ | 0.631 | 0.640 | 0.314 | 0.315 | 0.128 | 0.128 | 0.481 | 0.484 | -0.609 | -0.615 | 0.803 | 0.816 | 0.669 | 0.679 | 0.054 | 0.049
$M_{7}$ | -1.241 | -1.251 | -0.336 | -0.339 | 0.419 | 0.422 | -0.080 | -0.083 | 0.948 | 0.951 | -1.752 | -1.762 | -1.346 | -1.353 | 0.187 | 0.190
$M^{\prime}_{\eta}$ | 2.496 | 2.520 | 0.804 | 0.810 | -0.537 | -0.542 | 0.557 | 0.564 | -2.021 | -2.032 | 3.446 | 3.474 | 2.694 | 2.715 | -0.242 | -0.251
Table 4: The same table as Table.3 but with $g_{A}=0.75g_{A0}$ |  |  76Ge | 82Se |  96Zr |  100Mo |  116Cd |  128Te |  130Te |  136Xe
---|---|---|---|---|---|---|---|---|---
|  | a | b | a | b | a | b | a | b | a | b | a | b | a | b | a | b
$M_{1}$ | AA | 0.634 | 0.634 | 0.786 | 0.785 | 0.027 | 0.027 | 0.193 | 0.193 | 0.200 | 0.200 | 0.699 | 0.698 | 0.644 | 0.644 | 0.333 | 0.332
AP | 0.476 | 0.476 | 0.648 | 0.647 | 0.045 | 0.045 | 0.137 | 0.137 | -0.058 | -0.057 | 0.424 | 0.424 | 0.448 | 0.447 | 0.563 | 0.561
PP | -0.086 | -0.086 | -0.149 | -0.149 | -0.022 | -0.022 | -0.022 | -0.022 | 0.067 | 0.067 | -0.043 | -0.043 | -0.070 | -0.070 | -0.206 | -0.205
MM | -0.033 | -0.033 | -0.098 | -0.098 | -0.027 | -0.028 | -0.007 | -0.006 | 0.071 | 0.072 | 0.018 | 0.019 | -0.016 | -0.015 | -0.181 | -0.183
tot. | 0.991 | 0.991 | 1.187 | 1.185 | 0.022 | 0.022 | 0.302 | 0.302 | 0.280 | 0.281 | 1.097 | 1.098 | 1.007 | 1.006 | 0.508 | 0.505
$M_{2}$ | AA | -0.237 | -0.240 | -0.005 | -0.006 | 0.240 | 0.241 | -0.119 | -0.121 | -0.962 | -0.966 | -0.272 | -0.273 | -0.172 | -0.173 | 0.404 | 0.407
AP | 0.204 | 0.207 | -0.027 | -0.025 | -0.226 | -0.227 | 0.277 | 0.279 | 0.553 | 0.558 | 0.144 | 0.145 | 0.059 | 0.060 | -0.414 | -0.418
PP | -0.033 | -0.034 | -0.003 | -0.003 | 0.029 | 0.029 | -0.043 | -0.044 | -0.066 | -0.066 | -0.022 | -0.023 | -0.012 | -0.012 | 0.044 | 0.045
MM | -0.021 | -0.022 | -0.005 | -0.006 | 0.013 | 0.014 | -0.022 | -0.023 | -0.035 | -0.037 | -0.014 | -0.015 | -0.008 | -0.009 | 0.022 | 0.023
tot. | -0.087 | -0.089 | -0.039 | -0.040 | 0.056 | 0.057 | 0.092 | 0.091 | -0.509 | -0.511 | -0.165 | -0.165 | -0.134 | -0.134 | 0.056 | 0.058
$M_{3}$ | AA | -0.049 | -0.049 | 0.118 | 0.118 | -0.007 | -0.007 | 0.149 | 0.149 | 0.018 | 0.018 | -0.099 | -0.100 | -0.010 | -0.010 | 0.407 | 0.406
AP | -0.031 | -0.031 | -0.059 | -0.058 | -0.005 | -0.005 | -0.034 | -0.034 | -0.054 | -0.054 | -0.039 | -0.038 | -0.044 | -0.044 | -0.056 | -0.056
PP | 0.004 | 0.004 | -0.000 | -0.000 | -0.001 | -0.001 | -0.001 | -0.001 | 0.001 | 0.001 | 0.006 | 0.006 | 0.003 | 0.003 | -0.009 | -0.009
MM | 0.001 | 0.001 | -0.003 | -0.003 | -0.001 | -0.001 | -0.002 | -0.002 | -0.002 | -0.002 | -0.001 | -0.001 | -0.001 | -0.001 | -0.005 | -0.005
tot. | -0.074 | -0.074 | 0.056 | 0.056 | -0.014 | -0.014 | 0.112 | 0.112 | -0.038 | -0.038 | -0.133 | -0.133 | -0.053 | -0.053 | 0.336 | 0.336
$M^{\prime}_{3}$ | PP | -0.003 | -0.003 | -0.023 | -0.023 | 0.004 | 0.004 | -0.021 | -0.021 | -0.013 | -0.012 | -0.019 | -0.019 | -0.019 | -0.019 | -0.019 | -0.019
MM | -0.001 | -0.001 | -0.012 | -0.012 | 0.001 | 0.001 | -0.010 | -0.010 | -0.007 | -0.007 | -0.010 | -0.010 | -0.010 | -0.010 | -0.009 | -0.009
tot. | -0.004 | -0.004 | -0.035 | -0.035 | 0.005 | 0.005 | -0.031 | -0.031 | -0.020 | -0.020 | -0.029 | -0.029 | -0.030 | -0.030 | -0.028 | -0.028
$M_{4}$ | -0.146 | -0.146 | -0.123 | -0.122 | -0.071 | -0.071 | 0.120 | 0.120 | 0.143 | 0.143 | -0.187 | -0.188 | -0.140 | -0.140 | 0.119 | 0.119
$M_{5}$ | -0.061 | -0.062 | -0.214 | -0.214 | -0.258 | -0.257 | 0.048 | 0.048 | -0.124 | -0.124 | -0.128 | -0.128 | -0.161 | -0.161 | -0.280 | -0.279
$M_{\lambda}$ | 0.161 | 0.162 | 0.489 | 0.489 | 0.194 | 0.193 | 0.170 | 0.171 | 0.568 | 0.569 | 0.133 | 0.133 | 0.280 | 0.280 | 0.943 | 0.941
$M_{\eta}$ | 0.391 | 0.392 | 0.535 | 0.535 | 0.032 | 0.031 | 0.027 | 0.028 | 0.170 | 0.172 | 0.395 | 0.395 | 0.411 | 0.411 | 0.427 | 0.425
$M_{6}$ | 0.632 | 0.641 | 0.312 | 0.312 | 0.128 | 0.127 | 0.479 | 0.482 | -0.612 | -0.618 | 0.799 | 0.813 | 0.671 | 0.681 | 0.052 | 0.047
$M_{7}$ | -1.242 | -1.252 | -0.333 | -0.336 | 0.422 | 0.424 | -0.094 | -0.097 | 0.948 | 0.952 | -1.748 | -1.758 | -1.347 | -1.355 | 0.187 | 0.190
$M^{\prime}_{\eta}$ | 1.874 | 1.893 | 0.598 | 0.602 | -0.406 | -0.410 | 0.432 | 0.437 | -1.518 | -1.526 | 2.577 | 2.598 | 2.024 | 2.040 | -0.183 | -0.190

While PSF can be fairly well determined with decent accuracy within about
several percents, the NME’s have larger uncertainties, greater than a factor
of two YMN21 . In current work we adopt QRPA for the many-body calculations.
Compared to the LSSM approach, QRPA has the advantage of small computation
requirement with the price of less accuracy.

For the QRPA method with realistic forces several renormalization parameters
are needed to reproduce experimental data. We fit the pairing parameters by
the experimental pairing gaps. And we renormalize the residual forces by
$g_{pp}^{T=1}$ to eliminate $M_{F}^{2\nu}$ and $g_{pp}^{T=0}$ to reproduce the
measured $2\nu\beta\beta$ NME.

All our results for the NME’s are presented in Table 3 and Table 4. The two
tables correspond to two possible values of the axial vector coupling constant
$g_{A}$ commonly used in the literature. In these tables we also give the
results with two different src’s (short range correlations) marked as $a$ and
$b$ for AV-18 and CD-Bonn SFR09 , respectively and we find these different
src’s barely change the results. As we are aware, different choices of src’s
will change the NME’s for $0\nu\beta\beta(0^{+})$ by about $5-10\%$.
Therefore, $0\nu\beta\beta(2^{+})$ is less sensitive to the choice of of the
src’s. This suggests that the two nucleons involved in this decay are more far
away from each other, as we shall see in a quantitive analysis later.

As shown in eq. (26), except for $M_{1}$, the PP and MM terms are suppressed
by a factor of about $1/5\sim 1/10$. Therefore, since all the terms related to
spherical Bessel function $j_{3}$ in $M_{3}$ are contained in PP and MM terms,
they are heavily suppressed. Even not suppressed by this small factor, the MM
terms are generally smaller for $M_{1}$. In the contrary, for most cases, the
AP term plays an important role and they are supposed to be the LO
contribution CDV17 .

We first discuss the bare $g_{A}$ case. For $M_{1}$, besides the LO AA
contribution, the major correction of about 2/3 comes from AP term as
mentioned above. For most nuclei, this term leads to an enhancement for the
final results but for 116Cd, it gives an 30% reduction. For 136Xe, it gives
the largest contribution even larger than AA, this makes it the dominant term
in $M_{1}$ for this nucleus. Meanwhile for most nuclei, the PP term gives the
reductions from 10%$\sim$20% relative to AA term. Two exceptions are the 116Cd
and 136Xe cases: for the former, this term gives an enhancement about 30%
relative to the AA term which nearly cancels the reduction from the AP term as
we discussed above; for the latter, PP gives a 2/3 reduction relative to AA,
this term and the AP term together gives nearly the same contribution as AA.
The MM terms are supposed to be NLO contributions CDV17 . For $M_{1}$, it
gives negligible contributions for the nuclei 100Mo, 128Te and 130Te, while
corrections for the other nuclei ranges from 5% to 50%. For the extreme case
of 96Zr, the MM term has an equal size as AA due to the smallness of later.

For 96Zr, $M_{1}$ is heavily suppressed, it is one order of magnitude smaller
than other nuclei, and a careful check suggests that this smallness comes from
the cancellation from different intermediate states, this differs from the
case of 76Ge, where different intermediate states add up together (see fig.1
of FF21  ). For all other nuclei, the total $M_{1}$ is generally within one
order of magnitude, e.g. for 76Ge, 82Se, 128Te and 130Te, its value is around
$1$.

The case for $M_{2}$ is quite different from $M_{1}$. For the AA component,
different nuclei differ by several orders of magnitude, especially for 82Se,
it is heavily suppressed by the cancellations from different intermediate
states. This cancellation also reduces the NME for 100Mo. For other nuclei it
has a value of about $0.2\sim 0.3$, but for the two cases with semi magic
nuclei involved in the transition process, the results are somehow enhanced.
Further investigation is needed for the possible relation between the magicity
and this enhancement.

The correction from the pseudo-scalar current is much more pronounced for
$M_{2}$. For $M_{1}$ the AP term is about one order of magnitude larger than
the PP term. Unlike for $M_{1}$, AP terms here mostly appear for
cancellations, this causes an overall smallness of $M_{2}$ for most nuclei
compared to $M_{1}$. The semi magic nuclei - as already mentioned above - have
generally large AA terms, however for 136Xe, the AP term exactly cancels the
AA term. This leads to a suppressed $M_{2}$; meanwhile for 116Cd, the
cancellation from AP is smaller compared to AA. Thus we have the largest
$M_{2}$ for this nucleus. Like for 136Xe, a nearly exact cancellation between
AA and AP terms also happens for 76Ge and 96Zr. For 100Mo, we find that the
contribution from AP is about 5 times larger than that from AA. While for Te
isotopes, AP reduces the results by about 50%. The PP term behaves like NLO
contributions by the suppression from a small coefficient, their contributions
are generally within a magnitude of 10%. So does the NLO MM term which
contributes with about 10% for most cases.

Thus for most nuclei, $M_{2}$ is around 0.1 due to the cancellations between
AA and AP terms, with the exception of 82Se where all the components are
small. 116Cd is the one with the largest $M_{2}$ about $0.5$. As a result, our
calculations suggest that $M_{2}$ is generally one order of magnitude smaller
than $M_{1}$.

For most nuclei, $M_{3}$ is generally smaller than $M_{1}$ but similar to
$M_{2}$. However, unlike the smallness caused by the cancellation between AA
and AP terms for $M_{2}$, the magnitude of each component of $M_{3}$ is
generally much smaller than the counterparts in $M_{1}$ or $M_{2}$. For 100Mo
and 136Xe, the situation is a bit different: AA terms for them are close to
that of $M_{1}$, but on the other hand, there is no large enhancement from AP
term for $M_{3}$, therefore their $M_{3}$ are generally smaller. And we find
for $M_{3}$, PP and MM terms are negligible. In general, except 136Xe,
$M_{3}$’s are around or smaller than $0.1$.

We also find that the induced $M^{\prime}_{3}$ is much smaller and they barely
give visible contributions. But for 82Se and 130Te, $M^{\prime}_{3}$ is
comparable to $M_{3}$, and gives cancellation for the former nucleus and
enhancement for the latter.

In general, for the space-space components of hadronic currents, the inclusion
of pseudo-scalar and weak-magnetism terms especially the AP term will change
the NME drastically and these changes are usually not negligible. Whether
these contributions are positive or negative is nucleus dependent.

$M_{4}$ comes from the time-time component of hadronic currents and the
induced current will not contribute to this component. As in FF21 , they are
generally with the magnitude of $0.1$ and the difference for different nuclei
are generally within a factor of two. The largest value is found for 128Te and
the smallest from 96Zr, differed nearly by a factor of 3.

$M_{5}$ is the space-time component of hadronic currents. Their values are
generally close to $M_{4}$ with an magnitude around 1, but differ in details,
such as the phases. The relative magnitude of these two NME’s is nucleus
dependent.

$M_{\lambda}$ and $M_{\eta}$ are induced from the coupling of $\vec{q}$ and
the relative coordinate of the two decaying nucleons $\vec{r}$. The difference
of $M_{\lambda}$ and $M_{\eta}$ comes from the space-time and time-time
components of the nuclear currents, it can be defined as $\delta
M=2M_{4}-\sqrt{3/2}g_{A}M_{5}$. This suggests that for nuclei with $M_{4}$ and
$M_{5}$ close to each other, such as 82Se or 130Te, this difference is small.
On the other hand, this difference is also related to $g_{A}$. For all nuclei
in our calculation, $M_{\lambda}$ and $M_{\eta}$ are with the same phase, this
means they will cancel or add up with each other depending on the relative
sign of $\langle\eta\rangle$ and $\langle\lambda\rangle$. For 76Ge, 82Se,
128Te and 130Te, $M_{\lambda}$ is smaller than $M_{\eta}$, while for other
nuclei, $M_{\lambda}$ is larger. In general, except for 136Xe, $M_{\lambda}$
and $M_{\eta}$ is about one order of magnitude smaller than $M_{0\nu}(0^{+})$
(To compare our results with those in various literatures, we need to divide
current results by $g_{A0}^{2}$). 136Xe does have $M_{\lambda}$ larger than 1
in our convention, but it is still less than half of the values for
$M_{0\nu}(0^{+})$.

$M^{\prime}_{\eta}$ are actually induced by the coupling of the $\vec{q}$ term
in the neutrino propagator and the COM (Center Of Mass) coordinate of the two
decaying nucleons $\vec{r}_{+}$. The differences of $M_{6}$ and $M_{7}$
originate from the coupled angular momenta of the COM and relative coordinate,
thus are closely related to the orbitals of the decaying nucleons. In our
calculation, we find for most nuclei, these two terms give coherent
contribution but for 96Zr and 136Xe, they cancel with each other leading to
suppressed $M^{\prime}_{\eta}$’s. For most nuclei, we have also larger
$M_{7}$, but for 96Zr, $M_{7}$ is suppressed and is much smaller than $M_{6}$.
Also in most cases, $M^{\prime}_{\eta}$ is much larger than $M_{\lambda}$ and
$M_{\eta}$. Instead of one order of magnitude smaller than corresponding
$M^{0\nu}$($0^{+}$), they are generally smaller within a factor of one half.
This differs with previous studies Tom88  where these NMEs are supposed to be
several order of magnitude smaller than $M^{0\nu}(0^{+})$ and hence can be
safely neglected. We will study the consequence of this dominance for the role
of determining new physics parameters in the following section.

As presented in Table 3, different src’s produce actually negligibly smaller
deviations like for the light neutrino mass mechanism of
$0\nu\beta\beta(0^{+})$. On the other hand, a comparison of Table 3 and Table
4 suggests that quenching of $g_{A}$ affects the final results drastically. In
QRPA calculations the quenching affects the final results two-fold:
First we fit our model parameter $g_{pp}^{T=0}$ by the $2\nu\beta\beta$ NME
which is related to $g_{A}$. Therefore, quenching of $g_{A}$ will affect the
choice of this parameter and subsequently the individual NMEs.
Second the quenching changes the coefficients $C$’s in Table. 1 which are than
multiplied with the individual partial NME’s to get the final total NME’s.
For the current calculation, we find that the change of the parameter
$g_{pp}^{T=0}$ induced by quenching will not largely change the individual
NMEs, their changes are within 10%. Hence the change of the final results can
be directly connected to the quenched value of $g_{A}$ in the coefficients
$C$.

For $M_{\lambda}$, only the coefficient of $M_{4}$ are independent of $g_{A}$,
and all other coefficients have a dependence linear or squared on $g_{A}$.
Therefore, the quenching generally reduce $M_{\lambda}$ by a magnitude from
30% to 50% for different nuclei. The largest reduction in percentage is found
for 128Te, where the NME for the quenched case is only one third of that of
the bare $g_{A}$ case. The reduction for $M_{\eta}$ is similar to that of
$M_{\lambda}$, we now have no $g_{A}$ linear dependent terms any more. The
general reduction is from 30% to 70%, the largest reduction is found for 100Mo
with more than 70%. While for 96Zr an enhancement is observed, but for this
nucleus, all the individual NME’s are small and they cancel each other
significantly, the presence of quenching largely reduce the cancellation of
$M_{1}\sim M_{3}$ to $M_{4}$ which then leads to the enhancement. All these
behaviors for various nuclei stem from the interplay between $M_{4}$ and other
individual NME’s which are reduced by $g_{A}$. Drastic reductions are also
found for $M^{\prime}_{\eta}$ with a nearly fixed factor $q=g_{A}/g_{A0}$
since the two terms are both proportional to $g_{A}$ and the $g_{pp}^{T=0}$
dependence is weak for $M_{6}$ and $M_{7}$. Therefore, as for the decay to the
ground states, the role of $g_{A}$ is important. We need to understand the
origin of the quenching and the precise value of $g_{A}$ for the
$\beta\beta$-decay, we also need to determine wether the quenching is operator
dependent.

###  III.3 Constraints on L-R symmetric model

Figure 1: (Color online) Constraints on $\langle\lambda\rangle$ and
$\langle\eta\rangle$ from $0\nu\beta\beta(2^{+})$-decay. The region outside
the curves are excluded by the current limits obtained from the measurements
and our calculations.

Unlike for the decay to ground state, where the neutrino mass mechanism is
supposed to be dominant, the neutrino mass mechanism for
$0\nu\beta\beta(2^{+})$ can only be triggered by nuclear recoil Tom99  and
hence is suppressed by a factor of $Q/2M$ compared to $\langle\lambda\rangle$
and $\langle\eta\rangle$ terms. In Tom99 , an estimation suggests that their
NMEs are about several orders of magnitude smaller than that of LR mechanism.
Since this NME for the mass mechanism in the decay to excited state is
extremely hindered, we do not expect to to be able to observe this decay mode
with the dominance of the mass mechanism.

Therefore measurements of $0\nu\beta\beta(2^{+})$ can be used as perfect
constraints on new physics parameters such as $\langle\lambda\rangle$ or
$\langle\eta\rangle$. In this section, with the calculated NME and PSF as well
as the experimental limit presented in Table.2, we perform a simple analysis
on several extreme cases. Currently the most stringent constrain for half-
lives is that of 136Xe, which pushes the limit to $10^{25}$y while others are
two or three orders of magnitude shorter. For 96Zr, the lower limit is
$10^{20}$ years, which is much shorter than all others. While SDS15  provides
a way of probing different mechanism from comparisons of different nuclei,
current method would help distinguish the underlying new physics within one
isotope. This is done by comparing the decay rates to the ground and excited
state as we shall show.

####  III.3.1 $|\lambda|\gg|\eta|$

In the $\lambda$ dominant case (See Table.2), eq.(22) becomes
$\Gamma=G_{1}M_{\lambda}^{2}|\langle\lambda\rangle|^{2}$ and the most
stringent constraints for $\lambda$ are obtained from 136Xe. With our
calculated NME’s, we obtain the constraint $|\lambda|<2.43\times 10^{-5}$.
Other nuclei yield generally the same magnitude except 96Zr and 100Mo which
sets limits one order of magnitude larger. Since all these nuclei except 136Xe
have basically similar $M_{\lambda}$, these larger limits generally are due to
shorter half-life limits. In Table 2, we present constraints obtained in SDS15
from decay to ground states for two nuclei from QRPA calculations without
considering the hadronic induced currents. With a later calculation SSD17 ,
one finds that the NME for $0\nu\beta\beta(0^{+})$ is slightly enhanced, but
this will not change the general magnitude of the constraints. Therefore, we
can make a direct comparison. Our results suggest that for the $\lambda$
dominant case, the constraints for 76Ge differs by two orders of magnitude
while for 136Xe this difference is less than one order of magnitude while
comparing our results with those from SDS15 . The half-life lower limits used
for $0\nu\beta\beta(0^{+})$ for both nuclei are around $10^{25}$ yr, while for
$0\nu\beta\beta(2^{+})$, 76Ge has a half-life lower-limit one order of
magnitude shorter. However, 136Xe has the similar half-life lower limits to
the $2^{+}$ and to the ground state. Combining these results, we find that for
certain nuclei (136Xe here), for the case of $\lambda$ dominance, the half-
lives for the two modes are within a difference of two orders of magnitude,
this could perhaps be used in the future for the identification of the decay
mechanism if both modes are observed.

####  III.3.2 $|\lambda|\ll|\eta|$

For this case, eq.(22) becomes $\Gamma=(G_{1}M_{\eta}^{2}+G_{2}M^{\prime
2}_{\eta})|\langle\eta\rangle|^{2}$ and the most stringent constraints come
also from 136Xe. However, 76Ge and 130Te yield also constraints close to
136Xe. They all require that $|\eta|$ is smaller than about $10^{-6}$.
However, these constraints are far looser than that of
$0\nu\beta\beta(0^{+})$, where a limit around $|\eta|<10^{-9}$ is obtained
assuming the dominance of the $\eta$ mechanism. If this case is true in
nature, it will not be possible to observe the $2^{+}$ decay in these nuclei,
since we have an suppression for about 6 orders of magnitude for
$0\nu\beta\beta(2^{+})$ compared to $0\nu\beta\beta(0^{+})$ for 136Xe and
4$\sim$5 orders of magnitude for other nuclei. If future measurements do push
the half-life limit of $0\nu\beta\beta(2^{+})$ to 4$\sim$5 orders of magnitude
longer, we can then rule out both $\lambda$ and $\eta$ mechanisms.

####  III.3.3 $|\lambda|\sim|\eta|$

More general cases are presented in fig.1. For these cases, since
$\langle\eta\rangle$ and $\langle\lambda\rangle$ share the same phase
$e^{i\psi}$, eq.(22) can be written explicitly as
$\Gamma=G_{1}M_{\lambda}^{2}|\langle\lambda\rangle|^{2}+(G_{1}M_{\eta}^{2}+G_{2}M_{\eta}^{2})|\langle\eta\rangle|^{2}-2G_{1}M_{\lambda}M_{\eta}\langle\lambda\rangle\langle\eta\rangle$.
In the current convention, $\lambda$ is always positive and the sign for
$\eta$ is not definite. Since different nuclei give constraints of $\lambda$
and $\eta$ on different orders of magnitude, we use a logarithm scale in the
graph. The interference term of $\eta$ and $\lambda$ rotates the long axis of
the ellipse to the first quadrant. From the figure, we can clearly see that,
like in the $\lambda$ or $\eta$ dominant cases, 136Xe sets the most stringent
constraints. While 76Ge and 130Te set both less tight constraints for
$\lambda$ and $\eta$. And it is obvious that the future slight improvement of
the measurement will surely set more stringent constraints for these two
nuclei. This is highly probable for 130Te, of which current half-life lower
limit is about two orders of magnitude smaller than that for 136Xe. And better
half-life lower limits for other nuclei will also improve our analysis. A
combined analysis with decay to ground states could give us more hints on the
relative magnitude of $\langle\eta\rangle$ and $\langle\lambda\rangle$.

In general, the measurement of decay to the $2^{+}$ excited states offers a
way for the discrimination of underlying mechanism. For example, for 136Xe,
within our calculations, a half-time lower limit for about 2 orders of
magnitude longer than the observed half-life of $0\nu\beta\beta(0^{+})$ will
rule out the possibility of the existence of the $\lambda$ mechanism
dominance. However, these conclusions rely heavily on the NME calculations,
and need further verifications from other many-body calculations.

##  IV Conclusions and Outlooks

In this work, we systematically investigate the decay rates of
$0\nu\beta\beta(2^{+})$ under the LR symmetric model, where both the $\lambda$
and $\eta$ mechanisms are involved. We incorporate the contributions from
induced currents and find that the pseudo-scalar current is important for NME
calculations. Our results suggest that, the NME is nucleus dependent and may
differ by up to one order of magnitude, also the NME for the two different
mechanisms in the same nucleus can be differed by more than one order of
magnitude. We also improve the phase space factor calculations by using
numerical electron wave functions. These results lead to different constraints
on different nuclei both from different NME’s and different current lower
limits of decay half-lives. With the comparison to $0\nu\beta\beta(0^{+})$, we
can constrain the underlying mechanisms within individual nucleus. These
comparisons suggest that for the special $\lambda$ dominant case, 136Xe has
the potential of comparable decay half-lives for decays to both ground state
and $2^{+}$ excited state. Future experiments could shift the lower limit up
and set more stringent constraints on the new physics parameters for other
nucleus, and even have the potential of discovering the possible existence of
right-handed gauge bosons. Meanwhile, to draw a more solid conclusion, more
nuclear many-body calculations are needed to make a better prediction of the
NMEs which are keys for these investigations.

## acknowledgement

This work is supported by the ”Light of West China” and ”From Zero to One”
programs from CAS and Guangdong Major Project of Basic as well as Applied
Basic Research No. 2020B0301030008 from Guangdong Province. The numerical
calculations in this paper have been carried out on the supercomputing system
in the Southern Nuclear Science Computing Center. This work is also supported
by the National Key Research and Development Program of China
(2021YFA1601300).

## References

  * (1) R. N. Mohapatra and G. Senjanovic, Phys. Rev. D 23, 165 (1981)
  * (2) E. Lisi, A. M. Rotunno and F. Simkovic, Phys. Rev. D 92, 093004 (2015)
  * (3) L. Gráf, M. Lindner and O. Scholer, Phys. Rev. D 106, 035022 (2022)
  * (4) M. Duerr, M. Lindner and K. Zuber, Phys. Rev. D 84, 093004 (2011)
  * (5) M. Doi, T. Kotani and E. Takasugi, Prog. Theor. Phys. Suppl. 83, 1 (1985)
  * (6) D. Stefanik, R. Dvornicky, F. Simkovic and P. Vogel, Phys. Rev. C 92, 055502 (2015)
  * (7) T. Tomoda, Nucl. Phys. A 484, 635-646 (1988)
  * (8) T. Tomoda, Phys. Lett. B 474, 245-250 (2000)
  * (9) D. L. Fang and A. Faessler, Phys. Rev. C 103, 045501 (2021)
  * (10) M. L. Goldberger and S. B. Treiman, Phys. Rev. 111, 354-361 (1958)
  * (11) V. Cirigliano, W. Dekens, J. de Vries, M. L. Graesser and E. Mereghetti, JHEP 12, 082 (2017)
  * (12) E. Caurier, F. Nowacki, A. Poves and J. Retamosa, Nucl. Phys. A 654, 973c (1999); J. Menendez, A. Poves, E. Caurier and F. Nowacki, Nucl. Phys. A 818, 139 (2009); E. Caurier, F. Nowacki and A. Poves, Phys. Lett. B 711, 62 (2012).
  * (13) J. Menéndez, J. Phys. G 45, 014003 (2018)
  * (14) J. Barea, J. Kotila and F. Iachello, Phys. Rev. C 87, 014315 (2013); ibid 91, 034304 (2015).
  * (15) T. R. Rodriguez and G. Martinez-Pinedo, Phys. Rev. Lett. 105, 252503 (2010)
  * (16) L. S. Song, J. M. Yao, P. Ring and J. Meng, Phys. Rev. C 90, 054309 (2014); J. M. Yao, L. S. Song, K. Hagino, P. Ring and J. Meng, Phys. Rev. C 91, 024316 (2015); L. S. Song, J. M. Yao, P. Ring and J. Meng, Phys. Rev. C 95, 024305 (2017).
  * (17) P. K. Rath, R. Chandra, K. Chaturvedi, P. Lohani, P. K. Raina and J. G. Hirsch, Phys. Rev. C 88, 064322 (2013)
  * (18) F. Simkovic, G. Pantis, J. D. Vergados and A. Faessler, Phys. Rev. C 60, 055502 (1999).
  * (19) F. Šimkovic, V. Rodin, A. Faessler and P. Vogel, Phys. Rev. C 87, 045501 (2013).
  * (20) J. Hyvärinen and J. Suhonen, Phys. Rev. C 91, 024613 (2015)
  * (21) D. L. Fang, A. Faessler and F. Simkovic, Phys. Rev. C 97, no.4, 045503 (2018)
  * (22) J. M. Yao, J. Meng, Y. F. Niu and P. Ring, Prog. Part. Nucl. Phys. 126, 103965 (2022)
  * (23) M. Agostini, G. Benato, J. A. Detwiler, J. Menéndez and F. Vissani, [arXiv:2202.01787 [hep-ex]].
  * (24) K. Muto, E. Bender and H. V. Klapdor, Z. Phys. A 334, 187-194 (1989)
  * (25) F. Šimkovic, D. Štefánik and R. Dvornický, Fron. in Phys. 5, 57 (2017)
  * (26) E. Caurier, F. Nowacki, A. Poves and J. Retamosa, Phys. Rev. Lett. 77, 1954-1957 (1996)
  * (27) M. Horoi and A. Neacsu, Phys. Rev. D 93, no.11, 113014 (2016)
  * (28) M. Horoi and A. Neacsu, Phys. Rev. C 98, 035502 (2018)
  * (29) F. Ahmed and M. Horoi, Phys. Rev. C 101, 035504 (2020)
  * (30) S. Sarkar, Y. Iwata and P. K. Raina, Phys. Rev. C 102, 034317 (2020)
  * (31) Y. Iwata and S. Sarkar, Front. Astron. Space Sci. 8, 727880 (2021)
  * (32) P. K. Rath, R. Chandra, K. Chaturvedi and P. K. Raina, Front. in Phys. 7, 64 (2019)
  * (33) V. Cirigliano, W. Dekens, J. de Vries, M. L. Graesser and E. Mereghetti, JHEP 12, 097 (2018)
  * (34) J. Kotila and F. Iachello, Phys. Rev. C 85, 034316 (2012)
  * (35) J. Suhonen and O. Civitarese, Phys. Lett. B 308, 212-215 (1993)
  * (36) D. L. Fang and A. Faessler, Chin. Phys. C 44, no.8, 084104 (2020)
  * (37) F. Salvat, J. Fernadez-Varea, W. Williamson, Jr. Comput.Phys.Commun. 90, 151 (1995)
  * (38) I. J. Arnquist et al. [MAJORANA], Phys. Rev. C 103, no.1, 015501 (2021)
  * (39) A. S. Barabash, AIP Conf. Proc. 1894, no.1, 020002 (2017)
  * (40) F. Simkovic, A. Faessler, H. Muther, V. Rodin and M. Stauf, Phys. Rev. C 79, 055501 (2009)


# Efficient and Interpretable Traffic Destination Prediction using Explainable
Boosting Machines

Yasin Yousif<EMAIL_ADDRESS>
Department of Informatics
Clausthal University of Technology
Clausthal-Zellerfeld, Germany Jörg P. Müller<EMAIL_ADDRESS>
Department of Informatics
Clausthal University of Technology
Clausthal-Zellerfeld, Germany

###### Abstract

Developing accurate models for traffic trajectory predictions is crucial for
achieving fully autonomous driving. Various deep neural network models have
been employed to address this challenge, but their black-box nature hinders
transparency and debugging capabilities in a deployed system. Glass-box models
offer a solution by providing full interpretability through methods like
Generalized Additive Models (GAM). In this study, we evaluate an efficient
additive model called Explainable Boosting Machines (EBM) for traffic
prediction on three popular mixed traffic datasets: Stanford drone dataset
(SDD), Intersection Drone Dataset (InD), and Argoverse. Our results show that
the EBM models perform competitively in predicting pedestrian destinations
within SDD and InD while providing modest predictions for vehicle-dominant
Argoverse dataset. Additionally, our transparent trained models allow us to
analyse feature importance and interactions, as well as provide qualitative
examples of predictions explanation. The full training code will be made
public upon publication.

Keywords: Glass-box model, Traffic trajectory prediction, Partial dependence
graphs

##  1 Introduction

Traffic trajectory prediction is the problem of predicting the future path of
a traffic entity for a specific time horizon based on its current state, and
it has been widely studied in public traffic benchmarks. Many previous works
Bertugli et al. (2021); Cheng et al. (2021); Chiara et al. (2022) have been
evaluated on these benchmarks under standard conditions to assess their
prediction accuracy. By accurately predicting traffic trajectories, we can
anticipate and avoid accidents or unsafe situations on the road, which is
particularly useful for self-driving cars navigating safely.

Early methods to address this problem relied on rule-based approaches that
applied general rules to determine future movements, such as social force
models Helbing and Molnar (1995) or cellular automata models Maerivoet and De
Moor (2005). However, these methods lack the ability to extract data-based
patterns from traffic data. Later, deep learning-based models emerged Alahi et
al. (2016); Gupta et al. (2018), which automatically extracted rules from the
datasets and outperformed other methodologies on well-known traffic benchmarks
such as Argo Chang et al. (2019) and InD Bock et al. (2020).

In recent years, some deep learning models Sadeghian et al. (2019); Yue et al.
(2022); Bansal et al. (2018) have focused on enhancing the explainability of
their models in a post-hoc manner. This is because a black-box model that
provide accurate predictions using hidden processes for feature extraction and
prediction inference can decrease trust in the model and make debugging any
errors difficult.

This motivated works like Sadeghian et al. (2019); Bansal et al. (2018) to
extract attention maps of the input images or works like Yue et al. (2022) to
use a more interpretable input features. However, these remain partial
solutions, as they still considered an approximation of the true features
importance Lou et al. (2012).

An alternative approach is to use glass-box models such as linear regression,
decision trees, or GAM Hastie (2017) that provide transparent explanations at
the cost of lower predictive ability. Among these methods, GAM has been found
to offer the best predictive results in some applications (e.g., Caruana et
al. (2015); Lou et al. (2012)). Recent variants of GAM, such as Neural
Additive Models (NAM) or boosted trees (EBM Nori et al. (2019)), have also
shown competitive prediction results while maintaining transparency. In this
work, we use an EBM model to predict traffic trajectories and evaluate our
results on three well-known datasets in the field. We also extract feature
importance values and partial dependence graphs that provide a clear
description of how the model arrives at its predictions. These graphs show
each input feature’s contribution to the output variable, offering an exact
representation of the model’s decision-making process

In GAM, multiple sub-models are used, with each accessing one or two input
features. These sub-models were previously shown to perform competitively when
implemented using boosted trees, known as EBM. The implementation of EBM in
this work is based on the Interpretml Python library Nori et al. (2019).
However, there are limitations to GAM’s current formulation such as the
inability to process images as input since they need to be processed
collectively for useful information extraction. Additionally, having only a
single output for the GAM model is another restriction that can be addressed
by using full models for each output variable to predict. This approach
involves having two models for each $(x,y)$ coordinate of the predicted
trajectory in this work.

Due to the difficulty in predicting a single deterministic possible output
trajectory (known as mode), prediction is done for multiple modes with one EBM
for each. The clustering and separation of these modes is an independent
problem, addressed here by clustering the target variables

This work has two main contributions: 1) The evaluation of EBM as a glass-box
model for the first time by using it as a traffic destination prediction model
on three popular traffic datasets, namely SDD, InD, and Argoverse. This
provides insight into how transparent and interpretable models like EBM can
perform in real-world datasets. 2) We calculate partial dependence graphs and
feature importance for all input features with respect to the output variable.
These results offer insights into the interaction of features and their
significance for traffic prediction, as well as show some selected examples of
local explanations. The following sections outline GAM formula, the
methodology containing the preprocessing steps, experimental evaluation, and
finally, the conclusions.

##  2 GAM Formula

To define GAM and GA2M mathimtical formulas, let
$\mathcal{D}={(\mathbf{x}_{i},y_{i})}^{N}_{1}$ be the training dataset, and
$\mathbf{x}_{i}=(x_{i1},x_{i2},..,x_{ip})$ the input vector for the target
variable of $y_{i}$, and $x_{j}$ is the feature $j$ value. Therefore GA2M,
Caruana et al. (2015), will try to fit the dataset as the following
expression:

$g(E[y])=\beta_{0}+\sum_{j}f_{j}(x_{j})+\sum_{k\neq j}f_{kj}(x_{k},x_{j})$ (1)

Where $f_{j}$ is the model mapping function for feature $j$, $f_{kj}$ is for
the interaction between features $(j,k)$, $g$ is the link function and
$\beta_{0}$ is the bias parameter. This model will be completely explainable,
because one can draw the exact relationship between the input feature $x_{j}$
and the target $y$ as 2D plot or as heatmap for the interaction between two
features as in the last term.

With EBM, in addition to fitting each feature alone, a fast search for the
most useful feature interactions in the dataset is done as in Lou et al.
(2013), then the top $k$ features interactions are used in the model, where
$k$ is selected using cross validation.

List of the hyperparameters of EBM in this work are provided in appendix B.

##  3 Methodology

In our approach, we generate multiple outputs representing different modes and
assign probabilities for each mode. We detail the representation of multi-
modality and probability calculation in the first two subsections. Following
that is three subsections to focus on input and output shape specifications
for SDD, InD, and Argo datasets, respectively.

###  3.1 Multiple Modes Output

Most traffic prediction benchmarks require predicting multiple possible future
trajectories with a probability distribution over them, which is known as
multi-modal output Phan-Minh et al. (2020). This condition ease the prediction
task; however, in our case, EBM has a unimodal output, so we can train
separate models for each mode to obtain multi-modal output.

However, defining these modes is another challenge. Many approaches were used
in the literature, like using a set of rule-based modes Phan-Minh et al.
(2020), or with an implied clustering in multiple branches output Yousif and
Müller (2022).In this work, we manually cluster the target points using either
K-means clustering algorithm for SDD and InD datasets (see figure 1) or
uniform splitting over x-axis followed by y-axis to equally distribute samples
across modes for Argoverse dataset (see figure 2).

###  3.2 Outputs Probabilities

Each of the previous outputs should be assigned a probability normalized over
all the outputs. This will add flexibility to the prediction because it can
give an arbitrary number of outputs and then only the most probable $K$ output
will be taken. In this work, probability estimation of modes is done by
finding the sum of the loglikelihood of each of the mode-level models
$\mathcal{L}_{1}$ with the corresponding loglikelihood of its predictions
under the distribution of the main uni-modal model $\mathcal{L}_{2}$, as in
equation 2.

$\mathcal{L}_{mode_{x}}=\lambda_{1}\mathcal{L}_{1}(\pi_{mode_{x}}(x),\sigma_{mode_{x}})+\lambda_{2}\mathcal{L}_{2}(\pi_{mode_{x}}(x),\sigma_{all_{x}})$
(2)

After finding that for both of $x$ and $y$ axises, we take a weighted sum of
them for each mode in equation 3. Finally, the top 20 probable output will be
taken as final.

$\mathcal{L}_{mode}=\lambda_{3}\mathcal{L}_{mode_{x}}+\lambda_{4}\mathcal{L}_{mode_{y}}$
(3)

The weights $\lambda_{3}$ and $\lambda_{4}$ are determined based on relative
standard deviation between axes, while $\lambda_{1}$ and $\lambda_{2}$ can be
experimentally set

###  3.3 SDD

In the SDD experiments, our models predict the final point of each trajectory,
which the most challenging part to predict. This allows for comparison with
other methods in the SDD benchmark using the Final Displacement Error (FDE)
metric. In this work, we use pixels as the standard unit of measurement,
similar to previous works such as Phan-Minh et al. (2020); Chiara et al.
(2022); Deo and Trivedi (2020).

After loading data files from 60 scenes, they will be split according to a
common splitting standards as in Deo and Trivedi (2020), where only
pedestrians’ trajectories are considered. The training and validation sets are
combined and passed to the fitting function.

The input features of our model include the last 7 positions relative to the
8th position, with a timestep of 400 ms between each two positions. These
points are rotated by the latest heading angle, and we also consider the sum
of lost, occluded, and generated flags for the eight input points. Finally,
the width and height of the object at the 8th point are included as well.

The target variables consist of the $(x,y)$ coordinates of the object after
4.8 seconds, rotated and transferred in a manner similar to the input
trajectory case. This requires two EBM models for each coordinate.

Figure 1: Left: full training set of target variable. Right: target variable
set is split into multiple clusters.

Figure 1 on the left shows the target variable for both training and
validation sets of pedestrians combined. The data is concentrated towards the
centre, but there are some points that are far away due to noise in SDD as
mentioned earlier (Becker et al. (2018)). To address this issue, we removed
any point outside a specific range for both $x$ and $y$, which represents only
0.005% of the total data.

In previous works like Mangalam et al. (2021); Guo et al. (2022); Deo and
Trivedi (2020), evaluations typically involved using 20 modes as outputs.
However, since we started by clustering the data into 36 groups using k-means
algorithm (as shown in figure 1 on the right), and then trained two separate
models for each $x$ and $y$ within these clusters, only the most probable 20
pairs of predictions will be selected. Furthermore, we also trained global
uni-modal models to help identify the best models later in the process.

###  3.4 InD

For each scene, the data from its CSV file is read, and then we formulate the
input similarly to SDD, as in previous studies such as Bertugli et al. (2021);
Chiara et al. (2022). This includes the last 8 positions of the past 3.2
seconds, which consist of their position, velocity, acceleration, heading,
width, and height. The target variable is the traffic entity’s position after
4.8 seconds; however, only pedestrian predictions are considered in this case.
Similar to SDD, we split the target variable data into 50 modes.

###  3.5 Argoverse

####  3.5.1 Preprocessing Steps

For Argoverse Chang et al. (2019), we aim to represent all available input
data, including the road network and other nearby traffic entities. However,
before feeding this information into our EBM model, some preprocessing is
required to make it easier for the model while also maintaining
interpretability of the features.

We then split the dataset based on the type of ego traffic entity present in
each scene. This results in different sets of models being trained for various
types of entities (cars, bicycles, buses, pedestrians, and motorcyclists),
with approximately 88% of the data representing cars while the rest is
distributed among these other categories.

####  3.5.2 Road Network Data

The road network is initially represented as a binary image, where the
drivable area has high value. Next, different sub-areas representing various
modes are created based on the target variable (shown in figure 2). These mask
parts of the drivable area and find the geometric centre of each mode
rectangle. However, a tuple of zeros will be returned if there is no drivable
area available.

Figure 2: How the road network is represented (Left). Each mode rectangle
(Right) return the geometric center of the drivable area underneath it. Red
and green colors represent importance for x and y respectively

Using these same sub-area modes from figure 2 is an appropriate choice because
they represent the most relevant part of the road for its corresponding
prediction. Finally, the point of collision (POC) based on Post encroachment
time (PET) (Nasernejad et al. (2023)) is included as input features to capture
traffic movement in the area. A detailed explanation of this calculation can
be found in appendix A.

####  3.5.3 Defining Modes

Similar to SDD and InD, we rotate all points and move them to the last input
directions and positions for of their respective trajectory. For the challenge
conditions, 5 seconds of trajectory history are required as input, while a
prediction is made for 6 seconds into the future. To achieve this, we define
24 modes on the target variable (shown in figure 2 for vehicles) by splitting
uniformly along the x-axis and manually dividing along the y-axis to create
large central modes. Probabilities are calculated similarly to SDD and InD,
with the top 6 mode predictions being selected as final.

##  4 Experimental Evaluation

After training all models, we evaluate their performance using final
displacement errors displayed in table 1. The additive models show competitive
results on SDD and are comparable to state-of-the-art (SoTA) on InD without
road map input. However, they do not perform as well on Argoverse. One factor
for that is that cars predication is a harder problem than pedestrians’
predication and rely heavily on nearby elements and roadmap.

Table 1: Minimum final displacement errors for 20 modes on SDD and InD, and for 6 modes on Argoverse for EBM (ours) and a subset of SoTA previous methods Dataset | EBM | SEPT | ST-GAT | AC-VRNN | S-GAN | Y-Net | GOAL-SAR | P2T | PEC | TDOR
---|---|---|---|---|---|---|---|---|---|---
SDD (pixels) | 14.70 | - | - | - | 41.44 | 11.85 | 11.83 | 14.08 | 15.88 | 10.46
InD | 0.54 | - | 1 | 0.80 | 0.99 | 0.56 | 0.54 | - | - | -
Argoverse* | 3.88 | 1.09 | - | - | - | - | - | - | - | -
*Cars only predication. Updated leaderboard here
SEPT: Lan et al. (2023), ST-GAT: Huang et al. (2019),AC-VRNN: Bertugli et al.
(2021)
S-GAN: Gupta et al. (2018),Y-Net: Mangalam et al. (2021),GOAL-SAR: Chiara et
al. (2022)
P2T: Deo and Trivedi (2020),PEC: Mangalam et al. (2020),TDOR: Guo et al.
(2022)

Figures 3, 4, 5, 6, 7 and 8 show global feature importance bars for each
dataset using unimodal models. Additionally, partial dependence graphs of the
six most important features on ($x$ and $y$) are shown. It is noted in general
that the last steps are the most influencing input in all of the datasets.
Another interesting remark is that the acceleration is very crucial in InD
dataset models. Examples of local explanations can be found in appendix C.

Figure 3: SDD: Global Feature Average Importance for X and Y

(a) X-axis

(b) Y-axis

Figure 4: SDD: Partial Dependence Graph on X and Y axis for the best 6
features Figure 5: InD: Global Feature Average Importance for X and Y

(a) X-axis

(b) Y-axis

Figure 6: InD: Partial Dependence Graph on X and Y axes for the best 6
features Figure 7: Argo: Global Feature Average Importance for X

(a) X-axis

(b) Y-axis

Figure 8: Argo: Partial Dependence Graph on X and Y axis for the best 6
features

##  5 Conclusion

In this comprehensive analysis of generalized additive models for traffic
prediction, we have seen their notable performance when it comes to predicting
pedestrian movements. In fact, they outperform state-of-the-art deep learning
models in terms of accuracy. Furthermore, as a glass-box model, they offer
full transparency and interpretability. However, one potential drawback is
their limited ability to predict multiple variables per model, which requires
the training of numerous individual models for each point along the
trajectory. Future research will focus on developing more effective image
processing techniques for the road network, allowing us to fully leverage the
environmental information available in the data.

## Reproducibility Statement and Supplementary Material

The full training code is split into three notebooks for the three datasets,
with the instruction of how to run them here:

https://anonymous.4open.science/r/GAM4Traffic-5785/Readme.md

The full training (and testing) can be reproduced fairly fast on modest
hardware. In our case, the training was done on (Intel(R) Core(TM) i7-10850H
CPU)

## References

  * Alahi et al. (2016) A. Alahi, K. Goel, V. Ramanathan, A. Robicquet, L. Fei-Fei, and S. Savarese.  Social lstm: Human trajectory prediction in crowded spaces.  In _Proceedings of the IEEE conference on computer vision and pattern recognition_ , pages 961–971, 2016.
  * Bansal et al. (2018) M. Bansal, A. Krizhevsky, and A. Ogale.  Chauffeurnet: Learning to drive by imitating the best and synthesizing the worst.  _arXiv preprint arXiv:1812.03079_ , 2018.
  * Becker et al. (2018) S. Becker, R. Hug, W. Hübner, and M. Arens.  An evaluation of trajectory prediction approaches and notes on the trajnet benchmark.  _arXiv preprint arXiv:1805.07663_ , 2018.
  * Bertugli et al. (2021) A. Bertugli, S. Calderara, P. Coscia, L. Ballan, and R. Cucchiara.  Ac-vrnn: Attentive conditional-vrnn for multi-future trajectory prediction.  _Computer Vision and Image Understanding_ , 210:103245, 2021\.
  * Bock et al. (2020) J. Bock, R. Krajewski, T. Moers, S. Runde, L. Vater, and L. Eckstein.  The ind dataset: A drone dataset of naturalistic road user trajectories at german intersections.  In _2020 IEEE Intelligent Vehicles Symposium (IV)_ , pages 1929–1934. IEEE, 2020.
  * Caruana et al. (2015) R. Caruana, Y. Lou, J. Gehrke, P. Koch, M. Sturm, and N. Elhadad.  Intelligible models for healthcare: Predicting pneumonia risk and hospital 30-day readmission.  In _Proceedings of the 21th ACM SIGKDD international conference on knowledge discovery and data mining_ , pages 1721–1730, 2015.
  * Chang et al. (2019) M.-F. Chang, J. Lambert, P. Sangkloy, J. Singh, S. Bak, A. Hartnett, D. Wang, P. Carr, S. Lucey, D. Ramanan, et al.  Argoverse: 3d tracking and forecasting with rich maps.  In _Proceedings of the IEEE/CVF conference on computer vision and pattern recognition_ , pages 8748–8757, 2019.
  * Cheng et al. (2021) H. Cheng, W. Liao, M. Y. Yang, B. Rosenhahn, and M. Sester.  Amenet: Attentive maps encoder network for trajectory prediction.  _ISPRS Journal of Photogrammetry and Remote Sensing_ , 172:253–266, 2021.
  * Chiara et al. (2022) L. F. Chiara, P. Coscia, S. Das, S. Calderara, R. Cucchiara, and L. Ballan.  Goal-driven self-attentive recurrent networks for trajectory prediction.  In _Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition_ , pages 2518–2527, 2022.
  * Deo and Trivedi (2020) N. Deo and M. M. Trivedi.  Trajectory forecasts in unknown environments conditioned on grid-based plans.  _arXiv preprint arXiv:2001.00735_ , 2020.
  * Guo et al. (2022) K. Guo, W. Liu, and J. Pan.  End-to-end trajectory distribution prediction based on occupancy grid maps.  In _Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition_ , pages 2242–2251, 2022.
  * Gupta et al. (2018) A. Gupta, J. Johnson, L. Fei-Fei, S. Savarese, and A. Alahi.  Social gan: Socially acceptable trajectories with generative adversarial networks.  In _Proceedings of the IEEE conference on computer vision and pattern recognition_ , pages 2255–2264, 2018.
  * Hastie (2017) T. J. Hastie.  Generalized additive models.  In _Statistical models in S_ , pages 249–307. Routledge, 2017.
  * Helbing and Molnar (1995) D. Helbing and P. Molnar.  Social force model for pedestrian dynamics.  _Physical review E_ , 51(5):4282, 1995.
  * Huang et al. (2019) Y. Huang, H. Bi, Z. Li, T. Mao, and Z. Wang.  Stgat: Modeling spatial-temporal interactions for human trajectory prediction.  In _Proceedings of the IEEE/CVF international conference on computer vision_ , pages 6272–6281, 2019.
  * Lan et al. (2023) Z. Lan, Y. Jiang, Y. Mu, C. Chen, S. E. Li, H. Zhao, and K. Li.  Sept: Towards efficient scene representation learning for motion prediction.  _arXiv preprint arXiv:2309.15289_ , 2023.
  * Lou et al. (2012) Y. Lou, R. Caruana, and J. Gehrke.  Intelligible models for classification and regression.  In _Proceedings of the 18th ACM SIGKDD international conference on Knowledge discovery and data mining_ , pages 150–158, 2012.
  * Lou et al. (2013) Y. Lou, R. Caruana, J. Gehrke, and G. Hooker.  Accurate intelligible models with pairwise interactions.  In _Proceedings of the 19th ACM SIGKDD international conference on Knowledge discovery and data mining_ , pages 623–631, 2013.
  * Maerivoet and De Moor (2005) S. Maerivoet and B. De Moor.  Cellular automata models of road traffic.  _Physics reports_ , 419(1):1–64, 2005.
  * Mangalam et al. (2020) K. Mangalam, H. Girase, S. Agarwal, K.-H. Lee, E. Adeli, J. Malik, and A. Gaidon.  It is not the journey but the destination: Endpoint conditioned trajectory prediction.  In _Computer Vision–ECCV 2020: 16th European Conference, Glasgow, UK, August 23–28, 2020, Proceedings, Part II 16_ , pages 759–776. Springer, 2020.
  * Mangalam et al. (2021) K. Mangalam, Y. An, H. Girase, and J. Malik.  From goals, waypoints & paths to long term human trajectory forecasting.  In _Proceedings of the IEEE/CVF International Conference on Computer Vision_ , pages 15233–15242, 2021.
  * Nasernejad et al. (2023) P. Nasernejad, T. Sayed, and R. Alsaleh.  Multiagent modeling of pedestrian-vehicle conflicts using adversarial inverse reinforcement learning.  _Transportmetrica A: transport science_ , 19(3):2061081, 2023.
  * Nori et al. (2019) H. Nori, S. Jenkins, P. Koch, and R. Caruana.  Interpretml: A unified framework for machine learning interpretability.  _arXiv preprint arXiv:1909.09223_ , 2019.
  * Phan-Minh et al. (2020) T. Phan-Minh, E. C. Grigore, F. A. Boulton, O. Beijbom, and E. M. Wolff.  Covernet: Multimodal behavior prediction using trajectory sets.  In _Proceedings of the IEEE/CVF conference on computer vision and pattern recognition_ , pages 14074–14083, 2020.
  * Sadeghian et al. (2019) A. Sadeghian, V. Kosaraju, A. Sadeghian, N. Hirose, H. Rezatofighi, and S. Savarese.  Sophie: An attentive gan for predicting paths compliant to social and physical constraints.  In _Proceedings of the IEEE/CVF conference on computer vision and pattern recognition_ , pages 1349–1358, 2019.
  * Yousif and Müller (2022) Y. M. Yousif and J. P. Müller.  Generating explanatory saliency maps for mixed traffic flow using a behaviour cloning model.  In _International Workshop on Multi-Agent Systems and Agent-Based Simulation_ , pages 107–120. Springer, 2022.
  * Yue et al. (2022) J. Yue, D. Manocha, and H. Wang.  Human trajectory prediction via neural social physics.  In _European Conference on Computer Vision_ , pages 376–394. Springer, 2022.

##  Appendix A Representing Nearby Traffic

The future trajectory is also influenced by nearby entities, which need to be
compressed for better training and model interpretability. One approach to
achieve this is through direct measures related to movement such as Time to
collision (TTC) or Post encroachment alignment (PET). However, both of these
methods rely on time values while we require position predictions. Therefore,
in this section, the point of collision (POC) will be used as a measure. POC
assumes that two objects have constant velocity movement and predicts their
expected point of collision. To account for all possible scenarios, four
different directions of movement are assumed for the ego object, as shown in
figure 9. The mean of other objects’ POC with those directions is calculated
to give one point per direction. If no point is found for a particular
direction, a default mid-point will be returned.

Lastly, it should be noted that since one coordinate is redundant given the
movement direction; only four values representing the x coordinates of these
four points are taken as input features.

Figure 9: Finding the points of collision for four directions

##  Appendix B EBM Hyperparameters

List of the hyperparameters of EBM in this work are shown in table 2, such as
maximum number of leafs, and maximum training round, learning rate, maximum
feature bins, and outer bags.

Table 2: EBM Hyperparameters Hyperparameter | Value
---|---
Maximum feature bins | 256
Maximum interaction bins | 32
Maximum training rounds | 5000
Learning rate | 0.01
Maximum number of leaf | 3
Outer bags | 8
Validation size | 15%

##  Appendix C Qualitative Examples of Local Explanation

Figure 10: Argo: Local Example of the best mode in one case. Road highlighted
according to X values (red increase it and blue decrease it) Figure 11: Argo:
Local Example of the best mode in one case. Road highlighted according to X
values (red increase it and blue decrease it) Figure 12: Argo: Local Example
of the best mode in one case. Road highlighted according to X values (red
increase it and blue decrease it) Figure 13: SDD: Local Example for one point
prediction Figure 14: SDD: Local Example for one point prediction Figure 15:
SDD: Local Example for one point prediction Figure 16: InD: Local Example for
one point prediction

  *[GAM]: Generalized Additive Models
  *[EBM]: Explainable Boosting Machines
  *[SDD]: Stanford drone dataset
  *[InD]: Intersection Drone Dataset
  *[NAM]: Neural Additive Models


# First measurements with monolithic active pixel test structures produced in
a 65 nm CMOS process

M. Buckland11footnotetext: Now at STFC Daresbury Laboratory, Keckwick Ln,
Warrington WA4 4AD, UK.

###### Abstract

The Inner Tracking System (ITS) of the ALICE experiment at CERN will undergo
an upgrade during the LHC long shutdown 3, in which the three innermost
tracking layers will be replaced. This upgrade, named the Inner Tracking
System 3 (ITS3), employs stitched wafer-scale Monolithic Active Pixel Sensors
fabricated in a 65 CMOS process. The sensors are 260 in length and thinned to
less than 50 then bent to form truly half-cylindrical half-barrels. The
feasibility of this process for the ITS3 was explored with the first test
production run (MLR1) in 2021, whose goal was to evaluate the charged particle
detection efficiency and the sensor performance under non-ionising and
ionising radiation up to the expected levels for ALICE ITS3 of e13 (NIEL) and
10 (TID). Three sensor flavours were produced to investigate this process:
Analog Pixel Test Structure (APTS), Circuit Exploratoire 65 (CE65) and Digital
Pixel Test Structure (DPTS).

This contribution gives an overview of the MLR1 submission and test results,
describing the different sensor flavours and presenting the results of the
performance measurements done with particle beams for various chip variants
and irradiation levels.

##  1 Introduction

The importance of Monolithic Active Pixel Sensors (MAPS) for use in high
energy physics experiment vertex and tracking detectors has been established
in the last decade. The most recent implementation of a large-area MAPS
detector was for the ALICE Inner Tracking System 2 (ITS2) [1] at CERN, which
used the ALPIDE chip [2, 3]. This chip was fabricated in the TowerJazz 180
CMOS process and showed excellent performance in terms of detection efficiency
(≫99) and spatial resolution (about 5).

During the LHC Long Shutdown 3 (2026-2028), the ITS2 will undergo an upgrade
called the ITS3, where the three innermost tracking layers will be replaced.
The ITS3 employs wafer-scale MAPS with a length of 260 that are thinned to <50
and bent to radii of 18, 24, and 30 to form cylindrical half-barrels. To
obtain sensors of this length, a process called stitching is employed. In this
process the reticles of the CMOS imaging process are joined together to
produce a larger single sensor, removing the need for electrical services in
the detector. By utilising the natural stiffness of the cylindrical geometry,
the majority of the mechanical support can also be removed, using only carbon
foam spacers as support structures. Finally, the very low power consumption
(<20^2) of the chip will allow the detector to be cooled by air. These
reductions in the material within the sensitive volume will enable the ITS3 to
have a very low material budget of <0.05 X0 per layer. Altogether, these
upgrades will provide exceptional tracking and vertexing capabilities leading
to an improvement in the pointing resolution of the current detector by a
factor of two.

Due to the challenging design that the ITS3 poses, the Tower Partners
Semiconductor Co. (TPSCo) $65\text{\,}\mathrm{nm}$ CMOS imaging process [4]
was chosen as the starting point. The first submission in the 65 CMOS process,
in conjunction with the CERN EP R&D on monolithic pixel sensors [5], was
called MLR1 and contains many different test structures to fully explore the
CMOS process.

##  2 MAPS in the 65 CMOS imaging process

By using a smaller CMOS process node compared to the $180\text{\,}\mathrm{nm}$
CMOS process of the ITS2, more possibilities will be open, such as:

  * •

the ability to have complete coverage along the beam axis (z-direction) in the
detector with a single stitched sensor produced on 300 wafers;

  * •

lower power consumption when moving to deeper submicron processes.

However, moving to a new process also provides challenges, such as optimising
the design for sensor yield, the charge collection, and testing and verifying
the radiation harness.

Three main sensor flavours produced are the Analog Pixel Test Structure
(APTS), Circuit Exploratoire 65 (CE65), and Digital Pixel Test Structure
(DPTS), each measuring 1.5x1.5 in size and highlighted in Fig. 1.

Figure 1: The MLR1 reticle floor plan highlighting the APTS (orange), CE65
(green) and DPTS (blue) chips.

In addition to three sensor flavours, three process options were also explored
to investigate the charge collection properties of the CMOS process. These
processes, shown in Fig. 2, are called standard, modified, and modified-with-
gap and are similar to those used in the $180\text{\,}\mathrm{nm}$ CMOS
process [6]. In the standard process, the epitaxial layer is only partially
depleted, so some of the charges will undergo diffusion, and this process is
expected to have the largest charge sharing. For the modified process, a low-
dose n-type implant is added across the length of the pixel. This enables the
epitaxial layer to be fully depleted and increases the lateral electric field
to the collection diode. Thus, it better collects the signal charge and
accelerates the carriers towards the collection diode. This lateral electric
field is further increased in the modified-with-gap process where a gap in the
low-dose n-implant is added at the edges of the pixel.

Figure 2: The three process options implemented in the MLR1 chips.

The main goals of the MLR1 submission were to verify that the detection
efficiency was 99 and that the radiation hardness could reach the expected
levels for the ITS3 (e13 and 10).

##  3 Analog Pixel Test Structure

The APTS incorporates a 6x6 pixel matrix with direct analogue readout on the
central 4x4 pixels. Two versions of the output buffer were implemented: a
source-follower (APTS-SF) and a fast operational amplifier (APTS-OA) whose
focus was on time resolution. In addition, the sensor was produced in four
different pixel pitches ranging from 10 to 25. The goal of the APTS was to
explore the different sensor designs and processes.

The in-beam measurements of the APTS-SF shown in Fig. 3 (left) demonstrate the
impact of the different process types on the detection efficiency. While all
three can reach the desired 99, it is clear that the standard process has a
reduced performance at larger thresholds compared to the other two due to the
improved charge collection in the modified processes. The modified-with-gap
process shows the largest detection efficiency over the whole measured
threshold range. Comparing the detection efficiency for different pitches,
Fig. 3 (right), it can be seen that below a threshold of 200, there is minimal
difference among the pitches. However, above this value, larger pitches result
in higher efficiencies. Timing measurements were also performed with in-beam
measurements of two APTS-OA, resulting in a timing resolution of ($77\pm 5$)
$\mathrm{ps}$, Fig. 4.

Figure 3: Comparison of the detection efficiency vs. threshold for an APTS-SF
sensor for the three processes (left, see text for details) and different
pixel pitches (right). Figure 4: The time residual distribution of two APTS-OA
sensors fitted with a Gaussian function to extract the timing resolution.

##  4 Circuit Exploratoire 65

The CE65 is a “large"-area chip consisting of an analogue rolling shutter
readout with an integration time of 50. One type of the pixel matrix consists
of 64x32 pixels implemented with a pixel pitch of 15 and split into three
subvariants that differ by their in-pixel amplifier: AC, DC or SF. The other
type of pixel matrix contains 48x32 pixels implemented with a pitch of 25. The
goal of the CE65 was to study the pixel matrix uniformity.

Figure 5 shows the seed pixel distributions obtained from in-beam measurements
of the different CE65 variants. There is a clear distinction between the
standard and the modified-with-gap processes, with the latter having a larger
most probable value (MPV), signifying a larger charge collection depth. The
difference between the amplifiers is minimal for the modified-with-gap
process. However, for the standard process, it can be seen that the AC-coupled
amplifier has a larger MPV.

Figure 5: The in-beam seed pixel distribution comparing the response of the
different CE65 variants.

##  5 Digital Pixel Test Structure

The DPTS features a 32x32 pixel matrix with a pitch of 15 implemented in the
modified-with-gap process and contains a full digital front-end with
asynchronous readout [7]. The sensor is controlled by a set of external
reference currents and voltages and read out via a current mode logic (CML)
output [8, 9]. All the pixels are read out simultaneously via a differential
digital output that time encodes the pixel position and Time-over-Threshold
(ToT). The goal of the DPTS was to study the in-pixel full-digital front end.

The performance of DPTS chips for various irradiation levels taken at a
temperature of +$20\text{\,}\mathrm{\SIUnitSymbolCelsius}$ was evaluated using
in-beam measurements, the results of which are shown in Fig. 6. For all
irradiation levels, the sensor shows an excellent detection efficiency of 99
and a spatial resolution below the binary resolution (pixel pitch /
$\sqrt{12}$) while preserving a fake-hit rate below $10$ pixel-1 s-1. For the
detection efficiency, it can be seen that non-ionising irradiation leads to a
decrease in detection efficiency while ionising irradiation leads to an
increase in the fake-hit rate. For the spatial resolution, there is a
negligible impact of the irradiation on the performance of the tested
irradiation levels. Whereas the average cluster size shows a slight decrease
with increasing non-ionising dose.

Figure 6: The detection efficiency and fake-hit rate vs. threshold (top) and
the spatial resolution and average cluster size vs. threshold (bottom) for
DPTS chips irradiated to various levels [7].

To investigate the cause of the detection efficiency loss in the sensor, the
detection efficiency was studied as a function of the in-pixel position as
shown in Fig. 7 for a sensor irradiated to e15. It can be seen that the
further the track is away from the collection diode, in the centre of the
pixel, the smaller the detection efficiency is, an effect becoming
particularly acute in the corners of the pixel.

Figure 7: The in-pixel detection efficiency for a DPTS sensor irradiated to
e15. The grey circle represents the tracking resolution of the telescope used
to reconstruct the tracks [7].

##  6 Summary

The performance of the MLR1 chips was evaluated through extensive
characterisation in the laboratory and with in-beam measurement. The
measurements show that the MLR1 was a success thanks to the large number of
operational prototypes that allow the parameter space of the CMOS process to
be mapped out. Furthermore, the MLR1 structures exhibit excellent performance
in terms of detection efficiency (>99) and spatial resolution (3-4.5
$\mathrm{\SIUnitSymbolMicro m}$) from the in-beam measurements for all three
sensor flavours considered. The radiation hardness is demonstrated by the
sensors maintaining a detection efficiency of 99 for chips irradiated with a
dose at the expected ITS3 levels, e13 (NIEL) and 10 (TID). The radiation
hardness actually exceeds this goal as the desired performance is maintained
even for sensors irradiated up to e15 and operated at
+$20\text{\,}\mathrm{\SIUnitSymbolCelsius}$. In addition, the APTS-OA has
demonstrated a time resolution of ($77\pm 5$) $\mathrm{ps}$.

From the results of the MLR1, the detection efficiency and radiation hardness
have been validated and represent an important milestone in the R&D for ALICE
ITS3. The next step towards a wafer-scale bent sensor is the second submission
in the 65 process designated ER1, whose goal is the validation of stitching
and yield via the full-scale sensor prototypes.

## References

  * [1] ALICE Collaboration, _Alice upgrades during the lhc long shutdown 2_ , 2302.01238.
  * [2] M. Mager, _ALPIDE, the Monolithic Active Pixel Sensor for the ALICE ITS upgrade_ , _NIM-A_ 824 (2016) 434.
  * [3] G. Aglieri Rinella, _The ALPIDE pixel sensor chip for the upgrade of the ALICE Inner Tracking System_ , _NIM-A_ 845 (2017) 583.
  * [4] Tower Partners Semiconductor Co. _http://www.towersemi.com/_ .
  * [5] CERN EP R&D _https://ep-rnd.web.cern.ch/topic/monolithic-pixel-detectors_ .
  * [6] W. Snoeys et al., _A process modification for CMOS monolithic active pixel sensors for enhanced depletion, timing performance and radiation tolerance_ , _NIM-A_ 871 (2017) 90.
  * [7] G. Aglieri Rinella et al., _Digital pixel test structures implemented in a 65 nm cmos process_ , _NIM-A_ 1056 (2023) 168589.
  * [8] L. Cecconi et al., _Design and readout architecture of a monolithic binary active pixel sensor in TPSCo 65 nm CMOS imaging technology_ , _J. Inst._ 18 (2023) C02025.
  * [9] F. Piro et al., _A compact front-end circuit for a monolithic sensor in a 65 nm cmos imaging technology_ , _IEEE Transactions on Nuclear Science_ (2023) 1.