#  $0^{+}$ to $2^{+}$ neutrinoless double-$\beta$ decay of 76Ge, 82Se, 130Te
and 136Xe in the microscopic interacting boson model

J. Ferretti<EMAIL_ADDRESS>Department of Physics, University of
Jyväskylä, PO Box 35, FI-40014, Jyväskylä, Finland  R. Magaña Vsevolodovna
<EMAIL_ADDRESS>INFN, Sezione di Genova, via Dodecaneso 33, 16146
Genova (Italy)  J. Kotila<EMAIL_ADDRESS>Department of Physics,
University of Jyväskylä, PO Box 35, FI-40014, Jyväskylä, Finland  Finnish
Institute for Educational Research, University of Jyväskylä, P.O. Box 35,
FI-40014 Jyväskylä, Finland  Center for Theoretical Physics, Sloane Physics
Laboratory, Yale University, New Haven, Connecticut 06520-8120, USA  E.
Santopinto<EMAIL_ADDRESS>INFN, Sezione di Genova, via
Dodecaneso 33, 16146 Genova (Italy)

(August 30, 2024)

###### Abstract

Here, we study the neutrinoless double-$\beta$ ($0\nu\beta\beta$) decay
between the ground state and the first $2^{+}$ state of
${}^{76}\mbox{Ge}\rightarrow{}^{76}\mbox{Se}$,
${}^{82}\mbox{Se}\rightarrow{}^{82}\mbox{Kr}$,
${}^{130}\mbox{Te}\rightarrow{}^{130}\mbox{Xe}$ and
${}^{136}\mbox{Xe}\rightarrow{}^{136}\mbox{Ba}$ systems. The relevant nuclear
matrix elements (NMEs) involved in the process are calculated within the
formalism of the microscopic interacting boson model (IBM-2). The IBM-2 has
been widely used to obtain predictions for nuclear observables, such as the
spectrum, but also to explore the possible emergence of beyond-the-Standard
Model effects in the weak interactions of nuclei. Our calculations are carried
out by considering the exchange of a Majorana neutrino between two nucleons
($2N$-mechanism). In addition to NMEs, we calculate the associated leptonic
phase-space factors (PSFs) using electron radial wave functions, which are
obtained by solving numerically the Dirac equation of a screened Coulomb
potential that takes into account finite nuclear size. By combining our IBM-2
results for the NMEs with those for the PSFs along with experimental half-life
limits, we can set limits on the $\langle\lambda\rangle$ and
$\langle\eta\rangle$ couplings of left-right (L-R) models.

##  I INTRODUCTION

Neutrinos have a long story. Their existence was postulated by Pauli in 1930
to ensure the conservation of energy and angular momentum in $\beta$-decay
Pauli:1930pc . Fermi’s renowned theory of beta decay dates back to 1933
Fermi:1934hr . In 1956, neutrinos were first observed at Los Alamos by Cowan
and Reines via the study of inverse beta decay Cowan:1992xc . Several decades
later, neutrinos are still fascinating and mysterious particles.

Important questions regarding some of their main properties remain unsolved,
including the unknown mechanism that generates their masses and a complete
understanding of their mixing mechanism and mass hierarchy Giunti:2003qt .
Because of the lack in the standard model (SM) of a Yukawa coupling between
the Higgs boson and neutrinos, due to the absence of right-handed neutrinos,
the SM has to be extended to provide a neutrino mass term. Extensions of the
SM include the L-R symmetric Pati:1974yy ; Mohapatra:1974hk ;
Senjanovic:1975rk  and SUSY Mohapatra:1986su ; Vergados:1986td ; Hirsch:1995zi
; Babu:1995vh ; Hirsch:1995cg ; Faessler:1996ph  models.

Some important issues are directly related to the nature of neutrinos as
Fermi- or Majorana-type particles, a nature which could be directly assessed
via the experimental observation of neutrinoless double-beta
($0\nu\beta\beta$) decay process DellOro:2016tmg . However, despite of the
strenuous attempts by many experimental groups, e.g., Gando:2012zm ;
Agostini:2013mzu ; Alfonso:2015wka ; Auger:2012ar ; Arnold:2015wpy ;
Aalseth:2017btx ; Azzolini:2019yib , $0\nu\beta\beta$-decay has not yet been
observed.

Several theoretical investigations on $0\nu\beta\beta$ decay have been
published over the years (for a review see e.g., Refs. ejiri2019 ;
Agostini:2022zub ) in order to guide the experimentalists in their searches.
Owing to the low-energy character of $0\nu\beta\beta$ processes, these studies
necessarily involve elements of both particle and nuclear physics. In
particular, nuclear structure models are necessary in order to take care of
the nuclear matrix elements (NMEs) Barea:2013bz ; Simkovic:2007vu ;
Caurier:2007wq  entering the expression of the $0\nu\beta\beta$-decay half-
life.

Here, we show the results of a calculation of the $0^{+}$ to $2^{+}$
$0\nu\beta\beta$ decay of 76Ge, 82Se, 130Te and 136Xe, in which we consider
the exchange of a Majorana neutrino between two nucleons, the so-called
$2N$-mechanism, within an L-R symmetric model Doi:1985dx ; Tomoda:1990rs .
Specifically, in our study: I) we compute the relevant NMEs within the
microscopic interacting boson model (IBM-2) formalism Iachello:2006fqa  and
compare our results with previous calculations for the studied nuclei within
different nuclear structure models; II) we calculate the leptonic phase-space
factors (PSFs) by means of electron radial wave functions, obtained by solving
numerically the Dirac equation of a screened Coulomb potential that takes into
account finite nuclear size Kotila:2012zza ; and III) by combining the two
above elements, namely the results for the NMEs and the leptonic PSFs, with
the experimental limits on the half-life, we set limits on the
$\langle\lambda\rangle$ and $\langle\eta\rangle$ couplings of L-R models.
Experimental studies on this decay mode can be found e.g., in HM ;
PhysRevC.103.015501 ; LUCIFER:2015ozh ; CUPID:arxiv ; Arnaboldi:2002te ;
KamLAND-Zen:2015tnh . Previous calculations for the $0^{+}\rightarrow 2^{+}$
decay rate for the studied nuclei within different nuclear structure models
and via the $2N$-mechanism can be found in Refs. Tomoda:1990rs ; Doi:1985dx ;
Tomoda:1999zc ; Fang:2021tfc ; Tomoda:1988ew .

This article is organized as follows: In Sec. II the importance of
$0^{+}\rightarrow 2^{+}$ $0\nu\beta\beta$ decay is discussed and some details
on the calculation of $0\nu\beta\beta$ decay rates in L-R symmetric models are
provided. In Sec. III the calculation of IBM-2 wave functions is briefly
summarized, and in Sec. IV the decay operators needed for the description of
$0^{+}\rightarrow 2^{+}$ $0\nu\beta\beta$-decay are presented. In Sec. V the
numerical results for the ingredients needed for the calculation of the decay
rate are given and discussed. Finally, the conclusions are presented in Sec.
VI.

##  II $0^{+}\rightarrow 2^{+}$ $0\nu\beta\beta$ decay in the L-R symmetric
model

If we restrict ourselves only to long-range mechanisms for $0\nu\beta\beta$
decay, the most general effective Lagrangian is the Lorenz-invariant
combination of leptonic, $j_{\alpha}$, and hadronic, $J_{\alpha}$, currents
with definite tensor structure and chirality Kotila:2021xgw ; Ali:2007ec ;
Pas:1999fc ,

${\mathcal{L}}=\frac{G_{\rm F}\cos\theta_{\rm
c}}{\sqrt{2}}\left[U_{ei}\leavevmode\nobreak\
j_{V-A}^{\mu,i}J_{V-A,\mu}^{\dagger}+\displaystyle\sum_{\alpha,\beta}\epsilon_{\alpha,i}^{\beta}\leavevmode\nobreak\
j_{\beta}^{i}J_{\alpha}^{\dagger}+{\rm h.c.}\right]\mbox{ }.$ (1)

Here, $G_{\rm F}$ is the Fermi constant; $\theta_{\rm c}$ is the Cabibbo
angle; the hadronic and leptonic currents are defined as
$J_{\alpha}^{\dagger}=\bar{u}\mathcal{O}_{\alpha}d$ and
$j_{\beta}^{i}=\bar{e}\mathcal{O}_{\beta}\nu_{i}$, where the index $i$ spans
the neutrino mass eigenstates. The indices $\alpha,\beta$ are $V\mp A$, $S\mp
P$, $T\mp T_{5}$, where $S$, $P$, $T$ and $T_{5}$ stand for scalar, pseudo-
scalar, tensor, and pseudo-tensor, respectively. In Eq. (1) we have isolated
the standard model contribution proportional to $U_{ei}$, where $U_{ei}$ is
the PMNS mixing matrix element Pontecorvo:1957qd ; Maki:1962mu , from non-
standard contributions, which are those proportional to the couplings
$\epsilon_{\alpha,i}^{\beta}$.

By isolating the $V\mp A$ currents of the L-R models from those allowed in
other types of mechanisms, such as SUSY, and performing a non-relativistic
reduction of both the leptonic and hadronic currents, one can obtain the
expression of the $0\nu\beta\beta$ decay half-life for a $0^{+}\rightarrow
0^{+}$ transition in L-R models:

$\begin{array}[]{lll}\left[\tau_{1/2}^{0\nu}(0^{+}\rightarrow
0^{+})\right]^{-1}&=&C_{mm}^{(0)}\left(\frac{\left\langle
m_{\nu}\right\rangle}{m_{e}}\right)^{2}+C_{\lambda\lambda}^{(0)}\left\langle\lambda\right\rangle^{2}\\\
&+&C_{\eta\eta}^{(0)}\left\langle\eta\right\rangle^{2}+2C_{m\lambda}^{(0)}\frac{\left\langle
m_{\nu}\right\rangle}{m_{e}}\left\langle\lambda\right\rangle\\\
&+&2C_{m\eta}^{(0)}\frac{\left\langle
m_{\nu}\right\rangle}{m_{e}}\left\langle\eta\right\rangle+2C_{\lambda\eta}^{(0)}\left\langle\lambda\right\rangle\left\langle\eta\right\rangle.\end{array}$
(2)

The above equation is a complicated combination of the three parameters
$\langle m_{\nu}\rangle$, $\langle\lambda\rangle$, and $\langle\eta\rangle$
and their respective matrix elements and phase-space factors; for details on
the combinations $C_{ij}^{(0)},i,j=m,\lambda,\eta$ see e.g., Ref.
Kotila:2021xgw .

The new physics beyond the standard model is enclosed in the three parameters

$\langle m_{\nu}\rangle=\sum_{i}U_{ei}^{2}m_{i}\mbox{ },$ (3a)
$\langle\lambda\rangle=\lambda\sum_{j}U_{ej}V_{ej}\mbox{ },$ (3b) and
$\langle\eta\rangle=\eta\sum_{j}U_{ej}V_{ej}\mbox{ },$ (3c)

where $\langle m_{\nu}\rangle$ is the average neutrino mass obtained by
summing over mass $m_{i}$ of neutrino species $i$, and $\langle\lambda\rangle$
and $\langle\eta\rangle$ are the standard couplings of L-R models
Tomoda:1990rs . $U_{ej}$ and $V_{ej}$ are the mixing matrix elements of the
PMNS matrix Pontecorvo:1957qd ; Maki:1962mu  for the standard and non-standard
(L-R) mechanisms, respectively. Thus, the $0^{+}\rightarrow 0^{+}$ process can
occur because either of right- or left-handed leptonic currents.

By contrast, if $0\nu\beta\beta$-decay to a $2^{+}$ state is observed then in
addition to proving the Majorana character of neutrinos, the existence of
$V+A$ current would also be established, since as a first approximation, this
decay mode is triggered by right-handed leptonic currents only Tomoda:1988ew .
To be more specific, the combination of the lowest electron partial waves for
the $0^{+}\rightarrow 2^{+}$ transition is the $S$ and $P_{3/2}$ case since
the total angular momentum of the two-electron system should be 2. However, in
order to have a non-zero contribution of the L-L term that is proportional to
the average neutrino mass, $\langle m_{\nu}\rangle$, the leading term requires
the combination of $S$ and $D_{3/2}$ electron waves, making it negligible
compared with the contributions due to L-R terms, which are those proportional
to $\langle\lambda\rangle$ and $\langle\eta\rangle$). For more details see
(Doi:1985dx, , App. C).

