# PID Controller Optimization for Low-cost Line Follower Robots

Samet Oguten1, Bilal Kabas1 Email<EMAIL_ADDRESS><EMAIL_ADDRESS>1Department of Electrical-Electronics Engineering,
Abdullah Gül University
Kayseri, Turkey

###### Abstract

In this paper, modification of the classical PID controller and development of
open-loop control mechanisms to improve stability and robustness of a
differential wheeled robot are discussed. To deploy the algorithm, a test
platform has been constructed using low-cost and off-the-shelf components
including a microcontroller, reflectance sensor, and motor driver. This paper
describes the heuristic approach used in the identification of the system
specifications as well as the optimization of the controller. The PID
controller is analyzed in detail and the effect of each term is explained in
the context of stability. Lastly, the challenges encountered during the
development of controller and robot are discussed. Code is available at:
_https://github.com/sametoguten/STM32-Line-Follower-with-PID_.

###### Index Terms:

line tracking, PID control, control systems

##  I Introduction

Navigation is an essential part of Robotics applications. To complete a
specific task, it is crucial for a robot to follow a specific path. In
advanced applications, robots are equipped with sophisticated sensors such as
depth cameras, radars and lidars for enhanced perception capabilities.
Availability of such sensors enables the development of more sophisticated and
high-fidelity localization and mapping algorithms. However, these sensors are
expensive. Some certain tasks requiring repetitive motion in a relatively
fixed environment can be done by cost-effective applications. For instance, a
line tracking robot having only an array of reflectance sensors [1], [2], [3],
[4] can be used to carry a specific load from one point to another in a
factory. The aim of this research is to develop a controller algorithm for a
differential wheeled robot to follow a specific path. There are two main
requirements. The first requirement is that the robot should complete the
track. The second one is that the completion time should be as short as
possible. Throughout this research, only off-the-shelf and cost-efficient
components have been used.

PID control stand for proportional, integral, and derivative. PID controller
is a closed-loop control system widely used in industry. It tries to bring the
system to the state specified by the input by calculating the error. Error
calculation is done by comparing the reference input with the output of the
system. If the proportional gain in the PID control is increased, the response
of the system to disturbances becomes faster. However, after a certain level,
increasing the proportional gain can cause the system to show an oscillatory
unstable behavior. The integral control in the PID controller is used to keep
track of the past errors in the system. Thus, steady-state error is
eliminated. However, increasing integral gain can cause overshooting. The
resulting overshoot can cause the system to show an unstable behavior. The
main effect of the derivative controller is that it minimizes overshooting and
speeds up the transient response. It can be thought of as a mechanism to
detect future errors in the system. As the integral control minimizes the
steady-state error at the cost of overshooting, derivative control can be used
to minimize that overshooting to increase the robustness of the system.
However, a derivative controller may cause the system to be unstable in the
presence of noise as differentiation amplifies high frequency signals.

In this research, PID controller has been used to develop a high-speed and
high-fidelity line tracking controller. The PID controller calculates errors
according to the data coming from the reflectance sensor array. The main
reason why PID control is used during the development of the line tracking
controller is to prevent overshooting, minimize steady-state error, and
oscillation that may occur in the system [5].

##  II Test Platform

###  II-A Mechanical Components

Fig. 1 shows the two line-tracking robot used as a test platform throughout
this research. Two low-cost 6V DC gear motors have been used. These motors
spin at 250 RPM with 6 - 7 V DC input voltage. Since they have a gearbox, they
generate enough torque required for faster responses. One of the advantages of
these motors is that they are low-cost and a higher input voltage can be
applied to reach a higher RPM value for short operation times. The DC motors
have been mounted to the frame by using two plastic parts and two screws for
each.

(a) Top view (b) Side view (c) Rear view

Figure 1: Top, side and rear views of the line-tracking robots.

In line-following robots, wheels play an important role in the robot’s
handling of the path. High roadholding is needed to prevent the robot from
slipping, especially in sharp turns. The wheels have been selected to minimize
slipping and increase the robot’s tracking capabilities of the path. Wheels
with a diameter of 65 mm have been used to obtain high torque and minimize the
effect of disturbances caused by the track.

The frame enables the motors, sensors, motor driver used in the line-following
robot to stay as a whole. The frame is made of a hard plastic with 174 mm
length and 110 mm width. It is an off-the-shelf and low-cost component. The
center of mass directly affects stability of the robot. By adjusting the
location of the battery on the chassis, the center of mass can be manipulated
since the heaviest component is battery. As the center of mass gets closer to
the front of the robot, then slipping increases. On the other hand, as the
center of mass gets closer to the rear, then the front side goes up and down
in sharp turns. Hence, the center of mass has been centered relative to the
center of the frame. Additionally, the distance between the ground and the
chassis has been kept short which is another aspect affecting stability
especially in sharp turns. In this way, it has also become easier to place the
QTR-8 reflectance sensor closer to the ground.

###  II-B Electronic Components

####  II-B1 Motor Driver

To drive the two DC motors, a motor driver module having L298N dual full-
bridge driver has been used. This motor driver module is supplied with 13 V DC
voltage. It also has a 5 V DC voltage regulator which is used to supply STM32.
There are 4 input pins that PWM signals are applied to so that the direction
of rotation and the speed of the motors can be controlled. The output pins are
dedicated for DC motors and can supply maximum 2 A current for each. The motor
driver module has been mounted in front of the robot.

####  II-B2 Microcontroller

To read sensor values and drive the DC motors accordingly, a microcontroller
is needed. In this research STM32F103C8 microcontroller has been used. The
system clock has been set to 72 MHz. APB2 timer clocks work at 72 MHz.
Prescaler is 79 and counter period is 999 so that a 900 Hz PWM signal is
generated. The duty cycle is manipulated adaptively by the PID controller. 4
pins of the microcontroller have been dedicated for PWM generation.

####  II-B3 Battery

6 NiMH rechargeable batteries along with 4 AA alkaline batteries in series are
used to power the robot. In this way, a total of 13 V DC input voltage is
obtained. The positive terminal of the battery is connected to the 12 V input
pin of the motor driver. It has been observed that at least 7 V input voltage
is required to acquire 5 V output voltage from the internal regulator of the
motor driver so that both STM32 and QTR-8 reflectance sensor work efficiently.
The total value of current drawn by the system is approximately 1 Ampere. A
great portion of this current is consumed by DC motors.

Figure 2: Circuit schematic of the line tracking robot.

####  II-B4 Sensor

QTR-8 sensor is used for the robot to follow the line. There are eight
reflectance sensors on this module. These sensors consist of infrared light
emitters, MOSFETs switching in the presence of infrared light and capacitors.
MOSFETs control the ground connection of capacitors. The 8 pins which
correspond to a sensor are connected to STM32’s GPIO pins as shown in Fig. 2
for writing and reading operations. Before a reading operation, capacitors
should be charged. For this, a logic 1 is written to the GPIO ports in output
mode. If MOSFETs are exposed to infrared light, capacitors are connected to
the ground, so they are discharged. Thus, after a few milliseconds the GPIO
pins are switched to input mode. If the reading is logic 0 for a specific
sensor, it means that there is a reflective surface. If the reading is logic
1, so there is a non-reflective, i.e. dark in colour, surface. This mechanism
allows the microcontroller to understand where the line to be followed is. It
should be noted that the distance between sensors and the track should not
exceed 6 mm. Greater distance causes longer time taken for capacitors to
discharge consequently reducing the sampling rate from sensors.

##  III System Characteristics

Instead of using differential drive kinematics, a heuristic approach [1] has
been carried out to identify the system characteristics. Before implementing a
PID controller, a unity feedback system shown in Fig. 3 was experimented on
the robot. In paths where the line is straight, starting with non-zero error,
the robot oscillates left and right by an equal amount. This implies that
there is zero steady-state tracking error for a straight track. Consequently,
it can be concluded that this system has at least one pole at origin. Then,
the robot has been placed in the curved region of the track starting with non-
zero position error. It has been observed that the system oscillates more
towards the opposite side of the curvature. This implies that there is a
steady-state error in following the curved path [1]. In this case, curvature
can be thought of as a ramp or parabolic disturbance. This analysis shows that
using an integral controller will add one more pole to the system enhancing
its line-tracking capability in curved paths.

Figure 3: Block diagram of the unity feedback system used in system
identification.

##  IV Controller Design

The general structure of the controller algorithm is shown in Fig. 4. First, a
reading is obtained from sensors. Based on this reading, a positional
information is obtained so that the controller knows where the line is. Then,
error is calculated and fed to the PID controller which manipulates the error
term used for updating motor speeds. However, before feeding new PWM values,
an open-loop control checks if sensors miss the line. Open-loop control is
required to handle the corners in the test track. It also improves system
performance in curved paths.

###  IV-A Calculation of the Positional Error

Error is calculated based on sensor data. Each sensors’ reading corresponds to
a position information. This position information is calculated by using (1).
The reference position value is 4500 indicating the middle of the sensor.
Error is calculated by (2).

$x=\frac{\sum_{n=1}^{8}1000\cdot S_{n}}{\sum_{n=1}^{8}S_{n}}$ (1)

$e=4500-\frac{\sum_{n=1}^{8}1000\cdot S_{n}}{\sum_{n=1}^{8}S_{n}}$ (2)

If the error value is negative, this indicates that the robot is sliding to
the right. If the error value is positive, the robot deviates from the line to
the left. PID controller manipulates the speed of the motors. If the
calculated error is negative, the speed of the motor on the left is reduced,
the speed of the motor on the right is increased, and by this, the robot goes
to the left side and eliminates the error. If the error value is positive, the
speed of the motor on the left is increased and the speed of the motor on the
right is decreased. Then the robot goes to the right side and the error is
eliminated.

###  IV-B Proportional Control

Equation (3) is the mathematical expression for proportional control. It has
been observed that when only P control is used, the robot responds faster to
errors occurring during line tracking. But when the _k P_ gain goes above a
certain value, oscillation has been found to be quite high. When the _k P_
gain is decreased, the oscillation decreases, but the robot response to errors
has been found to be slower. For this reason, derivative control has been
added to the system.

$u(t)=k_{P}e(t)$ (3)

###  IV-C Integral Control

When proportional control is used alone, it has been observed that a steady-
state error occurs. The mathematical expression for the integral control is
given by (4). Integral term represents the summation of the past errors.
Consequently, combining with _k I_ integral gain, integral feedback minimizes
steady-state tracking error by introducing a pole at origin [5].

$u(t)=k_{I}\int_{t_{0}}^{t}e(\tau)\,d\tau$ (4)

As investigated in the system identification stage, the unity feedback
implementation without PID controller revealed that this system already has at
least one pole at origin. However, the system is not able to follow the curved
paths without steady-state error. To minimize the steady-state tracking error
in curved paths, the integral term is added to the controller. However, it has
been ineffective to sum all the past errors beginning from the start since the
shape of the path is not always straight. To improve the response of the
integral control for changes in the path’s curvature, only past five errors
are considered [1]. So, the classical integral control has been modified as
shown in (5).

$u(n)=k_{I}\times\sum_{n=0}^{4}e(n)$ (5)

###  IV-D Derivative Control

Equation (6) is the mathematical expression for the classical derivative
control. Derivative control has been added to eliminate oscillation and
overshooting caused by proportional and integral terms.

$u(t)=k_{D}\times\frac{de(t)}{dt}$ (6)

D control prevents overshooting in the system and improves the system’s
transient response. However, it should be noted that the D control is
sensitive to noise [5]. The reason why D control is sensitive to noise is that
it acts like a high pass filter amplifying instant changes in error. The
effect of noise has not been seen much during the development of the robot.
So, the classical derivative has been modified as in (7) so the calculation is
done by finding the difference between the current and previous error.

$u(n)=k_{D}\times(e(n)-e(n-1))$ (7)

###  IV-E Modified PID Controller

Equation (8) expresses the classical PID control. With the aforementioned
modifications, the overall mathematical expression for the PID controller is
given by (9). In this controller, proportional term is used to minimize
positional error, integral term is used to minimize steady-state error
occurring in curved paths and derivative term is used to minimize oscillatory
behavior and overshooting.

$u(n)=k_{P}e(t)+k_{I}\int_{t_{0}}^{t}e(\tau)\,d\tau+k_{D}\times\frac{de(t)}{dt}$
(8)

$u(n)=k_{P}e(n)+k_{I}\sum_{n=0}^{4}e(n)+k_{D}(e(n)-e(n-1))$ (9)

Figure 4: Flow chart of the line-tracking controller.

###  IV-F Open-loop Controller

The PID controller shows good performance in straight lines as well as curved
paths. However, for handling sharp turns, it shows a poor performance, that is
the robot goes off the track. Especially, in 45 degree turns, the robot was
turning opposite to the path when only the PID controller was used. A thorough
analysis has been conducted to understand why the PID controller fails in
sharp turns. The analysis has revealed that in fact, the PID controller
initially responds to the change in the path correctly for a short time.
However, when all sensors miss the line, last errors come from the opposite
side of the path due to the initial response of the PID controller causing the
robot to miss the track. Therefore, an open-loop control mechanism has been
developed for handling sharp turns. In this control, when all the sensors miss
the line, only right and left end sensors, i.e. sensor 1 and sensor 8
respectively, are considered to understand which direction the robot should
turn. If the last error reading comes from sensor 1, the robot should turn
right, if it comes from sensor 8, the robot should turn left by rotating the
two wheels in opposite directions. The amount of speed at which each motor
will rotate is determined by another control mechanism developed to handle
both curved paths and corners as fast as possible. If the sensors miss the
line for at least 25 reading which is the case when handling corners, then a
sharp turn is performed, i.e. inner wheel rotates backwards at 53% of the full
speed and outer wheel rotates forward at 70% of the full speed. This enables
the robot to quickly come back to the line in 90 and 45 degree turns. If the
sensors miss the line for a few times, i.e. not more than 25, which is the
case when handling curved paths, then the inner motor rotates backward at 20%
of the full speed and outer wheel rotates forward at full speed. This causes
the robot to turn slightly and go faster where the path is curved. Even though
integral control shows a good performance in handling curvatures, this open-
loop control improves the performance when handling slight turns.

