corrauthor _correspondence:_<EMAIL_ADDRESS>decdeceased alsoMEPHI
_also at:_ NRNU MEPhI, Moscow, Russia nowKIT _present address:_ Institut für
Experimentelle Teilchenphysik, Karlsruher Institut für Technologie, Karlsruhe
, Germany alsoLev _also at:_ Moscow Inst. of Physics and Technology, Russia
nowBosch _present address:_ Robert Bosch GmbH, Stuttgart, Germany alsoDubna
_also at:_ Dubna State University, Dubna, Russia nowBerkeley _present
address:_ Nuclear Science Division, Berkeley, USA nowTUM _present address:_
Physik Department, Technische Universität München, Germany 11institutetext:
INFN Laboratori Nazionali del Gran Sasso and Gran Sasso Science Institute,
Assergi, Italy 22institutetext: INFN Laboratori Nazionali del Gran Sasso and
Università degli Studi dell’Aquila, L’Aquila, Italy 33institutetext: INFN
Laboratori Nazionali del Sud, Catania, Italy 44institutetext: Institute of
Physics, Jagiellonian University, Cracow, Poland 55institutetext: Institut für
Kern- und Teilchenphysik, Technische Universität Dresden, Dresden, Germany
66institutetext: Joint Institute for Nuclear Research, Dubna, Russia
77institutetext: European Commission, JRC-Geel, Geel, Belgium 88institutetext:
Max-Planck-Institut für Kernphysik, Heidelberg, Germany 99institutetext:
Department of Physics and Astronomy, University College London, London, UK
1010institutetext: Dipartimento di Fisica, Università Milano Bicocca, Milan,
Italy 1111institutetext: INFN Milano Bicocca, Milan, Italy 1212institutetext:
Dipartimento di Fisica, Università degli Studi di Milano and INFN Milano,
Milan, Italy 1313institutetext: Institute for Nuclear Research of the Russian
Academy of Sciences, Moscow, Russia 1414institutetext: Institute for
Theoretical and Experimental Physics, NRC “Kurchatov Institute”, Moscow,
Russia 1515institutetext: National Research Centre “Kurchatov Institute”,
Moscow, Russia 1616institutetext: Max-Planck-Institut für Physik, Munich,
Germany 1717institutetext: Physik Department, Technische Universität München,
Germany 1818institutetext: Dipartimento di Fisica e Astronomia, Università
degli Studi di Padova, Padua, Italy 1919institutetext: INFN Padova, Padua,
Italy 2020institutetext: Physikalisches Institut, Eberhard Karls Universität
Tübingen, Tübingen, Germany 2121institutetext: Physik-Institut, Universität
Zürich, Zurich, Switzerland

# Pulse shape analysis in Gerda Phase II

The Gerda collaborationALNGS,corrauthor
M. AgostiniUCL,TUM  G. AraujoUZH  A.M. BakalyarovKU  M. BalataALNGS  I.
BarabanovINRM  L. BaudisUZH  C. BauerHD  E. BellottiMIBF,MIBINFN,dec  S.
BelogurovITEP,INRM,alsoMEPHI  A. BettiniPDUNI,PDINFN  L. BezrukovINRM  V.
BiancacciPDUNI,PDINFN  E. BossioTUM  V. BotheHD  V. BrudaninJINR  R.
BrugneraPDUNI,PDINFN  A. CaldwellMPIP  C. CattadoriMIBINFN  A.
ChernogorovITEP,KU  T. ComellatoTUM  V. D’AndreaAQU  E.V. DemidovaITEP  N. Di
MarcoALNGS  E. DoroshkevichINRM  F. FischerMPIP  M. FominaJINR  A.
GangapshevINRM,HD  A. GarfagniniPDUNI,PDINFN  C. GoochMPIP  P. GrabmayrTUE  V.
GurentsovINRM  K. GusevJINR,KU,TUM  J. HakenmüllerHD  S. HemmerPDINFN  R.
HillerUZH,nowKIT  W. HofmannHD  J. HuangUZH  M. HultGEEL  L.V.
InzhechikINRM,alsoLev  J. Janicskó CsáthyTUM  J. JochumTUE  M. JunkerALNGS  V.
KazalovINRM  Y. KermaïdicHD  H. KhushbakhtTUE  T. KihmHD  K. KilgusTUE  A.
KirschHD,nowBosch  I.V. KirpichnikovITEP  A. KlimenkoHD,JINR,alsoDubna  K.T.
KnöpfleHD  O. KochetovJINR  V.N. KornoukhovITEP,INRM  P. KrauseTUM  V.V.
KuzminovINRM  M. LaubensteinALNGS  A. LazzaroTUM  M. LindnerHD  I. LippiPDINFN
A. LubashevskiyJINR  B. LubsandorzhievINRM  G. LutterGEEL  C. MacolinoAQU  B.
MajorovitsMPIP  W. ManeschgHD  L. ManzanillasMPIP  M. MiloradovicUZH  R.
MingazhevaUZH  M. MisiaszekCR  Y. MüllerUZH  I. NemchenokJINR,alsoDubna  K.
PanasCR  L. PandolaCAT  K. PelczarGEEL  L. PertoldiTUM,PDINFN  P.
PiseriMILUINFN  A. PulliaMILUINFN  C. RansomUZH  L. RauscherTUE  M.
RedchukPDUNI,PDINFN  S. RiboldiMILUINFN  N. RumyantsevaKU,JINR  C.
SadaPDUNI,PDINFN  F. SalamidaAQU  S. SchönertTUM  J. SchreinerHD  M. SchüttHD
A.-K. SchützTUE,nowBerkeley  O. SchulzMPIP  M. SchwarzTUM  B. SchwingenheuerHD
O. SelivanenkoINRM  E. ShevchikJINR  M. ShirchenkoJINR  L. ShtembariMPIP  H.
SimgenHD  A. SmolnikovHD,JINR  D. StukovKU  A.A. VasenkoITEP  A.
VeresnikovaINRM  C. VignoliALNGS  K. von SturmPDUNI,PDINFN  V. WagnerHD,nowTUM
T. WesterDD  C. WiesingerTUM  M. WojcikCR  E. YanovichINRM  B. ZatschlerDD  I.
ZhitnikovJINR  S.V. ZhukovKU  D. ZinatulinaJINR  A. ZschockeTUE  A.J.
ZsigmondMPIP  K. ZuberDD  and G. ZuzelCR

(Received: date / Accepted: date)

###### Abstract

The GERmanium Detector Array (Gerda) collaboration searched for neutrinoless
double-$\beta$ decay in 76Ge using isotopically enriched high purity germanium
detectors at the Laboratori Nazionali del Gran Sasso of INFN. After Phase I
(2011-2013), the experiment benefited from several upgrades, including an
additional active veto based on LAr instrumentation and a significant increase
of mass by point-contact germanium detectors that improved the half-life
sensitivity of Phase II (2015-2019) by an order of magnitude. At the core of
the background mitigation strategy, the analysis of the time profile of
individual pulses provides a powerful topological discrimination of signal-
like and background-like events. Data from regular 228Th calibrations and
physics data were both considered in the evaluation of the pulse shape
discrimination performance. In this work, we describe the various methods
applied to the data collected in Gerda Phase II corresponding to an exposure
of 103.7 kg$\cdot$yr. These methods suppress the background by a factor of
about 5 in the region of interest around $Q_{\beta\beta}=2039$ keV, while
preserving $(81\pm 3)$% of the signal. In addition, an exhaustive list of
parameters is provided which were used in the final data analysis.

###### Keywords:

neutrinoless double-$\beta$ decay germanium detectors enriched 76Ge pulse
shape analysis

###### pacs:

23.40.-s $\beta$ decay; double-$\beta$ decay; electron and muon capture
27.50.+e mass $59\leq A\leq 89$ 29.30.Kv X- and $\gamma$-ray spectroscopy

††journal: Eur. Phys. J. C

##  1 Introduction

Neutrinoless double-$\beta$ ($0\nu\beta\beta$) decay is a hypothetical process
in which two neutrons in a nucleus are transformed simultaneously into two
protons with the emission of only two electrons. Such a process violates
lepton number conservation and requires the neutrino to be its own
antiparticle (Majorana particle). In combination with cosmological
observations and direct neutrino mass measurement, a non-zero $0\nu\beta\beta$
decay rate would highly constrain the standard light, left-handed neutrino
exchange mechanism via the effective Majorana neutrino mass or shed light on
alternative processes Cirigliano:2018yza .

The highest half-life sensitivity to $0\nu\beta\beta$ decay requires the
experiments to achieve large target mass, high detection efficiency, good
energy resolution and most complete elimination of background at the $Q$-value
of the decay ($Q_{\beta\beta}$). The goal of the GERmanium Detector Array
(Gerda) experiment was to realize a background-free111Number of expected
background events at full exposure in the region of interest below 1.
experiment for the first time. Gerda was located at the underground Laboratori
Nazionali del Gran Sasso (LNGS) of INFN, Italy. Gerda used up to 43 kg of high
purity germanium (HPGe) detectors enriched in the candidate isotope 76Ge up to
88%. They ensure high detection efficiency, low intrinsic background and
excellent energy resolution. The bare HPGe detectors were operated in liquid
argon (LAr), which served as cooling medium and as active shield against
environmental backgrounds at $Q_{\beta\beta}=2039$ keV. The details of the
experimental setup and its upgrade from Phase I to Phase II can be found
elsewhere Ackermann:2012xja ; Agostini:2017hit .

The Phase II data taking took place between December 2015 and November 2019
with an upgrade of the detector array and the surrounding LAr instrumentation
in 2018. Gerda operated three types of enriched HPGe detectors arranged in an
array of 7 strings: 7 semi-coaxial detectors, referred to as coaxial detectors
for brevity, from Phase I with a total mass of 15.6 kg; 30 Broad Energy
Germanium (BEGe) detectors (20.0 kg) Agostini:2014hra ; Agostini:2019mwn ; and
5 inverted coaxial (IC) detectors (9.6 kg) Cooper:201125 ; Agostini:2021wzn ,
which were installed in summer 2018. The accumulated exposure, product of
total detector mass and respective livetime, amounts to 103.7 kg$\cdot$yr for
Phase II. In order to avoid bias in the event selection criteria, Gerda
followed a strict blinding strategy, where events within a $Q_{\beta\beta}\pm
25$ keV energy window were processed only after the analysis had been
finalized.

In $0\nu\beta\beta$ decay, the two electrons deposit their energy in a small
volume (about 1 mm3 Abt:2007iy ) of a germanium detector producing single-site
events (SSEs). On the other hand, background events consisting of $\gamma$
rays from natural radioactivity interact mostly via Compton scattering
producing events with multiple separated energy depositions (multi-site
events, MSEs). Therefore, events with energy depositions in multiple germanium
detectors or in the LAr volume around the detector array are discarded as
background events. The time structure of the germanium signals is used to
identify MSEs in a single detector and additionally recognize events close to
the detector surface, together referred to as pulse shape discrimination
(PSD). In this paper, we present the PSD techniques applied to the Gerda Phase
II data that allow for a background-free operation of the experiment in
combination with the LAr veto Agostini:2017iyd . The PSD methods have been
improved and extended compared to Phase I Agostini:2013jta . They were applied
consistently to the complete Phase II dataset which was split in two parts,
one before and the other after the 2018 upgrade, to account for changes in the
readout electronics (cables and cross-talk). Alternative methods, e.g. applied
to parts of the data in intermediate releases Agostini:2017iyd ;
Agostini:2019hzm , are summarized in the Appendix.

##  2 Signal formation, readout and processing in germanium detectors

Gerda used p-type HPGe diodes of three different geometries called coaxial,
BEGe and inverted coaxial (Fig. 1). They all have a relatively thick
([0.8-2.6] mm) n+ electrode, formed by lithium diffusion, and a thin
($\sim$300 nm) p+ electrode created by boron implantation. The p+ and n+
electrodes are separated by a groove with non-conducting surface. Because of
their small p+ electrode the BEGe and IC detectors belong to the class of
point-contact detectors Luke:1989pcd  which exhibits intrinsic performance
advantages with respect to energy resolution and pulse shape discrimination.

The Ge detectors are operated under reverse bias voltage such that almost the
entire volume is depleted of free charge carriers. An interaction in the
active volume creates a number of electron-hole pairs proportional to the
deposited energy. The charge carriers drift according to the electric field
created by both the positive bias voltage applied to the n+ contact and the
volume charge density due to the net bulk impurity concentration. The
electrons are collected at the n+ contact, the holes at the p+ contact which
is used for readout. The n+ layer covers most of the crystal surface and
features a large Li concentration. Exhibiting zero electric field beyond the
p-n junction till the outer surface, a charge created in this layer will only
experience thermal diffusion with two possible outcomes: recombine (loss) or
reach the depleted volume (collection). As a result, two generic regions can
be identified, the dead layer with high probability of no charge collection
and a transition layer with partial charge collection.