Therefore, for the $0^{+}\rightarrow 2^{+}$ case the half-life can safely be
written without the dependence on the average neutrino mass

$\begin{array}[]{lll}\left[\tau_{1/2}^{0\nu\beta\beta}(0^{+}\rightarrow
2^{+})\right]^{-1}&=&g_{\rm
A}^{4}\Big{[}G_{1}\left(M_{\lambda}\langle\lambda\rangle\right.-\left.M_{\eta}\langle\eta\rangle\right)^{2}\\\
&+&G_{2}\left(M^{\prime}_{\eta}\langle\eta\rangle\right)^{2}\Big{]}\end{array}\mbox{
},$ (4)

where one can factorize the leptonic phase-space factor (PSF) Tomoda:1990rs ;
Kotila:2012zza , $G$, the nuclear matrix elements (NMEs), $M$, and the axial
vector coupling constant, $g_{\rm A}$. Neutrinoless double-beta decay to a
$2^{+}$ state thus provides information that is different from that one could
gather from the study of $0^{+}\rightarrow 0^{+}$ $0\nu\beta\beta$ processes.
Moreover, the observation of a $0^{+}\rightarrow 2^{+}$ $0\nu\beta\beta$-decay
may also possibly rule out those non-standard mechanisms in which no right-
handed gauge bosons or fermions are present.

Nucleus | $\epsilon_{d_{\nu}}$ | $\epsilon_{d_{\pi}}$ | $\kappa$ | $\chi_{\nu}$ | $\chi_{\pi}$ | $\xi_{1}$ | $\xi_{2}$ | $\xi_{3}$ | $c_{\nu}^{(0)}$ | $c_{\nu}^{(2)}$ | $c_{\pi}^{(0)}$ | $c_{\pi}^{(2)}$ | $c_{\pi}^{(4)}$
---|---|---|---|---|---|---|---|---|---|---|---|---|---
${}^{76}\mbox{Ge}$ Duval:1983tit  | 1.20 | 1.20 | -0.21 | 1.00 | -1.20 | -0.05 | 0.10 | -0.05 |  |  |  |  |
${}^{76}\mbox{Se}$ Kaup83  | 0.96 | 0.96 | -0.16 | 0.50 | -0.90 |  |  | -0.10 |  |  |  |  |
${}^{82}\mbox{Se}$ Kaup83  | 1.00 | 1.00 | -0.28 | 1.14 | -0.90 |  |  | -0.10 |  |  |  |  |
${}^{82}\mbox{Kr}$ Kaup79  | 1.15 | 1.15 | -0.19 | 0.93 | -1.13 | -0.10 |  | -0.10 |  |  |  |  |
${}^{130}\mbox{Te}$ Sambataro82  | 1.05 | 1.05 | -0.20 | 0.90 | -1.20 | -0.18 | 0.24 | -0.18 | 0.30 | 0.22 |  |  |
${}^{130}\mbox{Xe}$ Puddu80  | 0.76 | 0.76 | -0.19 | 0.50 | -0.80 | -0.18 | 0.24 | -0.18 | 0.30 | 0.22 |  |  |
${}^{136}\mbox{Xe}$111GS parameters fitted to reproduce the spectroscopic data of the low-lying energy states. |  | 1.31 |  |  |  |  |  |  |  |  | -0.04 | 0.01 | -0.02
${}^{136}\mbox{Ba}$ Puddu80  | 1.03 | 1.03 | -0.23 | 1.00 | -0.90 | -0.18 | 0.24 | -0.18 | 0.30 | 0.10 |  |  |
Table 1: Hamiltonian parameters employed in the IBM-2 calculation, along with
their references. All the values are in MeV, with the exception of those of
$\chi_{\pi}$ and $\chi_{\nu}$, which are dimensionless. The IBM-2 parameters
not shown here are set to zero.

##  III IBM-2 nuclear wave functions

IBM-2 is a nuclear structure model and was originally introduced as a
phenomenological approach to describing collective excitations in nuclei
Arima:1976ky . Soon afterwards, its relation with the shell model was
established Arima:1977vie ; Otsuka:1978zza ; Otsuka:1978zz . The starting
point of IBM-2 calculations of any nuclear observable, including weak decay
rates, is to obtain the nuclear wave functions of the nuclei of interest.
Realistic nuclear wave functions are obtained by fitting the IBM-2 parameters
in order to reproduce the experimental energy levels and other nuclear
properties, such as electromagnetic transition rates, quadrupole, and magnetic
moments etc. Iachello:2006fqa ; NPBOS ; Duval:1983tit ; Kaup79 ; Kaup83 ;
Sambataro82 ; Puddu80 , and the relevant two-body operators are derived in the
IBM-2 formalism Arima:1977vie ; Otsuka:1978zz ; Barea:2009zza .

The IBM-2 Hamiltonian describing the spectra of even-even nuclei, which is
often used in literature, and which is general enough for the phenomenological
studies, reads Otsuka:1978zz ; Iachello:2006fqa

$\begin{array}[]{rcl}H_{\rm
B}&=&\epsilon_{d}(\hat{n}_{d_{\pi}}+\hat{n}_{d_{\nu}})-\kappa\left(Q_{\nu}^{\rm
B}\cdot Q_{\pi}^{\rm B}\right)\\\
&+&\frac{1}{2}\xi_{2}\left[\left(d_{\nu}^{\dagger}s_{\pi}^{\dagger}-d_{\pi}^{\dagger}s_{\nu}^{\dagger}\right)\cdot\left(\tilde{d}_{\nu}s_{\pi}-\tilde{d}_{\pi}s_{\nu}\right)\right]\\\
&+&\displaystyle\frac{1}{2}\sum_{\rho}\sum_{K=1,3}\xi_{K}\left[d_{\nu}^{\dagger}\times
d_{\pi}^{\dagger}\right]^{(K)}\cdot\left[\tilde{d}_{\pi}\times\tilde{d}_{\nu}\right]^{(K)}\\\
&+&\displaystyle\sum_{L=0,2,4}C^{\rho}_{L}\Big{(}[d^{\dagger}_{\rho}d^{\dagger}_{\rho}]^{(L)}\cdot[\tilde{d}_{\rho}\tilde{d}_{\rho}]^{(L)}\Big{)}\mbox{
}.\end{array}$ (5)

i | $M_{i}$ | $C_{\lambda i}$ | $C_{\eta i}$ | $C^{\prime}_{\eta i}$
---|---|---|---|---
1 | ${\bm{\sigma}}_{1}\cdot{\bm{\sigma}}_{2}[\hat{\bf r}_{12}\times\hat{\bf r}_{12}]^{(2)}h(r_{12})$ | $\frac{1}{3}$ | $\frac{1}{3}$ | –
2 | $[{\bm{\sigma}}_{1}\otimes{\bm{\sigma}}_{2}]^{(2)}h(r_{12})$ | $-\frac{2}{3}$ | $-\frac{2}{3}$ | –
3 | $\left[[{\bm{\sigma}}_{1}\times{\bm{\sigma}}_{2}]^{(2)}\times[\hat{\bf r}_{12}\otimes\hat{\bf r}_{12}]^{(2)}\right]^{(2)}h(r_{12})$ | $\sqrt{\frac{7}{3}}$ | $\sqrt{\frac{7}{3}}$ | –
4 | $[\hat{\bf r}_{12}\times\hat{\bf r}_{12}]^{(2)}h(r_{12})$ | $(g_{\rm V}/g_{\rm A})^{2}$ | $-(g_{\rm V}/g_{\rm A})^{2}$ | –
5 | $\left[({\bm{\sigma}}_{1}+{\bm{\sigma}}_{2})\times[\hat{\bf r}_{12}\times\hat{\bf r}_{12}]^{(2)}\right]^{(2)}h(r_{12})$ | $-\sqrt{\frac{3}{2}}(g_{\rm V}/g_{\rm A})$ | – | –
6 | $\left[({\bm{\sigma}}_{1}-{\bm{\sigma}}_{2})\times[\hat{\bf r}_{12}\times\hat{\bf r}_{+12}]^{(1)}\right]^{(2)}\frac{r_{+12}}{r_{12}}h(r_{12})$ | – | – | $\sqrt{\frac{1}{2}}(g_{\rm V}/g_{\rm A})$
7 | $\left[({\bm{\sigma}}_{1}-{\bm{\sigma}}_{2})\times[\hat{\bf r}_{12}\times\hat{\bf r}_{+12}]^{(2)}\right]^{(2)}\frac{r_{+12}}{r_{12}}h(r_{12})$ | – | – | $-\sqrt{\frac{3}{2}}(g_{\rm V}/g_{\rm A})$
Table 2: 0+ $\rightarrow$ 2+ $0\nu\beta\beta$ decay via the $2N$-mechanism.
Here, we enlist the two-body operators, $M_{i}$, giving the dominant
contribution, as well as the coefficients $C_{\lambda i}$, $C_{\eta i}$ and
$C^{\prime}_{\eta i}$ of the two-body operators Tomoda:1988ew . The
coefficients which are not given explicitly are null.

In the previous expression, $\hat{n}_{d_{\rho}}=d_{\rho}^{\dagger}d_{\rho}$
and

$Q_{\rho}^{\rm
B}=d_{\rho}^{\dagger}s_{\rho}+s_{\rho}^{\dagger}\tilde{d}_{\rho}+\chi_{\rho}[d_{\rho}^{\dagger}\times\tilde{d}_{\rho}]^{(2)}$
(6)

represent the $d$-boson number operators and the boson quadrupole operators
for the proton ($\rho=\pi$) and neutron ($\rho=\nu$) pairs, respectively;
$s_{\rho}^{\dagger}$ and $d_{\rho}^{\dagger}$ are $s_{\rho}$\- and
$d_{\rho}$-boson creation operators, and the modified $d_{\rho}$-boson
annihilation operator satisfies $\tilde{d}_{\rho,m}=(-1)^{m}d_{\rho,-m}$. The
third term on RHS of Eq. (5) is the so-called Majorana term, which is relevant
to the proton-neutron mixed symmetry, and has been considered, e.g., in the
context of the isovector collective motion of valence shells. The last term on
the RHS of Eq. (5) corresponds to the interaction between like bosons, and
consists of $L=0,2$ and 4 components, respectively.

A detailed description of the IBM-2 Hamiltonian is given in Refs.
Arima:1977vie  and NPBOS . The Hamiltonian parameters are taken from the
literature Duval:1983tit ; Kaup83 ; Kaup79 ; Sambataro82 ; Puddu80 ;
Barea:2013bz . The values of the Hamiltonian parameters, together with the
references from which they are taken, are reported in Table 1.

##  IV 0$\nu\beta\beta$ decay operators for 0${}^{+}\rightarrow 2^{+}$
transitions in the IBM-2

In the present study, we focus on the $2N$-mechanism discussed by Tomoda
Tomoda:1988ew  in the context of L-R models Pati:1974yy ; Mohapatra:1974hk ;
Senjanovic:1975rk ; Hirsch:1996qw ; Doi:1985dx . To do so, we need to consider
a specific set of $V^{(2)}_{i}=M_{i}C_{\kappa i}$ operators, where the index
$\kappa$ refers to the $\lambda$, $\eta$ and $\eta^{\prime}$ contributions of
L-R models. The coefficients $C_{\kappa i}$ and the corresponding two-body
operators $M_{i}$ are enlisted in Table 2.

The seven operators in this table can be written as a combination of three
parts: the relative coordinate, $O_{i}(\hat{\bf r}_{12})h(r_{12})$, the
center-of-mass coordinate, $O_{i}(\hat{\bf r}_{+12})f(r_{+12})$, and the spin
part, $O_{i}({\bm{\sigma}}_{1},{\bm{\sigma}}_{2})$, where the following
notation for the coordinates is used: ${\bf r}_{12}={\bf r}_{2}-{\bf r}_{1}$,
${\bf r}_{+12}={\bf r}_{2}+{\bf r}_{1}$, $\hat{\bf r}={\bf r}/|{\bf r}|$.

