NUHEP-TH/21-17, FERMILAB-PUB-21-560-T, IPPP/21/46

#  $pp$ Solar Neutrinos at DARWIN

André de Gouvêa  Northwestern University, Department of Physics & Astronomy,
2145 Sheridan Road, Evanston, IL 60208, USA  Emma McGinness  University of
California Berkeley, Department of Physics, 366 Physics North, Berkeley, CA
94720, USA  University of California Berkeley, Department of Astronomy, 501
Campbell Hall, Berkeley, CA 94720, USA  Ivan Martinez-Soler  Northwestern
University, Department of Physics & Astronomy, 2145 Sheridan Road, Evanston,
IL 60208, USA  Colegio de Física Fundamental e Interdisciplinaria de las
Américas (COFI)
254 Norzagaray street, San Juan, Puerto Rico 00901  Theoretical Physics
Department, Fermilab, P.O. Box 500, Batavia, IL 60510, USA  Department of
Physics & Laboratory for Particle Physics and Cosmology, Harvard University,
Cambridge, MA 02138, USA  Yuber F. Perez-Gonzalez  Northwestern University,
Department of Physics & Astronomy, 2145 Sheridan Road, Evanston, IL 60208, USA
Colegio de Física Fundamental e Interdisciplinaria de las Américas (COFI)
254 Norzagaray street, San Juan, Puerto Rico 00901  Theoretical Physics
Department, Fermilab, P.O. Box 500, Batavia, IL 60510, USA  Institute for
Particle Physics Phenomenology, Durham University, South Road, Durham, United
Kingdom.

###### Abstract

The DARWIN collaboration recently argued that DARWIN (DARk matter WImp search
with liquid xenoN) can collect, via neutrino–electron scattering, a large,
useful sample of solar $pp$-neutrinos, and measure their survival probability
with sub-percent precision. We explore the physics potential of such a sample
in more detail. We estimate that, with 300 ton-years of data, DARWIN can also
measure, with the help of current solar neutrino data, the value of
$\sin^{2}\theta_{13}$, with the potential to exclude $\sin^{2}\theta_{13}=0$
close to the three-sigma level. We explore in some detail how well DARWIN can
constrain the existence of a new neutrino mass-eigenstate $\nu_{4}$ that is
quasi-mass-degenerate with $\nu_{1}$ and find that DARWIN’s sensitivity
supersedes that of all current and near-future searches for new, very light
neutrinos. In particular, DARWIN can test the hypothesis that $\nu_{1}$ is a
pseudo-Dirac fermion as long as the induced mass-squared difference is larger
than $10^{-13}$ eV2, one order of magnitude more sensitive than existing
constraints. Throughout, we allowed for the hypotheses that DARWIN is filled
with natural xenon or 136Xe-depleted xenon.

##  I Introduction

Multi-ton-scale, next-generation dark matter experiments are expected to
collect significant statistics of atmospheric and solar neutrinos. The DARWIN
collaboration recently argued that DARWIN (DARk matter WImp search with liquid
xenoN) can collect a large, useful sample of solar $pp$-neutrinos, measured
via elastic neutrino–electron scattering Aalbers et al. (2020). There, they
argued that the survivial probability of $pp$-neutrinos can be measured with
sub-percent precision and that one can measure the Weinberg angle at low
momentum transfers with 10% precision, independent from the values of the
neutrino oscillation parameters. Here, we explore other neutrino-physics-
related information one can obtain from a high-statistics, high-precision
measurement of the $pp$-neutrino flux.

In a nutshell, $pp$-neutrinos are produced in the solar core via proton-proton
fusion: $p+p\to{\rm{}^{2}H}+e^{+}+\nu$. The vast majority of neutrinos
produced by the fusion cycle that powers our Sun are produced via proton-
proton fusion. $pp$-neutrinos have the lowest energy among all solar neutrino
“types” (other types include $pep$-neutrinos, ${}^{7}{\rm Be}$-neutrinos,
${}^{8}{\rm B}$-neutrinos, and CNO-neutrinos) and are characterized by a
continuous spectrum that peaks around 300 keV and terminates around 420 keV.
Theoretically, the $pp$-neutrino flux is known at better than the percent
level Bahcall et al. (2001) given they are created early in the $pp$-fusion
cycle – they are the first link in the chain – and their flux is highly
correlated with the photon flux, measured with exquisite precision. For the
sake of comparison, the flux of ${}^{7}{\rm Be}$-neutrinos and ${}^{8}{\rm
B}$-neutrinos, which provide virtually all information on the particle-physics
properties of solar neutrinos, can be computed at, approximately, the 6% and
12% level, respectively Serenelli et al. (2011); Vinyoles et al. (2017). The
$pp$-neutrino flux has been directly measured, independent from the other
flux-types, by the Borexino collaboration Bellini et al. (2014), with 10%
precision.

A percent-level measurement of the $pp$-neutrino flux is expected to be
sensitive to new-physics effects in neutrino physics that are also at the
percent level. This includes, for example, effects from the so-called reactor
angle $\theta_{13}$ – not new physics but very small for solar neutrinos – and
the presence of new neutrino states or neutrino interactions. Furthermore, the
fact that $pp$-neutrinos have energies that are significantly lower than those
of the other solar neutrino types renders them especially well-suited to
constrain (or discover) new, very long oscillation lengths associate to very
small new neutrino mass-squared differences. These searches are expected to
add significantly to our ability to test the hypothesis that the neutrinos are
pseudo-Dirac fermions Wolfenstein (1981); Petcov (1982); Bilenky and
Pontecorvo (1983) (for relevant recent discussions, see, for example, de
Gouvêa et al. (2009); Donini et al. (2011); de Gouvêa et al. (2020); Martinez-
Soler et al. (2021)). Here, as far as new-physics hypotheses are concerned, we
concentrate on the search for new, very light neutrino states.

In Sec. II, we review the relevant features of the proposed DARWIN experiment
and provide information on how we simulate and analyze DARWIN data on
$pp$-neutrinos. In Sec. III, we show that a percent-level measurement of the
$pp$-neutrino flux allows for a “solar-neutrinos-only” measurement of
$\sin^{2}\theta_{13}$. In Sec. IV, we compute the sensitivity of DARWIN to the
hypothesis that there is a fourth neutrino with a mass $m_{4}$ that is quasi-
degenerate with the mass of the first neutrino state, $m_{1}$ (in the
Appendix, we discuss how this can be generalized). We concentrate on the
region of parameter space where the new mass-squared difference is
$10^{-13}~{}{\rm eV^{2}}\lesssim|m_{4}^{2}-m_{1}^{2}|\lesssim 10^{-6}$ eV2. We
add some concluding remarks in Sec. V.

##  II DARWIN as a Low-Energy Solar Neutrino Experiment

DARWIN is projected to be a large – 40 tons fiducial volume – liquid xenon
time-projection chamber, aimed at searching for weakly interacting massive
particles (WIMP) in the GeV to TeV mass range Aalbers et al. (2016) via
elastic WIMP–nucleon scattering. It will inevitably be exposed to a large flux
of solar and atmospheric neutrinos and is large enough that solar-neutrino
scattering events will occur at an observable rate.

According to Aalbers et al. (2020), DARWIN is expected to collect a sample of
almost ten thousand $pp$-neutrinos per year via elastic neutrino–electron
scattering:

$\nu_{\alpha}+e^{-}\to\nu_{\alpha}+e^{-},$ (II.1)

where $\alpha=e,\mu,\tau$ is the flavor of the incoming neutrino. The flavor
of the outgoing neutrinos is, of course, never observed. For $pp$-neutrino
energies, the cross section for $\nu_{e}e$-scattering is around six times
larger than that of $\nu_{a}e$-scattering, $a=\tau,\mu$ and the differences
between the cross sections for $\nu_{\mu}e$-scattering and
$\nu_{\tau}e$-scattering are negligible. At leading order in the weak
interactions, the differential cross section in the rest frame of the electron
is

$\frac{{\rm d}\sigma}{{\rm
d}T}(\nu_{\alpha}+e^{-}\to\nu_{\alpha}+e^{-})=\frac{2G_{F}^{2}m_{e}}{\pi}\left[a_{\alpha}^{2}+b_{\alpha}^{2}\left(1-\frac{T}{E_{\nu}}\right)^{2}-a_{\alpha}b_{\alpha}\frac{T}{E_{\nu}}\right],$
(II.2)

where $T$ is the kinetic energy of the recoil electron, $E_{\nu}$ is the
incoming neutrino energy, $m_{e}$ is the electron mass and $G_{F}$ is the
Fermi constant. The dimensionless couplings $a_{\alpha},b_{\alpha}$ are

$a_{e}=-\frac{1}{2}-\sin^{2}\theta_{W},~{}~{}~{}b_{e}=-\sin^{2}\theta_{W};~{}~{}~{}a_{a}=\frac{1}{2}-\sin^{2}\theta_{W},~{}~{}~{}b_{a}=-\sin^{2}\theta_{W},$
(II.3)

where $\theta_{W}$ is the weak mixing angle. DARWIN measures the kinetic
energy spectrum of the recoil electrons.

If filled with natural xenon, one expects a large number of electron-events in
the energy range of interest from the double-beta decays of ${}^{136}\rm Xe$.
These events are a powerful source of background for solar-neutrino studies
and, according to Aalbers et al. (2020), obviate the study of solar neutrinos
with energies higher than 1 MeV. They are a powerful nuisance for measurements
of the ${}^{7}\rm Be$-neutrino flux and have a significant but not decisive
impact on the measurement of the $pp$-neutrinos (around a 30% decrease in the
precision with which the overall $pp$-neutrino flux can be measured Aalbers et
al. (2020)). The reason one can measure the $pp$-neutrino flux in spite of the
${}^{136}\rm Xe$ background is that the shape of this particular background is
well known and the experiment can detect events over a large range of recoil-
electron energies, effectively measuring it with excellent precision. There is
the possibility of filling DARWIN with liquid xenon depleted of the double-
beta-decaying ${}^{136}\rm Xe$ isotope. This would allow the study of higher
energy solar neutrinos. Here we consider these two different scenarios, i.e.,
the ${}^{136}\rm Xe$-depleted version of DARWIN and the one where the
abundance of ${}^{136}\rm Xe$ agrees with natural expectations.

Other than the background from ${}^{136}\rm Xe$, for $pp$-neutrinos, the
double electron capture decay of ${}^{124}\rm Xe$ leads to two narrow peaks at
37 keV and 10 keV Aalbers et al. (2020) and, at higher recoil energies,
radioactive backgrounds from the detector components and the liquid volume
supersede the $pp$-neutrino events for recoil kinetic energies above 200 keV
or so. When simulating DARWIN data, we restrict our sample to events with
recoil kinetic energies below 220 keV and assume that, in this energy range,
the only sources of background are those from ${}^{136}\rm Xe$ and
${}^{124}\rm Xe$. We simulate the backgrounds using the results published in
Aalbers et al. (2020). When analyzing the simulated data, we marginalize over
the normalization of the two ${}^{124}\rm Xe$ lines, which we treat as free
parameters, and the normalization of the ${}^{136}\rm Xe$ recoil spectrum,
which we assume is independently measured with 0.1% precision. We assume the
shape of the ${}^{136}\rm Xe$ recoil spectrum is known with infinite
precision. For the ${}^{136}\rm Xe$-depleted version of DARWIN, we assume the
${}^{136}\rm Xe$-background is $1\%$ of the background presented in Aalbers et
al. (2020). We organize the simulated data into recoil-kinetic-energy bins
with 10 keV width, consistent with the recoil-kinetic-energy resolution quoted
in Aalbers et al. (2020), starting at 1 keV. We use a simple $\chi^{2}$-test
in order to address questions associated to the sensitivity of DARWIN to
different parameters and in order to combine simulated DARWIN data with those
from other experiments.

##  III Testing The Three-Massive-Neutrinos Paradigm

In the absence of more new physics, existing data reveals that the neutrino
weak-interaction-eigenstates $\nu_{\alpha}$, $\alpha=e,\mu,\tau$, are linear
combinations of the neutrino mass-eigenstates $\nu_{i}$ with mass $m_{i}$,
$i=1,2,3$:

$\nu_{\alpha}=U_{\alpha i}\nu_{i},$ (III.1)

where the $U_{\alpha i}$, $\alpha=e,\mu,\tau$, $i=1,2,3$, define the elements
of a unitary matrix. Here, we are only interested in solar neutrinos so all
accessible observables are sensitive to $|U_{ei}|^{2}$, $i=1,2,3$. These, in
turn, are parameterized with two mixing angles, $\theta_{12}$ and
$\theta_{13}$. Following the parameterization of the Particle Data Group Zyla
et al. (2020),

$|U_{e2}|^{2}=\sin^{2}\theta_{12}\cos^{2}\theta_{13},~{}~{}~{}~{}~{}|U_{e3}|^{2}=\sin^{2}\theta_{13},$
(III.2)