Figure 1: Geometry and weighting potential of a typical BEGe, an inverted
coaxial and a coaxial detector. The p+ and n+ contacts are indicated
schematically.

During the drift of charge carriers, charge is induced on the readout contact
as described by the Shockley-Ramo theorem HE2001250 :

$Q(t)=-Q_{0}\left[\Phi^{w}(\mathbf{r}_{h}(t))-\Phi^{w}(\mathbf{r}_{e}(t))\right]$
(1)

where $Q_{0}$ is the total charge carried by the holes or electrons and
$\Phi^{w}(\mathbf{r}_{h/e}(t))$ is the weighting potential along the drift
path of holes or electrons. The weighting potential is shown in Fig. 1 for the
three different detector geometries used in Gerda also indicating the geometry
of the p+ and n+ surfaces. The weighting potentials have been calculated using
the AGATA Detector Library Bruyneel:2016zih  pulse shape simulation package.

Due to their small p+ contact, BEGe and IC detectors have a weighting
potential distribution that is very small in most of the volume and sizable
only close to the p+ contact. This results in similar waveforms $Q(t)$ from
interactions in a large part of the volume. Multiple energy deposits can be
treated as a superposition of single interactions. In Fig. 2, normalized
example pulses from a BEGe detector are shown for a SSE, a MSE, an event close
to the p+ contact and an event near the dead layer of the n+ contact with
incomplete charge collection. As shown, surface events produce characteristic
pulse shapes being fast close to the p+ contact due to the strong electric
field and slow near the n+ contact due to the weak electric field and the
transition layer.

Figure 2: Examples of different type of normalized charge pulses (red) along
with the derived current pulses (blue) from a BEGe detector: single-site
event, multi-site event, event near the p+ contact and event near the n+
contact with incomplete charge collection.

In coaxial detectors, both the hole and electron drift play a role in the
pulse formation, which result in different pulse shapes throughout the volume
of the detector. In Fig. 3, simulated example pulses from different parts of
the detector are shown. Similarly to BEGe detectors, coaxial detectors also
show special pulse shape characteristics in case of surface events. Indeed,
energy deposits near the groove or the bottom of the borehole cause faster
pulses.

Figure 3: Simulated charge pulses in indicated parts (bulk, groove,
bottom/side of borehole) of a typical coaxial detector. The electronic
response of the readout chain and noise are not included here.

The signals induced on the p+ contact of the Gerda detectors are read out by
charge sensitive amplifiers located in the LAr about 35 cm above the detector
array. The signals are digitized at 25 MHz for $160~{}\mu$s and at 100 MHz for
10 $\mu$s. Both traces are centered at the rising edge of the charge pulse
$Q(t)$. The offline analysis of the digitized signals follows the procedures
described in Agostini:2011mh ; Agostini:2011xe . The 25 MHz traces used for
the energy reconstruction ensure the excellent energy resolution, while the
100 MHz traces are used in the PSD methods presented in the following
sections. The energy estimator $E$ is reconstructed with a zero-area cusp
filter Agostini:2015pta , whose parameters are optimized for each detector and
calibration run.

##  3 Overview of event samples and discrimination methods

Weekly calibration runs with 228Th sources are performed to determine the
energy scale and resolution of the detectors Agostini:2021xyz  and to
calibrate and train the PSD techniques. Figure 4 shows a calibration spectrum
highlighting the different event samples used in pulse shape analysis. The
most prominent feature is the full energy peak (FEP) at 2615 keV from 208Tl
decay. Its double escape peak (DEP) at 1593 keV is used as a sample of SSEs as
the electron and positron from pair production deposit their energy in a small
volume and both annihilation $\gamma$ rays leave the detector. The FEP at 1621
keV from 212Bi is used as a sample of MSEs that is sufficiently near in energy
to the DEP in order to avoid noise dependent biases. To test the performance
of MSE rejection, the FEP and single escape peak (SEP) of 208Tl, mostly
featuring MSEs are used while the Compton continuum region around
$Q_{\beta\beta}\pm 35$ keV (CC($Q_{\beta\beta}$)) serves to estimate the
background rejection in the $0\nu\beta\beta$ decay signal region.

Figure 4: Calibration spectrum highlighting the different event samples used
in the pulse shape analysis.

In physics data, the standard neutrino accompanied double-$\beta$
($2\nu\beta\beta$) decay provides another sample of signal-like events that is
equally distributed in the whole detector volume and used for the
investigations in the following. After applying the LAr veto cut
Agostini:2017hit , about 97% of the events in the $1000-1300$ keV region
originate from $2\nu\beta\beta$ decays. Beside MSEs from $\gamma$ rays, the
physics data have a significant amount of surface events from $\alpha$ and
$\beta$ decays that can be discriminated thanks to their specific pulse shape.
The physics spectrum at low energies is dominated by $\beta$ decays of 39Ar up
to its $Q$-value of 565 keV. However, these events are not used in the pulse
shape analysis due to their relatively high noise. A prominent background
source at $Q_{\beta\beta}$ is the $\beta$ decay of 42K, which is produced as a
progeny of the long-lived 42Ar and has a $Q$-value of 3525 keV. Beta particles
deposit their energy in germanium within a few mm resulting in events partly
in the dead and transition layers of the n+ surface. Such n+ surface events
can induce slow pulse shapes with incomplete charge collection. Apart from
possible HPGe bulk contamination, that have been shown to be insignificant
Abramov:2019hhh ; Agostini:2017lim , $\alpha$ particles can only reach the
active volume of the detector at the thin p+ contact or at the non-conducting
groove producing pulse shapes with fast rise as shown in Figs. 2 and 3. A
clean sample of $\alpha$ surface events is found in physics data above the 42K
$Q$-value. The most prominent structure at these energies is a broad peak at
5304 keV, the $\alpha$ energy of the 210Po $\alpha$ decay (238U decay chain)
Abramov:2019hhh ; Agostini:2014 .

Due to their different geometries, BEGe and IC detectors are treated
separately from coaxial detectors in the pulse shape analysis. In the case of
the BEGe and IC detectors one parameter, $A/E$, is used to classify background
events, where $A$ is the maximum current amplitude as indicated in Fig. 2 and
$E$ is the energy. As MSEs and surface events at the n+ contact are
characterized by longer, i.e. wider current pulses, they feature a lower $A/E$
value compared to SSEs, while surface events at the p+ contact show a higher
$A/E$ value Budjas:2009zu . Therefore, rejecting events on both sides of the
$A/E$ distribution of SSEs enhances the signal to background ratio. The
details of the $A/E$ analysis are presented in Sec. 4. In the case of coaxial
detectors an artificial neural network (ANN TMVA ) is used to discriminate
SSEs from MSEs similar to the approach applied in Phase I Agostini:2013jta .
To discard events close to the p+ contact, a dedicated cut on the risetime of
the pulses is applied. The training and optimisation of the ANN and risetime
cuts for coaxial detectors is described in Sec. 5. An additional cut is
applied to all detectors to remove events with slow or incomplete charge
collection. These events arise from energy depositions in a non-depleted
volume (n+ layer or insulating groove). These events are identified through
the difference between two energy estimates performed using the same digital
filter but different shaping times as summarized in Sec. 6. The signal
efficiency of $0\nu\beta\beta$ decay is estimated using the survival fraction
of DEP and $2\nu\beta\beta$ decay events by taking into account the energy
dependence of the different PSD techniques. The details of this extrapolation
to $Q_{\beta\beta}$ including systematic uncertainties are found in Sec. 7 and
B.

##  4 The A/E method for BEGe and IC detectors

The amplitude of the current pulse $A$ is computed after applying 3 times a
moving window average (MWA) filter with 50 ns length and interpolating the
pulse down to 1 ns sampling time. This filtering procedure optimizes the high
frequency noise attenuation while preserving the pulse shape information. The
energy estimator $E$ is determined by a pseudo Gaussian filter with a shaping
time of 10 $\mu$s. $A$ is then divided by $E$ before calibration, providing
the raw $A/E$ for each pulse. The raw $A/E$ is then corrected for time
stability and energy dependence before a cut value is defined.

For each calibration run, the $A/E$ distribution of events in the $1000-1300$
keV region is fitted with a Gaussian (SSEs) and a low-side tail (MSEs) as
described in Agostini:2013jta . The position of the Gaussian $\mu_{A/E}$ from
each calibration of the four years of data taking is used to define stable
time periods where the raw $A/E$ changes by less than its $\sigma_{A/E}$
resolution. Instabilities are mostly related to hardware changes and a few
detectors show a small systematic drift of the raw $A/E$. Physics data between
the stable periods are removed from the analysis causing a few percent
exposure loss. After normalizing the raw $A/E$ by the average $\mu_{A/E}$
within a given time period, the data of all calibrations are merged, only
separating before and after the upgrade of 2018.

The $A/E$ of SSEs in the Compton continuum of the merged calibration data show
a small linear energy dependence of a few percent per MeV, due to the larger
charge cloud size at higher energies that broadens the current pulse. The
energy dependence of $\mu_{A/E}(E)$ and $\sigma_{A/E}(E)$ is described by a
linear and a $\sqrt{b+c/E^{2}}$ type of function, respectively, as shown in
Fig. 5 for a BEGe (GD61A) and an inverted coaxial (IC74A) detector as
examples. In addition to the correction for the energy dependence, $A/E$ is
normalized to the mean of the $A/E$ distribution of the DEP, which lies about
0.25% above the SSE band.

Figure 5: Distribution of $A/E$ as a function of energy from 228Th calibration
data for the BEGe detector GD61A (top) and the IC detector IC74A (bottom). The
dashed lines indicate the linear energy dependence and the 1-$\sigma$ width of
SSEs.

The cut values for each detector are determined on the energy-independent
$A/E$ classifier defined as
$\zeta=([A/E]\,/\,\mu_{A/E}(E)-1)\,/\,\sigma_{A/E}(E)$. Its distribution is
centered around zero and has a standard deviation of one for SSEs. The low-
side $A/E$ cut against MSEs and n+ surface events is chosen to yield a DEP
survival fraction of 90%. The resulting cut values range from $\zeta=-1.9$ to
$-1.2$ and from $\zeta=-1.9$ to $-1.7$ for BEGe and IC detectors,
respectively. The high-side $A/E$ cut against p+ surface events is chosen at
$\zeta=3.0$ for each detector, in order to reject all $\alpha$ events in
physics data above 3525 keV. It has been shown that the high-side $A/E$ cut
discards events, including degraded $\alpha$ events, from a small volume
around the p+ contact Comellato:2020ljj  causing the survival fraction of
events after the high-side $A/E$ cut to be proportional to the detector mass.

Fig. 6 shows for BEGe and IC detectors the energy distribution of calibration
data before and after the $A/E$ cuts described above as well the corresponding
survival fractions. As expected the survival fraction of events in the FEPs
and SEP is much lower than in the DEP. Events in the Compton continuum are
discarded with about the same probability independent of their energy but
depending on the overall detector size, and more generally speaking from the
detector type. IC detectors discard a higher fraction of events because of the
higher probability of multiple scattering of $\gamma$ rays due to their larger
size.

Figure 6: Energy distribution of events from calibration data before and after
the $A/E$ cut (top) and their ratio (bottom) for BEGe and IC detectors.

The survival fractions of events in the DEP, FEP and CC($Q_{\beta\beta}$) are
shown in Fig. 7 for each detector before and after the upgrade. The DEP
survival fractions are slightly smaller than 90% due to the high $A/E$ cut.
The rejection of MSEs shows a small dependence on the detector position in the
string because of different electronic noise conditions. This effect was
reduced after the upgrade. The IC detectors (detector numbers above 35) reject
MSEs more efficiently than BEGe detectors.

Figure 7: Survival fraction of events in the 208Tl DEP, 212Bi FEP and
CC($Q_{\beta\beta}$) for each detector (see Table 3 for detector numbers and
types). Open (filled) symbols show the calibration dataset before (after) the
upgrade. The dashed lines separate the detector strings in the array. The
uncertainties are only statistical and smaller than the markers.