The neutrino potential, which comes from the electron-Majorana neutrino
exchange, is given by Tomoda:1986yz ; Vergados:1988xp

$h(r_{12})=-r_{12}\frac{\partial}{\partial
r_{12}}\,H(r_{12},\overline{A})\mbox{ }.$ (7)

In the above equation,

$\overline{A}=\langle
E_{N}\rangle-E(0^{+})+m_{e}+\frac{1}{2}Q_{\beta\beta}(2^{+})$ (8)

is the closure energy, where $E(0^{+})$ is the energy of the initial $0^{+}$
state, $\langle E_{N}\rangle$ the average energy of the intermediate excited
state, and $Q_{\beta\beta}(2^{+})$ is the $Q$-value of the $0^{+}\rightarrow
2^{+}$ decay process. The neutrino propagation function in Eq. (7),
$H(r_{12},\overline{A})$, is given by Tomoda:1986yz ; Vergados:1988xp

$H(r_{12},\overline{A})=\frac{4\pi}{(2\pi)^{3}}\int d{\bf
p}_{12}\frac{\text{exp}(i{\bf p}_{12}\cdot{\bf
r}_{12})}{p_{12}(p_{12}+\overline{A})}\mbox{ },$ (9)

where ${\bf p}_{12}$ is the conjugate momentum to the ${\bf r}_{12}$
coordinate. More details on the neutrino potential can be found in e.g.,
Barea:2009zza .

On introducing a proton (neutron) creation (annihilation) operator
$\pi^{\dagger}_{nljm}(\tilde{\nu}_{nljm})$ that acts on the single-particle
state $|nljm\rangle$, the second quantized fermion operator can be written as

$\begin{array}[]{lll}V^{(\lambda)}_{i}&=&\displaystyle-\frac{1}{4}\sum_{j_{1},j_{2}}\sum_{j_{1^{\prime}},j_{2^{\prime}}}\sum_{J,J^{\prime}}(-1)^{J+J^{\prime}}\sqrt{1+(-1)^{J}\delta_{j_{1},j_{2}}}\sqrt{1+(-1)^{J^{\prime}}\delta_{j_{1^{\prime}},j_{2^{\prime}}}}\\\
&&\times{\mathcal{G}}_{i}(j_{1},j_{2},J,j_{1^{\prime}},j_{2^{\prime}},J^{\prime},\lambda)(\pi^{\dagger}_{j_{1}}\times\pi^{\dagger}_{j_{2}})^{(J)}(\tilde{\nu}_{j^{\prime}_{1}}\times\tilde{\nu}_{j^{\prime}_{2}})^{(J^{\prime})}\mbox{
},\end{array}$ (10)

with $J,J^{\prime}=0,2$ for the current study, and for $i=1-5$

$\begin{array}[]{lll}{\mathcal{G}}_{i}&=&\sqrt{\frac{2}{3}}\displaystyle\sum_{kk^{\prime}}\sum_{k_{1}k_{2}}i^{k_{1}-k_{2}+\lambda_{2}}\frac{\hat{k}_{1}^{2}\hat{k}_{2}^{2}}{\hat{\lambda}_{2}^{2}}\,\langle
k_{1}0k_{2}0\lambda_{2}0\rangle\,v^{k_{1},k_{2};\lambda_{2}}(r_{1},r_{2})\,\hat{k}\hat{k}^{\prime}\hat{\lambda}_{1}\hat{\lambda}_{2}\left\\{\begin{array}[]{cccc}s_{1}&k_{1}&k\\\
s_{2}&k_{2}&k^{\prime}\\\
\lambda_{1}&\lambda_{2}&\lambda\end{array}\right\\}\hat{J}\hat{\lambda}\hat{J}^{\prime}\left\\{\begin{array}[]{cccc}j_{1}&j_{2}&J\\\
j_{1}^{\prime}&j_{2}^{\prime}&J^{\prime}\\\
k&k^{\prime}&\lambda\end{array}\right\\}\\\
&&\times\hat{j}_{1}\hat{k}\hat{j}_{1}^{\prime}\left\\{\begin{array}[]{cccc}\frac{1}{2}&l_{1}&j_{1}\\\
\frac{1}{2}&l_{1}^{\prime}&j_{1}^{\prime}\\\
s_{1}&k_{1}&k\end{array}\right\\}\hat{j}_{2}\hat{k}^{\prime}\hat{j}_{2}^{\prime}\left\\{\begin{array}[]{cccc}\frac{1}{2}&l_{2}&j_{2}\\\
\frac{1}{2}&l_{2}^{\prime}&j_{2}^{\prime}\\\
s_{2}&k_{2}&k\end{array}\right\\}\langle\frac{1}{2}\|\Sigma^{(s_{1})}\|\frac{1}{2}\rangle\,(-1)^{k_{1}}\hat{l}_{1}\langle
l_{1}0k_{1}l_{1}^{\prime}0\rangle\langle\frac{1}{2}\|\Sigma^{(s_{2})}\|\frac{1}{2}\rangle(-1)^{-k_{2}}\\\
&&\times\hat{l}_{2}\langle l_{2}0k_{2}l_{2}^{\prime}0\rangle
R^{(k_{1}k_{2}\lambda_{2})}(n_{1},l_{1},n_{2},l_{2},n_{1}^{\prime},l_{1}^{\prime},n_{2}^{\prime},l_{2}^{\prime})\mbox{
}.\end{array}$ (11)

In Eq. (11), $\Sigma^{(1)}={\bm{\sigma}}$, $\Sigma^{0}=1$,
$\langle\frac{1}{2}\|\Sigma^{(s)}\|\frac{1}{2}\rangle=\sqrt{2(2s+1)}$,
$|l_{i}-l^{\prime}_{i}|\leq k_{i}\leq l_{i}+l^{\prime}_{i}$,
$|j_{1}-j^{\prime}_{1}|\leq k\leq j_{1}+j^{\prime}_{1}$, and
$|j_{2}-j^{\prime}_{2}|\leq k^{\prime}\leq j_{2}+j^{\prime}_{2}$. Also, an
additional factor of $-\sqrt{3}$ is needed for ${\mathcal{G}}_{1}$ and of
$-\sqrt{3/2}$ for ${\mathcal{G}}_{2}$. ${\mathcal{G}}_{5}$ is evaluated in two
parts by applying additional factor of $1/3\sqrt{1/6}$. In our current study
the protons and neutrons occupy the same major shell and thus the
contributions of ${\mathcal{G}}_{6}$ and ${\mathcal{G}}_{7}$ vanish. Regarding
the radial integrals, indicated as
$R^{(k_{1}k_{2}\lambda_{2})}(n_{1},l_{1},n_{2},l_{2},n_{1}^{\prime},l_{1}^{\prime},n_{2}^{\prime},l_{2}^{\prime})$
in Eq. (11), their evaluation follows the procedure discussed in Horie1961 .

The calculation of the nuclear matrix elements of $0\nu\beta\beta$ decay
would, in principle, proceed by going through all the virtual intermediate
states in the odd-odd nucleus. However, this is a demanding task, which can be
greatly simplified by treating the sum over the intermediate states in the
closure limit. This is a good approximation in the case of $0\nu\beta\beta$
decay, as the energy of the virtual neutrino exchange between nucleons is much
larger than the typical excitation energy of the intermediate states
Doi:1985dx ; Tomoda:1990rs . In the closure approximation, one is left with
the calculation of two-body NMEs between even-even initial and final states.

In order to obtain the bosonic image of the fermionic $0\nu\beta\beta$ decay
operator, we proceed in a similar way to Ref. Barea:2009zza . The following
step is then to define a mapping between the IBM-2 $J^{P}=0^{+}$ and $2^{+}$
boson creation operators, $s^{\dagger}$ and $d^{\dagger}$, and the shell-model
creation operators of collective nucleon pairs,
$S^{\dagger}=\left(\rho_{j}^{\dagger}\times\rho_{j}^{\dagger}\right)^{(0)}$
and
$D^{\dagger}=\left(\rho_{j}^{\dagger}\times\rho_{j^{\prime}}^{\dagger}\right)^{(2)}_{M}$,
where the fermion operators $\rho^{\dagger}_{j}$ create nucleons (either
neutrons, $\rho=\nu$, or protons, $\rho=\pi$) with angular momentum $j$
Arima:1977vie ; Otsuka:1978zz ; Barea:2009zza . This procedure enables us to
find a direct correspondence between the matrix elements between fermionic
states in the $SD$ shell-model subspace and the matrix elements in the $sd$
bosonic space of the IBM-2. One has

$\left(\rho_{j}^{\dagger}\times\rho_{j}^{\dagger}\right)^{(0)}\rightarrow
A_{\rho}(j)\,s_{\rho}^{\dagger}$ (12a) and
$\left(\rho_{j}^{\dagger}\times\rho_{j^{\prime}}^{\dagger}\right)^{(2)}_{M}\rightarrow
B_{\rho}(j,j^{\prime})\,d_{\rho,M}^{\dagger}\mbox{ },$ (12b)

where the mapping coefficients $A_{\rho}(j)$ and $B_{\rho}(j,j^{\prime})$ are
obtained by means of the OAI method Otsuka:1978zz  and depend on the specific
normalization of the nuclear structure coefficients that one considers. Our
choice is to use the conventions for the $A_{\rho}(j)$ and
$B_{\rho}(j,j^{\prime})$ coefficients reported in Refs. Mardones:2016wgy ;
Barea:2014lza , which are based on the procedure for diagonalizing the Surface
Delta Interaction (SDI) of Ref. Pittel:1982  and the use of the commutator
method of Refs. Frank:1982zz ; Lipas:1990rs ; Barea:2009zza .

##  V 0${}_{1}^{+}$ to 2${}_{1}^{+}$ $0\nu\beta\beta$ decay

###  V.1 Nuclear matrix elements

Here, we give results for the NMEs relevant to the 0${}_{1}^{+}$ to
2${}_{1}^{+}$ $0\nu\beta\beta$ decay processes of 76Ge, 82Se, 130Te and 136Xe.
The nuclear matrix elements of the operators $M_{i}$ (with $i=1,...,7$) in
Table 2 are computed in the IBM-2 formalism Arima:1977vie ; Otsuka:1978zz ;
Barea:2009zza . The calculations can be made more realistic by introducing
short-range correlation effects, i.e. the two-body operators in Table 2 need
to be multiplied by the short-range correlation function, $f(r)$, squared.
Following Ref. Barea:2009zza , we make use of the Jastrow function,

$f(r)=1-ce^{-ar^{2}}(1-br^{2})\mbox{ },$ (13)

with Argonne parametrization $a=1.59$ fm-2, $b=1.45$ fm-2 and $c=0.92$ argonne
.

The finite size of the nucleon is taken into into account by substituting the
coupling constants $g_{\rm A}$ and $g_{\rm V}$ with the form factors,

$g_{\rm V}(p_{12}^{2})=\frac{g_{\rm V}}{\left(1+\frac{p_{12}^{2}}{M_{\rm
V}^{2}}\right)^{2}}\mbox{ }$ (14a) and $g_{\rm A}(p_{12}^{2})=\frac{g_{\rm
A}}{\left(1+\frac{p_{12}^{2}}{M_{\rm A}^{2}}\right)^{2}}\mbox{ }.$ (14b)

In the above equations, the constant $M_{\rm V}^{2}=0.71$ GeV2 is fixed by the
electromagnetic form factor of the nucleon Iachello:1972nu ; Bijker:2004yu
and the value of $g_{\rm V}=1$ by the conserved vector current (CVC)
hypothesis; the value of $M_{\rm A}^{2}=1.09$ GeV2 is determined from neutrino
scattering data Schindler:2006jq  and that of $g_{\rm A}$ from neutron decay
Zyla:2020zbs .