and unitarity uniquely determines the third matrix-element-squared:
$|U_{e1}|^{2}=1-|U_{e2}|^{2}-|U_{e3}|^{2}$. Combined fits to the existing data
reveal that the two independent mass-squared differences are $\Delta
m^{2}_{21}\equiv m_{2}^{2}-m_{1}^{2}\sim 10^{-4}$ eV2 and $|\Delta
m^{2}_{31}|\equiv m_{3}^{2}-m_{1}^{2}\sim 10^{-3}$ eV2. For more precise
values see, for example, Esteban et al. (2020).***See also http://www.nu-
fit.org. While $\Delta m^{2}_{21}$ is defined to be positive, the sign of
$\Delta m^{2}_{31}$ is still unknown; for our purposes here, it turns out,
this is irrelevant. The two mixing parameters of interest have been measured
quite precisely. According to Esteban et al. (2020), at the one-sigma level,

$\sin^{2}\theta_{12}=0.304^{+0.013}_{-0.012},~{}~{}~{}~{}\sin^{2}\theta_{13}=0.02221^{+0.00068}_{-0.00062}.$
(III.3)

The experiments that contribute most to these two measurements are
qualitatively different. $\theta_{12}$ is best constrained by solar neutrino
experiments – and is often referred to as the “solar angle” – while
$\theta_{13}$ is best constrained by reactor antineutrino experiments – and is
often referred to as the “reactor angle.”

We are interested in the solar $pp$-neutrinos. These have a continuous energy
spectrum that peaks around 300 keV and terminates at around 420 keV. The
matter-potential $V=\sqrt{2}G_{F}N_{e}$, where $G_{F}$ is the Fermi constant
and $N_{e}$ is the electron number-density, inside the Sun is
$V_{\odot}<2\times 10^{-5}$ (eV2/MeV) so, for neutrino energies $E<0.420$ MeV,
$|\Delta m^{2}_{21}|/(2E),|\Delta m^{2}_{31}|/(2E)\gg V_{\odot}$.†††For
example, at the center of the Sun, for neutrino energies less than 420 keV,
the “matter equivalent” of $\sin^{2}2\theta_{12}$ differs from its vacuum
counterpart by less than one percent. This, in turn, implies that, given what
is known about the neutrino mass-squared differences, matter effects can be
neglected. Including the fact that, for all practical purposes, solar
neutrinos lose flavor coherence as they find their way from the Sun to the
Earth, it is trivial to show that the $\nu_{e}$ survival probability is energy
independent and given by

$P_{ee}=|U_{e1}|^{4}+|U_{e2}|^{4}+|U_{e3}|^{4}.$ (III.4)

On the other hand, solar neutrino experiments cannot distinguish $\nu_{\mu}$
from $\nu_{\tau}$ – the neutrino energies are too small – but are potentially
sensitive to the combination $P_{ea}\equiv P_{e\mu}+P_{e\tau}$. In the three-
massive-neutrinos paradigm

$P_{ea}=1-P_{ee}.$ (III.5)

Given our current knowledge of mixing parameters, for $pp$-neutrinos, we can
indirectly infer that $P_{ee}=0.552\pm 0.025$, naively combining the
uncertainties in Eq. (III.3) in quadrature.

According to Aalbers et al. (2020), after 20 ton-years of exposure, DARWIN can
measure $P_{ee}$ with better than 1% accuracy. Assuming the three-massive-
neutrinos paradigm, this can be converted into a measurement of the relevant
mixing parameters. Fig. 1(top,left) depicts the allowed region of the
$\sin^{2}\theta_{12}\times\sin^{2}\theta_{13}$ parameter space assuming DARWIN
can measure $P_{ee}$ for $pp$-neutrinos at the 1% level, and assuming the
best-fit value is $P_{ee}=0.552$. There is very strong degeneracy between
different values of $\sin^{2}\theta_{12}$ and $\sin^{2}\theta_{13}$, for
obvious reasons. The degeneracies can be lifted by including constraints from
other neutrino experiments.

Figure 1: Top: One-, two- and three-sigma allowed regions of the
$\sin^{2}\theta_{12}\times\sin^{2}\theta_{13}$ parameter space assuming DARWIN
can measure $P_{ee}=0.552$ at the one percent level, excluding (left) and
including (right) external constraints on the neutrino-mixing parameters from
other solar experiments. The tiny empty ellipse in the right-hand panel
indicates the best-fit point. The open regions bound by dashed lines (left-
hand panel) represent one-, two- and three-sigma results from current 8B
neutrino experiments, as discussed in the text. Bottom: Solar-only $\chi^{2}$
as a function of $\sin^{2}\theta_{13}$, marginalized over
$\sin^{2}\theta_{12}$, assuming 300-ton-years of simulated DARWIN data. The
full line corresponds to the assumption of a depleted background, while the
dashed line is obtained including the expected natural ${}^{136}{\rm Xe}$
background.

It is interesting to investigate how well one can constrain neutrino-mixing
parameters using only solar-neutrino data. In order to estimate that, we add
to the hypothetical $pp$-neutrino measurement from DARWIN current information
from 8B neutrinos, mostly from the Super-Kamiokande and SNO experiments, see
Aharmim et al. (2013); Abe et al. (2016) and references therein. These provide
the strongest constraints on $\sin^{2}\theta_{12}$. Here, we address this in a
simplified but accurate way Nunokawa et al. (2006), postulating that 8B
experiments measure

$(P_{ee})_{{}^{8}\rm
B,~{}average}=(1-|U_{e3}|^{2})\left[0.9|U_{e2}|^{2}+0.1|U_{e1}|^{2}\right]+|U_{e3}|^{4},$
(III.6)

with 4% accuracy, consistent with the current uncertainty on
$\sin^{2}\theta_{12}$, mostly constrained by high-energy solar neutrino data.
$(1-|U_{e3}|^{2})\times 0.9$ (or $(1-|U_{e3}|^{2})\times 0.1$) is the average
probability that a 8B neutrino arrives at the surface of the Earth as a
$\nu_{2}$ (or a $\nu_{1}$). When the 8B data is treated as outlined above, it
translates into the open regions bound by dashed lines in Fig. 1(top,left).
Strong matter effects lead to the boomerang-shaped allowed region of the
parameter space and restrict the parameter space to values of
$\sin^{2}\theta_{12}\lesssim 0.5$. The results of the joint $pp-^{8}$B
analysis are depicted in Fig. 1(top,right). All degeneracies present in the
$pp$-neutrino data are lifted and one is constrained to small values of
$\sin^{2}\theta_{13}$ and $\sin^{2}\theta_{12}<0.5$.

The combined 8B and DARWIN data can rule out $\sin^{2}\theta_{13}=0$ with some
precision. This is important; it implies that a hypothetical DARWIN
measurement of the $pp$-neutrino flux, combined with the current 8B solar
neutrino data, can measure $\sin^{2}\theta_{13}$ in a way that is independent
from all non-solar measurements. The marginalized $\chi^{2}$ as a function of
$\sin^{2}\theta_{13}$ is depicted in Fig. 1(bottom), for 300 ton-years of
simulated DARWIN data and the current 8B sollar neutrino data. On average, if
the $pp$-neutrino flux can be measured at the percent level, we expect to
measure $\sin^{2}\theta_{13}$ at the 35% level and rule out
$\sin^{2}\theta_{13}=0$ at almost the three-sigma level. Here we consider the
two scenarios outlined earlier, one with natural xenon (dashed line), the
other with 136Xe-depleted xenon (solid line).

The precision on $\sin^{2}\theta_{13}$ obtained above is not comparable to
that of the current measurement of $\sin^{2}\theta_{13}$, Eq. (III.3).
However, these measurements are qualitatively different. The most precise
measurements of $\sin^{2}\theta_{13}$ come from reactor antineutrino
experiments and a baseline of order 1 km Bak et al. (2018); Adey et al.
(2018); de Kerret et al. (2020). The estimate discussed above is a “solar
only” measurement, i.e., it exclusively makes use of measurements of neutrinos
(and not antineutrinos) produced in the Sun. Current measurements of
$\sin^{2}\theta_{13}$ that make use of neutrinos (as opposed to
antineutrinos), from T2K and NOvA, are much less precise (at the 50%, see
Adamson et al. (2016); Abe et al. (2017)). Looking further into the future,
the DUNE experiment, for example, is expected to independently measure the
“neutrino-only” value of $\sin^{2}\theta_{13}$ at the 20% level de Gouvêa and
Kelly (2017) (or worse, depending on the assumptions made in the analysis).

##  IV Beyond the Three-Massive-Neutrinos Paradigm

The fact that the $pp$-neutrino flux can be computed with great precision,
combined with the sub-MeV $pp$-neutrino energies, allows a high-statistics
measurement of the $pp$-neutrino flux to meaningfully search for phenomena
beyond the thee-massive-neutrinos paradigm. Here we concentrate on testing the
hypothesis that the neutrinos produced in the Sun have a nonzero probability
of behaving as “sterile neutrinos” $\nu_{s}$, characterized by their lack of
participation in charged-current and neutral-current weak interactions.

We first discuss, in Sec. IV.1, the case where the oscillation probabilities
are energy-independent for the energies of interest, as in the case of the
thee-massive-neutrinos paradigm discussed in Sec. III. In particular, we test
the hypothesis that $P_{ee}+P_{ea}=1$ for $pp$-neutrinos. Then, in Sec. IV.2,
we compute DARWIN’s ability to constrain the hypothesis that there is a fourth
neutrino $\nu_{4}$ and that its mass is quasi-degenerate with $m_{1}$.

###  IV.1 Model-Independent Considerations

As discussed in Sec. II, we are interested in the shape and normalization of
the electron recoil-energy spectrum from neutrino–electron elastic scattering.
The differential cross-section for $\nu_{e}$ and $\nu_{a}$ scattering are
different, both in normalization and shape and hence, in principle, one can
obtain independent information on both $P_{ee}$ and $P_{ea}$.

We simulate and analyze 300 ton-years of DARWIN $pp$-data, as discussed in
Sec. II, and attempt to measure $P_{ee}$ and $P_{ea}$ independently. The
results are depicted in Fig. 2(left) for both the natural xenon (dashed) and
the 136Xe-depleted (solid) hypotheses. Strong departures from
$P_{ee}+P_{ea}=1$ are allowed and the “natural” data are not capable of ruling
out $P_{ea}=0$ at the three-sigma confidence level. The “depleted” data can
rule out $P_{ea}=0$ at the five-sigma confidence level. For both scenarios,
one can constrain the departure of $P_{ee}+P_{ea}$ from one, which we
interpret as the oscillation probability into sterile neutrinos $P_{es}\equiv
1-P_{ee}-P_{ea}$. The colorful diagonal lines in Fig. 2(left) correspond to
different constant values of $P_{es}$. Fig. 2(right) depicts $\chi^{2}$ as a
function of $P_{es}$, marginalized over $P_{ee}$ and restricing $P_{es}$ to
non-negative values for both scenarios. If DARWIN data are consistent with the
three-active-neutrinos paradigm, they will be capable of constraining
$P_{es}<0.35$ at the two-sigma confidence level even if DARWIN is filled with
natural xenon.

Figure 2: Left: One- and three-sigma allowed region of the $P_{ee}\times
P_{ea}$-plane, for 300 ton-years of simulated DARWIN data. The diagonal lines
correspond to constant $P_{es}\equiv 1-P_{ee}-P_{ea}$ values. The burgundy
line segment with positive slope corresponds to the values of
$(P_{ee},P_{ea})$ accessible via Eqs. (IV.9). Right: Marginalized $\chi^{2}$
as a function of $P_{es}$. The full line correspond to the assumption of a
depleted background, while the dashed line is obtained considering no cuts in
the ${}^{136}{\rm Xe}$ background.

###  IV.2 Fourth-Neutrino Hypothesis

We explore in more detail the scenario where there is one extra neutrino mass-
eigenstate $\nu_{4}$ with mass $m_{4}$. In this case, the interaction
eigenstates are, including $\nu_{s}$, related to the four mass-eigenstates via
a $4\times 4$ unitary matrix $U_{\alpha i}$, $\alpha=e,\mu,\tau,s$,
$i=1,2,3,4$. We will concentrate on the scenario where, among the four
$U_{si},$ only $U_{s1}$ and $U_{s4}$ are potentially nonzero.‡‡‡It is easy to
generalize this analysis assuming that only one of the $U_{sj}$, $j=1,2,3$,
and $U_{s4}$ are potentially nonzero. We spell this out in the Appendix. In
this case, we can parameterize the $|U_{ei}|^{2}$ entries of the mixing matrix
using three mixing angles $\theta_{12},\theta_{13},\theta_{14}$. Eqs. (III.2)
are still valid, along with

$|U_{e1}|^{2}=\cos^{2}\theta_{12}\cos^{2}\theta_{13}\cos^{2}\theta_{14},~{}~{}~{}~{}~{}|U_{e4}|^{2}=\cos^{2}\theta_{12}\cos^{2}\theta_{13}\sin^{2}\theta_{14}.$
(IV.1)

It is easy to check that $\sum_{i=1}^{4}|U_{ei}|^{2}=1$. The non-zero
“sterile” entries of the mixing matrix are

$|U_{s1}|^{2}=\sin^{2}\theta_{14},~{}~{}~{}~{}~{}|U_{s4}|^{2}=\cos^{2}\theta_{14}.$
(IV.2)