In order to check the validity of the $A/E$ corrections and cuts in the whole
dataset, survival fractions of the usual event samples are studied for each
calibration run. The average survival fraction from all BEGe detectors is
shown for each calibration in Fig. 8 as a function of time. The $A/E$ cut
shows, for both BEGe and IC detector types, a stable behaviour at the 3%
relative level during the whole data collection period when applied to Compton
continuum events at $Q_{\beta\beta}$. Residual instabilities stem on one side
from changes in the detector gain or leakage current and on the other side
from statistical fluctuations.

Figure 8: Average survival fractions of events in the 208Tl DEP, 212Bi FEP and
CC($Q_{\beta\beta}$) for BEGe detectors as a function of time. Each data point
represents a calibration run with its statistical uncertainty.

The raw $A/E$ is corrected in the same way for physics data and the same cut
values are applied as in calibration data. In order to check the validity of
the corrections and the cut values, Fig. 9 shows a comparison of the $A/E$
classifier between DEP events from calibration data and $2\nu\beta\beta$
events from physics data for BEGe and IC detectors. By construction, the $A/E$
classifier peaks at 0 and has a width of about 1 for these SSEs. The agreement
between physics and calibration data is satisfactory and confirms the applied
correction procedure. For each detector, the residual difference of the
average $A/E$, between calibration and physics data, is included in the
systematic uncertainties by shifting the cut value accordingly.

Figure 9: $A/E$ classifier distributions after LAr veto comparing events in
the DEP from calibration data and $2\nu\beta\beta$ decay events from physics
data from BEGe (left) and IC (right) detectors. The histograms are normalized
to their integrals.

The energy distributions of events before and after the $A/E$ cut are shown in
Fig. 10 for the whole physics data corresponding to 61.8 kg$\cdot$yr exposure
from BEGe (53.3 kg$\cdot$yr) and IC (8.6 kg$\cdot$yr) detectors. Events in the
$2\nu\beta\beta$ decay region survive the cut with a high probability, while
the 40K and 42K peaks at 1461 keV and 1525 keV, respectively, mostly featuring
MSEs, are reduced significantly. High energy events above 3525 keV coming from
p+ surface events are all discarded by the high $A/E$ cut by definition. In a
240 keV wide window around $Q_{\beta\beta}$, only 7 events in BEGe detectors
and 1 event in IC detectors survive the LAr veto and $A/E$ cuts. This results
in the corresponding unique background indices222The background index is
evaluated in the range between 1930 and 2190 keV without the two intervals
($2104\pm 5$) and ($2119\pm 5$) keV due to known $\gamma$ rays and without the
signal interval ($Q_{\beta\beta}\pm 5$) keV Agostini:2020prl . The quoted
uncertainties are statistical only. of Gerda in the region of interest of
$5.5_{-1.8}^{+2.4}\cdot 10^{-4}$ counts/(keV$\cdot$kg$\cdot$yr) and
$4.9_{-3.4}^{+7.3}\cdot 10^{-4}$ counts/(keV$\cdot$kg$\cdot$yr) for BEGe and
IC detectors, respectively.

Figure 10: Energy distributions of all Gerda Phase II physics data from the
BEGe and IC detectors before (grey) and after (blue) the $A/E$ cut. The
position of $Q_{\beta\beta}$ is indicated. The prominent 40K and 42K lines are
located at 1461 keV and 1525 keV, respectively.

##  5 The ANN and risetime methods for coaxial detectors

The 100 MHz waveform traces are used to compute the artificial neural network
input variables (IVs) and the risetime. They are first filtered with a MWA of
30 ns width three times. They are subsequently baseline subtracted and
normalized between 0 and 1 by the amplitude of the 25 MHz pole-zero corrected
traces. This amplitude is provided by a trapezoidal filter with a typical
precision of 0.2% at 2 MeV.

The ANN IVs are a list of the 50 times at which the resulting waveform reaches
[1%, 3%, …, 99%] of its amplitude (see Fig. 11). An interpolation of the 10 ns
wide gaps between data points allows for a more precise estimation. These IVs
are computed for all physics and calibration events found in Gerda. However,
given the degraded signal-to-noise ratio at low energy, only events above 1000
keV are considered in this pulse shape analysis to avoid loosely constrained
energy dependence correction. Calibration runs are used to optimize the
discrimination of SSEs from MSEs. The network is built on two hidden layers
with 50 and 51 neurons, using the TMVA - MLPBFGS algorithm TMVA .

Fig. 12 shows as an example the classifier distributions from the ANN training
of the ANG5 detector with the events from the indicated DEP and FEP peaks. The
distributions from the Compton events under the peaks are statistically
subtracted using the distributions of the events in the energy side-bands of
the peaks. The lower and upper side-bands are defined by selecting events
falling in the [$-9\sigma$,$-4.5\sigma$] and [$4.5\sigma$,$9\sigma$] energy
regions where $2.355\cdot\sigma$ is the full width at half maximum used to
quantify the energy resolution of the Ge detectors. The indicated ANN cut
keeps 90% of the events in the DEP peak.

Due to the significant change in hardware, the data taken before and after the
2018 upgrade periods have been trained separately. Finer splittings of the
data have yielded signal efficiencies and background rejection values which
agreed on the one-percent level. As a result of the limited statistics, a
minimal number of two datasets was preferred. Typically, about 10000 and 15000
events enter the signal (DEP) and background (212Bi FEP) samples,
respectively, for the ANN training of each detector. Similarly to the $A/E$
method, the ANN cut for each detector and each training period is set on the
208Tl DEP classifier distribution such that 90% of the events survive.

The risetime is estimated after interpolating the waveform with a 1 ns time
step. From studies on the rejection of $\alpha$ particles Lazzaro:2019pda ,
the [10%-90%] amplitude signal risetime was selected (see Fig. 11). This
parameter is used to reject $\alpha$ events on the p+ contact that develop
faster signals (see Fig. 3). The cut definition relies on the maximization of
the following figure of merit:

$f(x)=\varepsilon^{2}_{2\nu\beta\beta}(x)\cdot(1-\varepsilon_{\alpha}(x)),$
(2)

where $\varepsilon_{2\nu\beta\beta}(x)$ is the $2\nu\beta\beta$ survival
fraction at risetime cut $x$ and $\varepsilon_{\alpha}(x)$ is the
corresponding survival fraction of $\alpha$ events. Only physics data after
ANN-MSE and LAr veto, to increase purity of samples, are used for this figure
of merit that allows to reject most of the $\alpha$ particles while preserving
a high $2\nu\beta\beta$ decay signal survival fraction. Figure 13 depicts an
example of such a cut definition. On average, about 90% of the high energy
$\alpha$ particles are rejected.

Figure 11: Zoom on a typical normalized 100 MHz trace of a coaxial detector.
The 50 ANN input variables (red circles) and rise time estimates (dashed
green) are indicated. Figure 12: Trained ANN classifier values for events in
the 208Tl DEP (SSEs) and 212Bi FEP (MSEs) from the ANG5 detector. The
histograms are normalized to their integrals. Figure 13: Survival fractions
after the rise time cut for $2\nu\beta\beta$ decay and $\alpha$ particles
events from the Phase II data before the upgrade of the ANG4 detector. Also
shown are the figure of merit and the chosen cut value.

The calibration energy spectrum of all coaxial detectors before and after
applying the PSD cuts is shown in Fig. 14. As targeted, the ANN cut removes
preferentially the regions highly populated by MSEs (FEPs and SEP in
particular) and preserves 90% of the 208Tl DEP. On the contrary, the risetime
cut deployed to reject events with fast risetime is insensitive to these types
of events and hence has a high survival fraction for both SSEs and MSEs.

Figure 14: Energy distributions of calibration data events from the coaxial
detectors before and after the ANN and risetime cuts (top) and the
corresponding survival fractions (bottom). Figure 15: Survival fractions of
events in the 208Tl DEP, 212Bi FEP and CC($Q_{\beta\beta}$) for each coaxial
detector after ANN and risetime cuts. Open (filled) symbols show the
calibration dataset before (after) the 2018 upgrade. The uncertainties are
only statistical and smaller than the markers.

Various survival fractions for all coaxial detectors, after applying both ANN
and risetime cuts for the two data taking periods, are depicted in Fig. 15. In
addition, the stabilities of these cuts in three calibration energy regions
are plotted in Figs. 16 and 17. Apart from a slight improvement of the ANN
rejection of the Compton continuum events at $Q_{\beta\beta}$ after the 2018
upgrade due to an improved signal cable management, an overall 3% level
stability in PSD performance is observed over the course of Gerda Phase II.
Compared to Gerda Phase I Agostini:2013jta , on average a 7% relative
worsening is observed, mostly attributed to the different electronics scheme,
hence different noise. The risetime cut also shows a very stable behavior
during this period.

Figure 16: Average survival fractions of events after the ANN cut in the 208Tl
DEP, 212Bi FEP and CC($Q_{\beta\beta}$) for coaxial detectors as a function of
time. Each data point represents a calibration run with its statistical
uncertainty. Figure 17: Average survival fractions of events after the
risetime cut in the 208Tl DEP, 212Bi FEP and CC($Q_{\beta\beta}$) for coaxial
detectors as a function of time. Each data point represents a calibration run
with its statistical uncertainty.

The result of the pulse shape analysis applied on the full 41.8 kg$\cdot$yr
physics data exposure, after the LAr veto cut, is summarized in Fig. 18. The
ANN is preserving 80.5(3)% of the $2\nu\beta\beta$ decay event sample while
removing 62(1)% of the 40K and 42K lines. The reduction of the 214Bi FEPs at
1806 keV and 2204 keV is also visible. However, it is in general inefficient
at rejecting fast and degraded $\alpha$ events from about 3500 keV to 5500
keV, as 58.4(8)% of them survive. After applying also the risetime cut,
67.8(4)% of $2\nu\beta\beta$ decay events remain. The $\alpha$ background is
suppressed by a factor of $\sim$23 above 3500 keV, only 93 events out of 2169
survive the PSD cuts. The complementarity between the two cuts yields for the
coaxial detectors a dataset that remains background-free in the region of
interest with a background index of $5.0^{+2.6}_{-2.0}~{}\cdot~{}10^{-4}$
counts/(keV$\cdot$kg$\cdot$yr).

Figure 18: Top: ANN classifier of physics events from coaxial detectors as a
function of energy. Middle: Risetime of events surviving the ANN cut as a
function of energy. Bottom: Energy distribution of events from coaxial
detectors before and after the ANN and risetime cuts.

##  6 Events with incomplete charge collection

Events from the n+ layer or the groove featuring slow or incomplete charge
collection (see Fig. 2) have an uncertainty on energy reconstruction because
the ZAC filter Agostini:2015pta  is optimized for the FEP resolution with a
relatively short integration time. These particular events, especially for
coaxial detectors, can survive the ANN and risetime cuts. In order to discard
these events with uncertain ZAC energy reconstruction, an additional
classification based on energies reconstructed with pseudo Gaussian filters
with different integration times is performed. The energy is reconstructed
with an integration time of 4 $\mu$s ($E_{short}$) and 20 $\mu$s ($E_{long}$).
The ratio $E_{short}/E_{long}$ is then normalized to its average observed in
events from the 2615 keV line in calibration data. The classifier is defined
as:

$\delta E=\left(\frac{E_{short}/E_{long}}{\langle
E_{short}/E_{long}\rangle_{FEP}}-1\right)\cdot E$ (3)

where $E$ is the default ZAC-reconstructed energy for each event. With this
definition, the classifier $\delta E$ has an average of 0 keV.

The $E_{short}$/$E_{long}$ distribution normalization is performed time
dependently in order to account for possible instabilities of the readout
electronics: for each calibration run the mean of the fitted Gaussian to the
FEP distribution is used as a normalization factor for the following physics
events.

Figure 19 shows the resulting $\delta E$ distributions for calibration events
in the 208Tl FEP and in the CC($Q_{\beta\beta}$) as well as for physics events
in the $2\nu\beta\beta$ region for the ANG2 detector as an example. Events in
the Compton continuum have a higher fraction with large negative $\delta E$
values. This is due to the higher fraction of pulses with incomplete charge
collection that is not present in a peak where the whole energy has been
collected. Physics data in the $2\nu\beta\beta$ region shows a similar
Gaussian distribution as calibration data with no significant energy
dependence.

Figure 19: Distribution of the $\delta E$ classifier for calibration events in
the 208Tl FEP and in the Compton region around $Q_{\beta\beta}$ as well as for
physics events in the $2\nu\beta\beta$ region for the coaxial ANG2 detector.
The histograms are normalized to their integrals.