##  V Analysis

In this research, a heuristic approach [1] has been followed for system
identification. It has been found that this system has a pole at origin.
Furthermore, straight lines can be thought as step disturbance and curved
paths can be thought as ramp disturbance. Step and ramp response of an
arbitrary Type 1 transfer function has been simulated as shown in Fig. 5 and
analyzed to get a better understanding of the system under consideration. Fig.
6(a) shows the step response of an arbitrary Type 1 unity feedback system. As
expected, there is zero steady-state error for a step reference. However, Type
1 systems cannot track ramp references with zero steady-state error. Fig. 6(b)
shows the ramp response of the same Type 1 unity feedback system. Therefore,
integral control is added to the system to obtain zero-steady state error when
tracking ramp references. Fig. 6(c) shows the ramp response of the Type 1
system with integral control added.

Figure 5: A Type 1 unity feedback system.

(a) Step response of the Type 1 unity feedback system.

(b) Ramp response of the Type 1 unity feedback system.

(c) Ramp response when integral term is added.

Figure 6: Analyzing steady-state tracking error for different system
configurations.

##  VI Results

Fig. 7 shows the test track consisting of straight lines, curvatures, and
corners with 90 and 45 degree turns. The track is made of wooden plates having
a light color. A black electrical tape has been used to create the path so
that reflectance sensors can work efficiently.

Figure 7: A slightly modified version of the test track in [1].

Before implementing the open-loop controls, with 10.2 V battery voltage, the
robot could not handle 45 degree turns. After implementing the open-loop
control in which when all sensors miss the line, only leftmost and rightmost
sensors are considered, the robot completed the track in 9.2 seconds. To
further reduce the completion time, the battery voltage was increased up to 13
V. However, the robot started to fail to follow the line in sharp turns. For
handling sharp turns and improving system performance in curved paths, the
open-loop control has been modified. Finally, with 13 V battery voltage and
92% base PWM duty cycle, the robot successfully completed the track in 7.8
seconds.

##  VII Discussion

Considering the limitations of the DC motors and wheels, it has been
challenging to reduce the completion time below 7.5 seconds. It has been
observed that dust and dirt on the test track and wheels negatively affect
completion time and stability. When the track is dirty, it has been sometimes
observed that the robot slipped off the track and could not complete it.
Another limitation is that 6 V DC motors are slow. To address this issue, a
higher input voltage has been applied to the motor driver, i.e. 13 V. Yet
another difficulty was handling 45 degree turns. Initially, an integral
control based approach has been developed in which only a few past errors were
being considered after all sensors miss the line. However, this approach
suffered from the initial response of the PID controller to the 45 degree
turn. So, past errors were also coming from the opposite side of the path. To
eliminate this problem, a basic method has been developed in which only right
and left end sensors are considered to decide which direction to turn after
all sensors miss the line. Another limitation is posed by the distance between
each sensor on the QTR-8. The discrete structure of the sensor array causes
the robot to oscillate while following the line. If another module having
sensors placed more densely is used, it is expected that the amount of
oscillation is reduced since the sudden changes in error is mitigated. Yet
another limitation is that batteries are contributing to a great portion of
the robot’s weight. Also they cause the robot skid off the path in sharp
turns.

##  VIII Conclusion

In this paper, a PID controller scheme along with additional open-loop
controls has been developed for a line tracking robot consisting of low-cost,
off-the-shelf components such as STM32F103C8, QTR-8RC reflectance sensor,
L298N motor driver. During the construction of the robot, it has been ensured
that sensors are close to the ground, i.e. 3.5 mm, and the total weight is
approximately at the center of the chassis and close to the ground. Before
designing a controller, the system characteristics have been identified with a
simple unity feedback implementation. Then, a PID controller has been
implemented. However, to improve system stability and robustness, a number of
changes have been made. Since the PID controller was ineffective in handling
sharp turns, an open-loop control mechanism has been developed. To minimize
overshooting and oscillatory behavior, $k_{p}$, $k_{I}$ and $k_{D}$ gains have
been optimized through experiments. One of the main limitations affecting
stability of the system is the separation between each sensor in the QTR-8
module. This separation causes the error term to change suddenly consequently
the robot shows oscillatory behaviour.

## References

  * [1] V. Balaji, M. Balaji, M. Chandrasekaran, M. K. Khan, and I. Elamvazuthi, “Optimization of PID Control for High Speed Line Tracking Robots,” in _Procedia Computer Science_ , 2015.
  * [2] A. Latif, H. A. Widodo, R. Rahim, and K. Kunal, “Implementation of Line Follower Robot based Microcontroller ATMega32A,” _Journal of Robotics and Control (JRC)_ , 2020.
  * [3] K. M. Hasan, Abdullah-Al-Nahid, and A. Al Mamun, “Implementation of autonomous line follower robot,” in _2012 International Conference on Informatics, Electronics and Vision, ICIEV 2012_ , 2012.
  * [4] M. Pakdaman, M. M. Sanaatiyan, and M. R. Ghahroudi, “A line follower robot from design to implementation: Technical issues and problems,” in _2010 The 2nd International Conference on Computer and Automation Engineering, ICCAE 2010_ , 2010.
  * [5] G. Franklin, J. D. Powell, and A. Emami-Naeini, “Feedback control of dynamic systems, 3e,” in _American Society of Mechanical Engineers, Dynamic Systems and Control Division (Publication) DSC_ , 1994.


# Search for Light Exotic Fermions in Double-Beta Decays

Matteo Agostini<EMAIL_ADDRESS>Department of Physics and
Astronomy, University College London, Gower Street, London WC1E 6BT, UK.
Physik-Department, Technische Universität München, James-Franck-Straße, 85748
Garching, Germany  Elisabetta Bossio<EMAIL_ADDRESS>Alejandro
Ibarra<EMAIL_ADDRESS>Xabier Marcano<EMAIL_ADDRESS>Physik-Department,
Technische Universität München, James-Franck-Straße, 85748 Garching, Germany

###### Abstract

The Standard Model of Particle Physics predicts the double-$\beta$ decay of
certain nuclei with the emission of two active neutrinos. In this letter, we
argue that double-$\beta$ decay experiments could be used to probe models with
light exotic fermions through the search for spectral distortions in the
electron spectrum with respect to the Standard Model expectations. We consider
two concrete examples: models with light sterile neutrinos, singly produced in
the double-$\beta$ decay, and models with a light $Z_{2}$-odd fermion, pair
produced due to a $Z_{2}$ symmetry. We estimate the discovery potential of a
selection of double-$\beta$ decay experiment and find that future searches
will test for the first time a new part of the parameter space of interest at
the MeV-mass scale.

††preprint: TUM-HEP 1306/20

##  I Introduction

Many models of new physics contain new spin 1/2 particles, singlets under the
Standard Model gauge group, possibly related to the mechanism of neutrino mass
generation, or to the dark matter of the Universe. An archetype of such
fermions is the sterile neutrino, also called right-handed neutrino. In a
variant of this scenario, the singlet fermion is furnished with a $Z_{2}$
symmetry so it can only be produced in pairs. Unfortunately, in either of
these scenarios, the mass of the light exotic fermion and the coupling
strength to the Standard Model particles are free parameters of the model,
leaving a vast parameter space that must be probed in laboratory experiments,
or with astrophysical and cosmological observations.

From the phenomenological standpoint, sub-GeV sterile neutrinos are
particularly attractive, since they could be produced in charged lepton or in
hadron decays and possibly be discovered in a laboratory experiment. In fact,
there is an intense experimental program aiming to detect signals of sterile
neutrinos in the sub-GeV range, either in dedicated experiments, or as a by-
product of an experiment initially designed for a different purpose. The
current limits are fairly stringent below $\sim 100$ keV and above $\sim 100$
MeV (for an updated summary of the constraints, see e.g. Bolton  _et al._
(2020a)). However, the intermediate mass region remains relatively unexplored
and some of the current constraints date back to the 1990s Deutsch  _et al._
(1990); Schreckenbach  _et al._ (1983). The scenario with $Z_{2}$-odd singlet
fermions has been even less studied, and, to the best of our knowledge, there
is no laboratory constraint for the mass range $\sim 100$ keV $-100$ MeV.

In this letter, we explore the possibility of searching for the production of
light exotic fermions in double-$\beta$ decays. These decays are nuclear
transitions in which the atomic number $Z$ increases by two units while the
nucleon number $A$ remains constant. The process is allowed by the Standard
Model as long as two neutrinos and two electrons are emitted:
$(A,Z)\rightarrow(A,Z+2)+2e+2\bar{\nu}$ ($2\nu\beta\beta$ decay). If light
exotic fermions exist, they will also be produced in double-$\beta$ decays,
replacing one or both neutrinos in the final state. While neutrinos or exotic
fermions are challenging to detect, the energy distribution of the electrons
emitted in the process can be accurately measured. The shape of this
distribution carries information about the full kinematic of the process and
can be used to reconstruct which other particles have been emitted. The
fermion mass range accessible with this kind of search extends from a few
hundred keV to a few MeV, i.e., from the energy threshold of the experiments
till the Q-values of the decay.

In the next couple of decades, several experiments with the capability of
measuring the electron energy distribution in $2\nu\beta\beta$ decays will get
online Dolinski  _et al._ (2019). The primary physics goal of these
experiments is the discovery of double-$\beta$ decays with only the two
electrons in the final state: $(A,Z)\rightarrow(A,Z+2)+2e$. Observing this so-
called neutrinoless double-$\beta$ ($0\nu\beta\beta$) decay would have far-
reaching implications. It would prove that the lepton number is not a global
symmetry and that neutrinos have a Majorana mass component Furry (1939);
Schechter and Valle (1982). This physics goal already justifies major
investments for experimental infrastructure. Nevertheless, there are several
other valuable physics searches that can benefit from the low-background and
large target mass of these experiments Hirsch  _et al._ (1996); Díaz (2014);
Deppisch  _et al._ (2020a, b). Our proposed search for light exotic fermions
provides a valuable extension of their physics program.

We will discuss and estimate the discovery potential of double-$\beta$ decay
experiments under two illustrative scenarios that extend the Standard Model by
adding light exotic fermions. In the first one, we assume the existence of a
massive sterile neutrino $N$ that mixes with the electron neutrino. This opens
the possibility to observe a double-$\beta$ decay final state in which a
standard neutrino is replaced by a sterile neutrino ($\nu N\beta\beta$ decay):
$(A,Z)\rightarrow(A,Z+2)+2e+\nu+N$, with a modified spectrum due to the
sterile neutrino mass. This is indeed the same principle as the one used in
kink searches with single-$\beta$ decays Shrock (1980); Riis and Hannestad
(2011); Mertens  _et al._ (2015); Abada  _et al._ (2019); Aker  _et al._
(2020), which currently provide the strongest laboratory bounds in our mass
range of interest Holzschuh  _et al._ (2000); Deutsch  _et al._ (1990);
Schreckenbach  _et al._ (1983); Derbin  _et al._ (2018). As we will see, the
presence of light exotic fermions does not create a kink in the total
double-$\beta$ energy spectrum but still creates a continuous distortion that
is detectable by $0\nu\beta\beta$ decay experiments.

The second scenario that we will consider is characterized by the presence of
a new symmetry that forbids the production of a single exotic fermion in
double-$\beta$ decays. This is a typical scenario for models adding new exotic
fermions $\chi$ as dark matter candidates, which are charged under a $Z_{2}$
symmetry to make them stable. This kind of models cannot be tested through
single-$\beta$ decays, but would result in a new final state in double-$\beta$
decays ($\chi\chi\beta\beta$ decay). Thus, double-$\beta$ decay experiments
can provide the first laboratory-based constraints on several models.

In the following, we will review the phenomenology of the aforementioned
scenarios and present our calculations of how the electron energy distribution
is affected by the model parameters. The experimental identification of
distortions in the electron energy distribution requires a spectral fit whose
accuracy is limited by statistical and systematic uncertainties. We will
discuss how to set up a statistical analysis and evaluate the impact of
systematic uncertainties. Finally, we will derive sensitivity projections for
the most promising $0\nu\beta\beta$ experiments, assuming the systematic
uncertainties and experimental performance achieved by running experiments.

##  II  Double-$\boldsymbol{\beta}$ decay into sterile neutrinos

One of the simplest extensions of the Standard Model consists in adding to the
particle content one spin 1/2 particle, singlet under the gauge group, $N$.
The gauge symmetry allows a Yukawa coupling of $N$ to the Standard Model Higgs
doublet and lepton doublets, which leads upon electroweak symmetry breaking to
a Dirac mass term which couples $N$ to $\nu$, $m_{D}\overline{\nu}N$. For this
reason, $N$ is commonly denominated right-handed neutrino or sterile neutrino.
Furthermore, the gauge symmetry also allows a Majorana mass term for $N$,
$M\overline{N^{c}}N$ which we assume to be in the range 0.1-2 MeV. The two
parameters of the model, $m_{D}$ and $M_{N}$ are usually recast in terms of
the misalignment angle between the interaction and mass eigenstates,
$\sin\theta\simeq m_{D}/M$, and the mass of the heaviest eigenstate
$m_{N}\simeq M$.

If kinematically possible, sterile neutrinos could be produced in any decay
process involving Standard Model neutrinos, due to the active-sterile mixing
angle. The new decay channels lead to distortions in the energy distribution
for the visible particles compared to the Standard Model expectations, and
which then constitute a test for the sterile neutrino scenario. In this work
we concentrate on the double-$\beta$ decay. In the presence of sterile
neutrinos, and provided their mass is below the Q-value of the decay, the
decay channels $\nu N\beta\beta$ and $NN\beta\beta$ become kinematically
possible. The differential energy spectrum $\Gamma$ is given by the incoherent
superposition of three channels, and can be expressed as:

$\frac{d\Gamma}{dT}=\cos^{4}\theta\,\frac{d\Gamma_{\nu\nu}}{dT}\,\theta(T_{0}-T)+2\cos^{2}\theta\sin^{2}\theta\,\frac{d\Gamma_{\nu
N}}{dT}\,\theta(T_{0}-T-x_{N})+\sin^{4}\theta\,\frac{d\Gamma_{NN}}{dT}\,\theta(T_{0}-T-2x_{N})\,,$
(1)