| $M_{1}$ | $M_{2}$ | $M_{3}$ | $M_{4}$ | $M_{5}$ | $M_{6}$ | $M_{7}$ | $M_{\lambda}$ | $M_{\eta}$ | $M^{\prime}_{\eta}$
---|---|---|---|---|---|---|---|---|---|---
76Ge | 0.189 | -0.056 | -0.023 | -0.069 | -0.013 | 0 | 0 | 0.035 | 0.108 | 0
82Se | 0.003 | 0.080 | -0.003 | -0.007 | -0.011 | 0 | 0 | -0.051 | -0.053 | 0
130Te | 0.153 | -0.081 | -0.016 | -0.050 | -0.006 | 0 | 0 | 0.056 | 0.112 | 0
136Xe | 0.058 | -0.112 | -0.006 | -0.012 | -0.0001 | 0 | 0 | 0.077 | 0.092 | 0
Table 3: Nuclear matrix elements for 0+ $\rightarrow$ 2+ $0\nu\beta\beta$
decay via the $2N$-mechanism obtained by using IBM-2.

Our IBM-2 results for 76Ge, 82Se, 130Te, and 136Xe are given in Table 3. In
the last three columns, we give the combined NMEs

$\begin{array}[]{l}M_{\lambda}=\displaystyle\sum_{i=1}^{5}C_{\lambda
i}M_{i}\,,\mbox{ }\,M_{\eta}=\displaystyle\sum_{i=1}^{5}C_{\eta i}M_{i}\,,\\\
M^{\prime}_{\eta}=\displaystyle\sum_{i=6}^{7}C^{\prime}_{\eta
i}M_{i}\,,\end{array}$ (15)

where the values of the coefficients $C_{\lambda i}$, $C_{\eta i}$ and
$C^{\prime}_{\eta i}$ are given in Table 2. The values in Table 3 are
calculated by using unquenched values of $g_{V}=1$ and $g_{A}=1.269$, and
quenching can be implemented through coefficients $C$.

It is worth noting that, in the nuclei considered here, protons and neutrons
occupy the same major shell and thus the contributions of $M_{6}$ and $M_{7}$
vanish, leading to $M^{\prime}_{\eta}=0$. One also notices particularly small
$M_{1}$ value for the case of 82Se originating from the small bosonic
$d^{\dagger}_{\pi}\tilde{d}_{\nu}$ matrix element. The same also happens for
the first excited $0^{+}$ state in 82Se decay, leading to small IBM-2
${}^{82}\mbox{Se}(0^{+}_{1})\rightarrow{}^{82}\mbox{Kr}(0^{+}_{2})$
$0\nu\beta\beta$ NME, as can be seen from Ref. Barea:2015kwa .

Our results can be compared to those of the few existing studies on this
topic. Specifically, the $0^{+}\rightarrow 2^{+}$ $0\nu\beta\beta$ decay of
${}^{76}{\rm Ge}$ was studied in Refs. Tomoda:1988ew  and Fang:2021tfc  by
means of projected the Hartree Fock Bogoliubov (PHFB) method and quasiparticle
random-phase approximation (QRPA), respectively. We observe that the first
operator provides the largest contribution to both $M_{\lambda}$ and
$M_{\eta}$ in all of these three calculations. Our results,
$M_{\lambda}=0.035$ and $M_{\eta}=0.108$, stand in between the PHFB results
($M_{\lambda}=0.002$ and $M_{\eta}=0.061$) and QRPA results
($M_{\lambda}=0.008-0.228$ and $M_{\eta}=0.317-0.540$).

###  V.2 Leptonic phase-space integrals

The leptonic phase-space integrals, indicated as $G_{i}$ $(i=1,2)$ in Eq. (4),
are given by Tomoda:1988ew ; Tomoda:1990rs

$G_{i}=\frac{2}{\ln 2}\frac{(G_{\rm F}\cos\theta_{\rm
C})^{4}}{16\pi^{5}}\frac{m_{e}^{2}}{4R^{2}_{A}}\int^{Q_{\beta\beta}^{2^{+}}+m_{e}^{2}}_{m_{e}^{2}}f_{i}p_{1}p_{2}E_{1}E_{2}dE_{1},$
(16)

with $E_{2}=Q_{\beta\beta}^{2^{+}}+m_{e}^{2}-E_{1}$. In the above equation,
$G_{\rm F}=1.1663787(6)\cdot 10^{-5}$ GeV-2 is the Fermi coupling constant and
$\theta_{\rm C}$ the Cabibbo angle Zyla:2020zbs ; $E_{j}$ and
$p_{j}=\sqrt{E_{j}^{2}-m_{e}^{2}}$ are the energies and the asymptotic momenta
of the electrons, respectively;

$\begin{array}[]{rcl}f_{1}&=&\frac{3}{(m_{e}R_{A})^{2}}\,\left[\left|f^{-2-1}\right|^{2}+\left|f_{21}\right|^{2}+\left|f^{-1-2}\right|^{2}\right.\\\
&+&\left.\left|f_{12}\right|^{2}\right]\end{array}$ (17a) and
$\begin{array}[]{rcl}f_{2}&=&\frac{3}{(m_{e}R_{A})^{2}}\,\left[\left|f^{-2}{}_{1}\right|^{2}+\left|f^{-1}{}_{2}\right|^{2}+\left|f_{1}{}^{-2}\right|^{2}\right.\\\
&+&\left.\left|f_{2}{}^{-1}\right|^{2}\right]\end{array}$ (17b)

are combinations of electron wave functions as defined in Ref. (Tomoda:1990rs,
, Appendix 1).

In the calculation of $G_{1,2}$, we have used electron radial wave functions
obtained via a numerical solution of the Dirac equation with potential
Kotila:2012zza ; radial

$\displaystyle V(r)=\begin{cases}-\alpha Z_{\rm
F}\frac{3-(r/R_{A})^{2}}{2R_{A}}\times\varphi(r)\mbox{ },&r<R\mbox{ },\\\
-\frac{\alpha Z_{\rm F}}{r}\times\varphi(r)\mbox{ },&r\geq R\mbox{
},\end{cases}$ (18)

which includes finite size corrections to the Coulomb potential of the final
nucleus with charge $Z_{\rm F}$ and electron screening, due to the electronic
cloud described in the Thomas-Fermi approximation by the function
$\varphi(r)$. The thus obtained values of the phase-space integrals in Eq.
(16) are given in Table 4.

|  $G_{1}$ [$10^{-15}$ yr-1] |  $G_{2}$ [$10^{-15}$ yr-1] |  $Q$-value [keV]
---|---|---|---
76Ge | 1.669 | 1.157 | 1479.9
82Se | 12.357 | 9.159 | 2221.4
130Te | 18.464 | 14.462 | 1991.4
136Xe | 8.611 | 6.269 | 1639.3
Table 4: Phase-space factors $G_{1}$ and $G_{2}$ for 0+ $\rightarrow$ 2+ $0\nu\beta\beta$ decay. The $Q$-values for the above decays are reported in the last column. They are obtained by subtracting the 0${}_{1}^{+}$-2${}_{1}^{+}$ energy splitting in the levels of the daughter nuclei from the corresponding $Q$-values of the standard $0^{+}\rightarrow 0^{+}$ processes from Refs. Rahaman:2007ng ; Redshaw:2009zz ; Lincoln:2012fq ; CUORE:2017tlq ; CUPID:2018npf ; EXO:2017poz .  Decay | Collaboration |  $\tau_{1/2}^{0\nu,\rm exp}$ [y] | $|\langle\lambda\rangle|$ | $|\langle\eta\rangle|$
---|---|---|---|---
76Ge | Majorana PhysRevC.103.015501  | $>2.1\times 10^{24}$ | < $3.0\cdot 10^{-4}$ | $<9.7\cdot 10^{-5}$
82Se | CUPID-0 CUPID:arxiv  | $>3.0\times 10^{23}$ | $<2.0\cdot 10^{-4}$ | $<1.9\cdot 10^{-4}$
130Te | Gran Sasso Arnaboldi:2002te  | $>1.4\times 10^{23}$ | $<2.2\cdot 10^{-4}$ | $<1.1\cdot 10^{-4}$
136Xe | KamLAND-Zen KamLAND-Zen:2015tnh  | $>2.6\times 10^{25}$ | $<1.7\cdot 10^{-5}$ | $<1.4\cdot 10^{-5}$
Table 5: Calculated limits on the $\langle\eta\rangle$ and
$\langle\lambda\rangle$ L-R model couplings. These limits are estimated by
comparing our IBM-2 results for the 0+ $\rightarrow$ 2${}_{1}^{+}$ half-life,
computed by means of Eq. (4), with the experimental limits on the 76Ge, 82Se,
130Te and 136Xe half-life (90% C.L.).

Our results for PSFs are comparable to those of Ref. Tomoda:1990rs , where the
leptonic phase-space integrals were computed by making use of electron radial
wave functions approximated by their leading terms in a power series expansion
in $r$. As an example, converted to our notation, the resulting values for
76Ge read $G_{1}^{\rm pse}=1.865\times 10^{-15}$ yr-1 and $G_{2}^{\rm
pse}=1.296\times 10^{-15}$ yr-1. It is noteworthy that the values obtained
with approximate wavefunctions are slightly larger than the values of
$G_{1,2}$ reported in Table 4, as was also shown for the decays to $0^{+}$
states in Ref. Kotila:2012zza .

###  V.3 Limits on the $\langle\lambda$ and $\langle\eta\rangle$ couplings

By combining the calculated values of the leptonic PSFs in Table 4 with our
IBM-2 results for the NMEs in Table 3, we can use Eq. (4) to place limits on
the $\langle\lambda\rangle$ and $\langle\eta\rangle$ couplings in L-R models.

The upper limit on the value of the $\langle\lambda\rangle$ coupling in L-R
models is obtained by setting $\langle\eta\rangle$ to zero and equating Eq.
(4) to the experimental limit on the 0+ $\rightarrow$ 2${}_{1}^{+}$
$0\nu\beta\beta$ half-life of the mother nucleus. Analogously, by setting
$\langle\lambda\rangle=0$, one can implement the same procedure and obtain the
limit on the value of the $\langle\eta\rangle$ coupling. The experimental
results for the 0+ $\rightarrow$ 2${}_{1}^{+}$ $0\nu\beta\beta$ half-lives of
${}^{76}\mbox{Ge}$, ${}^{82}\mbox{Se}$, ${}^{130}\mbox{Te}$ and
${}^{136}\mbox{Xe}$ are extracted from Refs. PhysRevC.103.015501 ; CUPID:arxiv
; Arnaboldi:2002te ; KamLAND-Zen:2015tnh . Our upper limits on the absolute
values of the $\langle\eta\rangle$ and $\langle\lambda\rangle$ L-R model
parameters are reported in Table 5. The most stringent limits on these
parameters can be set by making use of the KamLAND-Zen experimental limits on
the 0+ $\rightarrow$ 2${}_{1}^{+}$ $0\nu\beta\beta$ half-life of 136Xe
KamLAND-Zen:2015tnh , leading to limits of $|\langle\lambda\rangle|<1.7\cdot
10^{-5}$ and $|\langle\eta\rangle|<1.4\cdot 10^{-5}$.

The dependence of the parameters $\langle\eta\rangle$ and
$\langle\lambda\rangle$ is also shown in Fig. 1, where limits on the
combination of the $\langle\lambda\rangle$ and $\langle\eta\rangle$ couplings
in L-R models are shown.

Figure 1: Limits on the combination of the $\langle\lambda\rangle$ and
$\langle\eta\rangle$ couplings in L-R models in the case of 76Ge (green dotted
line), 82Se (continuous blue line), 130Te (red dotted line), and 136Xe (purple
dot-dashed line). The parameter space outside of the parallel lines, the
shaded area, is excluded.

It is also very interesting to compare our results with those of the IBM-2
study of 0+ $\rightarrow$ 0+ $0\nu\beta\beta$ decays. While one can set even
more stringent limits on $\langle\eta\rangle$ and $\langle\lambda\rangle$ in
the case of 0+ $\rightarrow$ 0+ $0\nu\beta\beta$ transitions Kotila:2021xgw ,
namely $\sim 10^{-9}$ for $\langle\eta\rangle$ and $\sim 10^{-7}$ for
$\langle\lambda\rangle$, in this case one cannot disentangle the
$\langle\eta\rangle$ and $\langle\lambda\rangle$ dependencies from their
dependence on the $\langle m_{\nu}\rangle$ parameter. This is one of the
reasons why 0+ $\rightarrow$ 2${}_{1}^{+}$ $0\nu\beta\beta$ decay searches are
worth investigating and why 0+ $\rightarrow$ 2${}_{1}^{+}$ experimental
searches are conducted in parallel with those for standard 0+ $\rightarrow$ 0+
transitions.