A cut value is applied to the lower side of the $\delta E$ distribution and it
is set for each detector separately $3\sigma$ away from zero.333No event with
$\delta E$ classifier greater than $3\sigma$ was found in the background
analysis window before $A/E$ or ANN+risetime cuts, therefore no high $\delta
E$ cut was applied. The cut value is loose enough that more than 99% of signal
events are kept whereas those with uncertain energy (significant difference
between energies reconstructed by the ZAC and Gauss algorithms) are rejected.

For the BEGe and IC detectors, most of the $2\nu\beta\beta$ decay events
rejected by the $\delta E$ cut are also rejected by the $A/E$ cut. For BEGe
detectors out of 3477 events cut by either method only 7 are rejected by the
$\delta E$ cut (for energies above 1000 keV 12 events in total are cut by
$\delta E$). For the IC detectors none of the 521 $2\nu\beta\beta$ decay
events is removed by the $\delta E$ cut alone.

For the coaxial detectors the correlation between ANN, risetime and $\delta E$
cuts is weaker. Analyzing the $2\nu\beta\beta$ decay region for the full Phase
II dataset, one gets 4970 events (out of 15433) removed by either method, 1
event cut by all three methods and 57 by $\delta E$ only. 83 events from
coaxial detectors with $E>1000$ keV are cut by $\delta E$ only. Many of these
events show slow pulses that could originate from the detector surface but are
not cut by ANN or risetime cuts.

The survival fractions for different energy regions of physics and calibration
data have been studied for each detector separately. Survival fractions of
$2\nu\beta\beta$ decay events for the different detector types are presented
in Table 1 before and after the other PSD cuts. Values for DEP, FEP,
CC($Q_{\beta\beta}$) and $2\nu\beta\beta$ decay events for the two datasets
before and after the 2018 upgrade are very close to 100%. The Compton region
as well as the $2\nu\beta\beta$ region show lower acceptance because of the
contribution of slow pulses.

For BEGe and IC detectors the impact of the $\delta E$ cut is negligible,
while for coaxial detectors a small correction of the efficiency has to be
taken into account.

Table 1:  Survival fractions of $2\nu\beta\beta$ decay events for the $\delta E$ cut without and in combination with other PSD methods ($A/E$ or ANN+risetime). The uncertainties are statistical only. Detector | Survival fraction [%]
---|---
type | before PSD | after PSD
Coaxial | 99.57 $\pm$ 0.05 | 99.46 $\pm$ 0.07
BEGe | 98.47 $\pm$ 0.09 | 99.96 $\pm$ 0.02
IC | 98.58 $\pm$ 0.20 |  $>99.90$ (95% C.I.)

##  7 PSD detection efficiencies at Qββ

In the absence of signal proxies at $Q_{\beta\beta}$ in calibration and
physics data and due to a not sufficiently accurate modeling of the pulse
shape response in simulation, $\varepsilon_{\rm{PSD}}$ is estimated for point
contact and coaxial detectors from the extrapolation of the survival fractions
of the 208Tl DEP and of the $2\nu\beta\beta$ decay events at 1593 keV and 1150
keV to $Q_{\beta\beta}$ (see circles and squares in Fig. 20).

Table 2: Average signal detection efficiencies at $Q_{\beta\beta}$ of individual PSD methods for the different detector types and data taking periods. $\varepsilon_{\tiny{\textrm{PSD}}}$ is estimated via a MC sampling of individual values (see Table 3) and is thus different from the product of the individual PSD methods efficiencies reported here. | Dec 2015-May 2018 |  | July 2018-Nov 2019
---|---|---|---
| Coaxial | BEGe |  | Coaxial | BEGe | Inverted coaxial
Exposure $\mathcal{E}$ | 28.6 kg$\cdot$yr | 31.5 kg$\cdot$yr |  | 13.2 kg$\cdot$yr | 21.9 kg$\cdot$yr | 8.5 kg$\cdot$yr
$\varepsilon_{\tiny{\textrm{ANN}}}$ | 82.5% | - |  | 81.8% | - | -
$\varepsilon_{\tiny{\textrm{risetime}}}$ | 85.7% | - |  | 85.0% | - | -
$\varepsilon_{\tiny{\textrm{A/E}}}$ | - | 88.4% |  | - | 89.3% | 90.0%
$\varepsilon_{\tiny{\textrm{$\delta E$}}}$ | 99.5% | 100.0% |  | 99.7% | 100.0% | 99.7%
$\varepsilon_{\tiny{\textrm{PSD}}}$ | (69.1 $\pm$ 5.6)% | (88.2 $\pm$ 3.4)% |  | (68.8 $\pm$ 4.1)% | (89.0 $\pm$ 4.1)% | (90.0 $\pm$ 1.8)%
Figure 20: Extrapolation of reference survival fractions (circles) to
$Q_{\beta\beta}$ (squares) using the energy dependence deduced from indicated
two samples of down-scaled waveforms (see text). Examples are given for BEGe
(top), coaxial (middle) and IC detectors (bottom) for the $A/E$, ANN and
$\delta E$ methods, respectively. As to the BEGe - $A/E$ example on top: the
reference survival fraction is corrected for Compton events below the 208Tl
DEP while for the rescaled distribution (red) this correction is missing.

For that purpose, sets of down-scaled waveforms are used to evaluate the
energy dependence of the various PSD methods. This is based on the observation
that the energy dependence of PSD classifiers is dominated by the electronic
noise in Gerda. Two samples from calibration data have been considered, namely
208Tl DEP and CC($Q_{\beta\beta}$) events (see Fig. 4). The waveform
$\omega^{\prime}(k,E)$ at energy $E$ is produced by down-scaling the waveform
$\omega(k,E_{S})$, of true energy $E_{S}$, from one of these two samples and
superimposing it to a baseline sample $b(k)$ obtained from random triggers:

$\omega^{\prime}(k,E)=\omega(k,E_{S})\cdot\frac{E}{E_{S}}+b(k)\cdot\left(\sqrt{1-\left(\frac{E}{E_{S}}\right)^{2}}\right)$
(4)

where $k$ is the index running from 0 $\mu$s to 160 $\mu$s and from 0 $\mu$s
to 10 $\mu$s for the low and high frequency traces introduced in Sec. 2,
respectively. By applying the complete set of digital processing steps and
analysis methods described earlier, such a procedure allows to estimate the
contribution of the energy dependence, induced by the signal-to-noise ratio
variation, on the PSD methods results (for more details see Lazzaro:2019pda ).

Figure 20 shows the survival fractions of the 208Tl DEP and
CC($Q_{\beta\beta}$) datasets as a function of the down-scaled energy for the
$A/E$, ANN and $\delta E$ methods444The risetime method is not shown as its
survival fraction is insensitive to the signal-to-noise ratio at these
energies.. Each plot depicts the energy dependence of a given detector type,
as a result of a weighted average of the individual detector responses, to
illustrate the general trends. The reference efficiency is extrapolated to
$Q_{\beta\beta}$ using the average energy dependence of 208Tl DEP and
CC($Q_{\beta\beta}$) datasets555Energy dependencies are calculated from [1050,
1550] keV and [1550, 1900] keV for the 208Tl DEP and CC($Q_{\beta\beta}$)
datasets, respectively.. Table 2 shows how the different PSD procedures
contribute to the resulting average detection efficiencies
$\varepsilon_{\rm{PSD}}$ of the various detector types for each data taking
period. The efficiencies of the point contact detectors are almost 30% larger
than those of the coaxial detectors.

The analysis of the final Gerda results Agostini:2020prl  is based on time-
dependent and detector-wise datasets. Table 3 (Appendix A) reports the overall
PSD efficiency for each detector separately. Since the central values of the
individual detector efficiencies exhibit significant shifts, the average
overall PSD efficiencies $\varepsilon_{\rm{PSD}}$ of Table 2 have been
obtained, in fact, from an exposure weighted Monte-Carlo sampling, and not
from a simple average.

The systematic uncertainty of the extrapolation is estimated from the
difference of the slopes of the two down-scaled waveform samples. It is on
average about 0.5% and 1.3% for point contact and coaxial detectors,
respectively, the latter detector type being more sensitive to the noise due
to its larger p+ contact capacitance.

In addition, two other systematic effects have been taken into account: 1) the
difference between the calibration and physics data and 2) the difference
between signal proxies and $0\nu\beta\beta$ decay events. The former applies
only to point contact detectors as the 208Tl DEP from calibration data is used
as signal proxy. Shifting the $A/E$ cut, for each detector, by the $A/E$
distribution bias observed between the 208Tl DEP and $2\nu\beta\beta$ decay
(Fig. 9) leads to an average relative systematic uncertainty of 1.9%. The
latter makes use of pulse shape simulation (see B) to quantify the PSD
performance bias between $2\nu\beta\beta$ decay and 208Tl DEP events and those
coming from $0\nu\beta\beta$ decay. Indeed, the signal proxies feature much
lower energies hence different Bremsstrahlung contribution. In addition, 208Tl
DEP events have a higher probability to happen close to the detector surface
while $0\nu\beta\beta$ decays would occur homogeneously throughout the
detector bulk. A 2.3% and 1.5% absolute uncertainty has been estimated for
BEGe and IC detectors, respectively, while it amounts to 4% for coaxial
detectors. This last estimate is larger due to the difficulty to match the ANN
training performance in simulation with the data.

##  8 Summary

Nowadays, running a background-free $0\nu\beta\beta$ decay experiment is
essential to boost the $T^{0\nu}_{1/2}$ sensitivity on a reasonable time
scale. Over the past years, the Gerda collaboration demonstrated the
feasibility of such a program by upgrading its initial setup with additional
point contact detectors (BEGe and IC), a LAr veto instrumentation and lower
mass holders. As a consequence, the sensitivity linearly increased with the
exposure. The interplay between passive and active shielding techniques has
proven to be highly effective. In this paper, we focused on the ability
offered by germanium detectors to analyze the topological structure of the
recorded events. This topology, distinct for signal-like and background-like
events at a 100 ns time scale, is best scrutinized with the high frequency
based data acquisition system of Gerda. Using the Phase II dataset, we
confirmed the superior discriminating power of point contact detectors (BEGe
and IC) over the historical coaxial detectors thanks to the simple $A/E$
parameter. Their $0\nu\beta\beta$ decay PSD efficiency is 89% and 69%,
respectively for a similar background index of about $5\cdot 10^{-4}$
counts/(keV$\cdot$kg$\cdot$yr) after all cuts. We also demonstrated the high
and stable performance in LAr of the newly produced enriched inverted-coaxial
point contact germanium detectors. The LEGEND collaboration will deploy these
IC detectors for its future 76Ge $0\nu\beta\beta$ decay search program.

The total detection efficiencies for $0\nu\beta\beta$ decay are discussed and
listed for each detector in Appendix A.

###### Acknowledgements.

The Gerda experiment is supported financially by the German Federal Ministry
for Education and Research (BMBF), the German Research Foundation (DFG), the
Italian Istituto Nazionale di Fisica Nucleare (INFN), the Max Planck Society
(MPG), the Polish National Science Centre (NCN), the Foundation for Polish
Science (TEAM/2016-2/17), the Russian Foundation for Basic Research, and the
Swiss National Science Foundation (SNF). This project has received
funding/support from the European Union’s Horizon 2020 research and innovation
programme under the Marie Sklodowska-Curie grant agreements No 690575 and No
674896. This work was supported by the Science and Technology Facilities
Council, part of the U.K. Research and Innovation (Grant No. ST/T004169/1). J.
Huang and C. Ransom thank the UZH for the Postdoc and Candoc Forschungskredit
fellowships respectively. The institutions acknowledge also internal financial
support. The Gerda collaboration thanks the directors and the staff of LNGS
for their continuous strong support of the Gerda experiment.

##  Appendix A Efficiency tables per detector

The final results of Gerda on the search for $0\nu\beta\beta$ decay
Agostini:2020prl  have been obtained from data partitions, defined detector-
wise by stable data taking periods, that notably include the stability of PSD
performance. Table 3 lists for each germanium detector, before and after the
2018 upgrade, the PSD cut efficiencies $\varepsilon_{\rm{PSD}}$ as well as the
total efficiencies $\varepsilon_{0\nu\beta\beta}$ of detecting
$0\nu\beta\beta$ decays the latter entering Eq. (1) of Agostini:2020prl :

$\varepsilon_{0\nu\beta\beta}=f_{76}\cdot
f_{\rm{av}}\cdot\varepsilon_{\rm{fep}}\cdot\varepsilon_{\rm{PSD}}\cdot\varepsilon_{\rm{LAr}}.$
(5)