Given the quasi-two-flavors nature of these solar neutrino oscillations, to be
discussed momentarily, the entire physical parameter space is spanned by
either fixing $\Delta m^{2}_{41}>0$ and allowing $\sin^{2}\theta_{14}\in[0,1]$
or allowing both signs for $\Delta m^{2}_{41}$ and restricting
$\sin^{2}\theta_{14}\in[0,0.5]$ in such a way that $\nu_{4}$ is always “mostly
sterile.” Here, the former convention – to fix the sign of $\Delta
m^{2}_{41}>0$ – is most convenient. With this choice, when
$\sin^{2}\theta_{14}\in[0,0.5]$, the heaviest of the two quasi-degenerate
states (i.e., $\nu_{4}$) is mostly sterile, when
$\sin^{2}\theta_{14}\in[0.5,1]$, the lightest among the two quasi-degenerate
states (i.e., $\nu_{1}$) is mostly sterile. For historical reasons, we will
refer to $\sin^{2}\theta_{14}\in[0,0.5]$ as the light side of the parameter
space and $\sin^{2}\theta_{14}\in[0.5,1]$ as the dark side de Gouvêa et al.
(2000).

We are interested in the hypothesis that $\Delta m^{2}_{41}\ll\Delta
m^{2}_{21}$ and outside the reach of all current neutrino experiments. In this
case, the current neutrino oscillation data constrain the oscillation
parameters $\Delta m^{2}_{21}$, $\Delta m^{2}_{31}$, $\sin^{2}\theta_{12}$,
and $\sin^{2}\theta_{13}$ exactly as in the three-massive-neutrinos paradigm.
Furthermore, builiding on the discussion in Sec. III, it is easy to conclude
that the oscillation probabilities of interest $P_{e\alpha}$, $\alpha=e,a,s$,
are only functions of
$\sin^{2}\theta_{12},\sin^{2}\theta_{13},\sin^{2}\theta_{14}$, and $\Delta
m^{2}_{41}$. Further taking advantage of the fact that $|\Delta
m^{2}_{21}|/(2E),|\Delta m^{2}_{31}|/(2E)\gg V_{\odot}$, it is straightforward
to compute

$\displaystyle P_{ee}$ $\displaystyle=$
$\displaystyle|U_{e2}|^{4}+|U_{e3}|^{4}+(1-|U_{e2}|^{2}-|U_{e3}|^{2})^{2}P_{ee}^{2f}(\Delta
m^{2}_{41},\sin^{2}\theta_{14},V^{\rm eff}),$ (IV.3) $\displaystyle P_{es}$
$\displaystyle=$
$\displaystyle(1-|U_{e2}|^{2}-|U_{e3}|^{2})(1-P_{ee}^{2f}(\Delta
m^{2}_{41},\sin^{2}\theta_{14},V^{\rm eff})),$ (IV.4) $\displaystyle P_{ea}$
$\displaystyle=$ $\displaystyle 1-P_{ee}-P_{es},$ (IV.5)

where $P_{ee}^{2f}$ is the survival probability obtained in the scenario where
there are only two flavors, $\nu_{e}^{2f}$ and $\nu_{s}^{2f}$, characterized
by the mass-squared difference $\Delta m^{2}_{41}$ and the mixing angle
$\theta_{14}$, defined via
$\nu_{e}^{2f}=\cos\theta_{14}\nu_{1}+\sin\theta_{14}\nu_{4}$. Inside
$P_{ee}^{2f}$, the matter potential is replaced by an effective matter
potential $V^{\rm eff}$. It takes into account the neutral-current
contribution to the matter potential $V_{NC}=-\sqrt{2}/2G_{F}N_{n}$, where
$N_{n}$ is the neutron number density in the medium, while the charged-current
contribution is rescaled by
$(1-|U_{e2}^{2}-|U_{e3}|^{2})=\cos^{2}\theta_{13}\cos^{2}\theta_{12}$:

$V^{\rm
eff}=\sqrt{2}G_{F}\left(N_{e}\cos^{2}\theta_{13}\cos^{2}\theta_{12}-\frac{1}{2}N_{n}\right).$
(IV.6)

In the sun, the position-dependency of the electron and neutron number
densities are slightly different Bahcall et al. (2001).§§§There are relatively
more neutrons in the center of the sun relative to its edges. This is due to
the fact that most of the solar helium is concentrated in the core. In the
sun’s core, $N_{n}$ is around 50% of $N_{e}$ and $V^{\rm eff}$ is slightly
less than one half of the standard matter potential in the Sun.

Eq. (IV.3) allows us to estimate, in very general terms, the impact of the
sterile neutrinos. For $P^{2f}_{ee}=1$, we recover the three-active-neutrinos
result,
$P_{ee}=\cos^{4}\theta_{12}\cos^{4}\theta_{13}+\sin^{4}\theta_{12}\cos^{4}\theta_{13}+\sin^{4}\theta_{13}$,
Eq. (III.4). On the other hand, for $P^{2f}_{ee}=0$,
$P_{ee}=\sin^{4}\theta_{12}\cos^{4}\theta_{13}+\sin^{4}\theta_{13}$ such that,
given the current knowledge of oscillation parameters,

$P_{ee}\in[0.09,0.55].$ (IV.7)

$P_{es}$ values, on the other hand, are allowed to be as small as zero and as
large as 0.68.

In Sec. IV.1, we discussed that, very generically, DARWIN can rule out
$P_{es}<0.35$ at the two-sigma level. The situation here is more constrained
as $P_{ee},P_{ea},P_{es}$ are not only required to add up to one but depend on
the same oscillation parameters. We proceed to discuss the sensitivity of
DARWIN to the new oscillation parameters $\Delta
m^{2}_{41},\sin^{2}\theta_{14}$ by taking advantage of the fact that the
properties of $P_{ee}^{2f}$ are well known (see, for example, Giunti and Wook
(2007)).

For large-enough values of $\Delta m^{2}_{41}$, $P_{ee}^{2f}$ is well
approximated by averaged-out vacuum oscillations:

$P_{ee}^{2f,\rm ave}=1-\frac{1}{2}\sin^{2}2\theta_{14}.$ (IV.8)

This occurs, keeping in mind we are interested in energies below 420 keV, for
$\Delta m^{2}_{41}\gtrsim 10^{-5}$ eV2, when the solar matter effects can be
ignored. In this case,

$\displaystyle P_{ee}$ $\displaystyle=$ $\displaystyle
0.55-0.23\sin^{2}2\theta_{14},$ $\displaystyle P_{es}$ $\displaystyle=$
$\displaystyle 0.34\sin^{2}2\theta_{14},$ (IV.9) $\displaystyle P_{ea}$
$\displaystyle=$ $\displaystyle 0.45-0.11\sin^{2}2\theta_{14}.$

Here, it is impossible to distinguish the light from the dark side of the
parameter space since the oscillation probabilities are invariant under
$\sin^{2}\theta_{14}\leftrightarrow 1-\sin^{2}\theta_{14}$. Varying
$\sin^{2}2\theta_{14}\in[0,1]$, Eqs. (IV.9) define a line segment in the
$P_{ee}\times P_{ea}$-plane, depicted in Fig. 2(left) – burgundy line with
positive slope – keeping in mind the segment extends to $P_{ee}$ values below
0.4. Fig. 3 depicts $\chi^{2}$ as a function of $\sin^{2}\theta_{14}$ in the
regime where Eqs. (IV.9) are a good approximation, for 300 ton-years of
simulated DARWIN data, for both the natural (dashed) and 136Xe-depleted
(solid) scenarios. In this analysis, and in the upcoming analyses discussed
this subsection, we assume that $\sin^{2}\theta_{13}$ and
$\sin^{2}\theta_{12}$ are known with infinite precision. This is, currently, a
good approximation for $\sin^{2}\theta_{13}$ and will be a good approximation
for $\sin^{2}\theta_{12}$ once data from the JUNO experiment is analyzed An et
al. (2016). Similar results were recently presented and discussed in Goldhagen
et al. (2021). Where are assumptions agree, the estimated sensitivity also
agrees.

Figure 3: $\chi^{2}$ as a function of $\sin^{2}\theta_{14}$ for 300 ton-years
of simulated DARWIN data, in the regime when matter effects are not
significant and the new oscillatory effects driven by $\Delta m^{2}_{41}$
average out. The full line corresponds to the assumption of a depleted
background, while the dashed line is obtained considering no cuts in the
${}^{136}{\rm Xe}$ background.

For intermediate values of $\Delta m^{2}_{41}$, $P_{ee}^{2f}$ is well describe
by the strong MSW effect in the adiabatic regime. In this case, for a range of
energies,

$P_{ee}^{2f,\rm adiabatic}=\sin^{2}\theta_{14}.$ (IV.10)

For $pp$-neutrinos, this occurs for, very roughly, $\sin^{2}\theta_{14}\gtrsim
10^{-3}$ and $10^{-9}\lesssim\Delta m^{2}_{41}/(\rm eV^{2})\lesssim 10^{-6}$.
Under these conditions,

$\displaystyle P_{ee}$ $\displaystyle=$ $\displaystyle
0.09+0.46\sin^{2}\theta_{14},$ $\displaystyle P_{es}$ $\displaystyle=$
$\displaystyle 0.68-0.68\sin^{2}\theta_{14},$ (IV.11) $\displaystyle P_{ea}$
$\displaystyle=$ $\displaystyle 0.27+0.22\sin^{2}\theta_{14}.$

Here, oscillation probabilities are very different in the light and dark
sides. In particular, in the light side of the parameter space $P_{ee}$
($P_{es}$) is small (large) and increases (decreases) linearly with
$\sin^{2}\theta_{14}$. If DARWIN data are consistent with three-active
neutrinos, in this region of parameter space, small values of
$\sin^{2}\theta_{14}$ will be excluded while large values of
$\sin^{2}\theta_{14}$ are allowed.

For small-enough values of $\Delta m^{2}_{41}$, $P_{ee}^{2f}$ is well
described by the strong MSW effect in the very non-adiabatic regime and turns
out to be well approximated by vacuum oscillations,

$P_{ee}^{2f,\rm ave}=1-\sin^{2}2\theta_{14}\sin^{2}\left(\frac{\Delta
m^{2}_{41}L}{4E}\right).$ (IV.12)

This occurs, for $pp$-neutrinos, for $\Delta m^{2}_{41}\lesssim 10^{-9}$ eV2.
In this case,

$\displaystyle P_{ee}$ $\displaystyle=$ $\displaystyle
0.55-0.46\sin^{2}2\theta_{14}\sin^{2}\left(\frac{\Delta
m^{2}_{41}L}{4E}\right),$ $\displaystyle P_{es}$ $\displaystyle=$
$\displaystyle 0.68\sin^{2}2\theta_{14}\sin^{2}\left(\frac{\Delta
m^{2}_{41}L}{4E}\right),$ (IV.13) $\displaystyle P_{ea}$ $\displaystyle=$
$\displaystyle 0.45-0.22\sin^{2}2\theta_{14}\sin^{2}\left(\frac{\Delta
m^{2}_{41}L}{4E}\right).$

Here, again, it is impossible to distinguish the light from the dark side of
the parameter space. Given the average Earth–Sun distance $L=1.5\times
10^{12}$ m, the oscillation phase is

$\frac{\Delta m^{2}_{41}L}{4E}=4.8\left(\frac{\Delta
m^{2}_{41}}{10^{-12}~{}\rm eV^{2}}\right)\left(\frac{400~{}\rm
keV}{E}\right),$ (IV.14)

so we expect the vacuum oscillations to average out for $\Delta
m^{2}_{41}\gtrsim 10^{-11}$ eV2. This means that, for $10^{-11}\lesssim\Delta
m^{2}_{41}/(\rm eV^{2})\lesssim 10^{-9}$, the oscillation probabilities are
well described by Eqs. (IV.9).

Fig. 5 depicts contours of constant $P_{ee}$ in the $\Delta
m^{2}_{41}\times\sin^{2}\theta_{14}$-plane for $E_{\nu}=300$ keV. The other
parameters are fixed to their current best-fit values, Eq. (III.3). We assume
the matter potential is spherically symmetric and drops exponentially, $V^{\rm
eff}\propto e^{-r/r_{0}^{s}}$. We fit information from the prediction of the
B16-GS98 solar model Vinyoles et al. (2017) and obtain
$r_{0}^{s}=R_{\odot}/10.37$ where $R_{\odot}=6.96\times 10^{11}$ m is the
average radius of the Sun; see Fig. 4 for a comparison of the matter potential
in the standard case (left) and in the scenario of interest here (right,
labeled sterile neutrino).

Figure 4: Solar matter potential for active (left) and sterile (right)
neutrinos – the scenario of interest here – as function of the distance from
the center in units of the Solar radius, from the B16-GS98 Solar Model
Vinyoles et al. (2017). We also present our fitted exponential forms, where
$r_{0}^{e}=R_{\odot}/10.43$, and $r_{0}^{s}=R_{\odot}/10.37$, in dashed lines.

Under these circumstances, $P_{ee}^{2f}$ can be computed exactly Petcov
(1988). For simplified pedagogical discussions see, for example, de Gouvêa
(2004); Giunti and Wook (2007). We assume all solar neutrinos are produced in
the exact center of the Sun; we explicitly verified that the results we get
are very similar to the results we would have obtained by integrating over the
region where $pp$-neutrinos are produced. The region where matter effects are
strong and the adiabatic condition holds correspond to the vertical sides of
the constant $P_{ee}$ regions that form quasi-triangles. The “return” to
vacuum oscillations at low and high values of $\Delta m^{2}_{41}$ is
highlighted by the vertical, dark lines, which correspond to constant values
of the averaged-out vacuum oscillation probability. For a detailed discussion
of the boundary between the adiabatic and non-adiabatic transition, including
$L$ dependent effects, see Friedland (2001).