where $T=(E_{e_{1}}+E_{e_{2}}-2m_{e})/m_{e}$ is the sum of the kinetic
energies of the two electrons (normalized to the electron mass). This variable
is kinematically restricted to be $0\leq T\leq T_{0}$, $T_{0}-x_{N}$ and
$T_{0}-2x_{N}$, for $2\nu\beta\beta$, $\nu N\beta\beta$ and $NN\beta\beta$
respectively, with $x_{N}=m_{N}/m_{e}$ and $T_{0}=Q_{\beta\beta}/m_{e}$.
$Q_{\beta\beta}$ is the end-point energy assuming massless neutrinos and
depends on the particular nucleus.

The $2\nu\beta\beta$ decay has been studied in several works Doi  _et al._
(1981); Doi and others. (1981); Doi  _et al._ (1985); Haxton and Stephenson
(1984); Tomoda (1991); Kotila and Iachello (2013); Stefanik  _et al._ (2015);
Šimkovic  _et al._ (2018a), always assuming vanishing neutrino masses. In this
work we extend this analysis, leaving the masses of the invisible fermions as
free parameters. To this end, we follow and generalize Kotila and Iachello
(2013), and express the differential rate for a generic decay
$(A,Z)\rightarrow(A,Z+2)+2e+a+b$, with $a$ and $b$ being either $\nu$ or $N$,
as:

$\frac{d\Gamma_{ab}}{dT}=|\mathcal{M}_{2\nu}^{\rm
eff}|^{2}\,\frac{dG_{ab}^{(0)}}{dT}\,,\quad$ (2)

with $M_{2\nu}^{\rm eff}$ is the dimensionless nuclear matrix element (NME).
The phase space factor is given by

$G_{ab}^{(0)}=\frac{1}{g_{A}^{4}m_{e}^{2}}\frac{2\widetilde{A}^{2}}{3\log
2}\int_{m_{e}}^{E^{\rm max}_{e_{1}}}dE_{e_{1}}\int_{m_{e}}^{E^{\rm
max}_{e_{2}}}dE_{e_{2}}\int_{m_{a}}^{E_{a}^{\rm max}}dE_{a}\Big{(}\langle
K_{\mathcal{N}}\rangle^{2}+\langle L_{\mathcal{N}}\rangle^{2}+\langle
K_{\mathcal{N}}\rangle\langle
L_{\mathcal{N}}\rangle\Big{)}\,f_{11}^{(0)}\,\omega_{ab}\,,$ (3)

where

$\displaystyle\widetilde{A}$
$\displaystyle=\frac{1}{2}(Q_{\beta\beta}+2m_{e})+\langle
E_{\mathcal{N}}\rangle-E_{I}\,,$ (4) $\displaystyle\langle
K_{\mathcal{N}}\rangle$ $\displaystyle=\frac{1}{E_{e_{1}}+E_{a}+\langle
E_{\mathcal{N}}\rangle-E_{I}}+\frac{1}{E_{e_{2}}+E_{b}+\langle
E_{\mathcal{N}}\rangle-E_{I}}\,,$ (5) $\displaystyle\langle
L_{\mathcal{N}}\rangle$ $\displaystyle=\frac{1}{E_{e_{1}}+E_{b}+\langle
E_{\mathcal{N}}\rangle-E_{I}}+\frac{1}{E_{e_{2}}+E_{a}+\langle
E_{\mathcal{N}}\rangle-E_{I}}\,,$ (6)

with $E_{I}$ the energy of the initial nucleus, $\langle
E_{\mathcal{N}}\rangle$ a suitable excitation energy in the intermediate
nucleus and

$\omega_{ab}=\frac{g_{A}^{4}G_{\beta}^{4}}{64\pi^{7}}\,p_{a}E_{a}p_{b}E_{b}p_{e_{1}}E_{e_{1}}p_{e_{2}}E_{e_{2}}\,.$
(7)

Here, $E_{X}$ and $p_{X}=\sqrt{E_{X}^{2}-m_{X}^{2}}$ denote the energy and the
modulus of the 3-momentum of the particle $X=e_{1},e_{2},a,b$ with mass
$m_{X}$, and subject to the condition
$E_{b}=Q_{\beta\beta}+2m_{e}-E_{e_{1}}-E_{e_{2}}-E_{a}$ from energy
conservation. The integration limits read $E_{e_{1}}^{\rm
max}=Q_{\beta\beta}+m_{e}-m_{a}-m_{b}$, $E_{e_{2}}^{\rm
max}=Q_{\beta\beta}+2m_{e}-E_{e_{1}}-m_{a}-m_{b}$ and $E_{a}^{\rm
max}=Q_{\beta\beta}+2m_{e}-E_{e_{1}}-E_{e_{2}}-m_{b}$, which for
$m_{a},m_{b}\neq 0$ shifts the end-point of the spectrum to lower values.
Finally, the factor $f_{11}^{(0)}$ originates from the Coulomb interaction of
the electrons with the daughter nucleus, which we parameterize using the Fermi
function Doi  _et al._ (1985)

$f_{11}^{(0)}\simeq F_{0}(Z+2,E_{e_{1}})\,F_{0}(Z+2,E_{e_{2}})\,.$ (8)

In the limit $a,b=\nu$ and $m_{\nu}=0$ we recover the results of Kotila and
Iachello (2013).

The complete analysis of the spectrum requires a numerical evaluation of Eq.
(2). However, we also identified accurate analytical approximations for the
dominant $2\nu\beta\beta$ and $\nu N\beta\beta$ channels ($NN\beta\beta$ is
negligible for small $\sin\theta$), which assumes the Primakoff-Rosen
approximation Primakoff and Rosen (1959) for the Fermi function and neglects
the lepton energies in

$\langle K_{\mathcal{N}}\rangle\simeq\ \langle
L_{\mathcal{N}}\rangle\simeq\frac{4}{2\widetilde{A}-Q_{\beta\beta}-2m_{e}}\approx\frac{2}{\widetilde{A}}\,.$
(9)

Then, neglecting also the active neutrino mass but keeping the sterile
neutrino one, we obtain:

$\frac{dG_{ab}}{dT}\approx\frac{G_{\beta}^{4}m_{e}^{8}}{7200\pi^{7}\log
2}\bigg{[}\frac{2\pi\alpha(Z+2)}{1-\exp\big{(}-2\pi\alpha(Z+2)\big{)}}\bigg{]}^{2}\,F_{ab}(T)\,,$
(10)

where we have introduced the form factors

$\displaystyle F_{\nu\nu}(T)=$
$\displaystyle~{}T~{}(T^{4}+10T^{3}+40T^{2}+60T+30)~{}(T_{0}-T)^{5}\,,$ (11)
$\displaystyle F_{\nu N}(T)=$
$\displaystyle\,\frac{T}{2}(T^{4}+10T^{3}+40T^{2}+60T+30)\,\bigg{\\{}\Big{(}2(T_{0}-T)^{4}-9x_{N}^{2}(T_{0}-T)^{2}-8x_{N}^{4}\Big{)}\sqrt{(T_{0}-T)^{2}-x_{N}^{2}}$
$\displaystyle+15x_{N}^{4}(T_{0}-T)\Big{[}\log\left(T_{0}-T+\sqrt{(T_{0}-T)^{2}-x_{N}^{2}}\right)-\log
x_{N}\Big{]}\bigg{\\}}\,.$ (12)

Figure 1: Left: Phase space for different sterile neutrino mass hypotheses
computed using a numerical evaluation of Eq. (2) and the analytical
approximation of Eq. (10). Right: $2\nu\beta\beta$, $\nu N\beta\beta$ and
summed energy spectra for $m_{N}=500$ keV and an illustrative large mixing
$\sin^{2}\theta=0.2$.

The left panel of Fig. 1 displays the phase space resulting from our numerical
calculation applied to 76Ge ($Q_{\beta\beta}\simeq 2039$ keV), showing how the
presence of a final massive invisible fermion modifies the energy spectrum,
shifting the end-point as well as the peak to lower values as $m_{N}$
increases. We can also see that the analytical expression reproduces the full
result to a good approximation, up to a global factor of 1.67. This factor is
due to the approximation we considered, but it does not affect the analysis,
which is sensitive to the $G_{\nu N}/G_{\nu\nu}$ ratio. We stress that,
contrary to the single-$\beta$ decay spectrum, the end-point of the $\nu
N\beta\beta$ spectrum is smooth and therefore the total spectrum does not
manifest a kink at $Q_{\beta\beta}-m_{N}$. Nevertheless, the spectrum differs
from the standard one, allowing $0\nu\beta\beta$ experiments to probe this
hypothesis by measuring the $2\nu\beta\beta$ spectrum very accurately.

##  III Double-$\boldsymbol{\beta}$ decay into $\boldsymbol{Z_{2}}$-odd
singlet fermions

In general, the production of a pair of exotic particles is strongly
suppressed compared to the production of a single one and can be neglected.
However, there are scenarios in which the single production is forbidden and
only the double production can take place. These scenarios cannot be tested by
single-$\beta$ decay experiments, whereas $0\nu\beta\beta$ decay experiments
are still sensitive to the double production and have the unique opportunity
to explore these channels.

In this section we consider a variant of the previous model, in which the
symmetry group is extended by a discrete $Z_{2}$ symmetry, possibly related to
the dark matter sector, which is exact or mildly broken in the electroweak
vacuum. We assume that all Standard Model particles are even under the $Z_{2}$
symmetry, while the neutral singlet fermion is odd. We will denote the
$Z_{2}$-odd neutral singlet fermion as $\chi$ to differentiate it from the
sterile neutrino, since the $Z_{2}$ symmetry forbids the mass term
$m_{D}\overline{\nu}N$ characteristic of the latter. Correspondingly, the
$Z_{2}$ symmetry forbids the single-$\beta$ decay
$(A,Z)\rightarrow(A,Z+1)+e+\chi$, as well as the double-$\beta$ decay
$(A,Z)\rightarrow(A,Z+2)+2e+\nu+\chi$. However, the double-$\beta$ decay
$(A,Z)\rightarrow(A,Z+2)+2e+2\chi$ is possible ($\chi\chi\beta\beta$ decay).
The search for spectral distortions in the double-$\beta$ decay spectrum can
therefore probe scenarios with light $Z_{2}$-odd singlet fermions.

The total differential decay rate receives in this scenario two contributions:

$\displaystyle\frac{d\Gamma}{dT}=$
$\displaystyle~{}\frac{d\Gamma_{\nu\nu}}{dT}\,\theta(T_{0}-T)+\frac{d\Gamma_{\chi\chi}}{dT}\,\theta(T_{0}-T-2x_{\chi})\,,$
(13)

where $d\Gamma_{\nu\nu}/dT$ is the Standard Model contribution, defined in Eq.
(2), and $d\Gamma_{\chi\chi}/dT$ is the exotic contribution. Here, $T$ and
$T_{0}$ are defined as in Eq. (1), with $x_{\chi}=m_{\chi}/m_{e}$. The exotic
contribution is very model dependent. For definiteness, we will consider the
following effective interaction between the active neutrinos and $\chi$:

$\mathcal{L}=g_{\chi}\nu\nu\chi\chi\,,$ (14)

with a constant coupling $g_{\chi}$. We use the results of Deppisch  _et al._
(2020b), that considered a similar four-fermion scalar interaction as in Eq.
(14) but for neutrino self-interactions, to relate the decay rate for
$\chi\chi\beta\beta$ to that of $NN\beta\beta$. We obtain

$\frac{d\Gamma_{\chi\chi}}{dT}=\frac{g_{\chi}^{2}m_{e}^{2}}{8\pi^{2}R^{2}}|\mathcal{M}_{0\nu}|^{2}\frac{dG_{NN}^{(0)}}{dT}\,,$
(15)

where $\mathcal{M}_{0\nu}$ is the NME of the $0\nu\beta\beta$ process, $R$ is
the nuclear radius and the phase factor $G_{NN}^{(0)}$ is given in Eq. (3)
(replacing $m_{N}$ by $m_{\chi}$). The translation of experimental
measurements into constraints on $g_{\chi}$ will hence need to rely on NME
calculations.

If the $Z_{2}$ stabilizing symmetry is exact, $\chi$ could play a role in
cosmology. For example, if it thermalizes with the primeval plasma, it may
significantly alter the successful predictions of the standard Big Bang
nucleosynthesis scenario Sabti  _et al._ (2020). These strong conclusions
could be modified, for instance, if the $Z_{2}$ symmetry is only approximate,
so that $\chi$ is stable within the detectors, but is not cosmologically long-
lived. In this paper we will adopt a phenomenological approach, and we will
focus on the search of exotic light fermions produced in pairs in
double-$\beta$ decays, without addressing their possible implications on
cosmology that depend on additional parameters (such as the lifetime).

##  IV Discovery potential of future
$\boldsymbol{0\nu}\boldsymbol{\beta\beta}$ decay experiments

A large experimental program has been mounted to search for the
$0\nu\beta\beta$ decay of different isotopes and using different detection
techniques Dolinski  _et al._ (2019). The current-generation experiments set
the ground for the identification of the best techniques for further
investment. The next-generation experiments are currently being proposed,
designed, or in construction. Among the experiments proposed for the next
decade, we will focus on LEGEND Abgrall  _et al._ (2017), nEXO Kharusi  _et
al._ (2018) and CUPID Armstrong  _et al._ (2019). We selected these
experiments because their efficiencies and uncertainties in the search for
massive fermions can be inferred from the $2\nu\beta\beta$ analysis published
by their predecessors, i.e., GERDA Agostini  _et al._ (2015), EXO-200 Albert
_et al._ (2014) and CUPID-Mo Armengaud  _et al._ (2020). The parameters
assumed for each experiment are listed in Tab. 1.