For the coaxial and BEGe detectors, Table S1 of the Supplementary Material of
Agostini:2019hzm  provides the 76Ge enrichment fractions $f_{76}$ and the
active volume fractions $f_{\rm{av}}$. The electron containment efficiencies
$\varepsilon_{\rm{fep}}$ of the period before the upgrade are also provided
there while a new computation has been used for the period after the upgrade
that is reported in Table 3. Table 4 shows the active volume fractions for the
IC detectors. The efficiencies $\varepsilon_{\rm{LAr}}$ of the LAr veto cut
for the two data taking periods are 0.977(1) and 0.982(1), respectively
Agostini:2020prl . Equation 5 deliberately neglects the efficiencies of muon
veto and quality cuts both being larger than 0.999. The uncertainty of all
efficiencies are incorporated via a Monte-Carlo sampling from which we
retrieve the $\varepsilon_{0\nu\beta\beta}$ central value and its standard
deviation.

Table 3: Detector-wise PSD cut efficiencies $\varepsilon_{\tiny{\textrm{PSD}}}$ and the total detection efficiencies $\varepsilon_{0\nu\beta\beta}$ for $0\nu\beta\beta$ decays used in the final analysis of Agostini:2020prl  and corresponding exposures $\mathcal{E}$. For the period after the 2018 upgrade also the electron containment efficiencies $\varepsilon_{\rm{fep}}$ are listed (see text). Quoted uncertainties account for statistics and systematics. The channels with empty entries were used in anti-coincidence only. Dec 2015-May 2018 |  | July 2018-Nov 2019
---|---|---
# | Detector | $\mathcal{E}$ | $\varepsilon_{\tiny{\textrm{PSD}}}$ | $\varepsilon_{0\nu\beta\beta}$ |  | $\varepsilon_{\tiny{\textrm{fep}}}$ | $\mathcal{E}$ | $\varepsilon_{\tiny{\textrm{PSD}}}$ |  $\varepsilon_{0\nu\beta\beta}$
label | kg.yr |  | kg.yr
0 | GD91A | 1.167 | 0.890(29) | 0.610(25) |  | 0.895(2) | 0.736 | 0.891(23) | 0.612(22)
1 | GD35B | 1.455 | 0.895(27) | 0.634(24) |  | 0.899(2) | 0.880 | 0.891(30) | 0.632(26)
2 | GD02B | 1.181 | 0.878(22) | 0.598(21) |  | 0.891(2) | 0.733 | 0.885(22) | 0.603(21)
3 | GD00B | 1.193 | 0.859(27) | 0.583(24) |  | 0.893(2) | 0.753 | 0.885(21) | 0.601(21)
4 | GD61A | 1.340 | 0.876(21) | 0.606(21) |  | 0.899(2) | 0.858 | 0.901(22) | 0.624(21)
5 | GD89B | 0.639 | 0.869(21) | 0.572(21) |  | 0.885(2) | 0.707 | 0.894(22) | 0.587(22)
6 | GD02D | - | - | - |  | - | - | - | -
7 | GD91C | 0.234 | 0.867(23) | 0.592(22) |  | 0.892(2) | 0.736 | 0.894(29) | 0.611(25)
8 | ANG5 | 5.067 | 0.760(32) | 0.486(37) |  | - | - | - | -
9 | RG1 | 3.914 | 0.670(31) | 0.465(39) |  | 0.919(2) | 2.392 | 0.661(32) | 0.463(39)
10 | ANG3 | 4.557 | 0.653(29) | 0.449(39) |  | 0.920(2) | 2.711 | 0.681(32) | 0.472(40)
11 | GD02A | 1.016 | 0.904(23) | 0.621(21) |  | 0.889(2) | 0.640 | 0.886(24) | 0.609(22)
12 | GD32B | 1.257 | 0.909(95) | 0.620(67) |  | 0.890(2) | 0.780 | 0.90(20) | 0.62(13)
13 | GD32A | 0.489 | 0.865(22) | 0.583(23) |  | 0.883(2) | 0.529 | 0.886(34) | 0.597(29)
14 | GD32C | 1.377 | 0.878(23) | 0.609(21) |  | 0.898(2) | 0.872 | 0.898(57) | 0.624(42)
15 | GD89C | 1.056 | 0.870(22) | 0.581(23) |  | 0.880(2) | 0.698 | 0.894(27) | 0.597(25)
16 | GD61C | 1.072 | 0.874(28) | 0.594(24) |  | 0.889(2) | 0.734 | 0.883(24) | 0.602(22)
17 | GD76B | 0.698 | 0.835(21) | 0.538(20) |  | 0.879(2) | 0.426 | 0.858(22) | 0.553(21)
18 | GD00C | 1.479 | 0.892(21) | 0.618(22) |  | 0.890(2) | 0.956 | 0.902(22) | 0.618(22)
19 | GD35C | 1.202 | 0.916(39) | 0.635(30) |  | 0.889(2) | 0.739 | 0.895(44) | 0.620(34)
20 | GD76C | 1.371 | 0.902(43) | 0.614(32) |  | 0.899(2) | 0.964 | 0.911(47) | 0.621(35)
21 | GD89D | 0.945 | 0.879(32) | 0.576(26) |  | 0.880(2) | 0.617 | 0.891(51) | 0.584(37)
22 | GD00D | 1.550 | 0.889(21) | 0.613(21) |  | 0.898(2) | 0.954 | 0.902(31) | 0.622(26)
23 | GD79C | 1.206 | 0.877(21) | 0.596(20) |  | 0.896(2) | 0.953 | 0.905(21) | 0.615(20)
24 | GD35A | 1.460 | 0.882(20) | 0.618(21) |  | 0.901(2) | 0.901 | 0.896(22) | 0.630(22)
25 | GD91B | 0.420 | 0.876(50) | 0.601(37) |  | 0.893(2) | 0.734 | 0.879(35) | 0.603(28)
26 | GD61B | 1.172 | 0.896(22) | 0.615(21) |  | 0.895(2) | 0.808 | 0.899(21) | 0.616(21)
27 | ANG2 | 4.750 | 0.743(32) | 0.504(41) |  | 0.924(2) | 2.751 | 0.740(31) | 0.508(40)
28 | RG2 | 4.017 | 0.630(30) | 0.400(34) |  | 0.918(2) | 2.542 | 0.640(27) | 0.412(32)
29 | ANG4 | 4.521 | 0.683(32) | 0.477(40) |  | 0.920(2) | 2.783 | 0.707(31) | 0.498(39)
30 | GD00A | 0.945 | 0.885(22) | 0.598(22) |  | 0.884(2) | 0.582 | 0.864(24) | 0.584(22)
31 | GD02C | 1.462 | 0.899(22) | 0.618(21) |  | 0.898(2) | 0.925 | 0.904(21) | 0.623(21)
32 | GD79B | 0.861 | 0.896(39) | 0.608(31) |  | 0.894(2) | 0.469 | 0.900(24) | 0.612(23)
33 | GD91D | 1.003 | 0.858(44) | 0.588(34) |  | 0.895(2) | 0.813 | 0.894(44) | 0.614(34)
34 | GD32D | 1.232 | 0.874(23) | 0.617(22) |  | 0.897(2) | 0.775 | 0.873(25) | 0.618(23)
35 | GD89A | 0.976 | 0.864(21) | 0.585(21) |  | 0.888(2) | 0.591 | 0.883(52) | 0.597(39)
36 | ANG1 | 1.812 | 0.649(33) | 0.403(36) |  | - | - | - | -
37 | IC50B | - | - | - |  | 0.920(2) | 2.213 | 0.894(18) | 0.664(14)
38 | IC48A | - | - | - |  | 0.921(2) | 2.159 | 0.910(18) | 0.676(14)
39 | IC50A | - | - | - |  | 0.910(2) | 1.814 | 0.896(17) | 0.646(14)
40 | IC74A |  |  |  |  | 0.921(2) | 2.398 | 0.900(15) | 0.652(12)
Table 4: Active volume fractions $f_{\textrm{av}}$ of IC detectors used in the
final analysis of Agostini:2020prl . The $f_{\textrm{av}}$ values are pulled
from a preliminary analysis found in Miloradovic:2019phd . Quoted
uncertainties account for statistics and systematics. IC detectors
---
# | Detector |  $f_{\textrm{av}}$
label
36 | IC48B | 0.935(5)
37 | IC50B | 0.938(4)
38 | IC48A | 0.936(7)
39 | IC50A | 0.920(4)
40 | IC74A | 0.913(6)

##  Appendix B Pulse shape simulation

The analysis of Gerda is data-driven for all PSD methods thanks to the
selection of relevant energy regions for signal proxies. However, this comes
at the expense of neglecting the specific decay dynamics at $Q_{\beta\beta}$.
For instance, the amount of Bremsstrahlung at 2039 keV is larger than at 1593
keV, hence a potential $0\nu\beta\beta$ decay event has a different $A/E$
ratio compared to 208Tl DEP. Also, 208Tl DEP events do happen on average
closer to the detector surface while $0\nu\beta\beta$ decay are homogeneously
distributed all over the detector bulk (see Fig. 21). This results in a larger
event fraction of the latter population above the p+ contact with higher $A/E$
values (see Fig. 1).

Figure 21: Events distributions across the detector bulk from simulated 208Tl
DEP (left) and $0\nu\beta\beta$ decay (right) events. The position of an event
corresponds to the barycentre of individual energy depositions.

These effects are best studied with a Monte-Carlo simulation of the Gerda
experiment and subsequent pulse shape simulation (PSS) of individual events
occurring in the germanium detectors. For this purpose, the Geant4 based
simulation package (MaGe Boswell:2010mr ) and detector configuration used for
the background model of Gerda Abramov:2019hhh  was employed to generate 208Tl
and 212Bi decays in order to reproduce the calibration energy spectrum of each
germanium detector. In addition, $2\nu\beta\beta$ decays were simulated
homogeneously in each detector, to cross-check the PSD efficiencies in
simulation with the data, as well as $0\nu\beta\beta$ decays. The energy
spectrum of the two electrons emitted in the double-$\beta$ decay was sampled
according to the distribution given in reference TRETYAK199543  implemented in
Decay0 MC:double-beta .

The individual hit positions and deposited energies are used to calculate the
corresponding induced charge flow in the detector by means of electrostatic
simulation software (ADL Bruyneel:2016zih  and Fieldgen Radford:2015 ). Each
generated pulse is then convoluted by an optimized electronics response model
(bi-quad filter Panas:2018phd ) of the Gerda setup before taking into account
a realistic gain and offset. Subsequently, waveforms from baseline events,
recorded during the physics data acquisition to estimate the random
coincidence probability, are randomly picked-up and added to the convoluted
pulses to emulate the physical signal-to-noise ratio. Finally, the post-
processing is identical to the PSD methods described in Sec. 4 and 5 to
retrieve the $A/E$, ANN and risetime classifiers.

In Fig. 22, the PSD classifier distribution of the simulation is compared to
the data. For the particular case of coaxial detectors, two approaches have
been studied for the ANN analysis. First, the training and cut values from the
data were applied to the simulation. Second, the simulated pulses were fed
into the ANN for its training and the cuts were computed such that 90% of
simulated 208Tl DEP events survive. Due to the not accurate enough modeling of
the Gerda electronics response, the latter method was found to give results in
better agreement with the data.

Figure 22: Signal proxy and $0\nu\beta\beta$ decay events PSD classifier
distributions for an inverted coaxial (top) and a coaxial (bottom) detector.
The pulse shape simulation results (colored lines) are superimposed on the
data (grey area). The histograms are normalized to their integrals.

For IC detectors, the final systematic uncertainty from these data and
simulation comparison is estimated by computing the survival fraction
difference between 208Tl DEP and $0\nu\beta\beta$ decay for each individual
detector and then average it per detector type. For coaxial detectors, the
agreement of the $2\nu\beta\beta$ decay survival fractions between data and
simulation is used to quantify the systematic uncertainty on the
$0\nu\beta\beta$ signal efficiency. The corresponding BEGe detectors
systematic uncertainty was estimated in previous studies Wagner:2017phd .

##  Appendix C ANN for $\alpha$ event rejection