##  VI Conclusions

We have computed the neutrinoless double-$\beta$ ($0\nu\beta\beta$) decay
nuclear matrix elements between the ground state and the first excited $2^{+}$
state of ${}^{76}\mbox{Ge}\rightarrow{}^{76}\mbox{Se}$,
${}^{82}\mbox{Se}\rightarrow{}^{82}\mbox{Kr}$,
${}^{130}\mbox{Te}\rightarrow{}^{130}\mbox{Xe}$ and
${}^{136}\mbox{Xe}\rightarrow{}^{136}\mbox{Ba}$ within the framework of the
microscopic Interacting Boson Model (IBM-2) Arima:1977vie ; Otsuka:1978zz ;
Barea:2009zza ; Iachello:2006fqa  by considering the exchange of a Majorana
neutrino between two nucleons ($2N$-mechanism) Doi:1985dx ; Tomoda:1999zc .
The IBM-2 formalism was widely used in the past to obtain results for nuclear
observables, including the spectrum, the electromagnetic and the weak decays,
but also to explore the possible emergence of beyond-the-Standard Model
effects in the weak interactions of nuclei; see e.g., Refs. Barea:2009zza ;
Barea:2015kwa ; Santopinto:2018nyt ; Graf:2018ozy ; Deppisch:2020ztt . Our
results for the 76Ge(0${}_{1}^{+}$) $\rightarrow$ 76Se(2${}_{1}^{+}$) NMEs
stand in between the PHFB results and QRPA results from the literature. For
other reported cases, to our knowledge, this is the first calculation.

We have also calculated the relevant leptonic phase-space integrals
numerically by making use of exact Dirac wave functions with finite nuclear
size and electron screening Kotila:2012zza  in order to set some limits on the
standard couplings of L-R models, $|\langle\lambda\rangle|$ and
$|\langle\eta\rangle|$. As in the case of decays to $0^{+}$ states, the PSFs
are found to be slightly smaller than previous values obtained with
approximate wavefunctions.

The most stringent limits on the parameters $|\langle\lambda\rangle|$ and
$|\langle\eta\rangle|$ can be obtained from the 0+ $\rightarrow$ 2${}_{1}^{+}$
$0\nu\beta\beta$ half-life of 136Xe KamLAND-Zen:2015tnh , leading to
$|\langle\lambda\rangle|<1.7\cdot 10^{-5}$ and $|\langle\eta\rangle|<1.4\cdot
10^{-5}$. While one can set even more stringent limits on $\langle\eta\rangle$
and $\langle\lambda\rangle$ in the case of 0+ $\rightarrow$ 0+
$0\nu\beta\beta$ transitions, in this case one cannot disentangle the
$\langle\eta\rangle$ and $\langle\lambda\rangle$ dependencies from their
dependence on the $\langle m_{\nu}\rangle$ parameter, making 0+ $\rightarrow$
2${}_{1}^{+}$ $0\nu\beta\beta$ decay searches worth investigating further.

###### Acknowledgements.

This work was supported by the Academy of Finland, Grant No. 314733, 320062,
345869, and INFN, Italy.

## References

  * (1) W. Pauli, Phys. Today 31N9, 27 (1978).
  * (2) E. Fermi, Z. Phys. 88, 161 (1934).
  * (3) C. L. Cowan, F. Reines, F. B. Harrison, H. W. Kruse and A. D. McGuire, Science 124, 103 (1956).
  * (4) C. Giunti and M. Laveder, Neutrino mixing, hep-ph/0310238.
  * (5) J. C. Pati and A. Salam, Phys. Rev. D 10, 275 (1974) Erratum: [Phys. Rev. D 11, 703 (1975)].
  * (6) R. N. Mohapatra and J. C. Pati, Phys. Rev. D 11, 566 (1975).
  * (7) G. Senjanovic and R. N. Mohapatra, Phys. Rev. D 12, 1502 (1975).
  * (8) R. N. Mohapatra, Phys. Rev. D 34, 3457 (1986).
  * (9) J. D. Vergados, Phys. Lett. B 184, 55 (1987).
  * (10) M. Hirsch, H. V. Klapdor-Kleingrothaus and S. G. Kovalenko, Phys. Rev. Lett. 75, 17 (1995); Phys. Lett. B 352, 1 (1995); Phys. Lett. B 403, 291 (1997); Phys. Rev. D 57, 1947 (1998); Nucl. Phys. Proc. Suppl. 62, 224 (1998).
  * (11) K. S. Babu and R. N. Mohapatra, Phys. Rev. Lett. 75, 2276 (1995).
  * (12) M. Hirsch, H. V. Klapdor-Kleingrothaus and S. G. Kovalenko, Phys. Lett. B 372, 181 (1996) Erratum: [Phys. Lett. B 381, 488 (1996)]; H. Päs, M. Hirsch and H. V. Klapdor-Kleingrothaus, Phys. Lett. B 459, 450 (1999).
  * (13) A. Faessler, S. Kovalenko, F. Šimkovic and J. Schwieger, Phys. Rev. Lett. 78, 183 (1997).
  * (14) S. Dell’Oro, S. Marcocci, M. Viel and F. Vissani, Adv. High Energy Phys. 2016, 2162659 (2016).
  * (15) M. Auger et al. [EXO-200 Collaboration], Phys. Rev. Lett. 109, 032505 (2012); J. B. Albert et al. [EXO-200 Collaboration], Nature 510, 229 (2014).
  * (16) A. Gando et al. [KamLAND-Zen Collaboration], Phys. Rev. Lett. 110, 062502 (2013); Phys. Rev. Lett. 117, 082503 (2016). Addendum: [Phys. Rev. Lett. 117, 109903 (2016)].
  * (17) M. Agostini et al. [GERDA Collaboration], Phys. Rev. Lett. 111, 122503 (2013); Nature 544, 47 (2017).
  * (18) K. Alfonso et al. [CUORE Collaboration], Phys. Rev. Lett. 115, 102502 (2015); C. Alduino et al. [CUORE Collaboration], Phys. Rev. Lett. 120, 132501 (2018).
  * (19) R. Arnold et al. [NEMO-3 Collaboration], Phys. Rev. D 92, 072011 (2015).
  * (20) C. E. Aalseth et al. [Majorana Collaboration], Phys. Rev. Lett. 120, 132502 (2018).
  * (21) O. Azzolini et al., Phys. Rev. Lett. 123, 262501 (2019).
  * (22) H. Ejiri, J. Suhonen, and K. Zuber, Phys. Rep. 797, 1 (2019).
  * (23) M. Agostini, G. Benato, J. A. Detwiler, J. Menéndez and F. Vissani, Toward the discovery of matter creation with neutrinoless double-beta decay, arXiv:2202.01787.
  * (24) F. Simkovic, A. Faessler, V. Rodin, P. Vogel and J. Engel, Phys. Rev. C 77, 045503 (2008).
  * (25) E. Caurier, J. Menéndez, F. Nowacki and A. Poves, Phys. Rev. Lett. 100, 052503 (2008).
  * (26) J. Barea, J. Kotila and F. Iachello, Phys. Rev. C 87, 014315 (2013).
  * (27) M. Doi, T. Kotani and E. Takasugi, Prog. Theor. Phys. Suppl. 83, 1 (1985).
  * (28) T. Tomoda, Rept. Prog. Phys. 54, 53 (1991).
  * (29) F. Iachello and A. Arima, The Interacting Boson Model, Cambridge University Press (1987).
  * (30) J. Kotila and F. Iachello, Phys. Rev. C 85, 034316 (2012).
  * (31) B. Maier, Nucl. Phys. B (Proc. Suppl.) 35, 358 (1994).
  * (32) I. J. Arnquist et. al. [Majorana Collaboration], Phys. Rev. C 103, 015501 (2021).
  * (33) J. W. Beeman et al. [LUCIFER Collaboration], Eur. Phys. J. C 75, 591 (2015).
  * (34) O. Azzolini et al. [CUPID-0 Collaboration], Search of the neutrino-less double beta decay of 82Se into the excited states of 82Kr with CUPID-0, arXiv:2206.05130.
  * (35) C. Arnaboldi et al., Phys. Lett. B 557, 167 (2003).
  * (36) K. Asakura et al. [KamLAND-Zen Collaboration], Nucl. Phys. A 946, 171 (2016).
  * (37) T. Tomoda, Phys. Lett. B 474, 245 (2000).
  * (38) D. L. Fang and A. Faessler, Phys. Rev. C 103, 045501 (2021).
  * (39) T. Tomoda, Nucl. Phys. A 484, 635 (1988).
  * (40) H. Päs, M. Hirsch, H. V. Klapdor-Kleingrothaus and S. G. Kovalenko, Phys. Lett. B 453, 194 (1999).
  * (41) A. Ali, A. V. Borisov and D. V. Zhuridov, Phys. Rev. D 76, 093009 (2007).
  * (42) J. Kotila, J. Ferretti and F. Iachello, Long-range neutrinoless double beta decay mechanisms, arXiv:2110.09141.
  * (43) B. Pontecorvo, Sov. Phys. JETP 7, 172 (1958) [Zh. Eksp. Teor. Fiz. 34, 247 (1957)].
  * (44) Z. Maki, M. Nakagawa and S. Sakata, Prog. Theor. Phys. 28, 870 (1962).
  * (45) A. Arima and F. Iachello, Annals Phys. 99, 253 (1976); 111, 201 (1978); 123, 468 (1979).
  * (46) A. Arima, T. Ohtsuka, F. Iachello and I. Talmi, Phys. Lett. 66B, 205 (1977).
  * (47) T. Otsuka, A. Arima and F. Iachello, Nucl. Phys. A 309, 1 (1978).
  * (48) T. Otsuka, A. Arima, F. Iachello and I. Talmi, Phys. Lett. A 76, 139 (1978).
  * (49) P. D. Duval, D. Goutte and M. Vergnes, Phys. Lett. B 124, 297 (1983).
  * (50) U. Kaup and A. Gelberg, Z. Physik A 293, 311 (1979).
  * (51) G. Puddu, O. Scholten, and T. Otsuka, Nucl. Phys. A 348, 109 (1980).
  * (52) M. Sambataro, Nucl. Phys. A 380, 365 (1982).
  * (53) U. Kaup, C. Mönkemeyer, and P. V. Brentano, Z. Phys. A 310, 129 (1983).
  * (54) T. Otsuka and N. Yoshida, User’s manual of the program NPBOS, Report No. JAERI-M 85-094, 1985.
  * (55) J. Barea and F. Iachello, Phys. Rev. C 79, 044301 (2009).
  * (56) M. Hirsch, H. V. Klapdor-Kleingrothaus and O. Panella, Phys. Lett. B 374, 7 (1996).
  * (57) T. Tomoda, A. Faessler, K. W. Schmid and F. Grümmer, Nucl. Phys. A 452, 591 (1986).
  * (58) J. D. Vergados, A. Faessler and T. Tomoda, Nucl. Phys. A 490, 556 (1988).
  * (59) H. Horie and K. Sasaki, Prog. Theor. Phys. 25, 475 (1961).
  * (60) J. Barea, C. E. Alonso and J. M. Arias, Phys. Lett. B 737, 205 (2014).
  * (61) E. Mardones, J. Barea, C. E. Alonso and J. M. Arias, Phys. Rev. C 93, 034332 (2016).
  * (62) S. Pittel, P. D. Duval and B. R. Barrett, Ann. Phys. 144, 168 (1982).
  * (63) A. Frank and P. Van Isacker, Phys. Rev. C 26, 1661 (1982).
  * (64) P. Lipas, M. Koskinen, H. Harter, R. Nojarov and A. Faessler, Nucl. Phys. A 509, 509 (1990).
  * (65) F. Šimkovic, A. Faessler, H. Müther, V. Rodin, M. Stauf, Phys. Rev. C 79, 055501 (2009).
  * (66) F. Iachello, A. D. Jackson and A. Lande, Phys. Lett. 43B, 191 (1973).
  * (67) R. Bijker and F. Iachello, Phys. Rev. C 69, 068201 (2004).
  * (68) M. R. Schindler and S. Scherer, Eur. Phys. J. A 32, 429 (2007).
  * (69) P. A. Zyla et al. [Particle Data Group], PTEP 2020, 083C01 (2020).
  * (70) J. Barea, J. Kotila and F. Iachello, Phys. Rev. C 91, 034304 (2015).
  * (71) F. Salvat, J. Fernadez-Varea, and W. Williamson Jr., Comp. Phys. Comm. 90, 151 (1995).
  * (72) S. Rahaman et al., Phys. Lett. B 662, 111 (2008).
  * (73) M. Redshaw, B. J. Mount, E. G. Myers and F. T. Avignone III, Phys. Rev. Lett. 102, 212502 (2009).
  * (74) D. L. Lincoln, J. D. Holt, G. Bollen, M. Brodeur, S. Bustabad, J. Engel, S. J. Novario, M. Redshaw, R. Ringle and S. Schwarz, Phys. Rev. Lett. 110, 012501 (2013).
  * (75) J. B. Albert et al. [EXO Collaboration], Phys. Rev. Lett. 120, 072701 (2018).
  * (76) C. Alduino et al. [CUORE Collaboration], Phys. Rev. Lett. 120, 132501 (2018).
  * (77) O. Azzolini et al. [CUPID Collaboration], Eur. Phys. J. C 78, 888 (2018).
  * (78) E. Santopinto, H. García-Tecocoatzi, R. I. Magaña Vsevolodovna and J. Ferretti, Phys. Rev. C 98, 061601 (2018).
  * (79) L. Graf, F. F. Deppisch, F. Iachello and J. Kotila, Phys. Rev. D 98, 095023 (2018).
  * (80) F. F. Deppisch, L. Graf, F. Iachello and J. Kotila, Phys. Rev. D 102, 095016 (2020).