Isotope | Experiment |  $T_{1/2}^{2\nu\beta\beta}$ [yr] | Efficiency | Exposure [mol yr]
---|---|---|---|---
76Ge | GERDA / LEGEND-200 / LEGEND-1000 |  $2.0\cdot 10^{21}$ Agostini  _et al._ (2020a) | 79% Agostini  _et al._ (2015) |  $1.4\cdot{10^{3}}$ Agostini  _et al._ (2020b) / $1.4\cdot{10^{4}}$ Abgrall  _et al._ (2017) / $1.4\cdot{10^{5}}$ Abgrall  _et al._ (2017)
136Xe | EXO-200 / nEXO |  $2.2\cdot 10^{21}$ Albert  _et al._ (2014) | 85% Albert  _et al._ (2014) |  $1.7\cdot{10^{3}}$ Anton  _et al._ (2019) / $3.7\cdot{10^{5}}$ Kharusi  _et al._ (2018)
100Mo | CUPID-Mo / CUPID |  $7.1\cdot 10^{18}$ Armengaud  _et al._ (2020) | 91% Armengaud  _et al._ (2020) |  $6.5\cdot{10^{-1}}$ Armengaud  _et al._ (2020) / $2.7\cdot{10^{5}}$ Armstrong  _et al._ (2019)
Table 1: Parameters used for the sensitivity projections of each experiment.
The detection efficiency refers to fraction of $2\nu\beta\beta$ decay events
that populate the energy window of interest after all analysis cuts. The
exposure is given in number of moles of detector material per year of live
time.

While the target isotope and the backgrounds vary across these experiments,
the analysis to search for light exotic fermions is always conceptually the
same. The energy window of interest goes from the detector threshold to the
Q-value of the decay. In this window, the majority of the events is due to
$2\nu\beta\beta$ decays ($>95$% in the current-generation experiments). The
other events can be due to a multitude of processes, for instance, natural
radioactivity and cosmic rays. The most important parameters affecting the
sensitivity of an experiment to light exotic fermions are the exposure
$\mathcal{E}$, the background rate $R_{\text{B}}$ and the systematic
uncertainties due to the energy reconstruction $\sigma_{\text{sys}}$. The
exposure is given by the product of the number of observed nuclei and the
observation time. The background rate is primarily given by the
$2\nu\beta\beta$ decay rate – i.e. the inverse of the half-life – with a
subdominant contribution due to the other sources, i.e.
$R_{\text{B}}=R_{2\nu\beta\beta}+R_{\rm other}$. The systematic uncertainties
related to the event energy reconstruction can largely differ between
experiments, but in general, their impact can be parameterized through an
energy-dependent shape factor.

In the search for new physics, the upper limit on the rate of a new signal
$R_{\text{U.L.}}$ scales approximately as the sum in quadrature of the
statistical and systematic uncertainties. When the statistical uncertainty is
dominant, the upper limit on the number of signal counts will be proportional
to the fluctuations on the number of events in the analysis window, which is
dominated by the background contributions. By comparing the upper limit on the
number of signal counts ($R_{\text{U.L.}}\cdot\mathcal{E}$) with the Gaussian
fluctuations of the number of background counts
($\sqrt{R_{\text{B}}\cdot\mathcal{E}}$), we obtain that:
$R_{\text{U.L.}}\propto\sqrt{R_{\text{B}}/\mathcal{E}}$. When the statistical
uncertainty becomes comparable with the systematic one, the sensitivity starts
to saturate and, eventually, further increasing the exposure does not improve
the sensitivity.

To accurately quantify the sensitivity of the experiments, we have implemented
a comprehensive frequentist analysis framework. Distortions of the
double-$\beta$ decay energy distribution due to light exotic fermions are
searched through a binned maximum-likelihood fit based on a profile-likelihood
test statistic Cowan  _et al._ (2011). Each process possibly contributing to
the count rate in the energy window of interest should be added to the fit
through a probability distribution function (PDF) that describes its expected
event energy distribution. For this work, we use a PDF for the sought-after
signal and one for the dominant $2\nu\beta\beta$ decay. Both these PDFs are
based on the calculations described in the previous sections. We use a third
uniform probability distribution to account for other generic sub-dominant
background sources. The actual shape of this third PDF affects only marginally
our results. The parameters of the fit are typically the scaling factors of
each PDF, i.e. the number of events attributed to each process111If the
parameter of interest for the analysis is the ratio between two scaling
factors, the two scaling factors in the fit bases are replaced by a common
scaling factor and a parameter defining their ratio. This is the case for the
search for sterile neutrinos in which the parameter of interest is
$\sin^{2}\theta$. . The probability distributions of the test statistic are
computed from large ensembles of pseudo-data generated under different
hypotheses on the signal rate. This approach provides the right coverage by
construction, including when the signal rate is close to the physical border.

Systematic uncertainties can bias the event energy reconstruction. Many
detector-specific sources of bias should be considered. However, their overall
impact can be parameterized through a shape factor with the form
$f(E)=1+a\cdot E+b\cdot E^{2}+c/E$ where $a$, $b$ and $c$ are parameters that
are considered to be known with limited accuracy (i.e. $\sigma_{a}$,
$\sigma_{b}$ and $\sigma_{c}$). Current-generation experiments are typically
able to control the energy reconstruction bias at the per-cent level. To
incorporate this systematic uncertainty in our analysis, we randomize $a$, $b$
and $c$ during the generation of the pseudo-data and run the fit using not-
deformed PDFs Cousins and Highland (1992). The parameters are sampled from
normal distributions centred at 0 and with a sigma of
$10^{-3}\,\text{keV}^{-1}$, $10^{-6}\,\text{keV}^{-2}$ and
$10^{-3}\,\text{keV}$, respectively. We tested that this parameterization
covers the maximal distortions estimated by the experiments. The result of
this procedure is a broadening of the test statistic distribution and a
reduction of the power of the statistical analysis.

Figure 2 shows our projected sensitivity for a generic 76Ge experiment and a
sterile neutrino with a mass of 500 keV. The sensitivity for different
background rates and exposure values is given in terms of the median 90%
confidence level (C.L.) upper limit on $\sin^{2}\theta$ Cowan  _et al._
(2011). As mentioned above, the upper limit is proportional to the sum in
quadrature of the statistical and systematic uncertainties. The computed
sensitivities follow accurately the expected analytical expression:

$R_{\text{U.L.}}\propto\sqrt{\sigma_{\text{stat}}^{2}+\sigma_{\text{sys}}^{2}}~{}\sim\sqrt{k^{2}_{1}R_{\text{B}}/\mathcal{E}+k^{2}_{2}}\,.$
(16)

The parameters $k_{1}$ and $k_{2}$ cannot be computed without a full analysis.
They depend on the isotope specific parameters (e.g. molar mass, phase space
factor, Q-value and NME), the assumed mass of the sterile neutrino and the
systematic uncertainties. They also capture the capability of the fit to
differentiate between two spectral shapes.

Figure 2: Sensitivity of a 76Ge double-$\beta$ decay experiment searching for
massive sterile neutrinos with a mass of 500 keV, given as a function of
exposure and background level. The sensitivity is expressed in terms of the
median 90% C.L. upper limit on the squared of the mixing angle
($\sin^{2}\theta$), computed assuming no signal. The markers indicate the
value computed using a full frequentist analysis, while the lines show their
best fit with the analytical expression discussed in the text. The plot shows
the impact of additional backgrounds with a rate $R_{\rm other}$ which is
expressed w.r.t. the irreducible $2\nu\beta\beta$ decay rate
$R_{2\nu\beta\beta}$.

The vertical lines in the figure mark the exposure that has been collected by
GERDA, as well as the exposure planned for the two phases of LEGEND. LEGEND
has the potential to improve the sensitivity by an order of magnitude compared
to GERDA, assuming that the systematic uncertainties can be kept below the
statistical uncertainty. A significant improvement is expected even in the
conservative scenario in which the systematic uncertainties are not reduced
with respect to the current level.

After tuning the parameters $k_{1}$ and $k_{2}$, Eq. (16) describes accurately
the evolution of the sensitivity for all the considered experiments. It is
also valid for the search of the $Z_{2}$-odd fermions. In general, current-
generation experiments are dominated by the statistical uncertainty while
future experiments will need to fight against the systematic uncertainties. We
will report projections based on the conservative scenario in which the
systematic uncertainties will not be improved as well as the optimistic
scenario in which they will be sub-dominant. These two projections define the
ballpark for the sensitivity of future experiment and the results for
intermediate scenarios can be interpolated from these two cases.

Our sensitivity projections for the search of sterile neutrinos and of the
$Z_{2}$-odd fermions are shown in Fig. 3. In both cases, the sensitivity
evolution as a function of the fermion mass has a parabolic shape. Its minimum
depends on the Q-value of the decay and corresponds to the fermion mass for
which the experiment is most sensitive. The experiments quickly loose
sensitivity towards vanishing masses. This is because the smaller is the
fermion mass, the smaller is the spectral distortion. A similar loss of
sensitivity occurs at larger masses, where the fermion mass approaches the
maximum energy available in the decay and the phase space shrinks222Fermions
with masses up to the Q-value of the decay could be produced, although with
very suppressed phase space. When the simultaneous production of two fermions
is considered, the maximum mass is equal to a half of the Q-value.. The
exposure, isotope properties, and experimental parameters define the offset of
the curves. These projections are done assuming that each experiment can
extend its analysis window down to arbitrarily small energies. If this is not
the case, the upper edge of the curve will be lowered by the value of the
energy threshold and the offset will slightly increase because of the
reduction in detection efficiency. The sensitivity on the coupling constant
$g_{\chi}$ is computed from the decay rate of $\chi\chi\beta\beta$ using Eq.
(15). In addition to the systematic uncertainty considered for the sterile
neutrino search, our constraints on $g_{\chi}$ include also the uncertainties
due to the NME calculations.

Figure 3: Sensitivity to light exotic fermions for a selection of
double-$\beta$ decay experiments. The left panel shows the median 90% C.L.
upper limit on the squared of the sterile neutrino mixing angle as a function
of its mass. The experimental constraints are displayed through a band
covering a range of sensitivities that goes from the most optimistic scenario
(systematic uncertainty smaller than statistical uncertainty), till a
conservative scenario (systematic uncertainty at the level of past analyses).
Existing sterile neutrino constraints from single-$\beta$ decay experiments
Holzschuh  _et al._ (2000); Deutsch  _et al._ (1990); Schreckenbach  _et al._
(1983); Derbin  _et al._ (2018) and solar neutrinos Bellini  _et al._ (2013)
are shown in the background. The right panel shows the median upper limit on
the coupling constant between $Z_{2}$-odd fermion and the neutrinos assuming
the effective interaction in Eq. (14). The spread of the bands account for the
systematic uncertainties as in the case of the sterile neutrinos, but
additionally it covers also the full range of possible NME values found in
literature Barea  _et al._ (2015); Menéndez (2018); Horoi and Neacsu (2016);
López Vaquero  _et al._ (2013); Rodriguez and Martinez-Pinedo (2010); Song
_et al._ (2017); Hyvärinen and Suhonen (2015); Mustonen and Engel (2013); Fang
_et al._ (2018); Šimkovic  _et al._ (2018b); Coraggio  _et al._ (2020).

##  V Results and Outlook

The main results of our analysis are displayed in Fig. 3. The left panel shows
the projected upper limits on the active-sterile mixing, parametrized as
$\sin^{2}\theta$, as a function of the sterile neutrino mass. We also show in
the plot the current limits from single-$\beta$ decay experiments Holzschuh
_et al._ (2000); Deutsch  _et al._ (1990); Schreckenbach  _et al._ (1983);
Derbin  _et al._ (2018) and from solar neutrinos Bellini  _et al._ (2013). As
can be seen from the plot, the sensitivity of current double-$\beta$ decay
experiments is weaker than the existing limits, but only by a factor of a few.
The larger exposure of future double-$\beta$ decay experiments encourages a
dedicated search for these exotic decay channels.

Indeed, our sensitivity study demonstrates the potential of future
double-$\beta$ decay experiments to improve the current limits on the active-
sterile mixing angle in this mass range. Fig. 3 shows our projected
sensitivities for LEGEND-1000, nEXO and CUPID, assuming no improvement in the
systematic uncertainties with respect to the current-generation experiments
(upper side of the bands), or assuming that the systematic uncertainties will
be reduced below the statistical ones (lower side of the bands). In the most
conservative scenario, the sensitivity of future searches will be comparable
to the current limits. However, even with a modest reduction of the systematic
uncertainties, next-generation experiments will explore uncharted regions of
the sterile neutrino parameter space, even reaching $\sin^{2}\theta\sim
10^{-3}-10^{-4}$ for $m_{N}\sim(100-1600)$ keV.

Double-$\beta$ decay experiments have also the capability of probing models in
which only the double production of light exotic fermions is allowed. This
previously overlooked opportunity can lead to the first constraints on this
kind of models. This is shown in the right panel of Fig. 3 for the $Z_{2}$-odd
singlet $\chi$ introduced in Sec. III. The sensitivities of current
experiments for the effective coupling $g_{\chi}$ lie between $10^{-2}$ and
$10^{-3}~{}{\rm MeV}^{-2}$, which could be improved up to $10^{-4}~{}{\rm
MeV}^{-2}$ with a favorable NME value and negligible systematic uncertainties.

Assuming that this effective interaction is originated at some scale $\Lambda$
with $\mathcal{O}(1)$ dimensionless Wilson coefficient, this would imply that
$\Lambda\gtrsim 100$ MeV. This is comparable to the typical momentum transfer
in double-$\beta$ decay experiments, and therefore the effective field theory
approach employed in this work to recast the limits on the decay rate into
limits on the coupling strength should be taken with a grain of salt. Yet, the
search for distortions in the double-$\beta$ decay spectrum due to the
emission of two light exotic particles (fermions or scalars) is well motivated
theoretically and deserves further investigation.

To conclude, in this letter we have proposed and explored the possibility of
searching for the production of light exotic fermions using double-$\beta$
decay experiments. We have derived the double-$\beta$ decay spectrum when one
or two light exotic fermions are emitted in the final state. We have also
studied how these new channels can be searched for through a fit, and we have
estimated the projected sensitivity for various future experiments considering
the impact of both statistical and systematic uncertainties. We have found
that double-$\beta$ decay experiments constitute an promising avenue (and for
some scenarios unique) to search for new physics. We therefore encourage the
experimental collaborations to perform dedicated searches for this class of
scenarios, properly including all the systematic and statistical
uncertainties.