In order to mitigate the high energy $\alpha$ background, an analysis based on
a second ANN was developed at the beginning of Gerda Phase II Agostini:2017iyd
; Agostini:2018tnm . Other than in the ANN analysis described in Sec. 5, the
neural network applied to $\alpha$ events was trained for different datasets
using the highest possible statistics available Kirsch:2014phd . For the
training, the $2\nu\beta\beta$ decay events were used as signal-like sample,
while $\alpha$-induced events did serve as a background-like sample. From the
training point of view, all detectors had poor statistics for high energy
events. Therefore, the training of the ANN-$\alpha$ has been performed with
the whole available background dataset to optimize the selection. Later, two
changes have been introduced into the $\alpha$ cut: i) in order to avoid any
conflict with the determination of the survival fraction for the
$2\nu\beta\beta$ decay events, calibration events at $Q_{\beta\beta}\pm 10$
keV, were provided as a signal-like sample instead, ii) input variables
between 10% and 90% of the original trace were used for the supervised
learning, what helped to reduce energy dependence of the ANN-$\alpha$ cut.

Even though the statistics in the training data sets was much smaller compared
to the MSE-based training, a very efficient – and even better – separation
from the signal event class has been achieved. The discriminating parameter
was set such that for physics data with $E>3500$ keV a survival fraction of
10% yielded, while retaining, on average, $\sim$95% of 208Tl DEP events in the
calibration data. Only two detectors, RG1 and ANG4, featured slightly smaller
208Tl DEP survival fractions of about 92%. The corresponding statistical
uncertainty on the 10% cut on $\alpha$ events was in the range of (2 – 3)%.

For reliability and technical reasons, the original TMVA algorithm used in
Kirsch:2014phd ; Agostini:2017iyd ; Agostini:2018tnm  called TMlpANN (own ROOT
ANN TMultiLayerPercetron class) was replaced by the recommended and supported
MLPBFGS (ANN Multilayer Perceptron class) in 2018, both using the same
Broyden-Fletcher-Goldfarb-Shanno algorithm. Two existing caveats persisted in
the analysis. First, the low statistics of the background sample ($\sim$50–100
events), was not sufficient to test the ANN from the beginning. The new
algorithm required to split the data in training and testing samples hence an
even lower background statistics, which led to non-converging training in some
detectors. Second, few events with very fast signals, i.e. short risetime,
were found to survive the ANN-$\alpha$ cut in the region of interest while
they were thought to be $\alpha$ events with high confidence level. This put
in the question of reliability and robustness of this analysis approach, hence
the decision to switch to the simpler risetime cut at the cost of a $\sim$10%
loss in signal efficiency.

##  Appendix D Projective likehood analysis

An alternative analysis strategy was pursued for the coaxial detectors using a
projective likelihood (PL) to remove both MSEs and $\alpha$ background at
once. In this analysis the high-frequency waveform from the preamplifier
output was smoothed by averaging every 5 adjacent samples, differentiated and
then the maximum current value was found along with its corresponding time
($t_{0}$). For the analysis, 15 amplitudes before and 15 after $t_{0}$ were
extracted from the original trace (31 in total). In order to further reduce
the number of input variables sums of 4 neighboring amplitudes were calculated
establishing 5 new input parameters ($\Sigma_{1}$ – $\Sigma_{5}$, $\Sigma_{5}$
is a sum of 3 amplitudes) – see details in Panas:2018phd . In order to avoid
energy dependence, after the baseline subtraction, the pulses were normalized
with respect to their energies.

For training and definition of the cuts, calibration data were used. The
Compton edge events of the 2615 keV peak (region between 2350 keV and 2375
keV) were chosen as SSEs (signal), and multiple Compton scattered events
(region between 2450 keV and 2550 keV) as MSEs (background). The cut was
defined requesting 80% survival fraction for the DEP events. This allows to
compare the present analysis with the corresponding one from the Gerda Phase I
data. No dedicated cut for high energy events was defined.

The stability of the input data was monitored using the time distribution of
$\Sigma_{2}$. If instabilities (abrupt change of the input variable
distribution) were observed, the data in the affected channels was divided
into different (stable) sub-periods and analyzed separately. As a consequence,
the cut for a given survival probability, was time dependent. Such a change
for ANG5 has been observed between some runs, otherwise the distribution of
$\Sigma_{2}$ was stable within less than 4%.

The survival probability of the $0\nu\beta\beta$ decay events estimated from
the simple cut based on the PL classifier is $(65.5\pm 13.3)\%$, while the
background in the region of interest is reduced by 56%. Although, the signal
efficiency is relatively low with conservatively estimated systematic
uncertainty, PL allows also to eliminate efficiently the high energy events,
about 87% of them are rejected. An advantage of PL is the training performed
only on $\gamma$-rays from the calibration runs.

This method was not retained for the final $0\nu\beta\beta$ decay search
analysis due to lower signal efficiency and higher background in the region of
interest.

## References

  * (1) V. Cirigliano et al.,  A neutrinoless double beta decay master formula from effective field theory,  JHEP 12, 097 (2018).
  * (2) K.H. Ackermann et al. (Gerda collaboration),  The GERDA experiment for the search of $0\nu\beta\beta$ decay in 76Ge,  Eur. Phys. J. C 73(3), 2330 (2013).
  * (3) M. Agostini et al. (Gerda collaboration),  Upgrade for Phase II of the GERDA experiment,  Eur. Phys. J. C 78(5), 388 (2018).
  * (4) M. Agostini et al. (Gerda collaboration),  Production, characterization and operation of 76Ge enriched BEGe detectors in GERDA,  Eur. Phys. J. C 75(2), 39 (2015).
  * (5) M. Agostini et al. (Gerda collaboration),  Characterization of 30 76Ge enriched Broad Energy Ge detectors for GERDA Phase II,  Eur. Phys. J. C 79(11), 978 (2019).
  * (6) R.J. Cooper, D.C. Radford, P.A. Hausladen and K. Lagergren,  A novel HPGe detector for gamma-ray tracking and imaging,  Nucl. Instrum. Meth. A 665, 25 (2011).
  * (7) M. Agostini et al. (Gerda collaboration),  Characterization of inverted coaxial 76Ge detectors in GERDA for future double-$\beta$ decay experiments,  Eur. Phys. J. C 81(6), 505 (2021).
  * (8) I. Abt et al.,  Background reduction in neutrinoless double beta decay experiments using segmented detectors: A Monte Carlo study for the GERDA setup,  Nucl. Instrum. Meth. A 570, 479–486 (2007).
  * (9) M. Agostini et al. (Gerda collaboration),  Background-free search for neutrinoless double-$\beta$ decay of 76Ge with GERDA,  Nature 544, 47 (2017).
  * (10) M. Agostini et al. (Gerda collaboration),  Pulse shape discrimination for GERDA Phase I data,  Eur. Phys. J. C 73(10), 2583 (2013).
  * (11) M. Agostini et al. (Gerda collaboration),  Probing Majorana neutrinos with double-$\beta$ decay,  Science 365, 1445 (2019).
  * (12) P. N. Luke, F. S. Goulding, N. W. Madden and R. H. Pehl,  Low capacitance large volume shaped-field germanium detector,  IEEE Trans. Nucl. Sci. 36, 926 (1989).
  * (13) Zhong He,  Review of the Shockley–Ramo theorem and its application in semiconductor gamma-ray detectors,  Nucl. Instrum. Meth. A 463, 250 (2001).
  * (14) B. Bruyneel, B. Birkenbach and P. Reiter,  Pulse shape analysis and position determination in segmented HPGe detectors: The AGATA detector library,  Eur. Phys. J. A 52(3), 70 (2016).
  * (15) M. Agostini, L. Pandola and P. Zavarise,  Off-line data processing and analysis for the GERDA experiment,  J. Phys. Conf. Ser. 368, 012047 (2012).
  * (16) M. Agostini, L. Pandola, P. Zavarise and O. Volynets,  GELATIO: A General framework for modular digital analysis of high-purity Ge detector signals,  JINST 6, P08013 (2011).
  * (17) M. Agostini et al. (Gerda collaboration),  Improvement of the energy resolution via an optimized digital signal processing in GERDA Phase I,  Eur. Phys. J. C 75(6), 255 (2015).
  * (18) M. Agostini et al. (Gerda collaboration),  Calibration of the GERDA experiment,  Eur. Phys. J. C 81(8), 682 (2021).
  * (19) M. Agostini et al. (Gerda collaboration),  Modeling of GERDA Phase II data,  JHEP 03, 139 (2020).
  * (20) M. Agostini et al. (Gerda collaboration),  Limits on uranium and thorium bulk content in GERDA Phase I detectors,  Astropart. Phys. 91, 15 (2017).
  * (21) M. Agostini et al. (Gerda collaboration),  The background in the $0\nu\beta\beta$ experiment GERDA,  The European Physical Journal C 74(4), 2764 (2014).
  * (22) D. Budjáš et al.,  Pulse shape discrimination studies with a Broad-Energy Germanium detector for signal identification and background suppression in the GERDA double beta decay experiment,  JINST 4, P10007 (2009).
  * (23) https://root.cern/manual/tmva/.
  * (24) T. Comellato, M. Agostini and S. Schönert,  Charge-carrier collective motion in germanium detectors for $\beta\beta$-decay searches,  Eur. Phys. J. C 81(1), 76 (2021).
  * (25) A. Lazzaro,  Signal processing and event classification for a background free neutrinoless double beta decay search with the GERDA experiment,  PhD thesis Technische Universität München (2019),  http://mediatum.ub.tum.de/node?id=1507626.
  * (26) M. Agostini et al. (Gerda collaboration),  Final Results of GERDA on the Search for Neutrinoless Double-$\beta$ Decay,  Phys. Rev. Lett. 125, 252502 (2020).
  * (27) M. Miloradovic,  Characterisation of Inverted Coaxial Detectors and Calibration Source Production for the GERDA Experiment,  PhD thesis Universität Zurich (2020),  https://www.physik.uzh.ch/dam/jcr:94fa6df0-e441-413a-a31d-59f89e375970/.
  * (28) M. Boswell et al.,  MaGe \- a Geant4-based Monte Carlo Application Framework for Low-background Germanium Experiments,  IEEE Trans. Nucl. Sci. 58, 1212 (2011).
  * (29) V.I. Tretyak and Yu.G. Zdesenko,  Tables of double beta decay data,  Atomic Data and Nuclear Data Tables 61(1), 43–90 (1995).
  * (30) O. A. Ponkratenko, V. I. Tretyak and Yu. G. Zdesenko,  Event generator DECAY4 for simulating double-beta processes and decays of radioactive nuclei,  Physics of Atomic Nuclei 63(7), 1282–1287 (2000).
  * (31) D. Radford,  MAJORANA Detector Signal Generator (MJD SIGGEN),  http://radware.phy.ornl.gov/MJ/ (2015).
  * (32) K. Panas,  Development of Pulse Shape Discrimination Methods as Tool for Background Suppression in High Purity Germanium Detectors used in the GERDA Experiment,  PhD thesis Jagiellonian University in Cracow (2018),  https://fais.uj.edu.pl/documents/41628/139367252/ KPanas_PhDThesis.pdf.
  * (33) V. Wagner,  Pulse Shape Analysis for the GERDA Experiment to set a new limit on the half-life of $0\nu\beta\beta$ decay of 76Ge,  PhD thesis Universität Heidelberg (2017),  http://www.ub.uni-heidelberg.de/archiv/22621.
  * (34) M. Agostini et al. (Gerda collaboration),  Improved Limit on Neutrinoless Double-$\beta$ Decay of 76Ge from GERDA Phase II,  Phys. Rev. Lett. 120, 132503 (2018).
  * (35) A. Kirsch,  Search for the neutrinoless double $\beta$-decay in GERDA Phase I using a Pulse Shape Discrimination technique,  PhD thesis Universität Heidelberg (2014),  http://www.ub.uni-heidelberg.de/archiv/17149.


where the jump $\Pijump$ can be identifed as

\begin{align}
\Pijump  &= \frac{\Dx_{L}+\Dx_{R}}{2}\left\{\rho \dx \phi\right\}
\Pijump u^* &= \frac{\Dx_{L}+\Dx_{R}}{2}\left\{\mom \dx \phi\right\}.
\end{align}

At this point, the functions
$\left\{\rho \dx \phi\right\}$ and \(\left\{\mom \dx \phi\right\}\)
are still yet to be specified. Let us consider the constraint 3: if it is satisfied then for a state that verifies (<ref>) the
jumps $\Pijump$ and $\Pijump u^*$ necessarily take the value
\(\Pijump = - (\Pi_R - \Pi_L)\) and $\Pijump u^*=0$ . A simple choice that fulfills this requirement is

\begin{align}
\{
\rho \dx\phi
\}
(\rho_L +\rho_R)
\frac{\phi_R-\phi_L}{\Delta x_L +\Delta x_R}
\{
\rho \dx\phi
\}
(\rho_L +\rho_R)
\frac{\phi_R-\phi_L}{\Delta x_L +\Delta x_R}
\label{eq: discrete source term definition}
\end{align}