Figure 5: Contours of constant $P_{ee}$ in the $\Delta
m^{2}_{41}\times\sin^{2}\theta_{14}$-plane for $E=300$ keV.
$\sin^{2}\theta_{12},\sin^{2}\theta_{13}$ are fixed to their best-fit values,
Eq. (III.3). The vertical lines correspond to constant values of the averaged-
out vacuum oscillation probability.

We simulate 300 ton-years of DARWIN data consistent with the three-massive-
neutrinos paradigm and assuming the true values of $\sin^{2}\theta_{12}$ and
$\sin^{2}\theta_{13}$ are the ones in Eq. (III.3). We restrict our discussion
to values of $\Delta m^{2}_{41}<10^{-6}$ eV2. Larger values are constrained by
measurements of higher-energy solar neutrinos; these constraints have been
explored in de Holanda and Smirnov (2004, 2011), along with a detailed
discussion of the oscillation probabilities. The expressions we derive here
are contained in the analyses of de Holanda and Smirnov (2004, 2011) if one
explores them in the appropriate regime.

Fig. 6 depicts the region of the $\tan^{2}\theta\times\Delta m^{2}_{41}$–plane
inside of which 300 ton-years of DARWIN data is sensitive, at the 90%
confidence level, to the fourth neutrino for both the natural-xenon scenario
(dashed line) and the depleted-136Xe scenario (solid). On a log-scale, the
contour is symmetric relative to $\tan^{2}\theta_{14}=1$ when one cannot
distinguish the light from the dark side of the parameter space de Gouvêa et
al. (2000), a feature one readily observes, as advertised, for small values of
$\Delta m^{2}_{41}$. The impact of nontrivial matter effects is also readily
observable. For larger values of the $\Delta m^{2}_{41}$, the sensitivity to
small mixing angles is expected to “shut-off” quickly – see Fig. 5 – and would
return to values similar to those around $\Delta m^{2}_{41}\sim 10^{-10}$ eV2,
minus the tiny wiggles.

Figure 6: Two- and three-sigma sensitivity of 300 ton-years of DARWIN data to
a light sterile neutrino in the $\Delta
m^{2}_{41}\times\tan^{2}\theta_{14}$–plane. The dashed lines correspond to the
assumption of a full 136Xe background, while the full lines are for the
depleted case.

The low energies of the $pp$-neutrinos combined with the long Earth–Sun
distance render DARWIN a specially powerful probe of the hypothesis that
neutrinos are pseudo-Dirac fermions. This is the hypothesis that there are
right-handed neutrinos coupled to the left-handed lepton doublets and the
Higgs doublet via a tiny Yukawa coupling $y$ and that lepton number is only
slightly violated. In these scenarios, each of the neutrino mass eigenstates
is “split” into two quasi-degenerate Majorana fermions, each a 50–50 mixture
of an active neutrino (from the lepton doublet) and a sterile neutrino (the
right-handed neutrino). The mass splitting is small enough that, for most
applications, the two quasi-degenerate state act as one Dirac fermion. Pseudo-
Dirac neutrinos reveal themselves via active–sterile oscillations associated
with very large mixing and very small mass-squared differences.

In the language introduced here, a pseudo-Dirac neutrino corresponds to
$\sin^{2}2\theta_{14}=1$ (maximal mixing) and the small mass splitting leads
to a nonzero $\Delta m^{2}_{41}=4\epsilon m_{1}$ where $m_{1}\pm\epsilon$ are
the masses of the two quasi-degenerate states (here $\nu_{1}$ and $\nu_{4}$),
$m_{1}$ is the Dirac mass, proportional to the neutrino Yukawa coupling and
$\epsilon$ characterizes the strength of the lepton-number violating physics.
Fig. 7 depicts $\chi^{2}$ as a function of $\Delta m^{2}_{41}$ for
$\sin^{2}2\theta_{14}=1$ associated with 300 ton-years of simulated DARWIN
data for both the natural xenon (dashed) and the 136Xe-depleted (solid)
scenarios, assuming the data are consistent with no new neutrino states.
Current solar neutrino data exclude $\Delta m^{2}_{41}$ values larger than
$10^{-12}$ eV2 de Gouvêa et al. (2009); Donini et al. (2011) so DARWIN can
extend the sensitivity to $\Delta m^{2}_{41}$ – and $\epsilon$ – by an order
of magnitude.

Figure 7: Darwin sensitivity to the quadratic mass difference $\Delta
m^{2}_{41}$ in the case of maximal mixing, $\sin^{2}\theta_{14}=0.5$. The full
line correspond to the assumption of a depleted background, while the dashed
line is obtained considering no cuts in the ${}^{136}{\rm Xe}$ background.

##  V Concluding Remarks

Next-generation WIMP-dark-matter-search experiments will be exposed to a
large-enough flux of solar neutrinos that neutrino-mediated events are
unavoidable. The DARWIN collaboration recently argued that DARWIN can collect
a large, useful sample of solar $pp$-neutrinos, detected via elastic
neutrino–electron scattering Aalbers et al. (2020), and measure the survival
probability of $pp$-neutrinos with sub-percent precision. Here we explored the
physics potential of such a sample in more detail, addressing other concrete
neutrino-physics questions and exploring whether one can also extract
information from a precise measurement of the shape of the differential
$pp$-neutrino flux.

We estimate that, with 300 ton-years of data, DARWIN can not only measure the
survival probability of $pp$-neutrinos with sub-percent precision but also
determine, with the help of current solar neutrino data, the value of
$\sin^{2}\theta_{13}$, with the potential to exclude $\sin^{2}\theta_{13}=0$
close to the three sigma level. Such a $pp$-neutrino sample would allow one to
perform a “neutrinos-only” (and solar-neutrinos-only) measurement of
$\sin^{2}\theta_{13}$ and $\sin^{2}\theta_{12}$. Such a measurement can be
compared with, for example, reactor-based “antineutrinos-only” measurements of
the same mixing parameters and allow for nontrivial tests of the CPT-theorem
and other new physics scenarios.

DARWIN can also test the hypothesis that $pp$-neutrinos are oscillating into a
combination of active and sterile neutrinos. We estimate that DARWIN data can
exclude the hypothesis that the $pp$-neutrinos are “disappearing” in an energy
independent way – assuming their data are consistent with the three-active-
neutrinos paradigm – especially if the experiment manages to fill their
detector with 136Xe-depleted xenon.

We explored in some detail how well DARWIN can constrain the existence of a
new neutrino mass-eigenstate $\nu_{4}$ (mass $m_{4}$) that is quasi-degenerate
and mixes with $\nu_{1}$, i.e, $\Delta m^{2}_{41}\ll\Delta m^{2}_{21}$,
$U_{s1},U_{s4}\neq 0$, $U_{s2}=U_{s3}=0$. Our estimated sensitivity is
depicted in Fig. 6. It supersedes that of all current and near-future searches
for new, very light neutrinos. In particular, DARWIN can test the hypothesis
that $\nu_{1}$ is a pseudo-Dirac fermion as long as the induced mass-squared
difference is larger than $10^{-13}$ eV2. This is one order of magnitude more
sensitive than existing constraints de Gouvêa et al. (2009); Donini et al.
(2011).

Throughout, we allowed for the hypotheses that DARWIN is filled with natural
xenon or 136Xe-depleted xenon. We find that while the sensitivity of the
experiment with natural xenon is outstanding, a 136Xe-depleted setup is
significantly more sensitive when it comes to the measurements and searches
discussed here. In our discussions, we did not include time-dependent effects
(the seasonal and day-night effects). These can impact the sensitivity to new
neutrino states within a subset of the parameter space explored here. They
would not significantly modify the results discussed here but provide extra
handles on the new physics.

##  Appendix A Other Fourth-Neutrino Scenarios

We restricted our fourth-neutrino analyses to one new neutrino mass-eigenstate
$\nu_{4}$ and allowed for only “1–4” sterile mixing, i.e.,
$U_{s1}=\sin\theta_{14}$, $U_{s4}=\cos\theta_{14}$ while $U_{s2}=U_{s3}=0$.
Here we discuss some simple generalizations.

The case of one new neutrino mass-eigenstate $\nu_{5}$ that is quasi-
degenerate with $\nu_{2}$ and only “2–5” sterile mixing would also be
parameterized by a mass-squared difference $\Delta m^{2}_{52}$ (positive-
definite), assumed to be much smaller than $\Delta m^{2}_{21}$, $\Delta
m^{2}_{51}$ and $|\Delta m^{2}_{53}|$, and one mixing angle $\theta_{25}$:

$\displaystyle U_{e1}^{2}$
$\displaystyle=\cos^{2}\theta_{12}\cos^{2}\theta_{13},$ $\displaystyle
U_{e2}^{2}$
$\displaystyle=\sin^{2}\theta_{12}\cos^{2}\theta_{13}\cos^{2}\theta_{25},$
$\displaystyle U_{e3}^{2}$ $\displaystyle=\sin^{2}\theta_{13},$ $\displaystyle
U_{e5}^{2}$
$\displaystyle=\sin^{2}\theta_{12}\cos^{2}\theta_{13}\sin^{2}\theta_{25},$
(A.1) $\displaystyle U_{s1}^{2}$ $\displaystyle=0,$ $\displaystyle U_{s2}^{2}$
$\displaystyle=\sin^{2}\theta_{25},$ $\displaystyle U_{s3}^{2}$
$\displaystyle=0,$ $\displaystyle U_{s5}^{2}$
$\displaystyle=\cos^{2}\theta_{25}.$ (A.2)

Similar to Eq. (IV.3), and the equivalent expressions for $P_{ea}$ and
$P_{es}$, here

$\displaystyle P_{ee}$ $\displaystyle=$
$\displaystyle|U_{e1}|^{4}+|U_{e3}|^{4}+(1-|U_{e1}|^{2}-|U_{e3}|^{2})^{2}P_{ee}^{2f}(\Delta
m^{2}_{52},\sin^{2}\theta_{25},V_{25}^{\rm eff}),$ (A.3) $\displaystyle
P_{es}$ $\displaystyle=$
$\displaystyle(1-|U_{e1}|^{2}-|U_{e3}|^{2})(1-P_{ee}^{2f}(\Delta
m^{2}_{52},\sin^{2}\theta_{25},V_{25}^{\rm eff})),$ (A.4) $\displaystyle
P_{ea}$ $\displaystyle=$ $\displaystyle 1-P_{ee}-P_{es},$ (A.5)

where $P_{ee}^{2f}$ is the survival probability obtained in the scenario where
there are only two flavors, $\nu_{e}^{2f}$ and $\nu_{s}^{2f}$, characterized
by the mass-squared difference $\Delta m^{2}_{52}$ and the mixing angle
$\theta_{25}$, defined via
$\nu_{e}^{2f}=\cos\theta_{25}\nu_{2}+\sin\theta_{25}\nu_{5}$. Here, the
effective matter potential is

$V_{25}^{\rm
eff}=\sqrt{2}G_{F}\left(N_{e}\cos^{2}\theta_{13}\sin^{2}\theta_{12}-\frac{1}{2}N_{n}\right).$
(A.6)

On the other hand, the case of one new neutrino mass-eigenstate $\nu_{6}$ that
is quasi-degenerate with $\nu_{3}$ and only “3–6” sterile mixing would also be
parameterized by a mass-squared difference $\Delta m^{2}_{63}$ (positive-
definite), assumed to be much smaller than $\Delta m^{2}_{21}$, $|\Delta
m^{2}_{61}|$ and $|\Delta m^{2}_{62}|$, and one mixing angle $\theta_{36}$:

$\displaystyle U_{e1}^{2}$
$\displaystyle=\cos^{2}\theta_{12}\cos^{2}\theta_{13},$ $\displaystyle
U_{e2}^{2}$ $\displaystyle=\sin^{2}\theta_{12}\cos^{2}\theta_{13},$
$\displaystyle U_{e3}^{2}$
$\displaystyle=\sin^{2}\theta_{13}\cos^{2}\theta_{36},$ $\displaystyle
U_{e5}^{2}$ $\displaystyle=\sin^{2}\theta_{13}\sin^{2}\theta_{36},$ (A.7)
$\displaystyle U_{s1}^{2}$ $\displaystyle=0,$ $\displaystyle U_{s2}^{2}$
$\displaystyle=0,$ $\displaystyle U_{s3}^{2}$
$\displaystyle=\sin^{2}\theta_{36},$ $\displaystyle U_{s5}^{2}$
$\displaystyle=\cos^{2}\theta_{36}.$ (A.8)

Similar to Eq. (IV.3), and the equivalent expressions for $P_{ea}$ and
$P_{es}$, here

$\displaystyle P_{ee}$ $\displaystyle=$
$\displaystyle|U_{e1}|^{4}+|U_{e2}|^{4}+(1-|U_{e1}|^{2}-|U_{e2}|^{2})^{2}P_{ee}^{2f}(\Delta
m^{2}_{63},\sin^{2}\theta_{36},V_{36}^{\rm eff}),$ (A.9) $\displaystyle
P_{es}$ $\displaystyle=$
$\displaystyle(1-|U_{e1}|^{2}-|U_{e2}|^{2})(1-P_{ee}^{2f}(\Delta
m^{2}_{63},\sin^{2}\theta_{36},V_{36}^{\rm eff})),$ (A.10) $\displaystyle
P_{ea}$ $\displaystyle=$ $\displaystyle 1-P_{ee}-P_{es},$ (A.11)