# The Recursive Gauss-Newton Filter

Roaldje Nadjiasngar<EMAIL_ADDRESS>Michael Inggs
<EMAIL_ADDRESS>
###### Abstract

This paper presents a compact, recursive, non-linear, filter, derived from the
Gauss-Newton (GNF), which is an algorithm that is based on weighted least
squares and the Newton method of local linearisation. The recursive form
(RGNF), which is then adapted to the Levenberg-Maquardt method is applicable
to linear / nonlinear of process state models, coupled with the linear /
nonlinear observation schemes. Simulation studies have demonstrated the
robustness of the RGNF, and a large reduction in the amount of computational
memory required, identified in the past as a major limitation on the use of
the GNF.

###### keywords:

Gauss-Newton, filter, tracking, recursive

##  1 Introduction

The minimum variance algorithm has been used to estimate parameters from
batches of observations, accumulated over a defined period of time. The most
popular version of the minimum variance methods is the weighted least squares,
which are at the heart of adaptive filtering [1] [2]. The recursive least
squares (RLS) methods are efficient versions of the least squares approach,
and are applicable to estimation of future states from scalar input data
streams. However, recent studies [3] have seen the development of state space
recursive least squares (SSRLS) methods that show robustness in the estimation
of linear state space models. For the estimation of non-linear state space
models, a non-recursive filter called the Gauss-Newton filter (GNF) was
developed and has been successfully used in many applications [4] [5] . The
GNF algorithm is a combination of the Newton method of local linearization and
the least squares-like version of the minimum variance method[4]. It is used
to estimate process states that are governed by non-linear, autonomous,
differential equations, coupled with linear or non-linear observation schemes.
The GNF algorithm, although robust, requires significant processing power,
i.e. the amount of memory required. To improve the computational efficiency of
the GNF, studies of the use of Field Programmable Gate Arrays (FPGA) and other
co-processor technology have been made [6, 7]. Memory requirements were
identified in these studies as being the major stumbling block in
implementations on both on FPGA (low power and parallelism) and coprocessor
(ease of use) technology. This paper obtains a recursive form of the GNF with
zero memory. We then adapt the recursive filter to the Levenberg-Maquardt
method, renown for its robustness [8, 9, 10, 11, 12], widely used in non
linear curve fitting problems and neural networks algorithms. The contribution
of this paper is the derivation of a compact recursive form of the GNF that is
applicable to four major scenarios:

Case 1 : linear process dynamic and linear observation scheme.

Case 2 : linear process dynamic and non-linear observation scheme.

Case 3 : non-linear process dynamic and linear observation scheme.

Case 4 : non-linear process dynamic and non linear observation scheme.

The paper begins with an exposition of a state space model based on non-
linear, differential equations. This is followed, in Section 3, by the
derivation of a recursive GNF. In Section 4 we describe the adaptation of the
recursive equations of the filter to the Levenberg-Maquardt method. A complete
filter algorithm is presented. In Section 5, the state space situations to
which we can apply this new recursive form are demonstrated, with a look at
stability. We then demonstrate the power of the new recursive GNF in an
application to range and bearing only tracking of a manoeuvring target
(Section 6), before concluding with a summary of results achieved.

##  2 State space model based on non-linear differential equations

Consider the following autonomous, non-linear differential equation (DE)
governing the process state:

$DX(t)=F(X(t))$ (1)

in which $F$ is a non linear vector function of the state vector $X$
describing a process, such as the position of a target in space. We assume the
observation scheme of the process is a non-linear function of the process
state with expression :

$Y(t)=G(X(t))+v(t)$ (2)

where $G$ is a non-linear function of $X$ and $v(t)$ is a random Gaussian
vector. The goal is to estimate the process state from the given non-linear
state models. For linear differential equations (DEs), the state transition
matrix could be easily obtained. This, however, is not the case with non-
linear DEs. Nevertheless, there is a procedure, based on local linearization,
that enables us to get around this obstacle, which we will now present.

###  2.1 The method of local linearisation

The solution of the DE gives rise to infinitely many trajectories that are
dependent on the initial condition. However there will be one trajectory whose
state vector the filter will attempt to identify from the observations. We
assume that there is a known nominal trajectory with state vector $\bar{X}(t)$
that has the following properties:

  * 1.

$\bar{X}(t)$ satisfies the same DE as $X(t)$

  * 2.

$\bar{X}(t)$ is close to $X(t)$

The above-mentioned properties result in the following expression:

$X(t)=\bar{X}(t)+\delta X(t)$ (3)

where $\delta X(t)$ is a vector of time-dependent functions that are small in
relation to the corresponding elements of either $\bar{X}(t)$ or $X(t)$ . The
vector $\delta X(t)$ is called the perturbation vector and is governed by the
following DE (the derivation is shown in Appendix A):

$D(\delta X(t))=A(\bar{X}(t))\delta X(t)$ (4)

where $A(\bar{X}(t))$ is called a sensitivity matrix defined as follows:

$A(\bar{X}(t))=\left.\frac{\partial
F(X(t))}{\partial(X(t))}\right|_{\bar{X}(t)}$ (5)

.

Equation is thus a linear DE, with a time varying coefficient and has a the
following transition equation:

$\delta X(t+\zeta)=\Phi(t_{n}+\zeta,t_{n},\bar{X})\delta X(t)$ (6)

in which $\Phi(t_{n}+\zeta,t_{n},\bar{X})$ is the transition matrix from time
$t_{n}$ to $t_{n}+\zeta$ (increment $\zeta$). The transition matrix is
governed by the following DE:

$\frac{\partial}{\partial\zeta}\Phi(t_{n+\zeta},t_{n},\bar{X})=A(\bar{X}(t_{n}+\zeta))\Phi(t_{n+\zeta},t_{n},\bar{X})$
(7)

$\Phi(t_{n},t_{n},\bar{X})=I$ (8)

The transition matrix is a function of $\bar{X}(t)$ and can be evaluated by
numerical integration and in order to fill the values of
$A(\bar{X}(t_{n}+\zeta))$, $\bar{X}(t)$ has to be integrated numerically. We
will soon present a recursive algorithm that will avoid the computation of the
transition matrix. We have shown in this section that we can estimate the true
state of process by estimating the perturbation vector, which is governed by a
linear differential equation. The next task is to obtain a linear perturbation
observation from the non-linear observation scheme.

###  2.2 The observation perturbation vector

In this section we will adopt the notation $X_{n}$ and $Y_{n}$ for $X(t_{n})$
and $Y(t_{n})$ respectively. We define a simulated noise free observation
vector $\bar{Y}_{n}$ as follows:

$\bar{Y}_{n}=G(\bar{X}_{n})$ (9)

Subtracting $\bar{Y}_{n}$ from the actual observation $Y_{n}$ gives the
observation perturbation vector:

$\delta Y_{n}=Y_{n}-\bar{Y}_{n}$ (10)

In appendix A we show that the observation perturbation vector is related to
the state perturbation vector as follows:

$\delta Y_{n}=M(\bar{X}_{n})\delta X_{n}+v_{n}$ (11)

where $M(\bar{X}_{n})$ is the Jacobean matrix of G, evaluated at
$\bar{X}_{n}$. The matrix is also called the observation sensitivity matrix
and is defined as follows:

$M(\bar{X}_{n})=\left.\frac{\partial
F(X_{n})}{\partial(X_{n})}\right|_{\bar{X}_{n}}$ (12)

We now examine the sequence of observations.

###  2.3 Sequence of observation

We assume that $L+1$ observation are obtained with time stamps
$t_{n},t_{n-1},...,t_{n-L}$. Theses observations are assembled as follows :

$\left[\begin{array}[]{c}\delta Y_{n}\\\ \delta Y_{n-1}\\\ .\\\ .\\\ .\\\
\delta Y_{n-L}\end{array}\right]=\left[\begin{array}[]{c}M(\bar{X}_{n})\delta
X_{n}\\\ M(\bar{X}_{n-1})\delta X_{n-1}\\\ .\\\ .\\\ .\\\
M(\bar{X}_{n-L})\delta
X_{n-L}\end{array}\right]+\left[\begin{array}[]{c}v_{n}\\\ v_{n-1}\\\ .\\\
.\\\ .\\\ v_{n-L}\end{array}\right]$ (13)

Using the relationship:

$\delta X_{m}=\Phi(t_{m},t_{n},\bar{X})\delta X_{n}$ (14)

then, substituting Equation 13 the observation sensitity equation can be
written as:

$\mathbf{\delta Y}_{n}=\mathbf{T}_{n}\delta X_{n}+\mathbf{V}_{n}$ (15)

in which $\mathbf{T}_{n}$, the total observation matrix is defined as follows:

$\mathbf{T}_{n}=\left[\begin{array}[]{c}M(\bar{X}_{n})\\\
M(\bar{X}_{n-1})\Phi(t_{n-1},t_{n};\bar{X})\\\ .\\\ .\\\ .\\\
M(\bar{X}_{n-L})\Phi(t_{n-L},t_{n};\bar{X})\end{array}\right]$ (16)

The vectors $\mathbf{\delta Y}_{n}$ and $\mathbf{V}_{n}$ are large. The
solution of the equation can be obtained from the minimum variance estimation
as follows:

$\delta\hat{X}_{n}=(\mathbf{T}_{n}^{T}\mathbf{R}_{n}^{-1}\mathbf{T}_{n})^{-1}\mathbf{T}_{n}^{T}\mathbf{R}_{n}^{-1}\mathbf{\delta
Y}_{n}$ (17)

The estimate $\delta\hat{X}_{n}$ has a covariance matrix:

$S_{n}=(\mathbf{T}_{n}^{T}\mathbf{R}_{n}^{-1}\mathbf{T}_{n})^{-1}$ (18)

where $\mathbf{R}_{n}^{-1}$ is a block diagonal weight matrix, also called the
least squares weight matrix, but, in fact, if we define $R_{n}$ as the
covariance matrix of the the error vector $v_{n}$. Then $\mathbf{R}_{n}^{-1}$
is expressed as:

$\mathbf{R}_{n}^{-1}=\left[\begin{array}[]{cccccc}R_{n}^{-1}&0&.&.&.&0\\\
0&R_{n-1}^{-1}&&&&.\\\ .&&.&&&.\\\ .&&&.&&.\\\ .&&&&.\\\
0&.&.&.&0&R_{n-L}^{-1}\end{array}\right]$ (19)

In this section we arrived at a form of filter that uses the minimum variance
estimation initiated by Gauss and the local linearisation technique championed
by Newton, to estimate the estate of the process from the non linear
observation scheme. This filter is called Gauss-Newton filter (GNF) and is
described in detail in Morrison’s work [4, 13]. The GNF has been successfully
implemented in some practical applications:

[5] showing strong stability. The memory nature of the filter has made it
unattractive to researchers in the past, and even now, challenging [7] .
However recent developments have presented recursive form of the linear least-
squares for state space model [3] . We derive a recursive form of GNF using a
similar approach to M. B. Malik [3]. However, before we derive a recursive
form of the GNF filter, we rewrite the expression of $\mathbf{T}_{n}$ using
the backward differentiation:

$\Phi(t_{n-L},t_{n},\bar{X})=A(\bar{X}_{n-L})^{-1}\Phi(t_{n-L+1},t_{n},\bar{X})$
(20)

The expression is thus:

$\mathbf{\delta Y}_{n}\mathbf{T}_{n}=\left[\begin{array}[]{c}M_{0}\\\
M_{1}A_{1}\\\ M_{2}A_{2}\\\ .\\\ .\\\ .\\\ M_{L}A_{L}\end{array}\right]$ (21)

where

$A_{L}=\prod_{i=1}^{L}A(\bar{X}_{n-i})^{-1}$ (22)

and

$M_{L}=M(\bar{X}_{n-L})$ (23)

with

$A_{0}=I$ (24)

We now move to derive the Recursive Gauss Newton Filter in the next section.

##  3 The Recursive Gauss-Newton filter

To obtain the recursive form, we use an approach similar to M. B. Malik in
[3]. Suppose that the observations start arriving at $n=0$ and that all
initial values of the filter are available. In in order to maintain the filter
adaptiveness, a weight matrix function using a fading parameter $\lambda<1$ is
adopted, and is defined as follows:

$\mathbf{R}_{n}^{-1}=\left[\begin{array}[]{cccccc}R^{-1}&0&.&.&.&0\\\
0&\lambda R^{-1}&&&&.\\\ .&&.&&&.\\\ .&&&.&&.\\\ .&&&&.\\\
0&.&.&.&0&\lambda^{n}R^{-1}\end{array}\right]$ (25)

The following, further definitions are adopted:

$\mathbf{W}_{n}=\mathbf{T}_{n}^{T}\mathbf{R}_{n}^{-1}\mathbf{T}_{n}$ (26)

$\mathbf{\mathbf{\xi}}_{n}=\mathbf{T}_{n}^{T}\mathbf{R}_{n}^{-1}\delta\mathbf{Y}$
(27)

Resulting in:

$\delta\hat{X}_{n}=\mathbf{W}_{n}^{-1}\xi_{n}$ (28)

.

In the next section, the recursive update of the perturbation vector is
demonstrated.

###  3.1 The recursive update of $\mathbf{W}_{n}$

Using equation (21) and the definitions in equations (26) and (25) we have:

$\displaystyle\mathbf{W}_{n}$ $\displaystyle{}=$
$\displaystyle{}\sum_{j=1}^{L}\lambda^{j}R^{-1}\prod_{i=1}^{j}A(\bar{X}_{n-i})^{-T}M(\bar{X}_{n-j})^{T}$
(29)
$\displaystyle\times{}M(\bar{X}_{n-j})\prod_{i=0}^{j}A(\bar{X}_{n-i})^{-1}$
$\displaystyle{+}\>M(\bar{X}_{n})^{T}R^{-1}M(\bar{X}_{n})$

and

$\displaystyle\mathbf{W}_{n-1}$ $\displaystyle{}=$
$\displaystyle{}\sum_{j=1}^{L-1}\lambda^{j}R^{-1}\prod_{i=1}^{j}A(\bar{X}_{n-1-i})^{-T}M(\bar{X}_{n-1-j})^{T}$
(30)
$\displaystyle\times{}^{-1}M(\bar{X}_{n-1-j})\prod_{i=0}^{j}A(\bar{X}_{n-1-i})^{-1}$
$\displaystyle{+}\>M(\bar{X}_{n-1})^{T}R^{-1}M(\bar{X}_{n-1})$

Comparing equations (29) and (30) the following recursive equation is
obtained:

$\mathbf{W}_{n}=\lambda
A(\bar{X}_{n-1})^{-T}\mathbf{W}_{n-1}A(\bar{X}_{n-1})^{-1}+M(\bar{X}_{n})^{T}R^{-1}M(\bar{X}_{n})$
(31)

which is the discrete, quadratic, Lyapunov, difference equation.

###  3.2 The recursive form of $\mathbf{\mathbf{\xi}}_{n}$

Using equations (21) (25) (27) $\mathbf{\mathbf{\xi}}_{n}$ can be expressed
as:

$\displaystyle\xi_{n}$ $\displaystyle{}=$
$\displaystyle{}\sum_{j=0}^{L}\lambda^{j}R^{-1}\prod_{i=1}^{j}A(\bar{X}_{n-i})^{-T}M(\bar{X}_{n-j})^{T}\delta
Y_{n-j}$ (32) $\displaystyle{+}\>M(\bar{X}_{n})^{T}R^{-1}\delta Y_{n}$

and

$\displaystyle\xi_{n-1}$ $\displaystyle{}=$
$\displaystyle{}\sum_{j=1}^{L-1}\lambda^{j}R^{-1}\prod_{i=1}^{j}A(\bar{X}_{n-1-i})^{-T}M(\bar{X}_{n-1-j})^{T}$
(33) $\displaystyle{\times}\>\delta Y_{n-1-j}+M(\bar{X}_{n-1})^{T}R^{-1}\delta
Y_{n-1}$

Comparing equations (32) and (33) the following recursive equation is
obtained:

$\xi_{n}=\lambda A(\bar{X}_{n-1})^{-T}\xi_{n-1}+M(\bar{X}_{n})^{T}R^{-1}\delta
Y_{n}$ (34)

##  4 Adaptation to Levenberg and Maquardt

In order to guarantee local convergence of the recursive filter and also to
avoid the singularity of $\mathbf{W}_{n}$. we replace it by
$\mathbf{W}_{n}+\mu I$ as suggested by Levenberg and Maquardt. The presence of
the damping factor $\mu$ will have two effects:

  * 1.

for large value of $\mu$ the algorithm behaves as a steepest descent which is
ideal when the current solution is far from the local minimum. The convergence
will be slow but however guaranteed. We therefore have

$\delta\hat{X}_{n}=\frac{1}{\mu}\xi_{n}$ (35)

.

  * 2.

for $\mu$ very small the algorithm will behave as gauss newton with faster
convergence. The current step will be

$\delta\hat{X}_{n}=\mathbf{W}_{n}^{-1}\xi_{n}$ (36)

.

###  4.1 The Gain Ratio

The $\mu$ can be updated by the so called gain ratio. We consider the
following cost function which is

$E(\delta X_{n})=(\mathbf{\delta Y}_{n}-\mathbf{T}_{n}\delta
X_{n})^{T}{R}^{-1}(\mathbf{\delta Y}_{n}-\mathbf{T}_{n}\delta X_{n})$ (37)

The denominator of gain ratio is :

$E(0)-E(\delta X_{n})=\delta X_{n}^{T}(\xi_{n}+\mu\delta X_{n})$ (38)

We define :

$F(\delta X_{n})=(Y_{n}-G(\bar{X}_{n}+\delta
X_{n}))^{T}R^{-1}(Y_{n}-G(\bar{X}_{n}+\delta X_{n}))$ (39)

The gain ratio is therefore:

$\varrho=\frac{F(0)-F(\delta X_{n})}{E(0)-E(\delta X_{n})}$ (40)

A large value of $\varrho$ indicates that $E(\delta X_{n})$ is a good
approximation of $\bar{Y}$, and $\mu$ can be decreased so that the next
Levenberg-Marquardt step is closer to the Gauss-Newton step. If $\varrho$ is
small or negative then $E(\delta X_{n})$ is a poor approximation, then $\mu$
should be increased to move closer to the steepest descent direction. The
complete filter algorithm adapted from [9] is presented in Algorithm 1

$k:=0$;$\nu:=2$;$\bar{X}_{n}:=X_{n/n-1}$;

$\delta{Y}_{n}:={Y}_{n}-G(\bar{X}_{n})$;

${W}_{temp}=M(\bar{X}_{n})^{T}R^{-1}M(\bar{X}_{n})$;

$\mathbf{W}_{n}=\mathbf{W}_{n-1/n}+{W}_{temp}$;

$\xi_{temp}=M(\bar{X}_{n})^{T}R^{-1}\delta Y_{n}$;

$\xi_{n}=\xi_{n/n-1}+\xi_{temp}$;

$stop:=false$;$\mu=\tau*max(diag(\mathbf{W}_{n/n-1}))$;

While (not stop) and ($k\leq k_{max}$)

$k:=k+1$;

repeat;

solve $(\mathbf{W}_{n}+\mu I)\delta\hat{X}_{n}=\xi_{n}$;

if ($||\delta\hat{X}_{n}||\leq\varepsilon||\bar{X}_{n}||$)

stop:=true;

else

$X_{new}:=\bar{X}_{n}+\delta\hat{X}_{n}$;

$F(\delta X)=Y_{n}-G(X_{new})$;$F(0)=\delta Y_{n}^{T}R^{-1}\delta Y_{n}$;

$E(0)-E(\delta X_{n})=\delta X_{n}^{T}(\xi_{n}+\mu\delta X_{n})$;

$\varrho=\frac{F(0)-F(\delta X_{n})}{E(0)-E(\delta X)}$;

if $\varrho>0$

$\bar{X}_{n}=X_{new}$;

$\delta{Y}_{n}:={Y}_{n}-G(\bar{X}_{n})$;

${W}_{temp}=M(\bar{X}_{n})^{T}R^{-1}M(\bar{X}_{n})$;

$\mathbf{W}_{n}=\mathbf{W}_{n/n-1}+{W}_{temp}$;

$\xi_{temp}=M(\bar{X}_{n})^{T}R^{-1}\delta Y_{n}$;

$\xi_{n}=\xi_{n/n-1}+\xi_{temp}$;

$\mu=\mu*max(1/3,1-(2\varrho+1)^{3})$;$\nu:=2$;

else

$\mu:=\nu*\mu$;

$\nu:=2*\nu$;ssm

endif

endif

until($\varrho>0$)or(stop);

endwhile

$X_{n/n}=X_{new}$;

$X_{n/n+1}=\Phi(s)X_{n/n}$;

$\mathbf{W}_{n/n+1}=\lambda A({X}_{n/n})^{-T}\mathbf{W}_{n}A({X}_{n/n})^{-1}$;

$\xi_{n/n+1}=\lambda A({X}_{n/n})^{-T}\xi_{n}$;

Algorithm 1 L-M algorithm for tracking system

##  5 State Space Models

We will present four possible models to which the recursive GNF can be
applied:

  * 1.

Model 1, with linear process dynamic and linear observation scheme. In this
model the recursive formulation is similar to the derived forms except the
estimation is made directly for $X_{n}$ and that the observed perturbation
vector $\delta Y_{n}$ is replaced by the actual observation vector $Y_{n}$.The
sensitivity matrices in this case become the measurement and transition
matrices of the process. In this case the LM algorithm is not required.

  * 2.

Model 2, with linear process dynamic and non-linear observation scheme. The
recursive model of the filter remains the same except the state sensitivity
matrix becomes a the transition matrix of the process. The state perturbation
is estimated to obtain the estimate of the process state.

  * 3.

Model 3, with non-linear process dynamic and linear observation scheme. The
measurement sensitivity matrix has become the measurement matrix.

  * 4.

Model 4, with a non-linear process dynamic and non linear observation scheme.
The derived recursive form without any further modification is applicable to
this case.

###  5.1 Stability of the Recursive GNF