Relations (<ref>) and (<ref>) give a complete definition of the approximate Riemann solver (<ref>).
This solver yields a definition for the conservative numerical flux
\(\acousticop_\Delta(
\bU_{L},\Pi_{L},\phi_{L}
\)
and a source term discretization (located at the interface)
\(\sourceterm_\Delta(
\bU_{L},\Pi_{L},\phi_{L}
\)
thanks to the consistency in the integral sense. We get

\begin{align}
\acousticop_\Delta(
\bU_{L},\Pi_{L},\phi_{L}
\frac{\acousticop(\bU_R,\Pi_R)+\acousticop(\bU_L,\Pi_L)}{2}
-\frac{a}{2\rho_L}(\bU_L^* - \bU_L)
-\frac{a}{2\rho_R}(\bU_R - \bU_R^*)
\sourceterm_\Delta(
\bU_{L},\Pi_{L},\phi_{L}
[0,-\{\rho\dx\phi\},-\{\rho u\dx\phi\}]^T
\end{align}

so that for two neighbouring states $(\bU_{j}^n,\Pi_{j}^n,\phi_{j})$
and $(\bU_{j+1}^n,\Pi_{j+1}^n,\phi_{j+1})$ across the cell interface
$j+1/2$ that separates the cell $j$ and the cell $j+1$, the numerical
conservative flux $(0,\Pi^*_{j+1/2},\Pi^*_{j+1/2} u^*_{j+1/2})$ is defined by

\begin{equation}
(0,\Pi^*_{j+1/2},\Pi^*_{j+1/2} u^*_{j+1/2}) =
\acousticop_\Delta(
\bU^n_{j},\Pi^n_{j},\phi_{j}
\label{pressure system flux definition}
\end{equation}

and the discrete souce term $\sourceterm_{j}$ within the cell $j$ is
given by

\begin{align}
\sourceterm_{j}
\frac{\Delta x_{j+1/2}}{2\Delta x_j} \sourceterm_{j+1/2}
\frac{\Delta x_{j-1/2}}{2\Delta x_j} \sourceterm_{j-1/2}
\sourceterm_{j+1/2}
\sourceterm_\Delta(
\bU^n_{j},\Pi^n_{j},\phi_{j}
\end{align}

Let us now give some properties of the approximate Riemann solver.
Let us note
e_k^* = E_k^* - (u_k^*)^2/2
the following lemma is a direct consequence of (<ref>) that exhibits
a reminiscent property associated with the Riemann invariants associated of
the system (<ref>) when $\sourceterm=\vb{0}$.

\begin{align}
e_k^* - \frac{(\Pi_k^*)^2}{2a^2}
e_k - \frac{(\Pi_k)^2}{2a^2}
\Tau_k^* + \frac{\Pi_k^*}{a}
\Tau_k + \frac{\Pi_k}{a}
\label{eq: riemann invariants across acoustic waves}
\end{align}

The following positivity result is a direct consequence of (<ref>).

* If $a$ is chosen large enough then $\Tau_L^*>0$ and $\Tau_R^*>0$.

* $\Tau_L^*>0$ and $\Tau_R^*>0$ is equivalent to
$u_L -a\Tau_L = u_L -a/\rho_L < u^* < u_R +a\Tau_R = u_R +a/\rho_R$.

Following the lines of [16], we first prove two preliminary stability-related results. The differences from Lemma 1 of [16] is that the Riemann states we are dealing with here depend on the $\Pijump$ terms and that the specific volume we use is $\Tau$ instead of ${1}/{\rho}$ (that are different in the sub-system framework).  However, the proof turns out to be almost identical.

Consider the intermediate states defined by (<ref>).

and noting $s_k=s^{\EOS}(\Tau_k,s_k)$, we have

\begin{equation}
e_{k}^{*}-e^{\EOS}\left(\Tau_{k}^{*}, s_{k}\right)-\frac{\left(p^{\EOS}\left(\Tau_{k}^{*}, s_{k}\right)-\Pi^{*}_k\right)^{2}}{2 a^{2}} \geq 0,
\label{lemme 1}
\end{equation}

with $e^*_k = E^*_k - \frac{{u^*_k}^2}{2}$.

We only describe the case $k=R$. Consider the function:

\begin{multline}
\chi(\Tau)=e^{\EOS}\left(\Tau, s_{R}\right)-\frac{p^{\EOS}\left(\Tau, s_{R}\right)^{2}}{2 a^{2}}-e^{\EOS}\left(\Tau_{R}^{*}, s_{R}\right) +\frac{p^{\EOS}\left(\Tau_{R}^{*}, s_{R}\right)^{2}}{2 a^{2}} \\
+p^{\EOS}\left(\Tau_{R}^{*}, s_{R}\right)\left(\Tau+\frac{p^{\EOS}\left(\Tau, s_{R}\right)}{a^{2}}-\Tau_{R}^{*}-\frac{p^{\EOS}\left(\Tau_{R}^{*}, s_{R}\right)}{a^{2}}\right) .
\end{multline}

One can check that $\chi'(\Tau)=\left(p^{\EOS}\left(\Tau_{R}^{*}, s_{R}\right)-p^{\EOS}\left(\Tau, s_{R}\right)\right)\left(1-\rho^{2} c^{2}\left(\Tau, s_{R}\right) / a^{2}\right)$. We have $\partial_\Tau p< 0$ from <ref>, we also assume that $a$ is large enough. We have two different cases:

\begin{equation}
\begin{aligned}
\Tau^*_R < \Tau < \Tau_R & \implies  \chi'(\Tau) > 0 & \implies  \chi(\Tau^*_R)<\chi(\Tau)<\chi(\Tau_R).\\
\Tau^*_R > \Tau > \Tau_R& \implies   \chi'(\Tau) < 0&\implies \chi(\Tau^*_R)<\chi(\Tau)<\chi(\Tau_R).
\end{aligned}
\end{equation}

As $\chi(\Tau_R^*) =0$, we have $\chi(\Tau_R)>0$ , in both cases.
Accounting for (<ref>), we get

\begin{align}
\chi(\Tau_R) &=  e^{\EOS}\left(\Tau_R, s_{R}\right)-\frac{p^{\EOS}\left(\Tau_R, s_{R}\right)^{2}}{2 a^{2}}-e^{\EOS}\left(\Tau_{R}^{*}, s_{R}\right) +\frac{p^{\EOS}\left(\Tau_{R}^{*}, s_{R}\right)^{2}}{2 a^{2}}
\nonumber
\\
+p^{\EOS}\left(\Tau_{R}^{*}, s_{R}\right)\left(\Tau_R+\frac{p^{\EOS}\left(\Tau_R, s_{R}\right)}{a^{2}}-\Tau_{R}^{*}-\frac{p^{\EOS}\left(\Tau_{R}^{*}, s_{R}\right)}{a^{2}}\right)
\nonumber
\\
e_R^*-\frac{\qty(\Pi_R^*)^{2}}{2 a^{2}}-e^{\EOS}\left(\Tau_{R}^{*}, s_{R}\right) +\frac{p^{\EOS}\left(\Tau_{R}^{*}, s_{R}\right)^{2}}{2 a^{2}}
+p^{\EOS}\left(\Tau_{R}^{*}, s_{R}\right)
\left(
\frac{\Pi_R^*}{a^{2}}-\frac{p^{\EOS}\left(\Tau_{R}^{*}, s_{R}\right)}{a^{2}}
\right)
\nonumber
\\
&=   e_{R}^{*}-e^{\EOS}\left(\Tau_{R}^{*}, s_{R}\right)-\frac{\left(p^{\EOS}\left(\Tau_{R}^{*}, s_{R}\right)-\Pi^{*}_R\right)^{2}}{2 a^{2}}
\end{align}

Similar lines can be used for $k=L$.

We present a result concerning the behavior of the numerical scheme in the low Mach regime defined in
section <ref>: we consider a one-dimensional smooth solution of
the pressure subsystem (<ref>)
such that \(\partial_{\tilde{x}}\tilde{p} + \tilde{\rho}\widetilde{(\dx\phi)} = O(\mach^2)\).
Then, we proceed as in [32] by evaluating the truncation
error (in the sense of the Finite Difference) obtained by substituting these low Mach flow parameters into
the finite volume update formula derived from the fluxes (<ref>).
We obtain the following results.

In the low Mach regime, the rescaled discretization of the pressure system is consistent with

\begin{align}
\partial_{\tilde{t}} \tilde{\rho}&=0, & \partial_{\tilde{t}} (\tilde{\rho}\tilde{u})+\frac{1}{\mach^{2}} (\partial_{\tilde{x}}\tilde{p} + \tilde{\rho}\widetilde{(\dx\phi)} )
O(\Delta \tilde{t})+O\left(\frac{\Delta \tilde{x}}{\mach}\right)
\partial_{\tilde{t}} (\tilde{\rho}\tilde{E})+\partial_{\tilde{x}}(\tilde{p} \tilde{u}) &=O(\Delta \tilde{t})+O(\mach\Delta \tilde{x}).
\label{pressure subsystem low mach equivalent eq}
\end{align}

If one performs a similar evaluation for the full FSLP scheme, one can see that the truncation error term $O\left(\frac{\Delta \tilde{x}}{\mach}\right)$ that appears in the momentum equation of
(<ref>) will be the only error term whose magnitude is not uniform with respect to $\mach$. Similarly as in [28, 32, 1, 33], this
truncation error term can be traced back to the non-centered part of $\Pi_{j+1/2}^*$. To tackle this issue, we adopt the modification used in [32, 1] by replacing
$\Pi_{j+1/2}^*$ with

\begin{equation}
\Pi_{j+1 / 2}^{*, \theta}=\frac{1}{2}\left(\Pi_{j}^{n}+\Pi_{j+1}^{n}\right)-\theta_{j+1 / 2} \frac{a_{j+1 / 2}}{2}\left(u_{j+1}^{n}-u_{j}^{n}\right),
\label{p^* corrigé}
\end{equation}

where $\theta_{j+1/2}\in[0,1]$. This results in the update relation (<ref>) that is a finite approximation of (<ref>) with the flux definition (<ref>). We will see in
<ref> how this resulting modified flux can still be associated with an Approximate Riemann solver.

§ ALL-REGIME APPROXIMATE RIEMANN SOLVER FOR THE PRESSURE SUBSYSTEM

Following similar lines as in [32]: although the modified pressure scheme (<ref>) is defined as a flux scheme, it is possible  to find
an approximate Riemann solver \((\appRP^\theta,\PiRP^\theta,\TauRP^\theta)\)
that enables to retrieve the numerical flux
\({\acousticop}_{j+1/2}^\theta = (0,\Pi^{*,\theta}_{j+1/2},\Pi^{*,\theta}_{j+1/2} u^*_{j+1/2})\). We suppose that  \((\appRP^\theta,\PiRP^\theta,\TauRP^\theta)\)
has the same structure as  \((\appRP,\PiRP,\TauRP)\), we consider

\begin{equation}
\qty(\frac{x-\bar{x}}{t} ; \bU_{L},\Pi_L,\Tau_L,\phi_L, \bU_{R},\Pi_R,\Tau_R,\phi_R)
\begin{cases}
& \text { if  \(\frac{x-\bar{x}}{t} \leq-\frac{a}{\rho_L}\),}
\\
(\bU_{L}^{*,\theta},\Pi_L^{*,\theta},\Tau_L^{*,\theta},\Phi_L), & \text { if \(-\frac{a}{\rho_L} <\frac{x-\bar{x}}{t} \leq 0\)},
\\
(\bU_{R}^{*,\theta},\Pi^{*,\theta}_R,\Tau_R^{*,\theta},\Phi_R), & \text { if  \(0<\frac{x-\bar{x}}{t} \leq \frac{a}{\rho_R}\)},
\\
(\bU_{R},\Pi_R,\Tau_R,\Phi_R), & \text { if  \(\frac{a}{\rho_R} <\frac{x-\bar{x}}{t}\)},
\end{cases}
\label{eq: modified approx riemann solver def}
\end{equation}

where $\Pi_k$, $\Tau_k$ and $\Phi_k$
verify (<ref>), $k=L,R$.
The states $(\bU_{k}^{*,\theta}, \Pi^{*,\theta}_k, \Tau_k^{*,\theta})$, $k=L,R$ are yet to be defined.
First, we impose that \((\appRP,\PiRP,\TauRP)\) verifies the consistency in the integral sense

\begin{equation}
\begin{bmatrix}
\frac{1}{\alpha}\acousticop\left(\bU_{R}\right)
\frac{1}{\alpha}\acousticop\left(\bU_{L}\right)
\\
\frac{1}{\alpha} a^2(u_R - u_L)
\\
-(\frac{1}{\alpha} u_R - \frac{1}{\alpha} u_L)
\end{bmatrix}
\begin{bmatrix}
\bU_{L}^{*,\theta}-\bU_{L}
\\
\\
\end{bmatrix}
\frac{a}{\alpha \rho_R}
\begin{bmatrix}
\bU_{R}-\bU_{R}^{*,\theta}
\\
\\
\end{bmatrix}
\frac{\Dx_{L}+\Dx_{R}}{2}
\begin{bmatrix}
\frac{1}{\alpha}\left\{\sourceterm\right\}
\\
\\
\end{bmatrix}.
\label{intg form theta}
\end{equation}

We then enforce that the numerical flux resulting from (<ref>)
is $\acousticop_\Delta^{\theta}$, which boils down to require that

\begin{equation}
\begin{bmatrix}
\frac{1}{\alpha}\acousticop_\Delta^\theta
\\
\frac{1}{\alpha}a^2 u_\Delta^\theta
\\
\end{bmatrix}	\!\!
\bU_{L},\Pi_{L},\phi_{L}
\! = \!
\begin{bmatrix}
\acousticop(\bU_R,\Pi_R)\!+\!\acousticop(\bU_L,\Pi_L)
\\
a^2u_R + a^2 u_L
\\
-(u_R + u_L)
\end{bmatrix}
\!\!
\begin{bmatrix}
\bU_L^{,\theta} - \bU_L
\\
(\rho\Pi)_L^{,\theta} - (\rho\Pi)_L
\\
(\rho\Tau)_L^{,\theta} - (\rho\Tau)_L
\end{bmatrix}
\!\!
\begin{bmatrix}
\bU_R - \bU_R^{,\theta}
\\
(\rho\Pi)_R - (\rho\Pi)_R^{,\theta}
\\
(\rho\Tau)_R - (\rho\Tau)_R^{,\theta}
\end{bmatrix}
\label{eq: flux eq theta}
\end{equation}

Choosing $\rho_k^{*,\theta} = \rho_k$, $k=L,R$, relation (<ref>)
and (<ref>) provide a linear system with respect to
$u_k^{*,\theta}$, $\Pi_k^{*,\theta}$, $\Tau_k^{*,\theta}$ and $E_k^{*,\theta}$, $k=1,2$
whose solution is

\begin{align}
\rho_L^*&=\rho_L
\rho_R^*&=\rho_R
E^{*,\theta}_L & = E_L^* -(1-\theta)\frac{u_R-u_L}{2}u^*
, &
E^{*,\theta}_R & = E_R^* +(1-\theta)\frac{u_R-u_L}{2}u^*
, \\
u^{*,\theta}_L &= u^* - (1-\theta)\frac{u_R-u_L}{2}
u^{*,\theta}_R &= u^* + (1-\theta)\frac{u_R-u_L}{2}
\Pi_{L}^{*,\theta} &=\Pi_{L}^{*}
, &
\Pi_{R}^{*,\theta} &=  \Pi_{R}^{*}
, \\
\Tau_{L}^{*,\theta} &=\Tau_L^*
, &
\Tau_{R}^{*,\theta} &=\Tau_R^*.
\end{align}

We now turn to positivity-preserving related properties. Let us note
$e_k^{*,\theta}= E_k^{*,\theta} - {u_k^{*,\theta}}^2/2$, we have the following result.

Assuming again that $a$ is large enough, we have

\begin{equation}
e_{k}^{*,\theta}-e^{\EOS}\left(\Tau_{k}^{*,\theta}, s_{k}\right)
-\frac{\left(p^{\EOS}\left(\Tau_{k}^{*,\theta}, s_{k}\right)-\Pi^{*,\theta}_k\right)^{2}}{2 a^{2}}
\frac{(1-\theta)^2(u_R-u_L)^2}{8}\geq 0.
\label{lemme 2}
\end{equation}

Let us consider the case $k=R$, by (<ref>)
we get

\begin{align}
e_R^{*,\theta}- e_R^{*} &=  E_R^{*,\theta}-E_R^{*} -\frac{1}{2}( {u_R^{*,\theta}}^2 - {u_R^{*}}^2 )
\nonumber
\\
&=(1-\theta)\frac{u_R-u_L}{2} u^*
\frac{1}{2}
\qty(
(u^{*})^2 +u^{*}(1-\theta)(u_R-u_L)
+(1-\theta)^2 \frac{(u_R-u_L)^2}{4}
- (u^{*})^2
\nonumber
\\
&=-\frac{1}{8}(1-\theta)^2 (u_R-u_L)^2.
\end{align}

Using (<ref>), we obtain

\begin{multline}
e_{R}^{*,\theta}-e^{\EOS}\left(\Tau_{R}^{*,\theta}, s_{R}\right)
e^{\EOS}\left(\Tau_{R}^{*,\theta}, s_{R}\right)
-\frac{1}{8}(1-\theta)^2 (u_R-u_L)^2
e^{\EOS}\left(\Tau_{R}^{*,\theta}, s_{R}\right)
\\
\geq		\nonumber
-\frac{1}{8}(1-\theta)^2 (u_R-u_L)^2
\frac{
\qty(
p^{\EOS}\left(\Tau_{k}^{*,\theta}, s_{k}\right)-\Pi^{*,\theta}_k
)^{2}}{2 a^{2}}.
\end{multline}

Similar lines can be used for the case $k=L$.

The relation (<ref>) highlights the role of the inequality

\begin{equation}
\frac{1}{2 a^{2}}\left(p^{\EOS}(\Tau_{k}^{*, \theta}, s_{k})-\Pi^{*}_k\right)^{2}
\frac{(1-\theta)^{2}\left(u_{R}-u_{L}\right)^{2}}{8} \geq 0, \quad k=L, R
\label{inegalité sur Co}
\end{equation}

in obtaining stability properties for the modified scheme. We have the following proposition.

Let us note: $s_k^{*,\theta} = s^\EOS(\Tau_k^{*,\theta},e_k^{*,\theta})$, if (<ref>) is satisfied, then

* the modified approximate Riemann solver (<ref>) preserves the positivity of the internal energy, that is to say: $e_k^{*,\theta}>0$, $k=R,L$,

* the modified approximate Riemann solver (<ref>) verifies
$s_k^{*,\theta} \geq s_k $, $k=R,L$,

* the modified approximate Riemann solver (<ref>) is entropy satisfying in the sense that

\begin{equation}
-a(s^{*,\theta}_L -s_L)
+a(s_R - s^{*,\theta}_R )
\geq 0
\label{eq: approx riemann solver entropy inequality}
\end{equation}

If (<ref>) is satisfied, then
$e_{k}^{*,\theta} \geq e^{\EOS}\left(\Tau_{k}^{*,\theta}, s_{k}\right)$. By the assumption on the EOS, we have that $e_{k}^{*,\theta}>0$. Now, considering a fixed $\overline{\Tau}>0$, by (<ref>) we know that $e'\mapsto s^\EOS(\overline{\Tau},e')$ is increasing, thus we deduce that
$ s^\EOS(\Tau_k^{*,\theta},e^\EOS(\Tau_k^{*,\theta},s_k^{*,\theta})) = s_k^{*,\theta}
\geq
s^\EOS(\Tau_k^{*,\theta},e^\EOS(\Tau_k^{*,\theta},s_k)) = s_k
$, $k=L,R$. This implies (<ref>).

§ EIGENSTRUCTURE OF THE OFF-EQUILIBRIUM (<REF>\(\RELAXPARAM\!=\!0\))

We propose in this section to study the eigenstructure of the
relaxation system (<ref>\(\relaxparam\!=\!0\)).

Let us first express the acoustic part of (<ref>\(\relaxparam\!=\!0\)) using
a change of variables: accounting for
$e^P = E^P - (u^P)^2/2$, the evolution equations for $E^P$, for $\Pi^P$ and $\Tau^P$
in (<ref>) yield

\begin{align}
\dt(\rho^P e^P) + 2\Pi^P \dx u^P
2\dx u^P
\dt(\rho^P \Pi^P / a^2)
\label{eq: acoustic e eq}
\end{align}

We thus obtain the stationary equations

\begin{align}
\dt\qty
[e^P - \frac{(\Pi^P)^2}{2a^2}]
\dt\qty[\Tau^P + \frac{\Pi^P}{a^2}] &= 0
\label{eq: acoustic stationnary eq}
\end{align}

So now the acoustic subsystem (<ref>) takes the simple form

\begin{align}
\dt \phi & \!=\!0
\dt \rho^P &\!\!=\!0
\dt\qty
[e^P \!-\! \frac{(\Pi^P)^2}{2a^2}]
\dt (\rho^P u^P) + 2\dx\Pi^P +2\rho^P\dx \phi^P&\!\!=\!0
\dt (\rho^P\Pi^P) + 2 a^2 \dx u^P &\!\!=\!0
\dt\qty[\Tau^P + \frac{\Pi^P}{a^2}] &= 0
\label{eq: acoustic part v1}
\end{align}

We now turn to the advection part of (<ref>\(\relaxparam\!=\!0\)): the subsystem (<ref>) takes the simple form

\begin{align} %[left=\empheqlbrace]{align}
\dt \rho^A + \dx (2\rho^A u^P)&=0,
\dt\qty[
\rho^A\Tau^A-\frac{\rho^P\Pi^P}{a^2}
\dt b^A + 2 u^P\dx b^A&=0
b^A\in\{u^A, E^A, \Pi^A \}
\label{eq: advection part v1}
\end{align}

Therefore if we set

\begin{equation}
\bW^T = \qty [u^P , \Pi^P , \rho^P, \phi ,
e^P-\frac{(\Pi^P)^2}{2 a^2},
\Tau^P+\frac{\Pi^P}{a^2},
\rho^A\Tau^A-\frac{\rho^P\Pi^P}{a^2},
u^A , \Pi^A , E^A, \rho^A],
\end{equation}

we can see that (<ref>\(\relaxparam\!=\!0\)) can be recast
into the following quasilinear system

\begin{align}
\dt\bW + \diagomatrix(\bW)\dx\bW
\diagomatrix(\bW)
\qty[
\begin{array}{ccccccccccc}
0                     & \dfrac{2}{\rho^P} & 0 & 2 & 0 & 0 & 0 & 0 & 0 & 0 & 0
\\
\dfrac{2 a^2}{\rho^P} &      0            & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0
\\
0                     &      0            & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0
\\
0                     &      0            & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0
\\
0                     &      0            & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0
\\
0                     &      0            & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0
\\
0                     &      0            & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0
\\
0                     &      0            & 0 & 0 & 0 & 0 & 0 & 2 u^P & 0 & 0 & 0
\\
0                     &      0            & 0 & 0 & 0 & 0 & 0 & 0 & 2 u^P & 0 & 0
\\
0 & 0 & 0 & 0 & 0 & 0 & 0 &  0 & 0 & 2 u^P & 0
\\
2\rho^A & 0 & 0 & 0 & 0 &  0 & 0 & 0 & 0 & 0 & 2 u^P
\end{array}
\label{eq: overall system non-conservative}
\end{align}

It is then straightforward to see that the eigenvalues of $\diagomatrix(\bW)$ are
$2 u^P$ (with an algebraic multiplicity 4), $0$ (with an algebraic multiplicity 5) and
$\pm 2 a/\rho^P$.

The eigenvectors
and $\br_{\pm}$ that are respectively associated with $0$, $2 u^P$ and $\pm 2 a/\rho^P$ are

\begin{align}
\br_0^{(1)}
\br_0^{(2)}
\qty[0,-\rho^P,0,1,0,0,0,0,0,0,0]^T
\br_0^{(3)}
\br_0^{(4)}
\\
\br_0^{(5)}
% 				,
\\
\br_{u^P}^{(1)}
\br_{u^P}^{(2)}
\br_{u^P}^{(3)}
\br_{u^P}^{(4)}
\br_{+}
\qty
[1,a,0,0,0,0,0,0,0,0,-\frac{\rho^A\rho^P}{\rho^P u^P-a}]^T
\br_{-}
\qty
[1,-a,0,0,0,0,0,0,0,0,-\frac{\rho^A\rho^P}{\rho^P u^P+a}]^T
\end{align}

so that (<ref>) is hyperbolic and only involves linearly degenerate fields.