where $P_{ee}^{2f}$ is the survival probability obtained in the scenario where
there are only two flavors, $\nu_{e}^{2f}$ and $\nu_{s}^{2f}$, characterized
by the mass-squared difference $\Delta m^{2}_{63}$ and the mixing angle
$\theta_{36}$, defined via
$\nu_{e}^{2f}=\cos\theta_{36}\nu_{3}+\sin\theta_{36}\nu_{6}$. Here, the
effective matter potential is

$V_{36}^{\rm
eff}=\sqrt{2}G_{F}\left(N_{e}\sin^{2}\theta_{13}-\frac{1}{2}N_{n}\right).$
(A.12)

Qualitatively, the three scenarios – 1–4, 2–5, 3–6 – are identical modulo
relabelings of the mixing parameters. Quantitatively, however, there are
significant differences. The coefficients of the $P_{ee}^{2f}$ term in Eqs.
(IV.3), (A.3), (A.9) are, respectively,
$(\cos^{2}\theta_{12}\cos^{2}\theta_{13})^{2}\sim 0.5$,
$(\sin^{2}\theta_{12}\cos^{2}\theta_{13})^{2}\sim 0.1$, and
$(\sin^{2}\theta_{13})^{2}\sim 0.0005$. These numbers define the maximum
deviation of $P_{ee}$ from expectations from the three-massive-neutrinos
paradigm, $P_{ee}\sim 0.55$. Hence, very generically, 1–4 effects can be very
strong, as discussed in the text, 2–5 effects are at most of order 20%, and
3–6 effects are at the permille level. On the other hand, the effective
potentials are also quantitatively very different. The charged-current
contribution to $V^{\rm eff}_{25}$ (Eq. (A.6)) is suppressed relative to the
neutral-current one by a factor $\sin^{2}\theta_{12}\cos^{2}\theta_{13}\sim
0.3$. Since $N_{n}/N_{e}$ varies between, roughly, 0.5 and less than 0.1
between the center of the Sun and its edge, $V^{\rm eff}_{25}$ is
significantly smaller than $V_{\odot}$, almost vanishing at the Sun’s core,
when the charged- and neutral-current contributions, accidentally, almost
cancel out one another. $V^{\rm eff}_{36}$ (Eq. (A.12)), instead, is solidly
dominated by the neutral-current matter potential since the charged-current
contribution is suppressed by $\sin^{2}\theta_{13}\sim 0.02$. Not only is it
smaller than $V_{\odot}$, it has the opposite sign, a fact that qualitatively
impact the behavior of $P_{ee}^{2f}$.

The scenario where all neutrinos are pseudo-Dirac fermions is equivalent to
the combination of the 1–4, 2–5, and 3–6 scenarios spelled out above (see, for
example, de Gouvêa et al. (2009)). Note that such a combination is straight
forward; the effects of the different contributions simply “add up” without
too much interference, as long as the new mass-squared differences are
“isolated enough,” i.e., the three new mass-squared differences $\Delta
m^{2}_{41}$, $\Delta m^{2}_{52}$, and $\Delta m^{2}_{63}$ are much smaller
than all other mass-squared differences. For example,

$\displaystyle P_{ee}$ $\displaystyle=$
$\displaystyle(\cos^{2}\theta_{12}\cos^{2}\theta_{13})^{2}P_{ee}^{2f}(\Delta
m^{2}_{41},\sin^{2}\theta_{14},V_{14}^{\rm
eff})+(\sin^{2}\theta_{12}\sin^{2}\theta_{13})^{2}P_{ee}^{2f}(\Delta
m^{2}_{52}\sin^{2}\theta_{25},V_{25}^{\rm eff})+$ (A.13) $\displaystyle+$
$\displaystyle(\sin^{2}\theta_{13})^{2}P_{ee}^{2f}(\Delta
m^{2}_{63},\sin^{2}\theta_{36},V_{36}^{\rm eff}),$

where $V_{14}^{\rm eff}$ is given by Eq. (IV.6).

## Acknowledgements

We thank Pedro Machado for discussions of potential uses of DARWIN data for
neutrino physics. This work was supported in part by the US Department of
Energy (DOE) grant #de-sc0010143 and in part by the NSF grant PHY-1630782. The
document was prepared using the resources of the Fermi National Accelerator
Laboratory (Fermilab), a DOE, Office of Science, HEP User Facility. Fermilab
is managed by Fermi Research Alliance, LLC (FRA), acting under Contract No.
DE-AC02-07CH11359. This material is based upon work supported by the NSF grant
AST-1757792, a Research Experience for Undergraduates grant awarded to the
Center for Interdisciplinary Exploration and Research in Astrophysics (CIERA)
at Northwestern University. IMS is supported by the Faculty of Arts and
Sciences of Harvard University.

## References

  * Aalbers et al. (2020) J. Aalbers et al. (DARWIN), “Solar neutrino detection sensitivity in DARWIN via electron scattering,” Eur. Phys. J. C 80, 1133 (2020), eprint 2006.03114.
  * Bahcall et al. (2001) J. N. Bahcall, M. H. Pinsonneault, and S. Basu, “Solar models: Current epoch and time dependences, neutrinos, and helioseismological properties,” Astrophys. J. 555, 990 (2001), eprint astro-ph/0010346.
  * Serenelli et al. (2011) A. M. Serenelli, W. C. Haxton, and C. Pena-Garay, “Solar models with accretion. I. Application to the solar abundance problem,” Astrophys. J. 743, 24 (2011), eprint 1104.1639.
  * Vinyoles et al. (2017) N. Vinyoles, A. M. Serenelli, F. L. Villante, S. Basu, J. Bergström, M. C. Gonzalez-Garcia, M. Maltoni, C. Peña Garay, and N. Song, “A new Generation of Standard Solar Models,” Astrophys. J. 835, 202 (2017), eprint 1611.09867.
  * Bellini et al. (2014) G. Bellini et al. (BOREXINO), “Neutrinos from the primary proton–proton fusion process in the Sun,” Nature 512, 383 (2014).
  * Wolfenstein (1981) L. Wolfenstein, “Different Varieties of Massive Dirac Neutrinos,” Nucl. Phys. B 186, 147 (1981).
  * Petcov (1982) S. T. Petcov, “On Pseudodirac Neutrinos, Neutrino Oscillations and Neutrinoless Double beta Decay,” Phys. Lett. B 110, 245 (1982).
  * Bilenky and Pontecorvo (1983) S. M. Bilenky and B. Pontecorvo, “Neutrino Oscillations With Large Oscillation Length in Spite of Large (Majorana) Neutrino Masses?,” Sov. J. Nucl. Phys. 38, 248 (1983).
  * de Gouvêa et al. (2009) A. de Gouvêa, W.-C. Huang, and J. Jenkins, “Pseudo-Dirac Neutrinos in the New Standard Model,” Phys. Rev. D 80, 073007 (2009), eprint 0906.1611.
  * Donini et al. (2011) A. Donini, P. Hernandez, J. Lopez-Pavon, and M. Maltoni, “Minimal models with light sterile neutrinos,” JHEP 07, 105 (2011), eprint 1106.0064.
  * de Gouvêa et al. (2020) A. de Gouvêa, I. Martinez-Soler, Y. F. Perez-Gonzalez, and M. Sen, “Fundamental physics with the diffuse supernova background neutrinos,” Phys. Rev. D 102, 123012 (2020), eprint 2007.13748.
  * Martinez-Soler et al. (2021) I. Martinez-Soler, Y. F. Perez-Gonzalez, and M. Sen, “SN1987A still shining: A Quest for Pseudo-Dirac Neutrinos,” (2021), eprint 2105.12736.
  * Aalbers et al. (2016) J. Aalbers et al. (DARWIN), “DARWIN: towards the ultimate dark matter detector,” JCAP 11, 017 (2016), eprint 1606.07001.
  * Zyla et al. (2020) P. A. Zyla et al. (Particle Data Group), “Review of Particle Physics,” PTEP 2020, 083C01 (2020).
  * Esteban et al. (2020) I. Esteban, M. C. Gonzalez-Garcia, M. Maltoni, T. Schwetz, and A. Zhou, “The fate of hints: updated global analysis of three-flavor neutrino oscillations,” JHEP 09, 178 (2020), eprint 2007.14792.
  * Aharmim et al. (2013) B. Aharmim et al. (SNO), “Combined Analysis of all Three Phases of Solar Neutrino Data from the Sudbury Neutrino Observatory,” Phys. Rev. C 88, 025501 (2013), eprint 1109.0763.
  * Abe et al. (2016) K. Abe et al. (Super-Kamiokande), “Solar Neutrino Measurements in Super-Kamiokande-IV,” Phys. Rev. D 94, 052010 (2016), eprint 1606.07538.
  * Nunokawa et al. (2006) H. Nunokawa, S. J. Parke, and R. Zukanovich Funchal, “What fraction of boron-8 solar neutrinos arrive at the earth as a nu(2) mass eigenstate?,” Phys. Rev. D 74, 013006 (2006), eprint hep-ph/0601198.
  * Bak et al. (2018) G. Bak et al. (RENO), “Measurement of Reactor Antineutrino Oscillation Amplitude and Frequency at RENO,” Phys. Rev. Lett. 121, 201801 (2018), eprint 1806.00248.
  * Adey et al. (2018) D. Adey et al. (Daya Bay), “Measurement of the Electron Antineutrino Oscillation with 1958 Days of Operation at Daya Bay,” Phys. Rev. Lett. 121, 241805 (2018), eprint 1809.02261.
  * de Kerret et al. (2020) H. de Kerret et al. (Double Chooz), “Double Chooz $\theta_{13}$ measurement via total neutron capture detection,” Nature Phys. 16, 558 (2020), eprint 1901.09445.
  * Adamson et al. (2016) P. Adamson et al. (NOvA), “First measurement of electron neutrino appearance in NOvA,” Phys. Rev. Lett. 116, 151806 (2016), eprint 1601.05022.
  * Abe et al. (2017) K. Abe et al. (T2K), “Measurement of neutrino and antineutrino oscillations by the T2K experiment including a new additional sample of $\nu_{e}$ interactions at the far detector,” Phys. Rev. D 96, 092006 (2017), [Erratum: Phys.Rev.D 98, 019902 (2018)], eprint 1707.01048.
  * de Gouvêa and Kelly (2017) A. de Gouvêa and K. J. Kelly, “Neutrino vs. Antineutrino Oscillation Parameters at DUNE and Hyper-Kamiokande,” Phys. Rev. D 96, 095018 (2017), eprint 1709.06090.
  * de Gouvêa et al. (2000) A. de Gouvêa, A. Friedland, and H. Murayama, “The Dark side of the solar neutrino parameter space,” Phys. Lett. B 490, 125 (2000), eprint hep-ph/0002064.
  * Giunti and Wook (2007) C. Giunti and K. C. Wook, _Fundamentals of Neutrino Physics and Astrophysics_ (Oxford Univ., Oxford, 2007), URL https://cds.cern.ch/record/1053706.
  * An et al. (2016) F. An et al. (JUNO), “Neutrino Physics with JUNO,” J. Phys. G 43, 030401 (2016), eprint 1507.05613.
  * Goldhagen et al. (2021) K. Goldhagen, M. Maltoni, S. Reichard, and T. Schwetz, “Testing sterile neutrino mixing with present and future solar neutrino data,” (2021), eprint 2109.14898.
  * Petcov (1988) S. T. Petcov, “Exact analytic description of two neutrino oscillations in matter with exponentially varying density,” Phys. Lett. B 200, 373 (1988).
  * de Gouvêa (2004) A. de Gouvêa, in _Theoretical Advanced Study Institute in Elementary Particle Physics: Physics in D $\geqq$ 4_ (2004), eprint hep-ph/0411274.
  * Friedland (2001) A. Friedland, “On the evolution of the neutrino state inside the sun,” Phys. Rev. D 64, 013008 (2001), eprint hep-ph/0010231.
  * de Holanda and Smirnov (2004) P. C. de Holanda and A. Y. Smirnov, “Homestake result, sterile neutrinos and low-energy solar neutrino experiments,” Phys. Rev. D 69, 113002 (2004), eprint hep-ph/0307266.
  * de Holanda and Smirnov (2011) P. C. de Holanda and A. Y. Smirnov, “Solar neutrino spectrum, sterile neutrinos and additional radiation in the Universe,” Phys. Rev. D 83, 113011 (2011), eprint 1012.5627.


# Receding-Constraint Model Predictive Control using a
Learned Approximate Control-Invariant Set

Gianni Lunardi1, Asia La Rocca1, Matteo Saveriano1 and Andrea Del Prete1 1The
authors are with the Industrial Engineering Department, University of Trento,
Via Sommarive 11, 38123, Trento, Italy<EMAIL_ADDRESS>work has
been partially supported by the PRIN project DOCEAT (CUP n. E63C22000410001)
and the European Union under the NextGenerationEU project iNest (ECS
00000043).

###### Abstract