Note added. During the final stages of this project, we received a preprint by
Bolton et al. Bolton  _et al._ (2020b), also exploring the possibility of
searching for sterile neutrinos with double-$\beta$ decay experiments. While
they focus exclusively on sterile neutrinos (left and right-handed models) we
considered generic single and double production of light exotic fermions.
Their procedure to derive the sensitivity projections and their treatment of
the systematic uncertainties differ from ours. Nevertheless, in the aspects
where our analysis overlap, the results are qualitatively similar.

###### Acknowledgements.

X.M. would like to thank Claudia García-García for useful discussions. This
work has been supported by the Collaborative Research Center SFB1258, by the
Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under
Germany’s Excellence Strategy - EXC-2094 - 390783311, and by the Science and
Technology Facilities Council [grant number ST/T004169/1]. X.M. is supported
by the Alexander von Humboldt Foundation.

## References

  * Bolton  _et al._ (2020a) P. D. Bolton, F. F. Deppisch, and P. Bhupal Dev, JHEP 03, 170 (2020a), arXiv:1912.03058 [hep-ph] .
  * Deutsch  _et al._ (1990) J. Deutsch, M. Lebrun, and R. Prieels, Nucl. Phys. A 518, 149 (1990).
  * Schreckenbach  _et al._ (1983) K. Schreckenbach, G. Colvin, and F. Von Feilitzsch, Phys. Lett. B 129, 265 (1983).
  * Dolinski  _et al._ (2019) M. J. Dolinski, A. W. Poon, and W. Rodejohann, Ann. Rev. Nucl. Part. Sci. 69, 219 (2019), arXiv:1902.04097 [nucl-ex] .
  * Furry (1939) W. Furry, Phys. Rev. 56, 1184 (1939).
  * Schechter and Valle (1982) J. Schechter and J. Valle, Phys. Rev. D 25, 2951 (1982).
  * Hirsch  _et al._ (1996) M. Hirsch, H. Klapdor-Kleingrothaus, S. Kovalenko, and H. Pas, Phys. Lett. B 372, 8 (1996), arXiv:hep-ph/9511227 .
  * Díaz (2014) J. S. Díaz, Physical Review D 89 (2014), 10.1103/physrevd.89.036002.
  * Deppisch  _et al._ (2020a) F. F. Deppisch, L. Graf, and F. Šimkovic, Phys. Rev. Lett. 125, 171801 (2020a), arXiv:2003.11836 [hep-ph] .
  * Deppisch  _et al._ (2020b) F. F. Deppisch, L. Graf, W. Rodejohann, and X.-J. Xu, Phys. Rev. D 102, 051701 (2020b), arXiv:2004.11919 [hep-ph] .
  * Shrock (1980) R. Shrock, Phys. Lett. B 96, 159 (1980).
  * Riis and Hannestad (2011) A. S. Riis and S. Hannestad, JCAP 02, 011 (2011), arXiv:1008.1495 [astro-ph.CO] .
  * Mertens  _et al._ (2015) S. Mertens _et al._ , JCAP 02, 020 (2015), arXiv:1409.0920 [physics.ins-det] .
  * Abada  _et al._ (2019) A. Abada, A. Hernández-Cabezudo, and X. Marcano, JHEP 01, 041 (2019), arXiv:1807.01331 [hep-ph] .
  * Aker  _et al._ (2020) M. Aker _et al._ (KATRIN), (2020), arXiv:2011.05087 [hep-ex] .
  * Holzschuh  _et al._ (2000) E. Holzschuh, L. Palermo, H. Stussi, and P. Wenk, Phys. Lett. B 482, 1 (2000).
  * Derbin  _et al._ (2018) A. Derbin _et al._ , JETP Lett. 108, 499 (2018).
  * Doi  _et al._ (1981) M. Doi _et al._ , Prog. Theor. Phys. 66, 1739 (1981), [Erratum: Prog.Theor.Phys. 68, 347 (1982)].
  * Doi and others. (1981) M. Doi and others., Prog. Theor. Phys. 66, 1765 (1981), [Erratum: Prog.Theor.Phys. 68, 348 (1982)].
  * Doi  _et al._ (1985) M. Doi, T. Kotani, and E. Takasugi, Prog. Theor. Phys. Suppl. 83, 1 (1985).
  * Haxton and Stephenson (1984) W. Haxton and G. Stephenson, Prog. Part. Nucl. Phys. 12, 409 (1984).
  * Tomoda (1991) T. Tomoda, Rept. Prog. Phys. 54, 53 (1991).
  * Kotila and Iachello (2013) J. Kotila and F. Iachello, Phys. Rev. C 87, 024313 (2013), arXiv:1303.4124 [nucl-th] .
  * Stefanik  _et al._ (2015) D. Stefanik, F. Simkovic, and A. Faessler, Phys. Rev. C 91, 064311 (2015), arXiv:1506.00835 [nucl-th] .
  * Šimkovic  _et al._ (2018a) F. Šimkovic, R. Dvornický, D. Stefánik, and A. Faessler, Phys. Rev. C 97, 034315 (2018a), arXiv:1804.04227 [nucl-th] .
  * Primakoff and Rosen (1959) H. Primakoff and S. Rosen, Rept. Prog. Phys. 22, 121 (1959).
  * Sabti  _et al._ (2020) N. Sabti, J. Alvey, M. Escudero, M. Fairbairn, and D. Blas, JCAP 01, 004 (2020), arXiv:1910.01649 [hep-ph] .
  * Abgrall  _et al._ (2017) N. Abgrall _et al._ (LEGEND), AIP Conf. Proc. 1894, 020027 (2017), arXiv:1709.01980 [physics.ins-det] .
  * Kharusi  _et al._ (2018) S. A. Kharusi _et al._ (nEXO), (2018), arXiv:1805.11142 [physics.ins-det] .
  * Armstrong  _et al._ (2019) W. Armstrong _et al._ (CUPID), (2019), arXiv:1907.09376 [physics.ins-det] .
  * Agostini  _et al._ (2015) M. Agostini _et al._ , Eur. Phys. J. C 75, 416 (2015), arXiv:1501.02345 [nucl-ex] .
  * Albert  _et al._ (2014) J. Albert _et al._ (EXO-200), Phys. Rev. C 89, 015502 (2014), arXiv:1306.6106 [nucl-ex] .
  * Armengaud  _et al._ (2020) E. Armengaud _et al._ , Eur. Phys. J. C 80, 674 (2020), arXiv:1912.07272 [nucl-ex] .
  * Agostini  _et al._ (2020a) M. Agostini _et al._ (GERDA), JHEP 03, 139 (2020a), arXiv:1909.02522 [nucl-ex] .
  * Agostini  _et al._ (2020b) M. Agostini _et al._ (GERDA), (2020b), arXiv:2009.06079 [nucl-ex] .
  * Anton  _et al._ (2019) G. Anton _et al._ (EXO-200), Phys. Rev. Lett. 123, 161802 (2019), arXiv:1906.02723 [hep-ex] .
  * Cowan  _et al._ (2011) G. Cowan, K. Cranmer, E. Gross, and O. Vitells, Eur. Phys. J. C 71, 1554 (2011), [Erratum: Eur.Phys.J.C 73, 2501 (2013)], arXiv:1007.1727 [physics.data-an] .
  * Cousins and Highland (1992) R. D. Cousins and V. L. Highland, Nucl. Instrum. Meth. A 320, 331 (1992).
  * Bellini  _et al._ (2013) G. Bellini _et al._ (Borexino), Phys. Rev. D 88, 072010 (2013), arXiv:1311.5347 [hep-ex] .
  * Barea  _et al._ (2015) J. Barea, J. Kotila, and F. Iachello, Phys. Rev. C 91, 034304 (2015), arXiv:1506.08530 [nucl-th] .
  * Menéndez (2018) J. Menéndez, J. Phys. G 45, 014003 (2018), arXiv:1804.02105 [nucl-th] .
  * Horoi and Neacsu (2016) M. Horoi and A. Neacsu, Phys. Rev. C 93, 024308 (2016), arXiv:1511.03711 [nucl-th] .
  * López Vaquero  _et al._ (2013) N. López Vaquero, T. R. Rodríguez, and J. L. Egido, Phys. Rev. Lett. 111, 142501 (2013), arXiv:1401.0650 [nucl-th] .
  * Rodriguez and Martinez-Pinedo (2010) T. R. Rodriguez and G. Martinez-Pinedo, Phys. Rev. Lett. 105, 252503 (2010), arXiv:1008.5260 [nucl-th] .
  * Song  _et al._ (2017) L. Song, J. Yao, P. Ring, and J. Meng, Phys. Rev. C 95, 024305 (2017), arXiv:1702.02448 [nucl-th] .
  * Hyvärinen and Suhonen (2015) J. Hyvärinen and J. Suhonen, Phys. Rev. C 91, 024613 (2015).
  * Mustonen and Engel (2013) M. Mustonen and J. Engel, Phys. Rev. C 87, 064302 (2013), arXiv:1301.6997 [nucl-th] .
  * Fang  _et al._ (2018) D.-L. Fang, A. Faessler, and F. Simkovic, Phys. Rev. C 97, 045503 (2018), arXiv:1803.09195 [nucl-th] .
  * Šimkovic  _et al._ (2018b) F. Šimkovic, A. Smetana, and P. Vogel, Phys. Rev. C 98, 064325 (2018b), arXiv:1808.05016 [nucl-th] .
  * Coraggio  _et al._ (2020) L. Coraggio _et al._ , Phys. Rev. C 101, 044315 (2020), arXiv:2001.00890 [nucl-th] .
  * Bolton  _et al._ (2020b) P. D. Bolton, F. F. Deppisch, L. Gráf, and F. Šimkovic, (2020b), arXiv:2011.13387 [hep-ph] .


# A note on the formulation of the Ensemble Adjustment Kalman Filter

Ian Grooms

(Department of Applied Mathematics, University of Colorado, Boulder, Colorado,
USA 80309)

###### Abstract

The ensemble adjustment Kalman filter (EAKF; Anderson, 2001) is one of the
earliest ensemble square root filters. This note clarifies the correct
formulation of the EAKF, which depends on a careful treatment of an eigen-
decomposition of one of the matrices involved in the formulation.

Ensemble Kalman filters (EnKFs) are widely used for high-dimensional,
nonlinear, non-Gaussian data assimilation. EnKFs use an ensemble to
approximate the forecast mean and covariance, and then update these in a
manner consistent with the classical Kalman Filter formulas (Evensen, 2009).
The ensemble adjustment Kalman filter (EAKF; Anderson, 2001) is one of the
earliest EnKFs. The EAKF provides a deterministic update of the forecast
ensemble, unlike the perturbed-observation version of the EnKF (Burgers et
al., 1998, Houtekamer and Mitchell, 1998); as such it is a member of the class
of ensemble square root filters (Tippett et al., 2003). The presentation of
the EAKF in Anderson (2001) and Tippett et al. (2003) suffers from a lack of
clarity potentially leaving the practitioner confused as to the correct
implementation of the EAKF as well as confused about whether the EAKF is in
fact consistent with the Kalman filter update formula. This note clarifies the
correct formulation of the EAKF, and presents a novel detail required to
ensure consistency of the update formula with the Kalman filter covariance
update. No clarification is needed in the serial assimilation algorithm of
Anderson (2003), which bypasses the matrix algebra of the EAKF.

We follow the notation from Tippett et al. (2003). Denote the forecast
ensemble by $\bm{x}_{1}^{f},\cdots,\bm{x}_{m}^{f}$ where $m$ is the ensemble
size. The ensemble forecast covariance matrix is $\text{\bf P}^{f}=\text{\bf
Z}^{f}\text{\bf Z}^{fT}$ where $\text{\bf Z}^{f}$ is the scaled ensemble
pertubation matrix

$\text{\bf
Z}^{f}=\frac{1}{\sqrt{m-1}}\left[\bm{x}_{1}^{f}-\bm{\mu}^{f},\cdots,\bm{x}_{m}^{f}-\bm{\mu}^{f}\right]$
(1)

and $\bm{\mu}^{f}$ is the ensemble forecast mean. Denote the observation
matrix by H and the observation error covariance matrix by R.

The key to the EAKF is the formula for updating the ensemble perturbations in
such a way that they are consistent with the KF covariance update formula

$\text{\bf P}^{a}=\left(\text{\bf I}-\text{\bf KH}\right)\text{\bf P}^{f}$ (2)

where K is the Kalman gain matrix

$\text{\bf K}=\text{\bf P}^{f}\text{\bf H}^{T}\left(\text{\bf HP}^{f}\text{\bf
H}^{T}+\text{\bf R}\right)^{-1}.$ (3)

Defining $\text{\bf V}=\left(\text{\bf HZ}^{f}\right)^{T}$ and inserting
$\text{\bf P}^{f}=\text{\bf Z}^{f}\text{\bf Z}^{fT}$, the KF update formula
becomes

$\text{\bf P}^{a}=\text{\bf Z}^{f}\left[\text{\bf I}-\text{\bf
V}\left(\text{\bf V}^{T}\text{\bf V}+\text{\bf R}\right)^{-1}\text{\bf
V}^{T}\right]\text{\bf Z}^{fT}.$ (4)

The Woodbury matrix identity can be applied to write this as

$\text{\bf P}^{a}=\text{\bf Z}^{f}\left[\text{\bf I}+\text{\bf V}\text{\bf
R}^{-1}\text{\bf V}^{T}\right]^{-1}\text{\bf Z}^{fT}.$ (5)

The EAKF uses an adjustment matrix A to update the scaled ensemble
perturbation matrix as $\text{\bf Z}^{a}=\text{\bf AZ}^{f}$ in such a way that
$\text{\bf P}^{a}=\text{\bf Z}^{a}\text{\bf Z}^{aT}$ exactly satisfies Eqn.
(5). The EAKF adjustment matrix presented by Tippett et al. (2003) is

$\text{\bf A}=\text{\bf Z}^{f}\text{\bf C}\left(\text{\bf
I}+\mathbf{\Gamma}\right)^{-1/2}\text{\bf G}^{-1}\text{\bf F}^{T}$ (6)