The matrix $\mathbf{W}_{n}$ is the inverse of of the covariance matrix of the
filter and is therefore positive definite. As a consequence the solution of
the derived discrete Ly5apunov equation in (31) is unique with the sensitivity
matrix being stable. The eigenvalues of the inverse of the sensitivity matrix
are within an open unit circle and therefore the stability of athe system is
ensured by having $\lambda<1$.

##  6 Simulation: Range and Bearing tracking

In these simulation studies, we consider an example of a vehicle executing
various manoeuvres. During turn manoeuvres of unknown constant turn rate, the
aircraft dynamic model is :

$X_{n}=\left[\begin{array}[]{ccccc}1&\frac{sin(\Omega
T)}{\Omega}&0&-(\frac{1-cos(\Omega T)}{\Omega})&0\\\ 0&cos(\Omega
T)&0&-sin(\Omega T)&0\\\ 0&\frac{1-cos(\Omega T)}{\Omega}&1&\frac{sin(\Omega
T)}{\Omega}&0\\\ 0&sin(\Omega T)&0&cos(\Omega T)&0\\\
0&0&0&0&1\end{array}\right]X_{n-1}+v_{n}$ (41)

where the state of the vehicle is $X_{n}=[x,\dot{x},y,\dot{y},\Omega]$, with
$x$,$y$ the position coordinates and $\dot{x}$,$\dot{y}$ their corresponding
velocity components.The process noise $v_{k}\sim\mathcal{N}(0,Q)$ with
covariance matrix
$Q=diag\left[\begin{array}[]{ccc}q1BB^{T}&q1BB^{T}&q2T\end{array}\right]$
where,

$BB^{T}=\left[\begin{array}[]{cc}\frac{T^{4}}{4}&\frac{T^{3}}{2}\\\
\frac{T^{3}}{2}&T^{2}\end{array}\right]$ (42)

When the vehicle moves at nearly constant velocity its dynamic model is:

$X_{n}=\left[\begin{array}[]{ccccc}1&T&0&0&0\\\ 0&1&0&0&0\\\ 0&0&1&T&0\\\
0&0&0&1&0\\\ 0&0&0&0&1\end{array}\right]X_{n-1}+v_{n}$ (43)

The vehicle is observed by a radar located at the origin of the plane, capable
of measuring the range $r$ and and the bearing angle $\theta$. The measurement
equation is therefore:

$\left[\begin{array}[]{c}r_{n}\\\
\theta_{n}\end{array}\right]=\left[\begin{array}[]{c}\sqrt{x^{2}+y^{2}}\\\
tan^{-1}(\frac{y}{x})\end{array}\right]+w_{n}$ (44)

where the measurement noise is $w_{k}\sim\mathcal{N}(0,R)$ with covariance

$R=diag\left[\begin{array}[]{cc}\sigma_{r}^{2}&\sigma_{\theta}^{2}\end{array}\right]$

The following constants were used for data generation: $T=1s$;
$\Omega=-3^{0}s^{-1}$; $q1=0.01$m${}^{2}\rm{s}^{-4}$; $q2=1.75\times
10^{-4}\rm{s}^{-4}$; $\sigma_{r}$=10m; $\sigma_{\theta}=\sqrt{0.1}$mrad.

The vehicle starts at true initial state $X_{n}$=[10m, 25ms-1,400m,
0ms-1,-3ms-1] and moves at nearly constant velocity for $100$s, Then it
executes a turn manoeuvre from time index $n=101$ to $n=150$. After the
manoeuvre, the vehicle’s velocity remains nearly constant from $n=151$ to
$n=250$. At $n=251$ it starts a new turn manoeuvre at rate $\Omega$=3ms-1
until $n=400$. Finally from $n=400$ to $n=500$ it moves at nearly constant
velocity. Figure [1] describes the complete trajectory of the vehicle.

The filter uses a single model of a constant velocity to track the entire
manoeuvre:

$A(X_{n})=\left[\begin{array}[]{ccccc}1&T&0&0&0\\\ 0&1&0&0&0\\\ 0&0&1&T&0\\\
0&0&0&1&0\\\ 0&0&0&0&1\end{array}\right]$ (45)

The initial value $\mathbf{W}_{-1/0}=10^{-2}I$, where $I$ is an identity
matrix. The filter parameters are the following $k_{max}=200$,
$\varepsilon=1\times 10^{-24}$, $\tau=1\times 10^{-3}$, $\lambda=0.4$ The
filter initial state is generated randomly and then ensuring that it has the
same sign as the true state. This procedure guarantees the local convergence
of the first estimate. The experiment was repeated for 250 Monte Carlo runs
and the root means squared error (RMSE) is used as a performance metric. The
position RMSE is computed using the following expression:

$RMSE=\sqrt{\frac{1}{N}\sum_{i=1}^{N}\left((x_{n}^{i}-\hat{x_{n}})^{2}+(y_{n}^{i}-\hat{y_{n}})^{2}\right)}$
(46)

where $(x_{n}^{i},y_{n}^{i})$ and $(\hat{x_{n}},\hat{y_{n}})$ true and
estimated position coordinates respectively. The velocity root mean square
error (RMSE) is computed similarly. Figures [2] and [3] show the RMSE of the
position and velocity respectively. The position RMSE is not affected by
different manoeuvres while the velocity RMSE shows variation from different
manoeuvre states. The average values of the damping factor after complete
cycles of iteration is presented in Figure [3]. The damping factor increases
rapidly at the transition between manoeuvres. The average number of iterations
$k$ at convergence from Figure [4] shows similar variations.

##  7 Conclusions

The GNF with memory combines the minimum variance estimation and the Newton
method of local linearisation to estimate the process true state. The
recursive form for the Gauss-Newton filter has been derived in one compact
form that can be applied to all the four state and observation linearity and
nonlinearity scenarios:

Case 1 : linear process dynamic and linear observation scheme.

Case 2 : linear process dynamic and non-linear observation scheme.

Case 3 : with non-linear process dynamic and linear observation scheme.

Case 4 : non-linear process dynamic and non linear observation scheme.

The Hessian matrix of the filter which is computed recursively is augmented by
a damping factor as suggested earlier by Levenberg-Maquardt for non linear
curve fitting problems. The new filter is therefore a combination of Newtons
steepest descent and the Gauss-newton, ensuring its robustness. The presence
of a forgetting factor in the filter equations renders it capable of tracking
manoeuvring targets with a single filter dynamic model.

## Appendix A

###  A.1 The differential equation governing $\delta X(t)$

Starting from:

$\delta X(t)=X(t)-\bar{X}(t)$ (47)

The differentiation rule is applied:

$D\delta X(t)=F(\bar{X}(t)+\delta X(t))-F(\bar{X}(t))$ (48)

Let $F$ be defined as follows :

$F=\left[\begin{array}[]{c}f_{1}\\\ .\\\ .\\\ .\\\ f_{n}\end{array}\right]$
(49)

Equation becomes:

$D\delta X(t)=\left[\begin{array}[]{c}f_{1}(\bar{X}(t)+\delta X(t))\\\ .\\\
.\\\ .\\\ f_{n}(\bar{X}(t)+\delta
X(t))\end{array}\right]-\left[\begin{array}[]{c}f_{1}(\bar{X}(t))\\\ .\\\ .\\\
.\\\ f_{n}(\bar{X}(t))\end{array}\right]$ (50)

The Taylor first order approximation is applied:

$\displaystyle D\delta X(t)$ $\displaystyle{}=$
$\displaystyle{}\left[\begin{array}[]{c}f_{1}(\bar{X}(t))\\\ .\\\ .\\\ .\\\
f_{n}(\bar{X}(t))\end{array}\right]+\left[\begin{array}[]{c}\nabla
f_{1}(\bar{X}(t))^{T}\\\ .\\\ .\\\ .\\\ \nabla
f_{n}(\bar{X}(t))^{T}\end{array}\right]\delta X(t)$ (67)
$\displaystyle{-}\>\left[\begin{array}[]{c}f_{1}(\bar{X}(t))\\\ .\\\ .\\\ .\\\
f_{n}(\bar{X}(t))\end{array}\right]$

The following relation is obtained :

$D\delta X(t)=A(\bar{X}(t))\delta X(t)$ (68)

Where:

$A(\bar{X}(t))=\left[\begin{array}[]{c}\nabla f_{1}(\bar{X}(t))^{T}\\\ .\\\
.\\\ .\\\ \nabla f_{n}(\bar{X}(t))^{T}\end{array}\right]=\left.\frac{\partial
F(X(t))}{\partial(X(t))}\right|_{\bar{X}(t)}$ (69)

###  A.2 The relation between $\delta X_{n}$ and $\delta Y_{n}$

$\delta Y_{n}=G(\bar{X}_{n}+\delta X_{n})-G(\bar{X}_{n})$ (70)

As direct consequence of A.1 the following relationship is obtained:

$\delta Y_{n}=M(\bar{X}_{n})\delta X_{n}+v_{n}$ (71)

##  Appendix B Figure captions list

Figure 1: Target complete trajectory with manoeuvres

Figure 2: The Position RMSE is unaffected by the manoeuvres.

Figure 3: The velocity RMSE varies with manoeuvres.

Figure 4: The damping factor shows sharp peaks at start of manoeuvres.

Figure 5: The number of iterations increases during manoeuvres.

The figure numbering appears in the same order as the figures in the pdf
document

## Acknowledgment

The authors would like to thank Dr Norman Morrison for his contribution during
the research that leads to obtaining a recursive form of GNF. Dr Morrison has
been working on the GNF throughout his career and even in his retirement is
enthusiastic in providing teaching and insights into the fundamentals of
filter Engineering.

## References

  * [1] Recursive Least Squares With Linear Constraints.
  * [2] T. Kailath, A. H. Sayed, B. Hassibi, Linear Estimation, 2nd Edition, Prentice Hall, New Jersey, USA, 2000.
  * [3] M. B. Malik, State-space recursive least-squares: Part i, Signal Processing 84 (2004) 1709–1718.
  * [4] N. Morrison, Filter engineering -a practical approach: The Gauss-Newton and polynomial filters, to be published.
  * [5] N. Morrison, R. T. Lord, M. R. Inggs, The Gauss-Newton algorithm applied to track-while-scan radar, in: Proceedings of the IET International Conference on Radar Systems (RADAR 2007), Institution for Engineering and Technology, 2007\.
  * [6] J.-P. da Conceicao, Accelerating Gauss-Newton filters on FPGAs, Master’s thesis, University of Cape Town (Dec. 2011).
  * [7] J. Milburn, Co-processor offloading applied to passive coherent location with doppler and bearing data, Master’s thesis, University of Cape Town - RRSG (February 2010).
  * [8] D. W. Marquardt, An algorithm for least-squares estimation of nonlinear parameters, Journal of the Society for Industrial and Applied Mathematics 11 (2) (1963) pp. 431–441.
  * [9] O. T. K. Madsen, H.B. Nielsen, Method of non-linear least squares problems, 2nd Edition, Informatics and Mathematical Modelling, Technical University of Denmark, 2004.
  * [10] B. G. Kermani, S. S. Schiffman, H. T. Nagle, Performance of the Levenberg-Marquardt neural network training method in electronic nose applications, Sensors and Actuators B: Chemical 110 (1) (2005) 13 – 22.  doi:10.1016/j.snb.2005.01.008.
URL http://www.sciencedirect.com/science/article/pii/S09254%00505000961

  * [11] E. Derya, beyli, Analysis of EEG signals by implementing eigenvector methods/recurrent neural networks, Digital Signal Processing 19 (1) (2009) 134 – 143.  doi:10.1016/j.dsp.2008.07.007.
URL http://www.sciencedirect.com/science/article/pii/S10512%00408001243

  * [12] V. Singh, I. Gupta, H. Gupta, ANN-based estimator for distillation using Levenberg-Marquardt approach, Engineering Applications of Artificial Intelligence 20 (2) (2007) 249 – 259.  doi:10.1016/j.engappai.2006.06.017.
URL http://www.sciencedirect.com/science/article/pii/S09521%9760600114X

  * [13] N. Morrison, Introduction to Sequential Smoothing and Prediction, McGraw-Hill Book Company, 1969.