In recent years, advanced model-based and data-driven control methods are
unlocking the potential of complex robotics systems, and we can expect this
trend to continue at an exponential rate in the near future. However, ensuring
safety with these advanced control methods remains a challenge. A well-known
tool to make controllers (either Model Predictive Controllers or Reinforcement
Learning policies) safe, is the so-called _control-invariant set_ (a.k.a. safe
set). Unfortunately, for nonlinear systems, such a set cannot be exactly
computed in general. Numerical algorithms exist for computing approximate
control-invariant sets, but classic theoretic control methods break down if
the set is not exact. This paper presents our recent efforts to address this
issue. We present a novel Model Predictive Control scheme that can guarantee
recursive feasibility and/or safety under weaker assumptions than classic
methods. In particular, recursive feasibility is guaranteed by making the
safe-set constraint move backward over the horizon, and assuming that such set
satisfies a condition that is weaker than control invariance. Safety is
instead guaranteed under an even weaker assumption on the safe set, triggering
a safe task-abortion strategy whenever a risk of constraint violation is
detected. We evaluated our approach on a simulated robot manipulator,
empirically demonstrating that it leads to less constraint violations than
state-of-the-art approaches, while retaining reasonable performance in terms
of tracking cost and number of completed tasks.

##  I Introduction

Ensuring safety is crucial in all robotics applications. However, this is more
and more difficult with the recently increasing complexity of control methods
and robotic platforms. Indeed, recent data-driven approaches, often relying on
Reinforcement Learning (RL) algorithms, typically produce black-box policies
that are inherently hard to certify as safe. Moreover, even model-based
control methods for constrained nonlinear systems in practice struggle to
guarantee safety, which consists in recursive constraint satisfaction (a.k.a.
recursive feasibility). This is because the classic approach to guaranteeing
safety, both for Model Predictive Control (MPC) and for Quadratic-Programming-
based control methods, relies on the assumption of knowing a so-called _safe
set_ (a.k.a. control-invariant set) [1, 2], or a Control Barrier Function
(CBF) [3, 4]. However, exactly computing safe sets (or CBFs) for nonlinear
systems is not feasible in general. Therefore, practitioners must rely on
numerical methods to compute approximate versions of such sets (or functions)
[5, 6, 7, 8, 9, 10, 11, 12]. Unfortunately, safety guarantees are lost if the
used safe set is not exact.

In this paper, we present a novel MPC scheme that ensures: i) safety, assuming
the safe set is a _conservative_ approximation of a specific backward
reachable set; ii) recursive feasibility, assuming the safe set is N-step
control invariant, which is a weaker assumption than classic control
invariance. We compared our approach with classic MPC schemes: the standard
formulation (without terminal constraints but a longer horizon), and a
formulation using the safe set to constrain the terminal state. Our method
could successfully avoid constraint violation in more tests than the others,
being able to trade off performance and safety depending on the
conservativeness of the used safe set.

##  II Preliminaries

###  II-A Notation

  * •

$\mathbb{N}$ denotes the set of natural numbers;

  * •

$\\{x_{i}\\}_{0}^{N}$ denotes a discrete-time trajectory given by the sequence
$(x_{0},\dots,x_{N})$;

  * •

$x_{i|k}$ denotes the state at time step $k+i$ predicted when solving the MPC
problem at time step $k$;

###  II-B Problem statement

Let us consider a discrete-time dynamical system with state and control
constraints:

$x_{i+1}=f(x_{i},u_{i}),\qquad x\in\mathcal{X},\qquad u\in\mathcal{U}.$ (1)

Our goal is to design a control algorithm to ensure _safety_ (i.e., constraint
satisfaction), while preserving performance (i.e., cost minimization) as much
as possible. Let us define $\mathcal{S}$ as the set containing all the
equilibrium states of our system:

$\mathcal{S}=\\{x\in\mathcal{X}\ |\,\exists\,u\in\mathcal{U}:x=f(x,u)\\}.$ (2)

To achieve our goal, we rely on the _Infinite-Time Backward-Reachable Set_ [1]
of $\mathcal{S}$, which we denote as $\mathcal{V}$. Mathematically, it is
defined as the subset of $\mathcal{X}$ starting from which it is possible to
reach $\mathcal{S}$ in finite time:

$\displaystyle\mathcal{V}\triangleq\\{x_{0}\in\mathcal{X}\,|\,\exists\\{u_{i}\\}_{0}^{k},k\in\mathbb{N}:$
$\displaystyle x_{k+1}\in\mathcal{S},x_{i}\in\mathcal{X},$ (3) $\displaystyle
u_{i}\in\mathcal{U},\forall\,i=0,\dots,k\\}.$

As all backward reachable sets of equilibrium states, the set $\mathcal{V}$ is
a control-invariant set [1]. This means that, starting from inside
$\mathcal{V}$, it is possible to remain inside $\mathcal{V}$ indefinitely. If
we knew $\mathcal{V}$ we could use it to construct a safe controller. However,
we cannot reasonably assume to know it in general, but we rely instead on a
more realistic assumption.

######  Assumption 1.

We know a conservative approximation of the set $\mathcal{V}$:

$\hat{\mathcal{V}}\subseteq\mathcal{V}$ (4)

Note that $\hat{}\mathcal{V}$ is not control invariant in general.

######  Assumption 2.

We know an upper bound on the number of time steps needed to safely drive the
system to an equilibrium from a state in $\hat{}\mathcal{V}$, which we refer
to as $\bar{N}$.

As discussed in Section I, numerical methods exists to compute approximations
of $\mathcal{V}$. Among the others, the method in [12] can be made
conservative by an appropriate choice of a safety margin and it also produces
an estimate of $\bar{N}$, satisfying Assumption 2. Therefore, we used [12] in
our evaluation. Now we discuss different approaches to exploit
$\hat{\mathcal{V}}$ in an MPC formulation to try to achieve safety.

###  II-C Model Predictive Control and Recursive Feasibility

Let us consider the following MPC problem:

$\displaystyle\operatorname*{minimize}_{\\{x_{i}\\}_{0}^{N},\\{u_{i}\\}_{0}^{N-1}}$
$\displaystyle\quad\sum_{i=0}^{N-1}\ell_{i}(x_{i},u_{i})+\ell_{N}(x_{N})$ (5a)
$\displaystyle\operatorname*{subject\,to}$ $\displaystyle\quad x_{0}=x_{init}$
(5b) $\displaystyle\quad x_{i+1}=f(x_{i},u_{i})\qquad i=0\dots N-1$ (5c)
$\displaystyle\quad x_{i}\in\mathcal{X},u_{i}\in\mathcal{U}\qquad\;\;i=0\dots
N-1$ (5d) $\displaystyle\quad x_{N}\in\mathcal{X}_{N},$ (5e)

where $\ell(\cdot)$/$\ell_{N}(\cdot)$ is the running/terminal cost, $x_{init}$
is the current state, and $\mathcal{X}_{N}\subseteq\mathcal{X}$ is the
terminal set [13].

Even though MPC is one of the most suited frameworks for controlling
constrained systems, ensuring safety (i.e., constraint satisfaction) remains
challenging when the dynamics or the constraints are nonlinear. The most
common approach to ensuring safety is based on _recursive feasibility_ (RF),
which guarantees that, under the assumption of no disturbances/modeling
errors, if an MPC problem is feasible at the first loop, it remains feasible
forever.

RF is guaranteed if the MPC horizon $N$ is _sufficiently_ long (see Section
8.2 of [2]). However, in general we cannot know how long $N$ should be.
Moreover, even if $N$ were known, it may be too long to result in acceptable
computation times.

Alternatively, RF can be guaranteed by using the terminal set
$\mathcal{X}_{N}$ to constrain the final state inside a _control-invariant_
set (see Section II-D). While theoretically elegant, the practical issue with
this approach is that control-invariant sets are extremely challenging (if not
impossible) to compute for nonlinear systems/constraints. A special case of
this approach is when an _equilibrium_ state (or a set of equilibria) is used
as terminal set. This solves the issue of computing control-invariant sets,
but at the price of (potentially drastically) reducing the _basin of
attraction_ of the MPC.

Other approaches to RF exist that rely on the optimality properties of the
solution and the stability of the closed loop (e.g., Section 8.3 of [2]).
However, these approaches require controllability and other conditions on
running and terminal costs. Therefore, they are not applicable to arbitrary
cost formulations as the methods discussed in this paper.

###  II-D Terminal Constraint

As discussed above, a common way to ensure recursive feasibility in MPC is to
constrain the final state inside a control-invariant set, such as
$\mathcal{V}$. Unfortunately, we do not know $\mathcal{V}$, but only
$\hat{\mathcal{V}}$, which is not control invariant in general. Therefore,
using $\hat{\mathcal{V}}$ as terminal set in our MPC does not ensure RF. This
means that our MPC problem could become unfeasible, and at that point classic
MPC theory does not tell us what to do. A common strategy to deal with
unfeasibility is to relax the terminal constraint with a slack variable, which
is heavily penalized in the cost function [14, 15]. In this way, when the
terminal constraint cannot be satisfied, we can still get a solution that
allows us to keep controlling the system, in the hope that eventually the
terminal constraint be satisfied again. However, this approach does not ensure
_safety_ , nor RF, because the soft constraint allows the state to leave
$\hat{\mathcal{V}}$, which eventually can lead to constraint violations.

##  III Safe Model Predictive Control