where $\text{\bf Z}^{f}=\text{\bf FGU}^{T}$ is the singular value
decomposition (SVD) of $\text{\bf Z}^{f}$ and $\text{\bf V}\text{\bf
R}^{-1}\text{\bf V}^{T}=\text{\bf C}\mathbf{\Gamma}\text{\bf C}^{T}$ is the
eigenvalue decomposition of $\text{\bf V}\text{\bf R}^{-1}\text{\bf V}^{T}$
(cf. Eq. (19) of Tippett et al. (2003)).

The use of the notation $\text{\bf G}^{-1}$ in (6) implies that the expression
uses the version of the SVD that produces a square, invertible matrix G. If
the state dimension is $n$ and the rank of $\text{\bf Z}^{f}$ is $r$ then
$\text{\bf Z}^{f}$ is $n\times m$, and the SVD produces G of size $r\times r$,
F of size $n\times r$, and U of size $m\times r$. This is precisely the
interpretation given in a footnote by Tippett et al. (2003). In this
interpretation the only possible way for (6) to be meaningful is for $r=m$,
because otherwise the matrix $\left(\text{\bf
I}+\mathbf{\Gamma}\right)^{-1/2}$, which is $m\times m$, can’t be multiplied
by $\text{\bf G}^{-1}$, which is $r\times r$. But the rank of $\text{\bf
Z}^{f}$ is $r\leq$ min$\\{n,m-1\\}$, so (6) is a mathematical impossibility,
and the practitioner is left uncertain the correct formulation of the EAKF
adjustment matrix.

The situation can be fixed by using the version of the SVD that produces G of
size $r\times m$, F of size $n\times r$, and U of size $m\times m$, and by
replacing $\text{\bf G}^{-1}$ in (6) by the Moore-Penrose pseudoinverse
$\text{\bf G}^{\dagger}$. Here $\text{\bf G}^{\dagger}$ is an $m\times r$
diagonal matrix whose diagonal elements are $1/\sigma_{1},\ldots,1/\sigma_{r}$
where $\sigma_{i}$ are the singular values of $\text{\bf Z}^{f}$ that appear
on the diagonal of G. With this substitution the ensemble adjustment matrix
becomes

A $\displaystyle=\text{\bf Z}^{f}\text{\bf C}\left(\text{\bf
I}+\mathbf{\Gamma}\right)^{-1/2}\text{\bf G}^{{\dagger}}\text{\bf F}^{T}.$ (7)

This formulation is implied by the discussion of rank deficiency in appendix A
of Anderson (2001), though the consistency of this adjustment matrix with the
KF covariance update formula is not fully demonstrated when G is not
invertible.

We next show that this update formula leads to an analysis ensemble covariance
matrix that exactly satisfies (2), as long as the matrices C and
$\mathbf{\Gamma}$ are constructed carefully. We begin by simplifying the
expression for the scaled analysis ensemble perturbation matrix

$\displaystyle\text{\bf Z}^{a}$ $\displaystyle=\text{\bf Z}^{f}\text{\bf
C}\left(\text{\bf I}+\mathbf{\Gamma}\right)^{-1/2}\text{\bf
G}^{{\dagger}}\text{\bf F}^{T}\text{\bf FGU}^{T}$ $\displaystyle=\text{\bf
Z}^{f}\text{\bf C}\left(\text{\bf I}+\mathbf{\Gamma}\right)^{-1/2}\text{\bf
G}^{\dagger}\text{\bf G}\text{\bf U}^{T}.$ (8)

The matrix $\text{\bf G}^{\dagger}\text{\bf G}$ in this expression is an
$m\times m$ diagonal matrix whose first $r$ diagonal entries are 1, and whose
remaining entries are all 0. The product $\left(\text{\bf
I}+\mathbf{\Gamma}\right)^{-1/2}\text{\bf G}^{\dagger}\text{\bf G}$ is a
diagonal matrix whose first $r$ entries are the same as those of
$\left(\text{\bf I}+\mathbf{\Gamma}\right)^{-1/2}$, and whose remaining
entries are all 0. Denoting this product by $\left(\text{\bf
I}+\mathbf{\Gamma}\right)_{r}^{-1/2}$ we plug in to the update formula to find

$\text{\bf P}^{a}=\text{\bf Z}^{a}\text{\bf Z}^{aT}=\text{\bf Z}^{f}\text{\bf
C}\left(\text{\bf I}+\mathbf{\Gamma}\right)_{r}^{-1}\text{\bf C}^{T}\text{\bf
Z}^{fT}.$ (9)

The correct expression differs only in the r subscript, meaning that in this
formulation of EAKF the last $m-r$ elements on the diagonal of
$\left(\text{\bf I}+\mathbf{\Gamma}\right)_{r}^{-1}$ are zero when they should
not be. This suggests that the EAKF produces a posterior covariance that is
under-dispersed compared to the correct posterior covariance. Nevertheless,
this version of the EAKF still does produce the correct posterior covariance
if implemented carefully, as we now show.

First note that $\text{\bf VR}^{-1}\text{\bf V}^{T}=\text{\bf Z}^{fT}\text{\bf
H}^{T}\text{\bf R}^{-1}\text{\bf HZ}^{f}$, so the null space of $\text{\bf
Z}^{f}$ is equal to the null space of $\text{\bf VR}^{-1}\text{\bf V}^{T}$.
This implies that there are $m-r$ columns of C (the columns of C are
orthonormal eigenvectors of $\text{\bf VR}^{-1}\text{\bf V}^{T}$) that are in
the null space of $\text{\bf Z}^{f}$. Multiplying any one of these columns by
$\text{\bf Z}^{f}$ will yield $\bm{0}$, so the matrix $\text{\bf
Z}^{f}\text{\bf C}$ that appears in (7) will have $m-r$ columns equal to
$\bm{0}$. If the eigenvalue decomposition of $\text{\bf VR}^{-1}\text{\bf
V}^{T}$ is constructed so that these columns of C are the last $m-r$ columns
of C, then we have

$\text{\bf Z}^{f}\text{\bf C}\left(\text{\bf
I}+\mathbf{\Gamma}\right)_{r}^{-1/2}=\text{\bf Z}^{f}\text{\bf
C}\left(\text{\bf I}+\mathbf{\Gamma}\right)^{-1/2}$ (10)

because the $m-r$ trailing zeros on the diagonal of $\left(\text{\bf
I}+\mathbf{\Gamma}\right)_{r}^{-1/2}$ end up multiplying columns of $\text{\bf
Z}^{f}\text{\bf C}$ that are already $\bm{0}$. If, on the other hand, the
$m-r$ eigenvectors in the null space of $\text{\bf VR}^{-1}\text{\bf V}^{T}$
are not arranged as the last $m-r$ columns of C, then the $m-r$ trailing zeros
on the diagonal of $\left(\text{\bf I}+\mathbf{\Gamma}\right)_{r}^{-1/2}$ end
up multiplying nonzero columns of $\text{\bf Z}^{f}\text{\bf C}$, causing
erroneous cancellations and leading to an under-dispersed posterior ensemble.
If the columns of C are arranged correctly, then we have

$\displaystyle\text{\bf P}^{a}$ $\displaystyle=\text{\bf Z}^{a}\text{\bf
Z}^{aT}=\text{\bf Z}^{t}\text{\bf C}\left(\text{\bf
I}+\mathbf{\Gamma}\right)_{r}^{-1}\text{\bf C}^{T}\text{\bf Z}^{fT}$
$\displaystyle=\text{\bf Z}^{t}\text{\bf C}\left(\text{\bf
I}+\mathbf{\Gamma}\right)^{-1}\text{\bf C}^{T}\text{\bf Z}^{fT}$
$\displaystyle=\text{\bf Z}^{f}\left[\text{\bf I}+\text{\bf V}\text{\bf
R}^{-1}\text{\bf V}^{T}\right]^{-1}\text{\bf Z}^{fT}$ (11)

i.e. the EAKF update exactly satisfies (2).

The presentation of the EAKF in appendix A of Anderson (2001) avoids this
potential problem by using the language ‘singular value decomposition’ for the
eigenvalue decomposition of $\text{\bf V}\text{\bf R}^{-1}\text{\bf V}^{T}$.
Since $\text{\bf V}\text{\bf R}^{-1}\text{\bf V}^{T}$ is a symmetric non-
negative definite matrix, the SVD is the same as the eigenvalue decomposition.
The crucial difference is that the convention for the SVD is to always arrange
vectors corresponding to the null space as the last columns of the eigenvector
matrix C. Eigenvalue solvers often do not follow this convention, and scatter
the null space eigenvectors randomly in the columns of C.

In summary, the correct EAKF update formula is

$\text{\bf Z}^{a}=\text{\bf A}\text{\bf Z}^{f},\;\text{\bf A}=\text{\bf
Z}^{f}\text{\bf C}\left(\text{\bf I}+\mathbf{\Gamma}\right)^{-1/2}\text{\bf
G}^{{\dagger}}\text{\bf F}^{T}$ (12)

where $\text{\bf Z}^{f}=\text{\bf FGU}^{T}$ is the singular value
decomposition (SVD) of $\text{\bf Z}^{f}$ and $\text{\bf V}\text{\bf
R}^{-1}\text{\bf V}^{T}=\text{\bf C}\mathbf{\Gamma}\text{\bf C}^{T}$ is the
eigenvalue decomposition of $\text{\bf V}\text{\bf R}^{-1}\text{\bf V}^{T}$,
with the specific requirement that eigenvectors of $\text{\bf V}\text{\bf
R}^{-1}\text{\bf V}^{T}$ that are in the null space of $\text{\bf V}\text{\bf
R}^{-1}\text{\bf V}^{T}$ must be arranged as the final columns of C. The
purpose of this note is to (i) fully demonstrate the consistency of the EAKF
of Anderson (2001), (ii) clarify the correct treatment of the pseudoinverse in
the discussion of Tippett et al. (2003), and (iii) point out a potential
pitfall in the implementation of the EAKF that could lead to inconsistency and
an under-dispersed analysis ensemble if an eigenvalue decomposition is not
treated carefully.

## References

  * Anderson (2001) J. L. Anderson.  An ensemble adjustment Kalman filter for data assimilation.  _Mon. Wea. Rev._ , 129(12):2884–2903, 2001\.
  * Anderson (2003) J. L. Anderson.  A local least squares framework for ensemble filtering.  _Mon. Wea. Rev._ , 131(4):634–642, 2003.
  * Burgers et al. (1998) G. Burgers, P. Jan van Leeuwen, and G. Evensen.  Analysis scheme in the ensemble Kalman filter.  _Mon. Wea. Rev._ , 126(6):1719–1724, 1998.
  * Evensen (2009) G. Evensen.  _Data Assimilation: The Ensemble Kalman Filter_.  Springer, 2009.
  * Houtekamer and Mitchell (1998) P. L. Houtekamer and H. L. Mitchell.  Data assimilation using an ensemble Kalman filter technique.  _Mon. Wea. Rev._ , 126(3):796–811, 1998.
  * Tippett et al. (2003) M. K. Tippett, J. L. Anderson, C. H. Bishop, T. M. Hamill, and J. S. Whitaker.  Ensemble square root filters.  _Mon. Wea. Rev._ , 131(7):1485–1490, 2003.


D}))\otimes_{\mathbb{C}}\widetilde{F}^{-1}(\pi^{-1}_{\mathbf{x}}(\mathcal{L}^{\textup{univ}})),$

where $p_{1}\colon X\times\mathbb{C}^{r}\to X$ is the first projection. Since
$F^{-1}(\pi^{-1}_{\mathbf{x}}(\mathcal{L}^{\textup{univ}}))$ is
$G$-equivariant and $G$ acts trivially on objects originating from $X$ (for
instance $\mathcal{M}$ and $p_{1}^{-1}(\textup{DR}(\mathcal{M}|_{X\setminus
D}))$), (5.2) is a $G$-equivariant isomorphism. Therefore, by Lemma 4.1

$\pi_{*}^{G}(\textup{DR}_{X\times\mathbb{C}^{r}/\mathbb{C}^{r}}(\mathcal{M}(*D_{F}^{(r)})))|_{(X\setminus
D)\times(\mathbb{C}^{*})^{r}}\simeq\widetilde{\textup{DR}(\mathcal{M})|_{X\setminus
D}\otimes_{\mathbb{C}}F^{-1}\psi^{\textup{univ}}}.$

By adjunction, the isomorphism above induces a morphism

(5.3)
$\pi_{*}^{G}(\textup{DR}_{X\times\mathbb{C}^{r}/\mathbb{C}^{r}}(\mathcal{M}(*D_{F}^{(r)})))\rightarrow\widetilde{Rj_{*}}(\textup{DR}(\mathcal{M})|_{X\setminus
D}\otimes_{\mathbb{C}}F^{-1}\psi^{\textup{univ}}).$

We then make derived pullback of the morphism above by
$\widetilde{i}_{{\boldsymbol{\lambda}}}$ for every
${\boldsymbol{\lambda}}\in(\mathbb{C}^{*})^{r}$. But by definition

$\mathbf{L}\widetilde{i}_{{\boldsymbol{\lambda}}}^{*}\pi_{*}^{G}(\textup{DR}_{X\times\mathbb{C}^{r}/\mathbb{C}^{r}}(\mathcal{M}(*D_{F}^{(r)})))\simeq\mathbf{L}\widetilde{i}_{{\boldsymbol{\alpha}}}^{*}\textup{DR}_{X\times\mathbb{C}^{r}/\mathbb{C}^{r}}(\mathcal{M}(*D_{F}^{(r)}))$

for any ${\boldsymbol{\alpha}}\in Exp^{-1}({\boldsymbol{\lambda}})$. Since
$\mathcal{M}$ is regular, by Proposition 3.17