This section describes our novel MPC scheme, which relies on two components: a
safe task-abortion strategy (Section III-A, and a receding-constraint MPC
formulation (Section III-B), which can be used together (Section III-C).

###  III-A Safe Task Abortion

Our key idea to ensure safety relies on Assumption 1 and 2, and on the
following two assumptions.

######  Assumption 3.

We have access to two computational units, which we refer to as unit A and
unit B.

######  Assumption 4.

We can solve the following OCP for any $x_{init}\in\hat{\mathcal{V}}$, in at
most $N-1$ time steps:

$\displaystyle\operatorname*{minimize}_{\\{x_{i}\\}_{0}^{\bar{N}},\\{u_{i}\\}_{0}^{\bar{N}-1}}$
$\displaystyle\quad\sum_{i=0}^{\bar{N}-1}\ell_{i}(x_{i},u_{i})+\ell_{\bar{N}}(x_{\bar{N}})$
(6) $\displaystyle\operatorname*{subject\,to}$
$\displaystyle\quad\eqref{eq:mpc_initial_conditions},\eqref{eq:mpc_dynamics},\eqref{eq:mpc_path_constraints}$
$\displaystyle\quad x_{\bar{N}}=x_{\bar{N}-1},$

The choice of the cost function is irrelevant, and can simply be used to help
the solver to converge faster.

OCP (6) can be used to find a feasible trajectory to reach an equilibrium
state from $x_{init}$. Now we can describe our strategy to safely abort the
task in case we detect a risk of constraint violation. Let us assume that we
are using a classic MPC formulation with terminal constraint
$x_{N}\in\hat{\mathcal{V}}$, and that at the MPC loop $k$ our problem becomes
unfeasible. In this situation, we can follow these steps to safely abort the
task:

  1. 1.

unit A uses the MPC solution computed at loop $k-1$ to reach the terminal
state $x_{N|k-1}\in\hat{\mathcal{V}}$;

  2. 2.

in parallel, unit B solves OCP (6), using $x_{N|k-1}$ as initial state;

  3. 3.

after reaching $x_{N|k-1}$, we follow the solution of OCP (6) to safely reach
an equilibrium state.

This strategy allows us to reach a safe equilibrium state, where a stabilizing
controller can be used to maintain the system still. Actually, we do not need
to abort the task as soon as one MPC problem becomes unfeasible. While we
follow the last feasible solution, we can keep trying to solve OCP (5). This
strategy is summarized in Alg. 1 and it can guarantee safety, as stated in the
following Lemma.

######  Lemma 1.

Under Assumptions 1 to 4, the hard terminal-constraint MPC with safe task
abortion described in Alg. 1 guarantees that constraints are never violated.

  * Proof.This proof is straightforward. OCP (6) is always feasible because, by Assumption 1 and 2, from any state in $\hat{\mathcal{V}}$ we can reach an equilibrium in at most $\bar{N}$ time steps. Assumption 4 ensures that, by dedicating a computational unit to solving OCP (6), we get a solution before reaching the terminal state of the last feasible MPC problem, $x_{N|k-1}$. After reaching $x_{N|k-1}$, we follow the solution of OCP (6) to reach an equilibrium state, in which we can stay forever without violating the constraints. ∎

Our most critical assumption is probably Assumption 4, which relies on the MPC
horizon $N$ to be sufficiently long, and on $\bar{N}$ not to be too large, to
allow for enough computation time to solve the OCP. This may be challenging
because we can expect $\bar{N}$ to be rather large, since it must be
sufficient to allow the system to reach an equilibrium from any state in
$\hat{\mathcal{V}}$. At the same time, $N$ cannot be set too large because it
is proportional to the computation time of the MPC problem. However, learning-
based warm-start techniques could be used to speed-up computation [16, 17].

Algorithm 1 Terminal-Constraint MPC with Safe Abortion

1:Number of time steps $T$, Initial state $x_{0}$, Initial guess
$\\{x^{g}_{i}\\}_{0}^{N},\\{u^{g}_{i}\\}_{0}^{N-1}$, OCP (5), Safe-abort OCP
(6)

2:$fails\leftarrow 0$ $\triangleright$ Counter for failed OCP’s

3:for $t=0\rightarrow T-1$ do

4:
$\\{x_{i}^{*}\\}_{0}^{N},\\{u_{i}^{*}\\}_{0}^{N-1},feas\leftarrow\textsc{OCP}(x_{t},\\{x^{g}_{i}\\}_{0}^{N},\\{u^{g}_{i}\\}_{0}^{N-1})$

5: if $feas$ = True then $\triangleright$ If OCP’s solution is feasible

6: $fails\leftarrow 0$ $\triangleright$ Reset counter

7: else

8: if $fails=0$ then $\triangleright$ Start solving (6) in Unit B

9: $\textsc{SolveSafeAbortOCP}(x_{N-1}^{g})$

10: if $fails=N-1$ then $\triangleright$ Abort task

11: return FollowSafeAbortTrajectory()

12: $fails\leftarrow fails+1$ $\triangleright$ Increment counter

13:
$\\{x^{*}_{i}\\}_{0}^{N},\\{u^{*}_{i}\\}_{0}^{N-1}\leftarrow\\{x_{i}^{g}\\}_{0}^{N},\\{u_{i}^{g}\\}_{0}^{N-1}$
$\triangleright$ Copy last feasible solution

14: $x_{t+1}\leftarrow f(x_{t},u_{0}^{*})$ $\triangleright$ Simulate system

15:
$\\{x^{g}_{i}\\}_{0}^{N-1},\\{u^{g}_{i}\\}_{0}^{N-2}\leftarrow\\{x_{i}^{*}\\}_{1}^{N},\\{u_{i}^{*}\\}_{1}^{N-1}$

16: $x^{g}_{N},u^{g}_{N-1}\leftarrow x^{g}_{N-1},u^{g}_{N-2}$

####  III-A1 Safe-Abort for Robot Manipulators

During our tests, we have noticed that the safe-abort OCP (6) was hard to
solve for our numerical solver. Therefore, we suggest an alternative
formulation, which is equivalent to (6) for the case of robot manipulators,
but leads to less numerical issues with the solver. Given
$x_{init}=(q_{init},\dot{q}_{init})\in\hat{\mathcal{V}}$, where $q$ are the
joint angles and $\dot{q}$ are the joint velocities, OCP (6) can be
substituted by:

$\displaystyle\operatorname*{maximize}_{\\{x_{i}\\}_{0}^{\bar{N}},\\{u_{i}\\}_{0}^{\bar{N}-1}}$
$\displaystyle\quad d^{\top}\dot{q}_{0}$ (7)
$\displaystyle\operatorname*{subject\,to}$ $\displaystyle\quad q_{0}=q_{init}$
$\displaystyle\quad(I-dd^{\top})\dot{q}_{0}=0$ $\displaystyle\quad
d^{\top}\dot{q}_{0}\leq||\dot{q}_{init}||$
$\displaystyle\quad\eqref{eq:mpc_dynamics},\eqref{eq:mpc_path_constraints},x_{\bar{N}}=x_{\bar{N}-1},$

where $d=\frac{\dot{q}_{init}}{||\dot{q}_{init}||}$ is the initial velocity
direction. OCP (7) is inspired by the VBOC method [12]. Rather than fixing the
initial state as in (6), we fix only the joint angles and the direction of the
joint velocity vector, while maximizing the joint velocity norm. In this way,
the problem is feasible for any $x_{init}\in\mathcal{X}$. In practice, our
solver was always able to solve this formulation, even for cases where
$x_{init}\notin\hat{\mathcal{V}}$, making the Safe Task Abortion more
reliable.

###  III-B Receding-Constraint MPC

Figure 1: Example of Receding-Constraint MPC with $N=4$. After the MPC loop 3,
the receding constraint slides forward because $x_{4|3}\in\hat{\mathcal{V}}$.

Instead of relying exclusively on the final state to ensure safety, we could
exploit the fact that, as long as at least one state
$x_{r}\in\hat{\mathcal{V}}$ (with $1\leq r\leq N$), we know that
$x_{1}\in\mathcal{V}$ because from $x_{1}$ we can reach $x_{r}$. This suggests
that a less conservative constraint to include in our OCP would be:

$(x_{1}\in\hat{\mathcal{V}})\,\lor\,(x_{2}\in\hat{\mathcal{V}})\lor\,\dots\,\lor\,(x_{N}\in\hat{\mathcal{V}})$
(8)

Unfortunately, OR constraints are extremely challenging for numerical solvers.
Even if this constraint cannot be used, we can find other ways to exploit this
insight.

We suggest to adapt online the time step at which we constrain the state in
$\hat{\mathcal{V}}$. For instance, if at the MPC loop $k-1$ we had
$x_{r|k-1}\in\hat{\mathcal{V}}$, at the loop $k$ we know that it is possible
to have $x_{r-1|k}\in\hat{\mathcal{V}}$ (assuming no disturbances and modeling
errors), therefore we can impose this constraint in a hard way. This is
sufficient to ensure safety for $r$ loops, during which this _receding
constraint_ would slide backward along the horizon. However, once the receding
constraint reaches time step 0, we can no longer rely on it to ensure safety.
Therefore, we suggest to maintain also a soft constraint to encourage the
terminal state to be in $\hat{\mathcal{V}}$. This MPC formulation can be
stated as:

$\displaystyle\operatorname*{minimize}_{\\{x_{i}\\}_{0}^{N},\\{u_{i}\\}_{0}^{N-1},s}$
$\displaystyle\quad\sum_{i=0}^{N-1}\ell_{i}(x_{i},u_{i})+\ell_{N}(x_{N})+w_{s}||s||^{2}$
(9) $\displaystyle\operatorname*{subject\,to}$
$\displaystyle\quad\eqref{eq:mpc_initial_conditions},\eqref{eq:mpc_dynamics},\eqref{eq:mpc_path_constraints}$
$\displaystyle\quad x_{r}\in\hat{\mathcal{V}}$ $\displaystyle\quad
x_{N}\in\hat{\mathcal{V}}\oplus s.$

After solving the MPC at loop $k-1$, we can check whether
$x_{N|k-1}\in\hat{\mathcal{V}}$; if that is the case, at loop $k$ we can move
the receding constraint forward on $x_{N-1|k}$, which ensures safety for other
$N-1$ loops. Actually, we can even check whether $x_{i}\in\hat{\mathcal{V}}$,
for any $i>r$, and if that is the case we can set $r=i-1$ at the next loop.
The resulting algorithm is summarized in Alg. 2 (with SafeAbortFlag set to
false), and a simple example is depicted in Fig. 1.

Algorithm 2 Receding-Constraint MPC with Task Abortion

1:Number of time steps $T$, Initial state $x_{0}$, Initial guess
$\\{x^{g}_{i}\\}_{0}^{N},\\{u^{g}_{i}\\}_{0}^{N-1}$, OCP (9), SafeAbortFlag,
Safe-abort OCP (6)

2:$r\leftarrow N$ $\triangleright$ Receding constraint index

3:for $t=0\rightarrow T-1$ do

4: if SafeAbortFlag and $r=0$ then $\triangleright$ Abort task

5: return FollowSafeAbortTrajectory()

6:
$\\{x_{i}^{*}\\}_{0}^{N},\\{u_{i}^{*}\\}_{0}^{N-1}\leftarrow\textsc{OCP}(r,x_{t},\\{x^{g}_{i}\\}_{0}^{N},\\{u^{g}_{i}\\}_{0}^{N-1})$

7: $r\leftarrow r-1$ $\triangleright$ Recede constraint

8: for $k=r+2\rightarrow N$ do $\triangleright$ Search last state in
$\hat{\mathcal{V}}$

9: if $x_{k}^{*}\in\hat{}\mathcal{V}$ then

10: $r\leftarrow k-1$

11: if SafeAbortFlag and $r=0$ then $\triangleright$ Cannot recede

12: $\textsc{SolveSafeAbortOCP}(x_{1}^{*})$ $\triangleright$ Solve (6) in Unit
B

13: $x_{t+1}\leftarrow f(x_{t},u_{0}^{*})$ $\triangleright$ Simulate dynamics

14:
$\\{x^{g}_{i}\\}_{0}^{N-1},\\{u^{g}_{i}\\}_{0}^{N-2}\leftarrow\\{x_{i}^{*}\\}_{1}^{N},\\{u_{i}^{*}\\}_{1}^{N-1}$

15: $x^{g}_{N},u^{g}_{N-1}\leftarrow x^{g}_{N-1},u^{g}_{N-2}$

To clarify the theoretical properties of this receding-constraint formulation,
we first need to introduce the concept of $N$-step control invariant set.

  * Definition

A set $\mathcal{A}\subseteq\mathcal{X}$ is _$N$ -step control invariant_ if,
starting from any state in $\mathcal{A}$, it is possible either to remain in
$\mathcal{A}$, or to leave $\mathcal{A}$ and come back to it within $N$ time
steps:

$\displaystyle\forall x_{0}\in\mathcal{A}:$
$\displaystyle\,\exists\,\\{u_{i}\\}_{0}^{k-1},\,1\leq k\leq N,$ (10)
$\displaystyle
x_{k}\in\mathcal{A},\,x_{i}\in\mathcal{X},\,u_{i}\in\mathcal{U},\,\forall
i=0,\dots,k-1$

This is an extension of the well-known control invariance, with 1-step control
invariance being equivalent to classic control invariance. Now we can state
under which conditions the receding-constraint MPC is recursively feasible.

######  Theorem 1.

Assuming $\hat{\mathcal{V}}$ is $N$-step control invariant and the penalty on
the soft terminal constraint $w_{s}$ is sufficiently large, the Receding-
Constraint MPC formulation described in Alg. 2 (with _SafeAbortFlag_ set to
false) is recursively feasible.

  * Proof.Assume the receding-constraint formulation is feasible at the first MPC loop $k=0$, which implies that $x_{N|0}\in\hat{\mathcal{V}}$. This guarantees recursive feasibility for $N$ loops, during which the receding constraint can slide backward along the horizon. However, since $\hat{\mathcal{V}}$ is $N$-step control invariant by assumption, we know that in one of those $N$ loops it will be possible to satisfy the soft terminal constraint $x_{N}\in\hat{\mathcal{V}}$. This is because at each loop $k$, the MPC solver tries to satisfy condition (10) for a fixed value of $k$. Since we know that (10) is feasible for some $k\in[1,N]$, we can infer that the soft terminal constraint $x_{N|k}\in\hat{\mathcal{V}}$ must be feasible for some MPC loop $k\in[1,N]$. Under the assumption that $w_{s}$ is sufficiently large, we can infer that at that loop $k$, the soft terminal constraint will be exactly satisfied. When this happens, the hard receding constraint moves to time step $N-1$, ensuring recursive feasibility for another $N-1$ loops. At this point the same reasoning can be applied to ensure recursive feasibility indefinetely. ∎

This theorem highlights how the proposed receding-constraint MPC guarantees
recursive feasibility even if the set $\hat{\mathcal{V}}$ is not control
invariant. We rely indeed on a weaker condition, which is _$N$ -step control
invariance_. Our condition is weaker because any $1$-step control invariant
set is also $N$-step control invariant, for any $N>1$, therefore our approach
guarantees recursive feasibility for a larger class of sets, which contains
the class of control-invariant sets. Unfortunately, computing exactly an
$N$-step control invariant set is currently as hard as computing a standard
control invariant set. However, in practice, it is more likely that a
numerical method for approximating control invariant sets produces a set that
is $N$-step control invariant, rather than control invariant. Therefore, as
empirically shown in our results, our approach has a higher probability of
being recursively feasible than a terminal-constraint MPC, even if in practice
we cannot guarantee the assumptions of Theorem 1 to be satisfied.

###  III-C Safe Task Abortion with Receding Constraint

Since in practice we cannot guarantee that $\hat{\mathcal{V}}$ be $N$-step
control invariant, we cannot guarantee that the receding constraint
formulation be recursively feasible. Therefore, we may need to use the task-
abortion strategy when the receding constraint has reached time step 0. The
problem is that at that point we have only one time step to solve OCP (6). In
this paper, we assume that this computation time is enough, and we describe in
Alg. 2 (with SafeAbortFlag set to true) the Receding-Constraint MPC with Task
Abortion.

If one time step were not sufficient to solve (6), several solutions could be
explored. We briefly discuss them in the following, but we leave their
implementation for future work. A possible way to reduce computation time is
to pre-compute a warm-start for OCP (6), before $r$ reaches 0. While we do not
know in which state the system will be at that time, we can use the trajectory
predicted by the MPC as a guess. If this warm-start is not enough to solve OCP
(6) in one time step, we could modify the receding-constraint formulation to
ensure that the pre-computed _safe-abort trajectory_ starts exactly at the
state of the system when the task abortion is initiated. To achieve this, we
must modify the receding constraint from $x_{j|k}\in\hat{\mathcal{V}}$ to the
more conservative $x_{j|k}=x_{j+1|k-1}$. In other words, we constrain the
predicted state in $\mathcal{V}$ not to change across the MPC loops. This is
bound to deteriorate performance, but it should still outperform the standard
Terminal-Constraint MPC.

##  IV Results

This section presents our results111Our code is available at
https://github.com/idra-lab/safe-mpc. comparing five MPC formulations:

  * •

_Naive_ : a classic formulation without terminal constraint, i.e., problem (5)
with $\mathcal{X}_{N}=\mathcal{X}$.

  * •

_Soft Terminal_ (ST): it uses a soft terminal constraint set
$\mathcal{X}_{N}=\hat{\mathcal{V}}$ with a penalty weight of $10^{8}$.

  * •

_Soft Terminal With Abort_ (STWA): as the previous one, but it triggers the
safe abort whenever $x_{N|k}\notin\hat{\mathcal{V}}$.

  * •

_Hard Terminal With Abort_ (HTWA): it uses a hard terminal constraint set
$\mathcal{X}_{N}=\hat{\mathcal{V}}$, and it triggers the safe abort whenever
the OCP is unfeasible (as in Alg. 1).

  * •

_Receding_ : the novel formulation (9) described by Alg. 2, using soft
constraints for both $x_{r}\in\hat{\mathcal{V}}$ (penalty weight of $10^{8}$)
and $x_{N}\in\hat{\mathcal{V}}$ ($w_{s}=10^{5}$).

For the simulations, we have considered a planar triple pendulum, thus
$n_{x}=6,\,n_{u}=3$. We have used CasADi [18] for the symbolic computation of
the dynamics, costs and constraints, and Acados [19] to solve the OCPs and
integrate the dynamics. The OCP is a tracking problem with respect to a static
state, purposely chosen near the joint limits, to test the safety of the
controllers:

$x^{\text{ref}}=(q^{\text{max}}-0.05,\bar{q},\bar{q},0,0,0),$ (11)

with $\bar{q}=(q^{\text{max}}+q^{\text{min}})/2$. We have used as running cost
a least-squares function, penalizing deviations from $x^{\text{ref}}$ and
control efforts:

$\displaystyle l(x,u)$
$\displaystyle=||x-x^{\text{ref}}||_{Q}^{2}+||u||_{R}^{2}$ (12) $\displaystyle
Q$ $\displaystyle=\text{diag}([500,10^{-4}I_{5}]),\quad R=10^{-4}I_{3},$

where $I_{k}$ is the identity matrix with size $k$. Set membership to
$\hat{\mathcal{V}}$ is verified with the constraint:

$(1-\alpha)\phi(x)-||\dot{q}||\geq 0,$ (13)

where $\phi(\cdot)$ is a Neural Network (NN) computing an upper bound on the
joint velocity norm [12], and $\alpha\in[0,1]$ is a safety margin that we
introduced to ensure that $\hat{}\mathcal{V}\subseteq\mathcal{V}$.

We have run 100 simulations for each MPC formulation, starting from the same
100 random joint positions $q_{0}$ with $\dot{q}_{0}=0$. The time step of the
MPCs was $dt=$5\text{\,}\mathrm{ms}$$. The horizon of _Naive_ has been fixed
to $N=36$, so that each MPC iteration takes less than $4\text{\,}\mathrm{ms}$
(leaving $1\text{\,}\mathrm{ms}$ for further operations, to mimic the timing
limitations of a real-time application). We used instead shorter horizons for
the other approaches ($N=35$ for the three terminal-constrained MPC’s, and
$N=34$ for _Receding_), since their MPC iterations take more time due to the
additional constraints.

TABLE I: Number of times each controller completed the task, safely aborted it, or violated a constraint (with $\alpha=2\%$). MPC | Completed | Aborted | Failed
---|---|---|---
|  | (6) / (7) | (6) / (7)
Naive | 68 | - | 32
ST | 83 | - | 17
STWA | 58 | 9 / 16 | 33 / 26
HTWA | 60 | 10 / 17 | 30 / 23
Receding | 72 | 11 / 20 | 17 / 8
TABLE II: Number of times each controller completed the task, safely aborted it, or violated a constraint (with $\alpha=10\%$). MPC | Completed | Aborted | Failed
---|---|---|---
|  | (6) / (7) | (6) / (7)
Naive | 67 | - | 33
ST | 80 | - | 20
STWA | 56 | 20 / 36 | 24 / 8
HTWA | 54 | 21 / 37 | 25 / 9
Receding | 65 | 21 / 33 | 14 / 2

Table I reports the number of tasks completed, safely aborted, or failed by
each controller, using a safety margin $\alpha=2\%$. For the safe abort, we
have tested both formulation (6) and (7). In terms of safety, _Naive_ violated
the constraints the most, while _Receding_ violated them the least (when using
(7) for safe abort). In terms of performance, ST completed more tasks than the
others, but at the price of a higher number of failed tasks than _Receding_.
STWA and HTWA performed strictly worse than ST, completing less tasks and
failing more times. The lower number of completed tasks is explained by the
trigger of the safe abort, while the relatively high number of failures could
be explained by the small safety margin $\alpha$, which is not enough to
ensure $\hat{\mathcal{V}}\subseteq\mathcal{V}$.

Table II reports a similar comparison, but with a higher safety margin
$\alpha=10\%$. The number of completed tasks is slightly smaller for all
approaches using $\hat{\mathcal{V}}$, but the number of failures is remarkably
smaller for STWA, HTWA, and especially for _Receding_ , which failed only 2
times. The number of failures remained large when using (6) for the safe
abort, demonstrating the benefit of formulation (7).

Figure 2: Comparison between ST (task completed) and STWA (task aborted). The
last plot shows the value of the terminal constraint (13). The vertical line
highlights the start of the safe-abort trajectory. Figure 3: Comparison
between ST (task failed) and STWA (task aborted). The last plot shows the
value of the terminal constraint (13). The vertical line highlights the start
of the safe-abort trajectory.

Fig. 2 and 3 highlight the different risk-aversion levels of ST and STWA by
showing the joint trajectories of two simulations with $\alpha=10\%$. In both
cases STWA aborted the task. ST instead completed the first task, while it
failed the second one. ST is willing to take risks, which sometimes leads to
completing the task (Fig. 2), but sometimes it leads to failure (Fig. 3). STWA
is instead risk-averse, and it triggers a safe abort as soon as a risk of
constraint violation is detected, which leads to less completed tasks, but
also less failures.

TABLE III: Mean tracking cost and computation times for the MPC real-time iteration (RTI) and safe abort OCP. MPC | Cost ($\cdot 10^{3}$) | RTI ($\text{\,}\mathrm{ms}$) | Safe Abort ($\text{\,}\mathrm{s}$)
---|---|---|---
(6) | (7)
Naive | $2.867$ | $3.80$ | - | -
ST | $2.868$ | $3.74$ | - | -
STWA | $3.012$ | $4.13$ | $0.17$ | $2.34$
HTWA | $3.012$ | $3.94$ | $0.18$ | $1.98$
Receding | $2.872$ | - | $0.07$ | $3.36$

In terms of cost, Table III shows that the average cost for the completed
tasks (with $\alpha=10\%$) is comparable for the different formulations, thus
the tracking performance is not degraded by the extra constraints using
$\hat{}\mathcal{V}$. The same table also reports the computation times. The
99-percentile for the real-time iteration scheme [20] is always below the time
step duration (5 ms). We do not report the RTI computation times for
_Receding_ because of a technical issue. Indeed, the Python interface of
Acados does not support time-varying constraints. Therefore our current
implementation of _Receding_ actually soft constrains the whole state
trajectory in $\hat{\mathcal{V}}$, but then sets to zero the penalty weights
for all time steps except for $r$ and $N$, resulting in a much higher
computation time than needed. The Safe Abort column reports the maximum
computation times for the Task Abortion with the two methods. As previously
stated, OCP (7) reports a higher number of successes (see Table I and II) at
the cost of large computation times, while (6) reports good computation times
(satisfying Assumption 4), but with a high number of failures.

##  V Conclusions

We have presented a novel Receding-Constraint MPC formulation, which provides
recursive feasibility guarantees under a weaker assumption on the used safe
set with respect to classic approaches. Moreover, we have presented a task-
abortion strategy that allows to reach an equilibrium state whenever a risk of
constraint violation is detected. Our results on a 3-joint manipulator show
the improved safety of the presented Receding-Constraint MPC with respect to
other state-of-the-art methods.

Future research will focus on finding a safe-abort method that achieves high
success rates as (7), but with reasonable computation times as (6). For this,
we plan to extend the method in [12] to learn both the set $\hat{\mathcal{V}}$
and a policy to drive the state to an equilibrium. While this work focused on
model-based control methods, our approach could be applied in the future to
_safety filters_ for making black-box RL policies safe.

## References

  * [1] F. Blanchini, “Set invariance in control,” _Automatica_ , vol. 35, pp. 1747–1767, 1999.
  * [2] L. Grüne, J. Pannek, L. Grüne, and J. Pannek, _Nonlinear model predictive control_. Springer, 2017.
  * [3] A. D. Ames, J. W. Grizzle, and P. Tabuada, “Control barrier function based quadratic programs with application to adaptive cruise control,” in _IEEE Conference on Decision and Control_ , 2014, pp. 6271–6278.
  * [4] Z. Wu, F. Albalawi, Z. Zhang, J. Zhang, H. Durand, and P. D. Christofides, “Control Lyapunov-Barrier function-based model predictive control of nonlinear systems,” _Automatica_ , vol. 109, p. 108508, 2019. [Online]. Available: https://doi.org/10.1016/j.automatica.2019.108508
  * [5] B. Djeridane and J. Lygeros, “Neural approximation of pde solutions: An application to reachability computations,” in _IEEE Conference on Decision and Control_ , 2006, pp. 3034–3039.
  * [6] P.-A. Coquelin, S. Martin, and R. Munos, “A dynamic programming approach to viability problems,” in _IEEE International Symposium on Approximate Dynamic Programming and Reinforcement Learning_ , 2007, pp. 178–184.
  * [7] F. Jiang, G. Chou, M. Chen, and C. J. Tomlin, “Using neural networks to compute approximate and guaranteed feasible hamilton-jacobi-bellman pde solutions,” 2016. [Online]. Available: https://www.arxiv.org/abs/1611.03158
  * [8] V. Rubies-Royo and C. Tomlin, “Recursive Regression with Neural Networks: Approximating the HJI PDE Solution,” in _International Conference on Learning Representations_ , 2017.
  * [9] K. C. Hsu, V. Rubies-Royo, C. J. Tomlin, and J. F. Fisac, “Safety and Liveness Guarantees through Reach-Avoid Reinforcement Learning,” _Robotics: Science and Systems_ , 2021.
  * [10] C. Dawson, S. Gao, and C. Fan, “Safe Control With Learned Certificates : A Survey of Neural Lyapunov , Barrier , and Contraction Methods for Robotics and Control,” _IEEE Transactions on Robotics_ , vol. 39, no. 3, pp. 1749–1767, 2023.
  * [11] Y. Zhou, D. Li, Y. Xi, and Y. Xu, “Data-driven approximation for feasible regions in nonlinear model predictive control,” 2020. [Online]. Available: https://arxiv.org/abs/2012.03428
  * [12] A. La Rocca, M. Saveriano, and A. Del Prete, “VBOC: Learning the Viability Boundary of a Robot Manipulator using Optimal Control,” _IEEE Robotics and Automation Letters_ , 2023.
  * [13] D. Q. Mayne, J. B. Rawlings, C. V. Rao, and P. O. Scokaert, “Constrained model predictive control: Stability and optimality,” _Automatica_ , vol. 36, no. 6, pp. 789–814, 2000.
  * [14] E. C. Kerrigan and J. M. Maciejowski, “Soft constraints and exact penalty functions in model predictive control,” in _UKACC International Conference (Control)_ , 2000.
  * [15] M. N. Zeilinger, M. Morari, and C. N. Jones, “Soft constrained model predictive control with robust stability guarantees,” _IEEE Transactions on Automatic Control_ , vol. 59, no. 5, pp. 1190–1202, 2014.
  * [16] N. Mansard, A. Del Prete, M. Geisert, S. Tonneau, and O. Stasse, “Using a Memory of Motion to Efficiently Warm-Start a Nonlinear Predictive Controller,” in _IEEE International Conference on Robotics and Automation_ , 2018, pp. 2986–2993.
  * [17] G. Grandesso, E. Alboni, G. P. Papini, P. M. Wensing, and A. Del Prete, “CACTO: Continuous Actor-Critic With Trajectory Optimization - Towards Global Optimality,” _IEEE Robotics and Automation Letters_ , vol. 8, no. 6, pp. 3318–3325, 2023.
  * [18] J. A. E. Andersson, J. Gillis, G. Horn, J. B. Rawlings, and M. Diehl, “CasADi – A software framework for nonlinear optimization and optimal control,” _Mathematical Programming Computation_ , vol. 11, no. 1, pp. 1–36, 2019.
  * [19] R. Verschueren, G. Frison, D. Kouzoupis, J. Frey, N. van Duijkeren, A. Zanelli, B. Novoselnik, T. Albin, R. Quirynen, and M. Diehl, “Acados: a modular open-source framework for fast embedded optimal control,” _Mathematical Programming Computation_ , vol. 14, pp. 147–183, 2019.
  * [20] M. Diehl, R. Findeisen, F. Allgöwer, H. G. Bock, and J. P. Schlöder, “Nominal stability of real-time iteration scheme for nonlinear model predictive control,” _IEE Proceedings-Control Theory and Applications_ , vol. 152, no. 3, pp. 296–308, 2005.