$\mathbf{L}\widetilde{i}_{{\boldsymbol{\alpha}}}^{*}(\textup{DR}_{X\times\mathbb{C}^{r}/\mathbb{C}^{r}}(\mathcal{M}(*D_{F}^{(r)})))\simeq\textup{DR}(\mathbf{L}\widetilde{i}_{{\boldsymbol{\alpha}}}^{*}(\mathcal{M}(*D_{F}^{(r)}))\simeq
Rj_{*}(\textup{DR}(\mathcal{M})|_{X\setminus D}\otimes
L_{{\boldsymbol{\lambda}}}).$

Moreover, by construction

$\mathbf{L}\widetilde{i}_{{\boldsymbol{\lambda}}}^{*}\widetilde{Rj_{*}}(\textup{DR}(\mathcal{M})|_{X\setminus
D}\otimes_{\mathbb{C}}F^{-1}\psi^{\textup{univ}})\simeq{Rj_{*}}(\textup{DR}(\mathcal{M})|_{X\setminus
D}\otimes_{\mathbb{C}}F^{-1}\psi^{\textup{univ}})\otimes^{\mathbf{L}}_{\mathbb{C}[G]}\mathbb{C}_{{\boldsymbol{\lambda}}},$

where $\mathbb{C}_{{\boldsymbol{\lambda}}}$ is the residue field of
${\boldsymbol{\lambda}}\in(\mathbb{C}^{*})^{r}$. Since $\psi^{\textup{univ}}$
is locally free over $\mathbb{C}[G]$ and hence over $\mathbb{C}$,

${Rj_{*}}(\textup{DR}(\mathcal{M})|_{X\setminus
D}\otimes_{\mathbb{C}}F^{-1}\psi^{\textup{univ}})\otimes_{\mathbb{C}[G]}\mathbb{C}_{{\boldsymbol{\lambda}}}\simeq{Rj_{*}}(\textup{DR}(\mathcal{M})|_{X\setminus
D}\otimes_{\mathbb{C}}F^{-1}\psi^{\textup{univ}}\otimes_{\mathbb{C}[G]}\mathbb{C}_{{\boldsymbol{\lambda}}}).$

By the definition of $\psi_{\textup{univ}}$,

$F^{-1}\psi^{\textup{univ}}\otimes_{\mathbb{C}[G]}\mathbb{C}_{{\boldsymbol{\lambda}}}\simeq
L_{\boldsymbol{\lambda}}.$

Therefore, the $\mathbf{L}\widetilde{i}_{{\boldsymbol{\lambda}}}^{*}$ pullback
of the morphism (5.3) is an isomorphism for each ${\boldsymbol{\lambda}}$. By
[MFS13, Proposition 2.2 and Theorem 3.7], the morphism (5.3) is isomorphic. ∎

Although Theorem 1.1 is a special case of Theorem 1.2 (by Corollary 3.11), we
first give a direct proof of Theorem 1.1 as a warm-up. We keep using notation
in the proof of Theorem 5.3. By Lemma 3.7,

$\Psi_{F}(\mathcal{M})\stackrel{{\scriptstyle
q.i.}}{{\simeq}}\textup{cone}(\mathcal{M}(!D_{F}^{(r)})\rightarrow\mathcal{M}(*D_{F}^{(r)})).$

By Proposition 3.17,

$\mathbf{L}\widetilde{i}_{{\boldsymbol{\alpha}}}^{*}\textup{cone}(\mathcal{M}(!D_{F}^{(r)})\rightarrow\mathcal{M}(!D_{F}^{(r)}))\stackrel{{\scriptstyle
q.i.}}{{\simeq}}\textup{cone}(\mathcal{M}(!D)\to\mathcal{M}(*D)).$

Taking the relative de Rham functor, we then have that

$\textup{DR}_{X\times\mathbb{C}^{r}/\mathbb{C}^{r}}\mathbf{L}\widetilde{i}_{{\boldsymbol{\alpha}}}^{*}\textup{cone}(\mathcal{M}(!D_{F}^{(r)})\rightarrow\mathcal{M}(!D_{F}^{(r)}))$

is quasi-isomorphic to

$\mathscr{B}^{\bullet}\coloneqq\textup{cone}(j_{!}(\textup{DR}(\mathcal{M})|_{X\setminus
D}\otimes L_{{\boldsymbol{\lambda}}})\to
Rj_{*}(\textup{DR}(\mathcal{M})|_{X\setminus D}\otimes
L_{{\boldsymbol{\lambda}}})).$

But

$\mathscr{B}^{\bullet}\stackrel{{\scriptstyle
q.i.}}{{\simeq}}i_{D,*}i^{-1}_{D}Rj_{*}(\textup{DR}(\mathcal{M})|_{X\setminus
D}\otimes L_{{\boldsymbol{\lambda}}}).$

Thanks to [MFS13, Proposition 2.2], we have a natural isomorphism

$\textup{DR}_{X\times\mathbb{C}^{r}/\mathbb{C}^{r}}(\Psi_{F}(\mathcal{M}))\stackrel{{\scriptstyle
q.i.}}{{\simeq}}\widetilde{i}_{D,*}\widetilde{i}_{D}^{-1}\textup{DR}_{X\times\mathbb{C}^{r}/\mathbb{C}^{r}}(\mathcal{M}(*D_{F}^{(r)}))$

where $\widetilde{i}_{D}=(i_{D},\textup{id})\colon
D\times\mathbb{C}^{r}\hookrightarrow X\times\mathbb{C}^{r}$ is the associated
morphism. Since $\pi_{*}^{G}$ and $\widetilde{i}_{D,*}\widetilde{i}_{D}^{-1}$
commute, the proof of Theorem 1.1 is accomplished by Theorem 5.3 and
functoriality of sheafification.

###### Proof of Theorem 1.2.

By Lemma 3.13 and Theorem 3.12, we have

$\displaystyle\pi_{*}^{G}\big{(}\textup{DR}_{X\times\mathbb{C}^{r}/\mathbb{C}^{r}}\mathbb{D}R\Gamma_{[D_{I}\times\mathbb{C}^{r}]}\big{(}(\mathbb{D}\mathcal{M})(!D^{(r)}_{F})\big{)}\big{)}$
$\displaystyle\simeq\textup{DR}_{X\times(\mathbb{C}^{*})^{r}/(\mathbb{C}^{*})^{r}}\mathbb{D}R\Gamma_{[D_{I}\times(\mathbb{C}^{*})^{r}]}\big{(}\pi_{*}^{G}\big{(}(\mathbb{D}\mathcal{M})(!D^{(r)}_{F})\big{)}\big{)}.$

By [MFS13, Theorem 3.11], we have

$\displaystyle\textup{DR}_{X\times(\mathbb{C}^{*})^{r}/(\mathbb{C}^{*})^{r}}\mathbb{D}R\Gamma_{[D_{I}\times(\mathbb{C}^{*})^{r}]}\big{(}\pi_{*}^{G}\big{(}(\mathbb{D}\mathcal{M})(!D^{(r)}_{F})\big{)}\big{)}$
$\displaystyle\simeq\mathbf{D}\textup{DR}_{X\times(\mathbb{C}^{*})^{r}/(\mathbb{C}^{*})^{r}}R\Gamma_{[D_{I}\times(\mathbb{C}^{*})^{r}]}\big{(}\pi_{*}^{G}\big{(}(\mathbb{D}\mathcal{M})(!D^{(r)}_{F})\big{)}\big{)},$

where $\mathbf{D}$ denotes the duality functor for relative constructible
complexes (cf. [MFS13, §2.6]). By Proposition 3.17,
$(\mathbb{D}\mathcal{M})(!D^{(r)}_{F})$ is relative regular holonomic (cf.
[MFS19, Definition 2.1]) and thus so is
$\pi_{*}^{G}\big{(}(\mathbb{D}\mathcal{M})(!D^{(r)}_{F})\big{)}$. By
regularity and Lemma 2.11, DR and $R\Gamma$ commute and thus we have

$\displaystyle\mathbf{D}\textup{DR}_{X\times(\mathbb{C}^{*})^{r}/(\mathbb{C}^{*})^{r}}R\Gamma_{[D_{I}\times(\mathbb{C}^{*})^{r}]}\big{(}\pi_{*}^{G}\big{(}(\mathbb{D}\mathcal{M})(!D^{(r)}_{F})\big{)}\big{)}$
$\displaystyle\simeq\mathbf{D}R\Gamma_{D_{I}\times(\mathbb{C}^{*})^{r}}\textup{DR}_{X\times(\mathbb{C}^{*})^{r}/(\mathbb{C}^{*})^{r}}\big{(}\pi_{*}^{G}\big{(}(\mathbb{D}\mathcal{M})(!D^{(r)}_{F})\big{)}\big{)}.$

By (3.11), we know
$(\mathbb{D}\mathcal{M})(!D^{(r)}_{F})\simeq\mathbb{D}(\mathcal{M}(*D_{F}^{(r)}))$
and thus

$\displaystyle\mathbf{D}R\Gamma_{D_{I}\times(\mathbb{C}^{*})^{r}}\textup{DR}_{X\times(\mathbb{C}^{*})^{r}/(\mathbb{C}^{*})^{r}}\big{(}\pi_{*}^{G}\big{(}(\mathbb{D}\mathcal{M})(!D^{(r)}_{F})\big{)}\big{)}$
$\displaystyle\simeq\mathbf{D}R\Gamma_{D_{I}\times(\mathbb{C}^{*})^{r}}\textup{DR}_{X\times(\mathbb{C}^{*})^{r}/(\mathbb{C}^{*})^{r}}\big{(}\pi_{*}^{G}\mathbb{D}\big{(}\mathcal{M}(*D_{F}^{(r)})\big{)}\big{)}$
$\displaystyle\simeq\mathbf{D}R\Gamma_{D_{I}\times(\mathbb{C}^{*})^{r}}\mathbf{D}\textup{DR}_{X\times(\mathbb{C}^{*})^{r}/(\mathbb{C}^{*})^{r}}\big{(}\pi_{*}^{G}\big{(}\mathcal{M}(*D_{F}^{(r)})\big{)}\big{)},$

where in the second isomorphism above follows from (3.13) in Remark 3.10(2).
By [KS13, Proposition 3.1.11] and the construction of $\mathbf{D}$, we have

$\mathbf{D}R\Gamma_{D_{I}\times(\mathbb{C}^{*})^{r}}\mathbf{D}\simeq\widetilde{i}_{D_{I},*}\mathbf{D}{\tilde{i}}_{D_{I}}^{!}\mathbf{D}$

where $\widetilde{i}_{D_{I}}=(i_{D_{I}},\textup{id})\colon
D_{I}\times\mathbb{C}^{r}\hookrightarrow X\times\mathbb{C}^{r}$ is the
associated morphism. Moreover, by [KS13, Proposition 3.1.13] (see also [MFS13,
Remark 2.24]), for relative constructible complexes we have a functorial
isomorphism

$\mathbf{D}{\tilde{i}}_{D_{I}}^{!}\mathbf{D}\simeq\widetilde{i}_{D_{I}}^{-1}.$

Therefore,

$\displaystyle\mathbf{D}R\Gamma_{D_{I}\times(\mathbb{C}^{*})^{r}}\mathbf{D}\textup{DR}_{X\times(\mathbb{C}^{*})^{r}/(\mathbb{C}^{*})^{r}}\big{(}\pi_{*}^{G}\big{(}\mathcal{M}(*D_{F}^{(r)})\big{)}\big{)}$
$\displaystyle\simeq\widetilde{i}_{D_{I},*}\widetilde{i}_{D_{I}}^{-1}\textup{DR}_{X\times(\mathbb{C}^{*})^{r}/(\mathbb{C}^{*})^{r}}\big{(}\pi_{*}^{G}\big{(}\mathcal{M}(*D_{F}^{(r)})\big{)}\big{)}.$

By Theorem 5.3 we finally have

$\widetilde{i}_{D_{I},*}\widetilde{i}_{D_{I}}^{-1}\textup{DR}_{X\times(\mathbb{C}^{*})^{r}/(\mathbb{C}^{*})^{r}}\big{(}\pi_{*}^{G}\big{(}\mathcal{M}(*D_{F}^{(r)})\big{)}\big{)}\simeq\widetilde{\psi}_{D_{I}}(\textup{DR}_{X}(\mathcal{M})).$

∎

###### Proof of Corollary 1.3.

By Lemma 2.20, we have

$\textup{supp}_{(\mathbb{C}^{*})^{r}}\psi_{F}(\textup{DR}\mathcal{M})=\\{{\boldsymbol{\lambda}}\in(\mathbb{C}^{*})^{r}\mid\mathbf{L}\widetilde{i}^{*}_{{\boldsymbol{\lambda}}}(\widetilde{\psi}_{F}(\textup{DR}\mathcal{M}))\not=0\\}.$

Since $\mathbf{L}\widetilde{i}^{*}_{{\boldsymbol{\alpha}}}$ and DR commute, by
the Riemann-Hilbert correspondence for regular holonomic
$\mathscr{D}_{X}$-modules we have

$\\{{\boldsymbol{\alpha}}\in\mathbb{C}^{r}\mid\mathbf{L}\widetilde{i}^{*}_{{\boldsymbol{\alpha}}}(\textup{DR}_{X\times\mathbb{C}^{r}/\mathbb{C}^{r}}\Psi_{F}(\mathcal{M}))\not=0\\}=\\{{\boldsymbol{\alpha}}\in\mathbb{C}^{r}\mid\mathbf{L}\widetilde{i}^{*}_{{\boldsymbol{\alpha}}}(\Psi_{F}(\mathcal{M}))\not=0\\}$

By Theorem 1.1,

$G\backslash\\{{\boldsymbol{\alpha}}\in\mathbb{C}^{r}\mid\mathbf{L}\widetilde{i}^{*}_{{\boldsymbol{\alpha}}}(\textup{DR}_{X\times\mathbb{C}^{r}/\mathbb{C}^{r}}\Psi_{F}(\mathcal{M}))\not=0\\}=\\{{\boldsymbol{\lambda}}\in(\mathbb{C}^{*})^{r}\mid\mathbf{L}\widetilde{i}^{*}_{{\boldsymbol{\lambda}}}(\widetilde{\psi}_{F}(\textup{DR}\mathcal{M}))\not=0\\}.$

Thus, it is enough to prove

$\textup{supp}_{\mathbb{C}^{r}}\Psi_{F}(\mathcal{M})=\\{{\boldsymbol{\alpha}}\in\mathbb{C}^{r}\mid\mathbf{L}\widetilde{i}^{*}_{{\boldsymbol{\alpha}}}(\Psi_{F}(\mathcal{M}))\not=0\\}.$

But the equality above is the analytification of [BVWZ21a, Proposition 3.4.3]
as $\Psi_{F}(\mathcal{M})$ is $(n+1)$-Cohen-Macaulay. We leave details for
interested readers; see also [BVWZ21a, §3.6]. ∎

###  5.3. Proof of the local index formula

In this subsection, we prove the local index formula (1.3) in §1.3. We keep
notations as in §1.3 and §3. We focus on a small neighborhood $W\subseteq X$
around $x$.

We write $\mathcal{F}^{\bullet}=\textup{DR}(\mathcal{M})$ for simplicity and
let $q$ be the generic point of an irreducible component of
$\textup{supp}_{(\mathbb{C}^{*})^{r}}\psi_{F}\mathcal{F}^{\bullet}$. By
Corollary 1.5, $\bar{q}$ is always a translated subtori of codimension-one. We
pick a general point ${\boldsymbol{\lambda}}\in\bar{q}$ and a point
${\boldsymbol{\alpha}}\in Exp^{-1}({\boldsymbol{\lambda}})$. Then for
$i\in\\{1,\dots,r\\}$ we consider the closed embedding

$\delta_{\hat{\boldsymbol{\lambda}}_{i}}\colon\mathbb{C}^{*}\hookrightarrow(\mathbb{C}^{*})^{r},\quad
t\mapsto(\lambda_{1},\dots,\lambda_{i-1},t,\lambda_{i+1},\dots,\lambda_{r}).$

By pullback the quasi-isomorphism in Theorem 1.1 and Corollary 3.18 we obtain:

(5.4)
$\widetilde{\delta}_{\hat{\boldsymbol{\lambda}}_{i}}^{*}\big{(}\pi_{*}^{G}\textup{DR}_{X\times\mathbb{C}/\mathbb{C}}\big{(}\Psi_{f_{i}}(\mathcal{M}_{\hat{\boldsymbol{\alpha}}_{i}})\big{)}\big{)}\simeq\widetilde{\psi}_{f_{i}}(\mathcal{F}_{\hat{\boldsymbol{\alpha}}_{i}}),$

with
$\mathcal{F}_{\hat{\boldsymbol{\alpha}}_{i}}=\textup{DR}(\mathcal{M}_{\hat{\boldsymbol{\alpha}}_{i}})$.
By the definition of $\mathcal{M}_{\hat{\boldsymbol{\alpha}}_{i}}$, (5.4) is
independent of the choices of ${\boldsymbol{\alpha}}\in
Exp^{-1}({\boldsymbol{\lambda}})$.

We pick a sufficient small neighborhood $V_{{\boldsymbol{\alpha}}}$ of
${\boldsymbol{\alpha}}\in\mathbb{C}^{r}$ such that
$Exp(V_{{\boldsymbol{\alpha}}})\simeq V_{{\boldsymbol{\alpha}}}$ (since $Exp$
is the universal covering) and write

$V_{\alpha_{i}}=\delta_{\hat{\boldsymbol{\alpha}}_{i}}^{-1}(V_{{\boldsymbol{\alpha}}}).$

We then take restriction of (5.4) on $W\times Exp(V_{\alpha_{i}})$. Since
$\Psi_{f_{i}}(\mathcal{M}_{\hat{\boldsymbol{\alpha}}_{i}})|_{W\times
V_{\alpha_{i}}}$ is only supported on $W\times\\{\alpha_{i}\\}$, we can treat
$\Psi_{f_{i}}(\mathcal{M}_{\hat{\boldsymbol{\alpha}}_{i}})|_{W\times
V_{\alpha_{i}}}$ as a regular holonomic $\mathscr{D}_{X}$-module (over W) and
hence we get

(5.5)
$\textup{DR}_{X}(\Psi_{f_{i}}(\mathcal{M}_{\hat{\boldsymbol{\alpha}}_{i}})|_{W\times
V_{\alpha_{i}}})\simeq\psi_{f_{i}}(\mathcal{F}_{\hat{\boldsymbol{\alpha}}_{i}})_{m_{\alpha_{i}}}$

where $m_{\alpha_{i}}$ is the maximal ideal of $\alpha_{i}\in\mathbb{C}$.

######  Lemma 5.4.

With notations as above, we have

$\chi_{x}(\psi_{F}(\mathcal{F}^{\bullet}),q)=\chi_{x}(\psi_{F}(\mathcal{F}^{\bullet}_{\hat{\boldsymbol{\alpha}}_{i}}),m_{\alpha_{i}}).$

###### Proof.

Since the length function is additive with respect to short exact sequence,
the required identity follows from the general choice of
${\boldsymbol{\alpha}}$ in $\bar{q}$. ∎

We now apply [Gin86, Index Theorem 8.2] to the regular holonomic
$\mathscr{D}_{X}$-module
$\Psi_{f_{i}}(\mathcal{M}_{\hat{\boldsymbol{\alpha}}_{i}})|_{V_{\alpha_{i}}}$.
In consequence, the formula (1.3) follows from (5.5), Lemma 5.4 and
Proposition 3.19.

###  5.4. A question about an algebraic relative Riemann-Hilbert
correspondence

Both Theorem 1.1 and Theorem 5.3 can be seen as some special (regular)
relative Riemann-Hilbert correspondence over arbitrary dimensional base
spaces. See [FFS21] for the general regular relative Riemann-Hilbert
correspondence over curves.

Let us wrap up this article by asking the following question, motivated by our
method in approaching Theorem 1.1:

######  Question 5.5.

For $\mathcal{F}^{\bullet}$ a constructible complex of sheaves of $R$-modules
on a complex manifold $X$ with $R$ a commutative $\mathbb{C}$-algebra of
finite type, is it possible to construct a (regular) holonomic
$\mathscr{D}_{X\times\mathbb{C}^{r}/\mathbb{C}^{r}}$-module
$\widetilde{\mathscr{M}}$ (or a complex with (regular) relative holonomic
cohomology sheaves) supported on $X\times\textup{Spec}^{\textup{an}}R$ such
that naturally

$\textup{DR}_{X\times\mathbb{C}^{r}/\mathbb{C}^{r}}(\widetilde{\mathscr{M}})\stackrel{{\scriptstyle
q.i.}}{{\simeq}}\widetilde{\mathcal{F}}^{\bullet}$

by fixing an embedding $\textup{Spec }R\hookrightarrow\mathbb{C}^{r}$?

## References

  * [BBDG18] Alexander Beilinson, Joseph Bernstein, Pierre Deligne, and Ofer Gabber, _Faisceaux pervers_ , Société mathématique de France, 2018.
  * [Bei87] Alexander Beilinson, _How to glue perverse sheaves_ , $K$-theory, arithmetic and geometry (Moscow, 1984–1986), Lecture Notes in Math., vol. 1289, Springer, Berlin, 1987, pp. 42–51.
  * [BG12] Alexander Beilinson and Dennis Gaitsgory, _A corollary of the b-function lemma_ , Selecta Mathematica 18 (2012), no. 2, 319–327.
  * [Bj93] Jan-Erik Björk, _Analytic ${\mathcal{D}}$-modules and applications_, Mathematics and its Applications, vol. 247, Kluwer Academic Publishers Group, Dordrecht, 1993.
  * [BL94] Joseph Bernstein and Valery Lunts, _Equivariant sheaves and functors_ , Lecture Notes in Mathematics, vol. 1578, Springer-Verlag, Berlin, 1994.
  * [Bry86] Jean-Luc Brylinski, _Transformations canoniques, dualité projective, théorie de Lefschetz, transformations de Fourier et sommes trigonométriques_ , Astérisque 140 (1986), no. 141, 3–134.
  * [Bud15] Nero Budur, _Bernstein-sato ideals and local systems_ , Annales de l’Institut Fourier 65 (2015), no. 2, 549–603 (en).
  * [BLSW17] Nero Budur, Yongqiang Liu, Luis Saumell, and Botong Wang, _Cohomology support loci of local systems_ , Michigan Mathematical Journal 66 (2017), 295–307.
  * [BVWZ21a] Nero Budur, Robin van der Veer, Lei Wu, and Peng Zhou, _Zero loci of Bernstein-Sato ideals_ , Inventiones mathematicae 225 (2021), no. 1, 45–72.
  * [BVWZ21b] by same author, _Zero loci of Bernstein-Sato ideals II_ , Sel. Math. New Ser. 27, 32 (2021), no. 3, 1022–1824.
  * [BW15] Nero Budur and Botong Wang, _Cohomology jump loci of quasi-projective varieties_ , Ann. Sci. Éc. Norm. Supér. (4) 48 (2015), no. 1, 227–236.
  * [BW17] by same author, _Local systems on analytic germ complements_ , Adv. Math. 306 (2017), 905–928.
  * [FF18] Luisa Fiorot and Teresa Monteiro Fernandes, _t-structures for relative D-modules and t-exactness of the de Rham functor_ , J. Algebra 509(2018), 419–444.
  * [FFS21] Luisa Fiorot, Teresa Monteiro Fernandes, and Claude Sabbah, _Relative regular Riemann-Hilbert correspondence_ , Proceedings of the London Mathematical Society 122 (2021), no. 3, 434–457.
  * [FS19] Teresa Monteiro Fernandes and Claude Sabbah, _Riemann-Hilbert correspondence for mixed twistor-modules_ , Journal of the Institute of Mathematics of Jussieu 18 (2019), no. 3, 629–672.
  * [Gin86] Victor Ginsburg, _Characteristic varieties and vanishing cycles_ , Inventiones Mathematicae 84 (1986), no. 2, 327–402.
  * [Gyo93] Akihiko Gyoja, _Bernstein-sato’s polynomial for several analytic functions_ , Journal of Mathematics of Kyoto University 33 (1993), no. 2, 399–411.
  * [HTT08] Ryoshi Hotta, Kiyoshi Takeuchi, and Toshiyuki Tanisaki, _$D$ -modules, perverse sheaves, and representation theory_, Progress in Mathematics, vol. 236, Birkhäuser Boston, Inc., Boston, MA, 2008, Translated from the 1995 Japanese edition by Takeuchi.
  * [Kas83] Masaki Kashiwara, _Vanishing cycle sheaves and holonomic systems of differential equations_ , Algebraic geometry (Tokyo/Kyoto, 1982), Lecture Notes in Math., vol. 1016, Springer, Berlin, 1983, pp. 134–142.
  * [Kas03] by same author, _$D$ -modules and microlocal calculus_, Translations of Mathematical Monographs, vol. 217, American Mathematical Society, Providence, RI, 2003, Translated from the 2000 Japanese original by Mutsumi Saito, Iwanami Series in Modern Mathematics. MR 1943036
  * [Kas77] by same author, _$B$ -functions and holonomic systems. Rationality of roots of $b$-functions_, Invent. Math. 38 (1976/77), no. 1, 33–53. MR 0430304
  * [KS13] Masaki Kashiwara and Pierre Schapira, _Sheaves on manifolds: With a short history.Les débuts de la théorie des faisceaux, by Christian Houzel_ , vol. 292, Springer Science & Business Media, 2013.
  * [Lau83] Gérard Laumon, _Sur la catégorie dérivée des ${\mathcal{D}}$-modules filtrés_, Algebraic geometry (Tokyo/Kyoto, 1982), Lecture Notes in Math., vol. 1016, Springer, Berlin, 1983, pp. 151–237. MR 726427
  * [Mai16a] Philippe Maisonobe, _Filtration Relative, l’Idéal de Bernstein et ses pentes_ , arXiv preprint arXiv:1610.03354 (2016).
  * [Mai16b] by same author, _L’Idéal de Bernstein d’un arrangement libre d’hyperplans linéaires_ , preprint arXiv:1610.03356 (2016).
  * [Mal83] B. Malgrange, _Polynômes de Bernstein-Sato et cohomologie évanescente_ , Analysis and topology on singular spaces, II, III (Luminy, 1981), Astérisque, vol. 101, Soc. Math. France, Paris, 1983, pp. 243–267.
  * [Mal04] B. Malgrange, _On irregular holonomic D-modules_ , in _Eléments de la théorie des systemes différentiels géométriques_ , Séminaires & Congrés, vol. 8 (2004), 391-410.
  * [MFS13] Teresa Monteiro Fernandes and Claude Sabbah, _On the de Rham complex of mixed twistor $\mathcal{D}$-modules_, Int. Math. Res. Not. IMRN (2013), no. 21, 4961–4984.
  * [MFS19] Teresa Monteiro Fernandes and Claude Sabbah, _Riemann-Hilbert correspondence for mixed twistor D-Modules_ , Journal of the Institute of Mathematics of Jussieu 18 (2013), no. 3, 629–672.
  * [NS96] Nitin Nitsure and Claude Sabbah, _Moduli of pre- $\mathcal{D}$-modules, perverse sheaves and the Riemann-Hilbert morphism-I_, Math. Annalen 306 (1996), 47–73.
  * [Sab87a] Claude Sabbah, _Proximité évanescente. I. La structure polaire d’un ${\mathcal{D}}$-module_, Compositio Math. 62 (1987), no. 3, 283–328.
  * [Sab87b] by same author, _Proximité évanescente. II. Équations fonctionnelles pour plusieurs fonctions analytiques_ , Compositio Math. 64 (1987), no. 2, 213–241.
  * [Sab90] by same author, _Modules d’Alexander et $\mathcal{D}$-modules_, Duke Mathematical Journal 60 (1990), no. 3, 729–814.
  * [Sch12] Pierre Schapira, _Microdifferential systems in the complex domain_ , vol. 269, Springer Science & Business Media, 2012.
  * [Sta18] The Stacks Project Authors, _Stacks Project_ , https://stacks.math.columbia.edu, 2018.
  * [Wu20] Lei Wu, _Characteristic cycles associated to holonomic $\mathcal{D}$-modules_, arXiv preprint arXiv:2011.14527 (2020).
  * [Wu21] by same author, _On the comparison of nearby cycles via $b$-functions_, J. Singul. 23 (2021), 92-106.
  * [WZ21] Lei Wu and Peng Zhou, _Log D-modules and index theorems_ , Forum of Mathematics, Sigma 9 (2021).