11institutetext: Max-Planck-Institut für Kernphysik, Saupfercheckweg 1, 69117
Heidelberg, Germany 22institutetext: PreussenElektra GmbH, Osterende,
Brokdorf, Germany
E-mail address<EMAIL_ADDRESS>
# Pulse shape discrimination for the CONUS experiment in the keV and sub-keV
regime

H. Bonet 11 A. Bonhomme 11 C. Buck 11 K. Fülber 22 J. Hakenmüller 11 J.
Hempfling 11 J. Henrichs 11 G. Heusser 11 M. Lindner 11 W. Maneschg 11 T. Rink
11 E. Sánchez García 11 J. Stauber 11 H. Strecker 11 R. Wink 22

###### Abstract

Point-contact p-type high-purity germanium detectors (PPC HPGe) are
particularly suited for detection of sub-keV nuclear recoils from coherent
elastic scattering of neutrinos or light dark matter particles. While these
particles are expected to interact homogeneously in the entire detector
volume, specific classes of external background radiation preferably deposit
their energy close to the semi-active detector surface, in which diffusion
processes dominate that subsequently lead to slower rising pulses compared to
the ones from the fully active bulk volume. Dedicated studies of their shape
are therefore highly beneficial for the understanding and the rejection of
these unwanted events. This article reports about the development of a data-
driven pulse shape discrimination (PSD) method for the four 1 kg size PPC HPGe
detectors of the Conus experiment in the keV and sub-keV regime down to 210
eVee. The impact of the electronic noise at such low energies is carefully
examined. It is shown that for an acceptance of 90 % of the faster signal-like
pulses from the bulk volume, approx. 50 % of the surface events can be
rejected at the energy threshold and that their contribution is fully
suppressed above 800 eVee. Applied to the Conus background data, such a PSD
rejection cut allows to achieve an overall $(15-25)\,\%$ reduction of the
total background budget. The new method allows to improve the sensitivity of
future Conus analyses and to refine the corresponding background model in the
sub-keV energy region.

###### pacs:

29.40.WkSolid-state detectors and 27.50.+e 59$\leq$A$\leq$ 89 and
13.15.+gNeutrino interactions and 07.85.NcX-ray and $\gamma$-ray spectrometers

##  1 Introduction

The development of sub-keV energy threshold detectors opened the possibility
to detect tiny nuclear recoils from coherent elastic neutrino-nucleus
scattering (CE$\nu$NS) or to search for light dark matter interaction. High-
purity germanium (HPGe) detectors are specifically well suited for these
searches and widely used Soma2016 ; Aalseth2013 ; Collar2021b ; Belov2015 ;
Agnolet2017 ; Augier2021 ; Arnaud2018 ; Yang2018 . The Conus experiment aims
at detecting the CE$\nu$NS signal from reactor antineutrino CE$\nu$NS at the
Brokdorf nuclear power plant (Germany) by deploying four 1 kg size p-type
point contact (PPC) germanium detectors, specifically designed to have an
ultra-low noise and a very low intrinsic background contaminationBonet2021d .
Thanks to the very low background level achieved and an energy threshold of
210 eVee (electron equivalent energy), Conus has already provided stringent
constraints on the Standard Model CE$\nu$NS interactionBonet2021  and neutrino
physics beyond it Bonet2021b ; Bonet2022 . The challenge of the detection of
CE$\nu$NS at a reactor site was recently underlined by our precise measurement
of the quenching factor in germanium in the sub-keV energy region Bonhomme2022
and calls for further background reduction efforts and deeper understanding of
the signals at energies below 500 eVee, i.e. in the region of interest for the
CE$\nu$NS events search. For the Conus experiment, the background was
thoroughly studied and supported by detailed Monte-Carlo (MC) simulation
Bonet2021_bkg . It was shown that about half of the background consists of
internal contamination inside the end cap that produces radiation in the
vicinity of the diode susceptible to interact at the surface of the detectors.
In PPC HPGe detectors, these surface interactions lead to anomalous signals,
which can fall in the region of interest and complicate data interpretation. A
further step in the background rejection was therefore initiated by recording
the pulse waveforms of the signals in order to further discriminate surface
from bulk events.

This article reports the first results on pulse shape studies performed on
Conus data and it is organized as follows. The origin of the different types
of pulse shapes is exposed in Sec. 2. The methods developed to study and
categorize the different pulses are detailed in Sec. 3. Their validation and
the improvements expected for the Conus data analyses are discussed in Sec. 4
and is followed by the conclusion in Sec. 5.

##  2 Signal formation, readout and processing

###  2.1 Origin of pulses

(a)

(b)

(c)

Figure 1: Scheme of a diode used for the Conus detectors (a) with its
simulated electric potential (b) and weighting potential (c).

(a)

(b)

Figure 2: (a): Simulated charge signal $Q(t)$ induced at the p+ readout
electrode for single energy depositions at different interaction vertices
along the vertical axis of the cylindrical diodes with height 62 mm –
represented in Fig. 1(a). Residuals of the pulses w.r.t. a reference signal in
the center of the diode volume (z=30 mm) are shown in the bottom panel.
Deviations between pulses become negligible ($\lesssim$10%) for interactions
vertices farther away than $\sim$ 20 mm from the p+ contact. (b): Measured
pulses for the Conus-1 detector (each averaged over ten waveforms of 10 keVee
events) for a generic bulk event interaction (black), a bulk interaction in
the vicinity of the p+ contact (blue) and a surface event interaction (red).
The decay following the rising edge is due to the response of the electronics
(see Sec. 2.2).

The Conus experiment uses four 1 kg sized p-type point contact high-purity
germanium (PPC HPGe) detectors to pursue the goal of the detection of sub-
keVee signals from CE$\nu$NS. The design of the cylindric Ge diodes (diameter
and height of 62 mm) is shown in Fig. 1(a). A full description of the Conus
detectors can be found in Bonet2021d . The detectors are electrically depleted
by applying a positive reverse bias voltage of a few 1000 V on the surrounding
n+ contact, fabricated by lithium diffusion. The corresponding layer has a
typical thickness of 1 mm and can be decomposed into a dead layer (DL) and a
semi-active transition layer (TL) subject to incomplete charge collection. In
particular, the TL is evaluated to be $(0.16-0.24)\,\text{mm}$ thick from
calibration measurementsBonet2021d . The boron-implanted p+ contact has a
thickness of 100-200 nm and a small diameter of $(2.5-3)\,\text{mm}$; it
guarantees a detector capacitance $\lesssim$ 1 pF, which is crucial for low-
noise applications. The bottom layer separating p+ and n+ contacts is fully
passivated. Ionizing events in the electrically fully depleted bulk
volume111In the case of the Conus detectors, it corresponds to 91-95% of the
total crystal volume Bonet2021d . create electron-hole pairs proportional to
the deposited energy which drift in opposite directions according to the
electric field created by both the applied positive bias voltage and the space
charge of the depleted HPGe diode. The characteristic associated fast pulses
(FP) have a rise time of about 300 ns. However, an energy deposition within
the TL experiences a very weak electric field and diffuses until it either
recombines or reaches the bulk volume. This leads to partial charge collection
and a larger collection time due to the diffusion process (typically 1
$\mu$s). Such events will hence be referred thereafter as slow pulses (SP).
These events have an erroneous energy reconstruction which especially enhances
the event rate at energies below 10 keV. While particles of interest such
neutrinos or dark matter particles interact in the entire detector volume,
i.e. surface and bulk, and lead therefore mainly to FPs, specific classes of
background radiation are absorbed preferably in the outer detector layers
leading to the formation of SPs. This difference can be used, to discriminate
background from signal-like events.

To better understand the formation of pulses inside the bulk volume, the
electric potential inside the diode was calculated using fieldgen, a part of
the siggen software packageRadford2014 ; alvis2019multisite . It is shown for
illustration in Fig. 1(b). The p+ contact – collecting holes – is used for the
readout. The simulation framework siggen for signal pulse formation in
germanium detectors consists of two parts: the field generation code
mjd_fieldgen and the signal generation code mjd_siggen. Using the known
crystal impurity concentrations as well as their gradients, the dead layer and
passivation layer thickness and the size of the point contact as detector
specific input parameters, the electric field within the diode is determined.
If for a given radiation source the distribution and location of the energy
depositions within the diode are known e.g. from Geant4 MC simulations, the
expected read-out signal is determined by siggen according to the physics
described in the following. The time-dependent charge signal $Q(t)$ induced at
the readout electrode by the movement of the charge carriers can be calculated
using the Shockley-Ramo theorem Zhong2001 :

$Q(t)=-q_{0}\left[\,W(\mathbf{r}_{h}(t))-W(\mathbf{r}_{e}(t))\,\right]\,,$ (1)

where $q_{0}$ represents the total charge carried by the holes and electrons,
$\mathbf{r}_{h/e}(t)$ their trajectory paths, and $W$ the weighting potential,
i.e. the electric potential for the readout electrode and the others
electrodes at unit and zero potential, respectively. The corresponding
dimensionless $W$ for one Conus diode is represented in Fig. 1(c). The induced
charge $Q(t)$ increases over time until all charges are collected at the
electrodes. Due to the smallness of the p+ contact, the weighting potential is
homogeneous over a large part of the diode volume and becomes large only
within a few mm around the p+ contact. The induced charge mainly results from
the drift of the charges in these high $W$ regions and therefore, the position
dependency of the signal is rather weak for interactions within the full
volume. Induced signals shapes only start to differentiate for interactions
very close to the p+ contact. This is illustrated in Fig. 2(a), where the
induced charge signal $Q(t)$ at the readout electrode was calculated using
siggen for different interaction points inside the diode volume. From these
simulations, it appears that pulses start to have very fast rise times (very
fast pulses, hereafter VFP) for a distance from the point contact smaller than
20 mm, about ten times the size of the p+ contact. Outside of this volume,
bulk events have a very similar shape. A pulse corresponding to a vertex far
from the p+ contact will then simply be delayed compared to a signal from an
interaction closer to it, since the charges first have to drift before
entering the high field region where they induce a significant charge signal.
This characteristic feature provides a very good sensitivity above $\sim\,1$
keVee to distinguish single-site interactions from multi-site (MS)
interactions, resulting for example from multiple Compton scattering
processes, whose energy depositions are clearly distinguishable in their
signal shape. Such shape analyses were successfully implemented via e.g. the
A/E method with low capacitance large HPGe detectors Agostini2022 . Exemplary
recorded waveforms from different interaction vertices inside the Conus-1
diode are shown in Fig.2(b) for illustration.

###  2.2 Data acquisition and processing

Charge signals at the p+ contact are read out by a custom-built transistor
reset preamplifier (TRP), chosen to fulfill both ultra-low noise and low
background specifications Bonet2021d . In order to tag spurious events
following the reset of the baselines after saturation of the dynamic range,
inhibit TRP signals are generated. Artificially generated charge signals can
be injected at the test input of the preamplifier allowing to monitor the
response of the whole electronic chain.

A data acquisition system (DAQ) system based on commercially available CAEN
digitizers offering pulse processing implemented at the FPGA level was
installed for the data collection of Conus Run-5, lasting from May 21, 2021
until December 5, 2022 along with the pre-existing Lynx DAQ system. The CAEN
DAQ is operated via the CoMPASS software from CAEN and data are further
preprocessed and analysed offline. The HPGe signals from the four Conus
detectors were sampled at a rate of 100 MHz rate by a CAEN V1782 multi-channel
analyzer. Inputs with a 10 $\mu$s AC coupling were used in order to get rid of
the typical TRP step-like rising baseline level. A combination of a slow and a
fast triangular discriminator was used for the trigger. This gave better
results in terms of detection efficiency of small signals compared to the
standard RC-CR2 algorithm. The fast and slow triangular discriminators have
shaping times of 0.8 $\mu$s and 2.4 $\mu$s, respectively. The discriminator
thresholds were set as low as possible to optimize the detection efficiency of
low signals without exceeding 1 % of dead-time. This results in noise trigger
rates of about 500 Hz. The trigger efficiency was determined by injecting
artificial signals produced by a pulse generator with the same rise time as of
the physical signals. Fine-grained scans allowed to precisely measure the
detector response as a function of the energy. They are shown for the four
detectors in Fig. 3. The Conus-1 detector exhibits a better trigger efficiency
due to its slightly faster signal rise time combined with a lower noise level
compared to the other three detectors. A suitable description of the
experimental trigger efficiency curves was obtained using the fit function:

$\varepsilon_{trig}=0.5\cdot\left(1+\text{erf}\left(\frac{E_{ee}-t_{1}}{t_{2}}\right)\right)\,,$
(2)

with the error function erf and $t_{1}$ being in
$(225-295)\,\text{eV}_{\text{ee}}$ and $t_{2}$ in
$(105-138)\,\text{eV}_{\text{ee}}$, respectively, depending on the detector.
Note that the measured trigger efficiencies with this DAQ are slightly worse
compared to the ones obtained with the pre-existing Lynx DAQ system Bonet2021d
but this configuration was favored since it was found to offer a better
compromise between accepted noise events and signal efficiency, resulting in a
better signal-over-background ratio in the region of interest for the
CE$\nu$NS search.

Figure 3: Measured trigger efficiencies for the four Conus detectors (points)
with the CAEN DAQ during Run-5 and associated best-fit description (dashed
lines) using the analytical description of Eq. (2).

The energy of each event was reconstructed by a trapezoidal shaping filter
with a shaping time of 12 $\mu$s and a trapezoidal flat top of
$(3-6)\,\mu\text{s}$ depending on the detector, optimized in terms of energy
resolution. Signals separated by less than the trapezoidal shaping time are
identified as pile-up events and are rejected from the analysis chain. Because
of the very low count rates of the Conus detectors, the typical encountered
pile-up rates during data collection were very low – about 2 Hz. The logic
signals from the $\mu$-veto system as well as the TRP reset signals from the
Ge detectors were acquired with a CAEN V1725 module. The two modules shared a
common synchronization clock. All triggering events were saved and the time
cuts ($\mu$ and TRP) were applied offline in the analysis chain. In addition,
the raw waveforms in 20 $\mu$s windows are saved for all events triggering the
Ge detectors and are used for the offline pulse shape studies presented in
this article.

##  3 Methods

###  3.1 Rise time fit

The physical information allowing to discriminate bulk from surface events is
contained in the $\mathcal{O}$(100 ns) rise time of the signals and was
already exploited for PPC HPGe by CoGENT Aalseth2013  and TEXONO Yang2018b .
This rise time can be extracted on an event-by-event basis by fitting the
waveforms with the following analytical description:

$f(t)=A_{0}\left[\text{tanh}(\frac{t-t_{0}}{\tau})+1\right]\text{exp}(-\tau_{c}(t-t_{0}))+P_{0}\,,$
(3)

where $A_{0}$ is the amplitude of the signal, $t_{0}$ the timing offset,
$P_{0}$ the baseline noise level and $\tau$ the parameter of interest related
to the rise time of the event. The exponential decay accounts for the AC
coupling of the DAQ with a characteristic time $\tau_{C}\,\sim\,(7.5-8)\mu$s.
This parameter, being a constant of the system, was fixed in the fit in order
to help the fit converging at low energy. Exemplary fits for FPs induced by
bulk events and SPs induced by surface events are shown in Fig. 4. For all
four Conus detectors the rise time $\tau$ equals $(120-180)\,\text{ns}$ for
the FP (equivalent to t${}_{10-90\,\%}$ of $(275-325)\,\text{ns}$). This rise
time results both from the charge collection time in the diode and from the
time response of the preamplifier. The difference between the detectors is
explained by different settings of the preamplifiers. Depending on the nature
of the energy deposition, SPs are observed with rise times up to $\sim$ 1
$\mu$s, justifying for visibility reasons the logarithmic scale used for most
of the rise time distributions shown in this article.

Complementary methods (two sub-types of the A/E method Budjas2009  and an
integral-ratio method) were tested but were shown to be less efficient below a
few keVee compared to the $\tau$ fit method. Alternative parametrisations were
also tested for the pulse description and the use of the hyperbolic tangent in
Eq. (3) was found to be the best compromise between computation time and
discrimination power Henrichs2021 .

(a)

(b)

(c)

(d)

Figure 4: Fit examples of pulses at high (top) and low (bottom) energies for a
FP (left) and a SP (right) of the Conus-1 detector.

###  3.2 Physical event samples

For the energy calibration and the study of the shape of the pulses, advantage
is taken of the complementarity of various event samples.

First, calibration data from a 228Th source emitting $\gamma$-rays up to 2.6
MeV are available. These measurements were repeated every 2–4 weeks during the
whole data collection period allowing to monitor the stability of the rise
time over time. The $\sim$15 kBq source is inserted inside the Conus shield
Bonet2021  at a 20 cm distance from the diodes allowing for reasonable
statistics in a few hours and provides a good mixture of single-site events
(SPs, FPs and VFPs), MS and pile-up events.

Second, the Conus background data were considered. Before the application of
any veto, the experiment is dominated by muon-induced backgrounds. Hence, data
recorded in reactor-OFF before the application of the $\mu$-veto cut provide a
high-statistics sample of $\mu$-induced events which can be split in an
electromagnetic and a neutron-induced component Bonet2021_bkg . In the region
of interest below 10 keVee, the first component is sub-dominant and is
expected to contribute more to SP events, while the second dominant component
consists mainly of neutron-induced recoils inducing FP bulk events given the
large mean free-path of the neutrons. The application of the $\mu$-veto cut
suppresses $\sim$ 97 % of this background and the remaining $\mu$-induced
components are found to be of the same order of magnitude as for the detector-
dependent intrinsic contamination of 210Pb, present in the soldering in the
vicinity the diodes. Both type of backgrounds can result in SP and FP energy
depositions. However and according to the MC background model, a significant
amount of the 210Pb decays are expected to occur close to and even on the
surface of the diode close to the point contact, which is covered by the
passivation layer (cf. Fig. 1). Therefore, these decays will produce almost no
SPs: as there is no transition layer, they will induce energy deposition
within the bulk volume in the close vicinity of the p+ contact without
experiencing partial charge collection.

Finally, neutron activation in the germanium itself is responsible of the
appearance of characteristic lines in the background spectrum. The K-shell
line at 10.37 keV and L-shell line at 1.30 keV originating from atomic de-
excitation after electron capture in 68,71Ge are in particular well visible in
the Conus background, with a counting rate of about 16 counts/d/kg (2
counts/d/kg) for the K-shell (L-shell) during Run-5. In situ activation of
71Ge is the main contribution due to the long storage of the detectors inside
the shieldBonet2022 . The corresponding events are particularly interesting
since they are expected to be homogeneously distributed within the full Ge
volume and to contribute mainly to FP events in contrast to the complex
vertices distribution of the background. A clean sample of these events can be
used to reproduce artificially FP events, see Sec. 3.3.

The rise time of background events recorded during reactor-OFF after
application of the $\mu$-veto is shown in Fig. 5 as a function of energy. The
fast and slow populations (around log${}_{10}\tau\sim 2.1$ and
log${}_{10}\tau\sim 2.9$, respectively) are clearly separated for the highest
energies and start to mix in the sub- keVee regime due to the growing
influence of electronic noise. The gap between the two main bands is populated
by MS or pile-up events which are too close to be discarded by the pile-up
rejection of the DAQ and whose pulse fit leads to artificially slow rise
times. The lower band gathers both the FPs (log${}_{10}\tau\gtrsim 2.1$)
resulting from interactions inside the bulk volume and VFPs
(log${}_{10}\tau\lesssim 2.1$) for interactions close to the p+ contact as
illustrated in Fig. 2(b). The 10.37 keVee K-shell atomic de-excitation line
following Ge EC decays is particularly visible within the FP part of the band
since for these homogeneously distributed events, only the small fraction of
the decays close to the p+ contact are expected to lead to VFPs.

To illustrate the complementary in terms of pulse shape between the
aforementioned physical samples, Fig. 6 presents their rise time distributions
at mid (left pannel) and low energy (right panel) for two exemplary Conus
detectors (similar distributions as for Conus-4 are obtained for the two
others). From this it can be seen that the discrimination power between the
two main populations (SPs and FPs) is detector dependent: Conus-1 exhibits a
better separation than the other Conus detectors thanks to its slightly faster
response. Information on the underlying background components can be gained by
comparing the FP distributions at high energy before and after the application
of the $\mu$-veto, in red and blue, respectively. It appears that the
remaining background with applied $\mu$-veto contains a VFP component –
negligible without $\mu$-veto – dominating the distribution for Conus-1 and of
equal weight as the FPs for Conus-4. This is consistent with the findings
about the background decomposition detailed in Bonet2021_bkg , showing that
the background of Conus-1 is dominated by the intrinsic 210Pb contamination
(responsible for the VFP component since the decays occur mainly in the p+
contact region according to the MC model) whereas the tendency is inverted for
the other detectors.

At low energy however, the two populations mix due to the growing influence of
electronic noise. The large difference in terms of fraction of SPs between the
228Th calibration data and the background data confirms the high-dependency of
vertex distributions.

Figure 5: Rise time of background events in dependence of energy after
application of the $\mu$-veto for the Conus-1 detector. Bulk events are
characterised by a fast rise time and populate the lower band. The 10.37 keV
Ge activation line is clearly visible within the bulk event population.
Surface events produce slow signals which are found in the upper band. The red
line shows exemplary a cut with a 97 % FP signal acceptance (SA).

(a)

(b)

(c)

(d)

Figure 6: Rise time distributions for different event populations and for two
detectors in two different energy ranges: (a): Conus-1 for energies in (2, 8)
keVee. (b): Conus-1 for energies in (0.5, 1) keVee. (c): Conus-4 for energies
in (2, 8) keVee. (d): Conus-4 for energies in (0.5, 1) keVee. Background data
collected in reactor-OFF before (after) application of the $\mu$-veto cut are
shown in red (blue), and calibration data from a 228Th radioactive source
irradiation are given in orange. The upper range used for normalisation is
indicated by the gray arrow.

###  3.3 Artificial generation of bulk signals

In order to determine the cut acceptance for neutrino events, the rise time
distribution of homogeneously distributed bulk events has to be known. As
mentioned above, no calibration source allows to have a pure homogeneous bulk
event sample over the whole energy range. To overcome this obstacle,
artificial signals from an electronic pulser mimicking bulk events were
generated and injected through the whole DAQ chain and treated in the same way
as measured physical signals.

The correct shape of these signals can be determined either via a validated
electric field simulation of the diode or via a data-driven method consisting
in deconvolving the measured bulk signals at higher energy from the DAQ chain
response, the latter was chosen in this work. The deconvolution was performed
in a semi-empirical way as follows: the impulse response function of the DAQ
chain was first measured independently for all detectors by injecting a fast
rectangular signal in the test input of the preamplifier. The input signal
mimicking a bulk signal was then determined iteratively as the signal whose
convolution with the impulse response function matches a real measured
physical output pulse. In practise, a measured signal was used as initial
pulse and was distorted in an iterative way to achieve a match to the output
signal. About ten iterations are sufficient to obtain a very good agreement
($\lesssim$ 1 ns difference in terms of rise time) between the convolution of
the calculated input pulse and the real pulses Stauber2022 . This procedure is
illustrated in Fig. 7. The obtained red signal is then used as description of
a pure bulk signal and injected in the DAQ chain of the detectors for the
determination of the rise time distribution. The energy dependence of the
distribution was studied by a simple scaling of this input signal. The
hypothesis of a constant shape independent of the energy was verified using
calibration data between 1 and 10 keVee.

Fig. 8 illustrates the evolution of the measured rise time distributions as a
function of the energy obtained from a scan over the amplitude of the
theoretical input signal mimicking bulk events for Conus-1. As the energy
decreases, the relative importance of noise increases and explains the larger
smearing of the distributions. Note that for energies $\lesssim$ 1 keVee, the
distributions also deviate from a Gaussian shape and exhibit an asymmetric
behavior. From these distributions, the signal acceptance (SA) of FPs can be
estimated for any rise time cut value. At very low energy ($\lesssim$ 0.4
keVee), a small fraction of the rise time fits do not converge and the
corresponding events are discarded from the analysis chain, leading to an
additional reduction in the SA. It amounts to about 1 % at 300
eV${}_{\text{ee}}$ and 3 % at 250 eV${}_{\text{ee}}$.

Figure 7: Illustration of the iterative method: starting from a high-energy
measured real pulse ($o_{meas}(t)$, gray dotted line), the unknown input pulse
corresponding to bulk events ($i_{th}(t)$, red line) is determined iteratively
such that its convolution with the detector response ($o_{meas}(t)$, blue
line) matches the real pulse. Figure 8: Rise time distributions from
artificially generated bulk events for different energies, from 300 eVee to
10.35 keVee. The electronic noise is responsible for the broadening of the
distributions and for the deviation from a Gaussian below 1 keVee.

##  4 Results and discussion

###  4.1 Validation and systematic uncertainties

The four physics samples described in Sec. 3.2 can be used to validate
qualitatively and quantitatively the correct classification of the pulses via
their rise time.

First, the rise time distribution of the 228Th calibration data were compared
to a MC simulation (see Bonet2022  for a more detailed description of the
228Th MC simulations) in order to validate the ratio of bulk events over
surface events. This ratio depends on the energy of the $\gamma$-rays, the
width of the detector-dependent transition layers and the geometry of the full
experimental setup. From a Geant4 MC simulation of such a calibration
measurement, the fraction of interactions happening in the transition layers
can be determined and compared to the amount of pulses identified as slow in
the data. The fraction of SPs is about 20 % above 15 keV and increases to 30 %
(45 %) towards low energy, in the $(6-12)\,\text{keV}$ and $(2-6)\,\text{keV}$
ranges, respectively, due to the larger attenuation. An agreement at the 5 %
level is achieved Henrichs2021 . This demonstrates on one hand the accurate
description of the detector transition layer and on the other hand the
consistency of the surface event identification over the energy, showing no
additional unexpected component within this population. A percent-level
comparison was achieved by coupling the vertex distributions from the MC with
a full pulse shape simulation, allowing to take also into account MS events
corresponding to multiple energy depositions within one event, that lead to an
artificially larger fitted rise time if they occur within a too small time
interval.

Second, as a validation of the procedure of the artificial generation of bulk
events described in Sec. 3.3, the measured rise time distributions resulting
from the injected theoretical bulk signal were compared with the sample of
real physical events from the 10.37 keVee line from Ge activation. As
mentioned above, these events are expected to be homogeneously distributed in
the diode volume and should in consequence be predominantly made of bulk
events. Their distribution is obtained statistically after subtraction of the
background under the line and are shown exemplary for Conus-1 and Conus-4 in
Fig. 9. The remaining sub-component with faster rise time values corresponds
to VFPs originating from bulk events close to the p+ contact. The fast
component of the distributions is compared to the rise time distribution of
the artificial bulk signals, also shown in gray. By construction, an agreement
in terms of mean rise time is achieved at the ns level. Most importantly, the
width of the distributions is found to be in agreement at the $(5-10)\,\%$
level, indicating that the influence of noise is mostly responsible for the
smearing of the distribution and that no additional process contributes
significantly to the variability in the shape of the physical bulk pulses. The
small remaining differences – both in terms of mean value and width of the
distributions as reported in Tab. 1 – are taken as systematic uncertainties,
their impact is discussed below.

Figure 9: Rise time distribution of events in the 10.37 keV line after background subtraction for Conus-1 (black) and Conus-4 (blue). The sub-component with faster rise time corresponds to the VFPs, which correspond to energy depositions in the bulk region close to the p+ contact. The rise time distribution of the artificially generated fast bulk signals is shown in gray and is in agreement with the measured distributions. The fraction of VFPs is extracted via a fit of the data with two Gaussian distributions modeling the two components. Detector |  $\Delta\mu$ (log10) |  $\Delta\mu$ (ns) |  $\Delta\sigma$ (%) | VFP fraction (%) | rp+ (mm)
---|---|---|---|---|---
Conus-1 | 0.001 | 0.5 | 5 | 8 $\pm$ 1 | 19.5 $\pm$ 0.8
Conus-2 | 0.002 | 0.8 | 4 | 13 $\pm$ 1 | 22.7 $\pm$ 0.6
Conus-3 | 0.004 | 1.7 | 10 | 12 $\pm$ 2 | 22.1 $\pm$ 1.2
Conus-4 | -0.001 | -0.3 | 5 | 11 $\pm$ 1 | 21.7 $\pm$ 0.2
Table 1: Observed deviation of mean ($\Delta\mu$) and width ($\Delta\sigma$)
of the rise time distributions of FPs events originating from the artificially
pulser-generated bulk signals and the homogeneous 10.37 keV data event sample.
The fraction of VFPs can also be extracted from these two data distributions.
They can be directly translated into an effective half-sphere volume around
the p+ contact of radius rp+ in which the charge collection happens very fast.

Third, the stability of the rise time of the signals was monitored over time
using the bulk event distribution of 228Th calibration data. Fig. 10 shows the
mean position of the rise time distributions over a large fraction of the
Run-5 data collection period for the four Conus detectors. A maximal deviation
of $\sim$ 3 ns is observed for the Conus-2 detector and is related to small
instabilities of the pulse shape observed for this particular detector from
228Th measurements. This value is taken as additional source of systematic
uncertainty.

Figure 10: Deviation of the mean rise time of the bulk population over time
for the four Conus detectors, extracted from bulk events of the 228Th
calibration data with energies within (5, 15) keVee. The shown data cover the
period from May 2021 until June 2022, corresponding to a large fraction of the
Run-5 data taking period.

Fourth, an elementary simulation of event samples was developed in order to
get a qualitative understanding of the broadening of the rise time
distribution due to electronic noise. It also allows to investigate the impact
of the small systematic uncertainties on the mean position and the width of
the rise time distribution of bulk events mentioned above. For this purpose,
the theoretical input pulses mimicking bulk events are convoluted with the
response function of the detector including an accurate model of the
electronic noise. The noise samples are generated for each detector from their
noise frequency spectrum, measured directly via dedicated baseline samples
with 1 ms acquisition windows acquired with an external arbitrary trigger. The
resulting simulated rise time distributions match very well the measured
distributions even at the lowest energies; an example is presented in Fig. 11.

The impact of the systematic uncertainties on the distributions was then
investigated by simulating input pulses with slightly varying rise times and
by adding an artificial variability in the signals covering the small
mismatches at high energy observed in the comparison with the 10.37 keV line
as well as the time variations of the rise time. As shown in Fig. 11, such
small changes do not affect the distributions which are dominated by the
impact of the noise. The systematic uncertainties discussed above can
therefore safely be neglected in this energy range.

Figure 11: Measured (black) and simulated (red) rise time distributions for
fast bulk events of 380 eVee. The impact of a 3 ns shift in the mean rise time
of the input pulse (orange) and an additional variability in the input pulses
(blue) do not affect the distribution, mostly determined by noise.

###  4.2 Slow pulse rejection power

A rise time cut can be set from the measured rise time distributions of
artificial bulk events, providing a good reference for homogeneously
distributed vertices in the diode.

In order to evaluate the SP rejection capability of the PSD cut, 228Th
calibration data were used as they provide a high statistics sample of SPs.
The rise time distribution of the SP population only is obtained by
subtraction of the FP distribution – determined from pulser scans as described
in the previous section – from the 228Th rise time distribution. In this way,
the SP rejection power can be quantified for each detector as a function of
the energy for a given FP SA cut. This is illustrated in Fig. 12(a) in the
region of interest for a constant FP SA of 90 %.

The uncertainties reported contain both statistic and systematic
contributions. Systematic uncertainties were estimated in a conservative way
by varying the integration window for the normalisation of the FP in the
subtraction procedure. Note that the SP rejection numbers reported here are
given for information only, helping to quantify the discrimination power of
the method and allowing a direct comparison with other experiments. They are
not used directly in any of the CONUS analyses since their values only affect
the amount of background entering in the region of interest, which is measured
directly during the reactor-off periods. However, any drifts in time causing a
time dependent SP contamination could affect the level of background entering
the ROI. This is in particular crucial for the lowest energy points having the
highest SP contamination. In the same way as for the FP population as shown in
Fig. 11, it was therefore cross-checked with the event generator that small
rise-time drifts have a negligible impact on the SP acceptance.

Close to the $\sim$ 210 eVee energy threshold, about half of the SPs can be
rejected and a perfect separation between SPs and FPs is reached between 0.7
and 0.8 keVee depending on the selected detector. Variability between
detectors is clearly visible, with the best rejection power obtained for the
Conus-1 detector. With the adoption of a more conservative cut – optimized
with SA of FPs of 97 % for the CE$\nu$NS analysis Stauber2022  – a rejection
of $\sim\,$15 % of the slow pulse component close to the energy threshold is
achieved as shown in Fig. 12(b).

Compared to previous similar studies with PPC HPGe detectors, this represents
a notable improvement, in particular for the lowest energies: the SA reported
both by TEXONO Li2014  and CDEX-1 Zhao2016  start to drop for energies
$\lesssim\,$0.8 keVee whereas for COGENT Aalseth2013 , a constant 90 % SA is
conserved at the cost of a less strict SP rejection (less than 50 % rejection
below 0.8 keVee). This indicates that the separation between the SP and FP
populations achieved in this work is better down to lower energies. It mainly
results from the very low noise level achieved for the CONUS detectors, a
factor two lower than the experiments quoted above with typical pulser FWHM in
the $(140-160)\,\text{eV}_{\text{ee}}$ range. As mentioned in 3.1, the use of
the rise-time fit was also shown to help gaining efficiency toward low
energies compared to alternative simple methods, less cost-intensive but whose
discriminating observables are too affected by noise for sub-keV energies
Henrichs2021 .

(a)

(b)

Figure 12: Rejection efficiency of SPs as a function of the energy. The rise
time distributions of the SP populations are obtained from the subtraction of
the fast component from the 228Th rise time distribution.
(a): Rejection efficiency for the four Conus detectors and a constant 90 % SA
of FPs.
(b): Average trends over the four Conus detectors for a SA of 97 % (red) and
90 % (blue) of SPs for comparison. The y errors bars correspond to the root-
mean-square value over the detectors.

###  4.3 PSD application to CONUS background data

The pulse shape discrimination power provided by the rise time observable can
be directly applied to reduce the background level of the Conus experiment.
Indeed, due to their very large mean free path, neutrinos are expected to
interact homogeneously in the entire detector and therefore mainly in the
active bulk volume, whereas certain types of background radiation are expected
to interact closer to the semi-active detector surface. Thanks to the very
careful mitigation measures in the design of the detectors, the background
reduction achieved for Conus is moderate since a large part of the background
after application of the $\mu$-veto consists of FPs (see e.g. blue
distributions in Fig. 6). The impact of a 90 % FP bulk-events SA cut is shown
exemplary for Conus-1 in Fig. 13. An overall $(15-25)\,\%$ reduction of the
background level is achieved, depending on the detector and on the considered
energy range. In the CE$\nu$NS region of interest below 1 keVee, a background
reduction of about $(15-20)\,\%$ is achieved. A more conservative selection of
FP-like signals at 97 % still leads to a non-negligible $(5-10)\,\%$
reduction. These numbers are consistent with the MC background model
expectation, predicting a contribution of SPs of about 40 % in this region.
This number has to be coupled with the SP rejections shown in Fig. 12(b) to
retrieve the typical background reductions found in the data. Searches beyond
the standard model carried out at higher energy ($\sim$1–10 keV) with the
Conus data will also benefit from this background reduction. For this region
where a complete SP rejection is achieved, an overall background reduction of
about 25 % is obtained. The gain in terms of shape of the background will also
be of particular relevance: the surface event contamination is indeed expected
to raise at low energy and can degenerate with most of the searched signals
(CE$\nu$NS or beyond Standard Model physics), also having a typical rise
towards low energies. Any background reduction flattening the background
energy spectrum is therefore of high importance to gain sensitivity in these
searches.

Figure 13: Background energy spectra with $\mu$-veto for Conus-1 for Run-5,
before (black) and after (red) the application of a PSD cut with 90 % SA of
FPs. The rejected event spectrum is shown is blue.

More generally, pulse shape studies allow to improve knowledge on the detector
and can help refining and constraining the understanding of the background. An
illustration of this is the study of the VFPs corresponding to interactions in
the vicinity of the p+ contact, also discussed in Li2022 . At high-energy,
they are distinguishable by their prompt rise time from the homogeneous bulk
events as discussed in Sec. 2.1 and are responsible for the sub-component
observed e.g. in Fig. 9. Using the rise time distributions of the
homogeneously distributed events from the 10.37 keV Ge activation line, it is
possible to get back to information on the weighting potential close to the p+
contact: their fraction should correspond to the volume fraction of this high
$W$ region over the whole diode volume. These VFPs amount to $(8-13)\,\%$ of
the bulk events, depending on the Conus detectors, as reported in Tab. 1.
These numbers are relevant for the validation of the full pulse shape
simulation, currently under development. This translates into a half-sphere of
$(19.5-22.7)\,\text{mm}$ radius around the p+ contact Henrichs2021 . These
relatively small fractions for Conus moderate the suggestion expressed in
Li2022  regarding the potential of additional background reduction by
discriminating this type of events. Moreover, at low energy, the two
populations are not distinguishable anymore due to the overwhelming noise.
Efforts in reducing the electronic noise would allow a better discrimination
of these population. Moreover, a reduced p+ contact size leading to lower
noise would also modify the electric potential close to the contact, which
would further reduce the fraction of these VFPs.

Using calibration sources of known energy, the rise time information is of
high-relevance for the validation of the charge collection inefficiency in the
transition layers as it can provide an additional and independent cross-check
to the method described in Bonet2021d  using 241Am sources.

Finally, coupled to the full MC simulation of the Conus background, the rise
time distributions are expected to provide further insights on the composition
and vertices distributions of the different background components. This will
allow to independently validate and complement the results already reported in
Bonet2021_bkg .

##  5 Conclusion

This article reports about the first pulse shape data analysis in the sub-keV
and keV energy regime conducted in the framework of the Conus experiment. It
exploits the capability of the used ultra-low threshold point-contact p-type
high-purity germanium detectors to distinguish and discriminate surface events
from bulk events by means of different signal rise time distributions. A
fitting procedure of the waveform has been selected and implemented that
allows to distinguish rise times down to the detector energy thresholds of
$\sim$ 210 eVee. A detailed and comprehensive understanding of the different
pulse shape populations was achieved via the study of physics data (Conus
background data and radioactive source calibration data), artificially
injected data (via a pulser generator) and pulse shape simulation data. It was
verified that the pulse shapes of bulk events are homogeneous and independent
of the interaction vertices except for about 10 % of the bulk events, which
interact very close to the p+ contact and exhibit a very fast rise time. A
data-driven method was developed to generate artificially signals mimicking
bulk events that can then be reliably used to characterize the bulk rise time
distributions at any energy. The sub-keVee region of interest was studied in
detail and the growing relative influence of the electronic noise towards the
noise threshold on the rise time distributions was modeled and understood.
About half of the slow pulses can be rejected at energies as low as $\sim$ 210
eVee while keeping a bulk event acceptance greater than 90 %. A complete slow
pulses suppression is achieved above $\sim$ 800 eVee. Thanks to these studies,
a total background reduction of about $(15-25)\,\%$ is achieved for the Conus
experiment with negligible systematic uncertainties, which agrees with our
Monte Carlo background model expectation predicting $\sim\,$40% of slow pulses
in the total background of the Conus detectors. This additional background
reduction and flattening will improve the sensitivity of the ongoing Conus
analyses. Beyond that, Conus plans further pulse shape studies to validate
pulse shape simulation codes and to refine background models by providing
independent information on different background components and thus reducing
associated systematic uncertainties.

Acknowledgements We thank all the technical and administrative staff who
helped building the experiment, in particular the MPIK workshops and Mirion
Technologies (Canberra) in Lingolsheim. We express our gratitude to the
Preussen Elektra GmbH for great support and for hosting the Conus experiment.
The Conus experiment is supported financially by the Max Planck Society (MPG).

## References

  * (1) A. K. Soma et al.  _Characterization and Performance of Germanium Detectors with sub-keV Sensitivities for Neutrino and Dark Matter Experiments_.  Nucl. Instrum. Meth. A, 836:67–82, 2016.  doi:10.1016/j.nima.2016.08.044.  arXiv:1411.4802.
  * (2) C. E. Aalseth et al.  _CoGeNT: A search for low-mass dark matter using $p$-type point contact germanium detectors_.  Phys. Rev. D, 88:012002, 2013.  doi:10.1103/PhysRevD.88.012002.
  * (3) J. Colaresi et al.  _First results from a search for coherent elastic neutrino-nucleus scattering at a reactor site_.  Phys. Rev. D, 104:072003, 2021.  doi:10.1103/PhysRevD.104.072003.
  * (4) V. Belov et al.  _The $\upnu$GeN experiment at the Kalinin Nuclear Power Plant_.  Journal of Instrumentation, 10(12):P12011–P12011, 2015.  doi:10.1088/1748-0221/10/12/p12011.
  * (5) G. Agnolet et al.  _Background studies for the MINER Coherent Neutrino Scattering reactor experiment_.  Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment, 853:53–60, 2017.  ISSN 0168-9002.  doi:https://doi.org/10.1016/j.nima.2017.02.024.
  * (6) C. Augier et al.  _Ricochet Progress and Status_.  In 19th International Workshop on Low Temperature Detectors. 2021.
  * (7) Q. Arnaud et al.  _Optimizing EDELWEISS detectors for low-mass WIMP searches_.  Phys. Rev. D, 97:022003, 2018.  doi:10.1103/PhysRevD.97.022003.
  * (8) L. T. Yang et al.  _Limits on light WIMPs with a 1 kg-scale germanium detector at 160 eVee physics threshold at the China Jinping Underground Laboratory_.  Chin. Phys. C, 42(2):023002, 2018.  doi:10.1088/1674-1137/42/2/023002.  arXiv:1710.06650.
  * (9) H. Bonet et al.  _Large-size sub-keV sensitive germanium detectors for the CONUS experiment_.  Eur. Phys. J. C, 81(3):267, 2021.  doi:10.1140/epjc/s10052-021-09038-3.  arXiv:2010.11241.
  * (10) H. Bonet et al.  _Constraints on Elastic Neutrino Nucleus Scattering in the Fully Coherent Regime from the CONUS Experiment_.  Phys. Rev. Lett., 126(4):041804, 2021.  doi:10.1103/PhysRevLett.126.041804.  arXiv:2011.00210.
  * (11) H. Bonet et al.  _Novel constraints on neutrino physics beyond the standard model from the CONUS experiment_.  JHEP, 05:085, 2022.  doi:10.1007/JHEP05(2022)085.  arXiv:2110.02174.
  * (12) H. Bonet et al.  _First upper limits on neutrino electromagnetic properties from the CONUS experiment_.  Eur. Phys. J. C, 82(9):813, 2022.  doi:10.1140/epjc/s10052-022-10722-1.  arXiv:2201.12257.
  * (13) A. Bonhomme et al.  _Direct measurement of the ionization quenching factor of nuclear recoils in germanium in the keV energy range_.  Eur. Phys. J. C, 82(9):815, 2022.  doi:10.1140/epjc/s10052-022-10768-1.  arXiv:2202.03754.
  * (14) H. Bonet et al.  _Full background decomposition of the CONUS experiment_.  Eur. Phys. J. C, 83(3):195, 2023.  doi:10.1140/epjc/s10052-023-11240-4.  arXiv:2112.09585.
  * (15) D. Radford.  _Siggen Simulation Software Package_.  https://github.com/radforddc/icpc_siggen, 2014.
  * (16) S. Alvis et al.  _Multisite event discrimination for the MAJORANA DEMONSTRATOR_.  Physical Review C, 99(6):065501, 2019.
  * (17) Z. He.  _Review of the Shockley–Ramo theorem and its application in semiconductor gamma-ray detectors_.  Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment, 463(1):250–267, 2001.  ISSN 0168-9002.  doi:https://doi.org/10.1016/S0168-9002(01)00223-6.
  * (18) M. Agostini et al.  _Pulse shape analysis in Gerda Phase II_.  The European Physical Journal C, 82(4), 2022.  doi:10.1140/epjc/s10052-022-10163-w.
  * (19) L. Yang et al.  _Bulk and surface event identification in p-type germanium detectors_.  Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment, 886:13–23, 2018.  ISSN 0168-9002.  doi:https://doi.org/10.1016/j.nima.2017.12.078.
  * (20) D. Budjáš, M. B. Heider, O. Chkvorets, N. Khanbekov and S. Schönert.  _Pulse shape discrimination studies with a Broad-Energy Germanium detector for signal identification and background suppression in the GERDA double beta decay experiment_.  Journal of Instrumentation, 4(10):P10007–P10007, 2009.  doi:10.1088/1748-0221/4/10/p10007.
  * (21) J. Henrichs.  Development of a Pulse Shape Analysis for the CONUS Experiment.  Master thesis, Ruprecht-Karls-Universität, Heidelberg, 2021.  http://hdl.handle.net/21.11116/0000-0009-C7E1-8.
  * (22) J. Stauber.  Pulse Shape Discrimination for the CONUS Experiment in the sub-keV Regime.  Master thesis, Ruprecht-Karls-Universität, Heidelberg, 2022.  http://hdl.handle.net/21.11116/0000-000A-EBAC-C.
  * (23) H. B. Li et al.  _Differentiation of Bulk and Surface Events in p-type Point-Contact Germanium Detectors for Light WIMP Searches_.  Astropart. Phys., 56:1–8, 2014.  doi:10.1016/j.astropartphys.2014.02.005.  arXiv:1311.5957.
  * (24) W. Zhao et al.  _Search of low-mass WIMPs with a $p$-type point contact germanium detector in the CDEX-1 experiment_.  Phys. Rev. D, 93:092003, 2016.  doi:10.1103/PhysRevD.93.092003.
  * (25) R.-M.-J. Li et al.  _Identification of anomalous fast bulk events in a $p$-type point-contact germanium detector_.  Nucl. Sci. Tech., 33(5):57, 2022.  doi:10.1007/s41365-022-01041-x.


# Optimal terminal sliding mode control for second-order motion systems

Michael Ruderman This work has received funding from the EUs H2020-MSCA-RISE
research and innovation programme under grant agreement No 734832.M. Ruderman
is with Faculty of Engineering and Science, University of Agder, Norway.
Email<EMAIL_ADDRESS>
###### Abstract

Terminal sliding mode (TSM) control algorithm and its non-singular refinement
have been elaborated for two decades and belong, since then, to a broader
class of the finite-time controllers, which are known to be robust against the
matched perturbations. While TSM manifold allows for different forms of the
sliding variable, which are satisfying the $q/p$ power ratio of the measurable
output state, we demonstrate that $q/p=0.5$ is the optimal one for the second-
order Newton’s motion dynamics with a bounded control action. The paper
analyzes the time-optimal sliding surface and, based thereupon, claims the
optimal TSM control for the second-order motion systems. It is stressed that
the optimal TSM control is fully inline with the Fuller’s problem of optimal
switching which minimizes the settling time, i.e. with time-optimal control of
an unperturbed double-integrator. Is is also shown that for the given plant
characteristics, i.e. the overall inertia and control bound, there is no need
for additional control parameters. The single surface design parameter might
(but not necessarily need to) be used for driving system on the boundary layer
of the twisting mode, or for forcing it to the robust terminal sliding mode.
Additional insight is given into the finite-time convergence of TSM and
robustness against the bounded perturbations. Numerical examples with
different upper-bounded perturbations are demonstrated.

##  I INTRODUCTION

Terminal sliding mode control belongs to the class of finite time controllers
[1], see e.g. [2] for survey, for which a finite-time convergence can be
guaranteed. The finite-time convergence is here valid for both, the _reaching
phase_ during which the state trajectories reach the predefined sliding
surface, and the _sliding phase_ of a convergence to the state equilibrium. In
other words, when the control task is formulated as a zero-reference and non-
zero initial value problem, the controlled system attains the origin in the
finite-time. Yet, for the reaching phase one refers to a finite-time
convergence of the sliding variable to zero, while for the second (sliding)
phase a finite-time convergence of the output tracking error to zero is meant.
For more details on the reaching and sliding phases and convergence of the
sliding mode control systems we refer to the seminal literature e.g. [3, 4,
5].

Based on [6, 7] the so-called terminal sliding mode control has been proposed
in [8] and, later, brought into a non-singular form in [9]. The difference
between both will be summarized in the preliminaries provided in Section II.
Remarkable is a fact that for the second-order systems, which are in focus of
this work, the terminal sliding mode surface, correspondingly control,
coincides with the Fuller’s problem of optimal relay switching [10]. The
latter has also been brought into context and analyzed in terms of designing
high-order sliding mode control algorithms in [11]. It should be recalled that
for the second-order Newtonian dynamics of a relative motion, it is always
possible to bring the system to an equilibrium in finite time, by using a
bounded control, i.e. through inherently bounded actuation forces. The well-
known approach to this problem leads to the co-called bang-bang control [12],
which is the time-optimal solution of an unperturbed double-integrator control
problem. We will also briefly sketch it in Section III for convenience of the
reader.

In the rest of the paper, we will show that an optimal terminal sliding mode
control can be designed for motion systems with the bounded perturbations and
control actions, based on the time-optimal switching of double-integrator. The
main results with analysis are given in Section IV. Section V demonstrates
three numerical examples of the control performance for different-type matched
perturbations. The main conclusions are briefly drawn in Section VI.

##  II PRELIMINARIES

The so-called terminal sliding mode control suggests the first-order terminal
sliding variable, cf. [9],

$s=x_{2}+\beta x_{1}^{q/p},$ (1)

where $\beta>0$ is a design parameter and $p>q$ are the positive odd constants
[9]. Recall that the terminal sliding mode controllers belong to the class of
second-order sliding modes with finite-time convergence, see e.g. [13]. For a
terminal sliding mode, the control of a dynamic system
$\ddot{x}_{1}=\dot{x}_{2}=u$, in the most simple case, is given by

$u=-\alpha\,\mathrm{sign}(s),$ (2)

where $\alpha>0$ is the control parameter. The requirement on $p,q$ to be odd
[9] can be relaxed, and a frequently encountered terminal sliding surface is
written as, cf. [13],

$s=x_{2}+\beta\sqrt{|x_{1}|}\mathrm{sign}(x_{1}).$ (3)

By taking the time derivative of the sliding surface one can show that, cf.
[14],

$\dot{s}=-\alpha\,\mathrm{sign}(s)+\beta
x_{2}\bigl{(}2\sqrt{|x_{1}|}\bigr{)}^{-1}.$ (4)

One can recognize that a singularity occurs if $x_{2}\neq 0$ and $x_{1}=0$.
This situation does not occur in an ideal sliding mode, i.e. $s=0$. However,
the terminal sliding mode control with surface (1) cannot guarantee an always
bounded control signal before the system is on the sliding manifold $s=0$,
i.e. during the reaching phase which depends on the initial conditions. Once
on the surface, the sliding variable dynamics (4) reduces to

$2\dot{s}=-2\alpha\,\mathrm{sign}(s)-\beta^{2}\,\mathrm{sign}(x_{1}).$ (5)

Thereupon, one can show that the existence condition $s\dot{s}<0$ of an ideal
sliding mode is fulfilled for $\beta^{2}<2\alpha$. For that case, the
trajectories of the system reach the surface (3) and remain there for all
times afterwards. This control system behavior is often denoted as a _terminal
mode_ , cf. [15, 14]. When the design parameters are assigned as
$\beta^{2}>2\alpha$, the trajectories of the system do not remain on the
sliding surface, while the control still switches at $s=0$, according to (2),
(3). That case the control system trajectories proceed in the so-called
_twisting mode_ , cf. [15, 14].

For overcoming the singularity of the sliding surface (1), the non-singular
terminal sliding mode control with

$s=x_{1}+\beta^{-1}x_{2}^{p/q},$ (6)

has been proposed in [9], while $\beta$, $p$, and $q$ parameters are kept the
same as before. This control works around the surface singularity problem
while, at the same time, (1) and (6) describe one and the same switching
surface, seen from the geometrical interpretation, in the $(x_{1},x_{2})$
coordinates.

##  III TIME OPTIMAL CONTROL OF SECOND-ORDER MOTION SYSTEMS

We consider a generalized111We mean here the generalized coordinates $x$ and
forces $u$ for not explicitly distinguishing between the translational and
rotational degrees of freedom (DOFs). As consequence, $m$ is used for denoting
equally both quantities, inertial mass and moment of inertia respectively.
second-order motion system

$\ddot{x}(t)=m^{-1}u(t)$ (7)

of an inertial body $m$, which is first not affected by any counteracting,
correspondingly disturbing, forces. For the sake of simplicity and without
loss of generality we focus on the 1DOF motion, so that the scalar state
variables of motion are $x_{1}=x$ and $x_{2}=\dot{x}$. The control input
force, here and for rest of the paper, is constrained as $u\in[-U,+U]$. Here
it is worth emphasizing that the bounded (often denoted as _saturated_)
control signal, correspondingly actuation force, is in focus of our optimal
terminal sliding mode control. Consequently, $U$ will appear as a system
parameter in the following analysis and control synthesis.

Assuming an initial state $[x_{1}(0),x_{2}(0)]^{T}$ and a final state
$[x_{1}(t_{f}),x_{2}(t_{f})]^{T}$ the time-optimal control [10, 12, 16] will
minimize the cost criterion

$J(u)=t_{f}=\int\limits_{0}^{t_{f}}dt,$ (8)

thus minimizing the overall convergence (or settling [10]) time $t_{f}$.
Following the Pontryagin’s minimum principle and solving the constrained (by
control limit $U$ and initial and final conditions) optimization problem (8)
leads to the Hamiltonian function

$H=\lambda_{0}+\lambda_{1}(t)x_{2}(t)+\lambda_{2}(t)u(t),$ (9)

with $\lambda_{0-2}$ to be the vector of Lagrange multipliers. Minimizing the
Hamiltonian function yields the time-optimal control of the form, cf. with
[16],

$u^{o}(t)=\left\\{\begin{array}[]{ll}+U,&\hbox{for }\lambda_{2}^{o}(t)<0,\\\
0,&\hbox{for }\lambda_{2}^{o}(t)=0,\\\ -U,&\hbox{for }\lambda_{2}^{o}(t)>0,\\\
\end{array}\right.$ (10)

where $\lambda_{2}^{o}(t)$ is the control-related optimal Lagrange multiplier.
Substituting this control law into $x_{2}$-dynamics results in a two-point
boundary value problem, see [16] for details, from which the solution of
optimal Lagrange multiplier gives

$\lambda_{2}^{o}(t)=-c_{1}^{o}\,t+c_{2}^{o}.$ (11)

Finally, the optimal constant values $(c_{1}^{o},c_{2}^{o})\neq(0,0)$ need to
be found such that the two-point boundary value problem is solved. With closer
look on (10), (11) one can recognize that: (i) the motion system (7) is always
fully either accelerated or decelerated, (ii) the calculated multiplier (11)
determines the time instant of exactly one switching between both control
actions, e.g. from $+U$ to $-U$ when $x_{1}(0)<x_{1}(t_{f})$. Such control
behavior is well known as the so-called “bang-bang” control, while the optimal
switching time, resulting from (11), is always depending on the initial and
final states. Here it is worth recalling that the time-optimal control (10),
solved in terms of a switching time instant (11), was formulated as the
Fuller’s problem in [10] and further analyzed for an optimal sliding mode
control in [11], in both cases in terms of the state variables,
correspondingly output tracking error.

For the controlled system (7), (10) we are next to analyze the decelerating
trajectory, i.e. the state trajectory after the time-optimal switching at
$t=\tau$, and that for $u(t)=-U\;\forall\;\tau\leq t\leq t_{f}$. For the sake
of simplicity and without loss of generality we will assume
$x_{1}(t_{f})=x_{2}(t_{f})=0$, that for the rest of the paper, meaning zero
equilibrium point. Consequently, all our following analysis and developments
reduce to the problem of convergence to a stable equilibrium in the origin,
while other control tasks with $\bigl{(}x_{1}(t_{f}),x_{2}(t_{f})\bigr{)}\neq
0$ can be converted to that one by an appropriate transformation of the state
coordinates. For the phase-plane analysis, the system (7) with $u=-U$ can be
rewritten as

$\dot{x}\,d\dot{x}=-Um^{-1}dx,$

and after integration of both sides and substitution of $x_{1}$ and $x_{2}$
one obtains the parabolic trajectory

$x_{1}=-0.5m\,U^{-1}x_{2}^{2}.$ (12)

Figure 1: Decelerating trajectory of time optimal control (10).

Obviously, the upper branch of the $(x_{1},x_{2})$-parabola constitutes the
motion trajectory until $t=t_{f}$. The motion trajectory is steered with the
maximal possible deceleration forced by $u(t)=-U$. An example of the
converging state trajectory, driven by the time-optimal control (10) during
the decelerating phase, is shown in Fig. 1.

##  IV OPTIMAL TERMINAL SLIDING MODE CONTROL

Now, with the preliminaries from Section I and time optimal control summarized
in Section II, we are in the position to formulate the optimal terminal
sliding mode control for the perturbed second-order motion systems.

The proposed optimal terminal sliding surface is given by

$s=x_{1}+\alpha mU^{-1}x_{2}^{2}\,\mathrm{sign}(x_{2}),$ (13)

with $\alpha>0$ to be the single (optional) design parameter. Note that the
proposed switching, correspondingly sliding, surface (13), cf. [10, eq. (53)],
has the same form as the non-singular terminal sliding mode algorithm, see
[14], with $q/p=0.5$ rate cf. Section II. Important to notice is that both
surface shaping factors, $U$ and $m$, do not represent any free parameters to
be tuned. They are rather the inherent physical quantities which are
characterizing the motion system plant (7). One can recognize that the
trajectory (12) (of a time-optimal controlled system) coincides with the
sliding surface (13) for $\alpha=0.5$. Further we will show that for that
case, the boundary layer of the so-called twisting mode appears, cf. [15, 14].
The control value is

$u(s)=-U\mathrm{sign}(s),$ (14)

in which the control gain factor is the maximal possible one, i.e. inherently
limited by the constrained actuator force.

The sufficient condition for existence of the terminal sliding mode on the
surface (13) can be shown by assuming the candidate Lyapunov function

$V=0.5s^{2},$ (15)

for which

$0.5\frac{d}{dt}s^{2}<0$ (16)

is, respectively, required. That means for the trajectory which has once
reached and crossed the surface $s=0$ at $t=t_{r}$, it remains on it
$\forall\>t>t_{r}$. The corresponding motion along $s=0$, denoted as terminal
mode, requires $s\dot{s}<0$ to be always fulfilled. For showing it, we
exemplary consider $x_{2}>0$, thus focusing, yet without loss of generality,
on the upper parabolic branch in the second quadrant of the phase-plane only.
Next we summarize $\alpha mU^{-1}=:C$, for the sake of simplicity, and with
these simplifications in mind the terminal sliding surface (13) can be written
as a parabolic equation $s=x_{1}+Cx_{2}^{2}$. Taking the time derivative and
substituting the control (14) into the dynamics (7) one obtains

$\dot{s}=x_{2}-2m^{-1}CUx_{2}\,\mathrm{sign}(s).$ (17)

Multiplying both sides of (17) with $s$, the existence condition (16) results
in

$x_{2}s-2m^{-1}CUx_{2}|s|<0.$ (18)

Substituting back the plant parameters, instead of $C$, and solving the
inequality (18) with respect to the single free factor $\alpha$, the existence
condition of the optimal terminal sliding mode results in

$\alpha>0.5.$ (19)

Important to notice is that if the design parameter is selected as
$0<\alpha\leq 0.5$, the existence condition (16) of the sliding mode becomes
violated, and the trajectories do not remain on $s=0$ upon crossing.
Nevertheless, the manifold (13) appears further on as a switching surface, and
the trajectories will reach asymptotically the origin according to (7) and
(14). This convergence behavior appears as a twisting mode, in which the
frequency of switching (upon crossing the surface (13) with $0<\alpha\leq
0.5$) increases towards infinity as the trajectory approaches the origin.
Recall that the trajectory converges towards origin while circulating around
it. The convergence of the twisting mode can be directly proved geometrically,
since $\|x_{1},x_{2}\|(t_{s2})<\|x_{1},x_{2}\|(t_{s1})$ for all pairs of two
consecutive switching (correspondingly $s$-crossing) at $t_{s2}>t_{s1}$. This
is valid as long as $\alpha>0$. When $\alpha\rightarrow 0$, the switching
surface $s\rightarrow x_{1}$ and the twisting occurs without convergence,
while representing the boundary case of a harmonic oscillator. The latter is
nothing but the double-integrator with a zero-crossing relay in feedback. The
stability of the twisting mode has also been explicitly addressed in [15] by
means of the Lyapunov’s method. It is also worth noting that for
$0.25\leq\alpha\leq 0.5$, the sliding surface still falls into class of the
Fuller’s problem, cf. [10, 11]. Also one can easily recognize that the
positive $\alpha$ appears, from a plant parameters viewpoint, as a gain factor
which captures the multiplicative uncertainties in $m$ and $U$, cf. (13).

The phase-portrait with trajectories of the system (7), (14) is exemplary
shown in Fig. 2 together with the corresponding surface (13), that for the
different $\alpha$-values.

Figure 2: Sliding surface and trajectories for different $\alpha$ values.

Here the system parameters are assigned as in the numerical examples provided
in Section V, while the initial values
$\bigl{(}x_{1}(t_{0}),x_{2}(t_{0})\bigr{)}=(-1,0)$ are assumed for all three
$\alpha$ values. One can directly notice that once the switching surface is
reached, the control with $\alpha=0.3$, which violates (19), does not stay
sliding on the surface, and the trajectory converges in a twisting mode. Note
that here the chattering, due to non-ideal switching and numerical
computational issues, does not appear until the states reach the origin,
meaning until
$\bigl{(}x_{1}(t_{f}),x_{2}(t_{f})\bigr{)}\rightarrow\mathbf{0}$.

In order to establish the reachability condition, see [17, 5], one consider
the candidate Lyapunov function (15) for which the following conditions should
be satisfied: (i) $V(s)$ is positive definite $\forall\,s$, (ii) $\dot{V}<0$
for $s\neq 0$. For showing the finite-time convergence and, therefore, global
stability the so-called $\eta$-reachability [17] condition

$s\dot{s}<-\eta|s|$ (20)

can be evaluated, in which $\eta>0$ is a small constant which ensures the
reachability condition is satisfied. One slightly modifies the above condition
(ii) to be

$\dot{V}\leq-\gamma\sqrt{V},\quad\gamma=\sqrt{2}\eta.$ (21)

Integrating both sides of inequality (21) one obtains, cf. [5],

$\sqrt{V}(t)\leq-0.5\gamma t+\sqrt{V(0)}.$ (22)

Since reaching the sliding surface at finite-time $t_{r}$ means $V(t_{r})=0$,
one can obtain, out from (22), the following

$t_{r}\leq\gamma^{-1}2\sqrt{V(0)}=\gamma^{-1}2\sqrt{0.5}|s(0)|.$ (23)

One can recognize that a bounded finite time can be guaranteed by (23) and it
depends on the initial value of the sliding manifold only. In order to
evaluate the $\eta$-reachability condition for the control system (7), (13),
(14), we first rewrite (20) as

$\mathrm{sign}(s)\dot{s}<-\eta.$ (24)

Taking the time derivative of (13) and substituting the control system
dynamics (7), (14), one obtains

$\dot{s}=x_{2}-2\alpha|x_{2}|\mathrm{sign}(s).$ (25)

Combining (24) and (25) and separating the variables we achieve the
$\eta$-reachability condition as

$x_{2}\mathrm{sign}(s)<-\eta+2\alpha|x_{2}|.$ (26)

One can see that it is always possible to find a small positive constant
$\eta$ such that the condition (26) is fulfilled $\forall\>x_{2}$.

Once the existence and reachability conditions, (19) and (26) respectively,
are satisfied it is indicative to analyze the control system robustness
against the perturbations. The external and internal perturbations222Here we
refer to the internal perturbations as well, since these may arise from the
plant uncertainties and internal dynamics which are not captured by the motion
system in form of the double-integrator (7). matched by the control will
directly affect the right-hand side of (7) and, in worst case, produce a
counteracting (disturbing) force during the reaching phase, for which
$|u|(t)=U$ during the time $0<t<t_{r}$. On the other hand, the sudden ‘fast’
perturbations, once they are occurring during the sliding mode, will diverge
the trajectories from the sliding manifold, thus giving rise (again) to the
reachability problem. Such transient loss of the sliding mode leads, in turn,
to $u(t)=\mathrm{const}\in\\{-U,U\\}$ during the time $t_{d}<t<t_{r}$, where
$t_{d}>0$ is the last instant of losing the sliding manifold $s$ due to
appearance of the perturbations. In order to guarantee that the optimal
terminal sliding mode control can compensate for the unknown perturbations
$\xi(t)$, the single boundedness condition

$|\xi|<U,$ (27)

is required for all times $t>0$ of the control operation.

In the following, we will demonstrate different examples of the optimal
terminal sliding mode control applied to the perturbed motion system. For that
purpose, the system plant is extended by the matched upper bounded
perturbation (27), thus resulting in the modified (7) which becomes

$\ddot{x}(t)=m^{-1}\bigl{(}u(t)+\xi(t)\bigr{)}.$ (28)

##  V NUMERICAL EXAMPLES

The assumed numerical parameters, here and for the rest of the paper, are
$m=0.1$, $U=1$, $\alpha=0.6$. All the controlled motion trajectories start at
the initial state $\bigl{(}x_{1}(t_{0}),x_{2}(t_{0})\bigr{)}=(-1,0)$, while
the control set point is placed in the origin, meaning
$\bigl{(}x_{1}(t_{f}),x_{2}(t_{f})\bigr{)}_{ref}=\mathbf{0}$. The implemented
and executed simulations are with 1kHz sampling rate and the (most simple)
first-order Euler numerical solver.

###  V-A Motion system with nonlinear friction

The motion system dynamics (7) is first considered as being affected by the
nonlinear Coulomb friction $f(x_{2},t)$ which, moreover, includes the
continuous presliding transitions, see [18] for details. It is worth
emphasizing that the modeled Coulomb friction force is not discontinuous at
zero velocity crossing, while a saturated friction force at steady-state
results in $F_{c}\,\mathrm{sign}(\dot{x})$, where $F_{c}>0$ is the Coulomb
friction constant. The friction force, which appears as an internal
perturbation $\xi=-f$, is assumed to be bounded by the half of the control
amplitude, i.e. $F_{c}=0.5$.

Figure 3: Phase portrait of state trajectory (a) and convergence of output
tracking variable $x_{1}(t)$ and down-scaled $x_{2}(t)$-state over
perturbation (b).

The controlled state trajectory is shown in Fig. 3 (a), while the region
around origin is additionally zoomed-in for the sake of visualization. One can
see that after the state convergence, in the sliding mode, a stable low-
amplitude ($\times 10^{-5}$) limit cycle around zero appears due to the by-
effects caused by the nonlinear presliding [18] friction. This is not
surprising since $F_{c}=0.5U$, and no explicit presliding friction
compensation is performed by the terminal sliding mode control which is only
switching between $u=\pm U$. The convergence of both state variables ($x_{2}$
is down-scaled by factor 0.1 for the sake of visualization) is shown in Fig. 3
(b) over the progress of frictional perturbation. It can also be side-noted
that the feasibility of the Coulomb friction compensation by means of a
classical (first-order) sliding mode control with discontinuous control action
was experimentally demonstrated, in combination with a time-optimal bang-bang
strategy, in a case study presented in [19].

###  V-B Motion system with harmonic perturbation

Figure 4: Phase portrait of state trajectory (a) and convergence of output
tracking variable $x_{1}(t)$ and down-scaled $x_{2}(t)$-state over
perturbation (b).

As next, we assume an internal (or external) perturbation to be a periodic
function of time, modeled by $\xi(t)=0.5\sin(20t)$. Note that such disturbing
harmonic forces (or torques) may occur in various types of the motors, or more
generally actuators, for example see [20]. The sinusoidal amplitude is
selected to be the half of the control value, and the angular frequency of 20
rad/s is assigned for providing a sufficient number of periods during both the
transient and steady-state response. The controlled state trajectory and the
converged state and perturbation values are shown in Fig. 4. Also the zoomed-
in steady-state output error is plotted.

###  V-C Motion system with random binary perturbation

Finally, we consider the motion system affected by an external stepwise
perturbation which is appearing in form of a random binary signal. The
perturbation amplitude is kept at the same level as in both previous examples
(cf. Sections V-A and V-B), so that the perturbed plant (28) has
$\xi(t)\in\\{-0.5,0.5\\}$ which is a random realization of the binary
variable.

Figure 5: Phase portrait of state trajectory (a) and convergence of output
tracking variable $x_{1}(t)$ and down-scaled $x_{2}(t)$-state over
perturbation (b).

The controlled state trajectory is shown in Fig. 5 (a), and, alike, in (b) the
states convergence and the perturbation value are plotted together as the time
series. One can recognize that the motion trajectory is highly perturbed. Yet,
both states converge to zero, piecewise in the sliding mode and piecewise in
the twisting-like mode. The latter appears transiently each time the binary
perturbation changes the sign. It can be stressed that from all three types of
the perturbation examples, this one constitutes the ‘worst’ case which is
largely affecting the control system. This is not surprising since the
differential changes of the $\xi$-amplitude are of the same magnitude as $U$.
The stepwise perturbation $\xi(t)$ does not allow system to continuously
staying in the sliding mode and forces it multiple times into the reaching
phase. This is directly reflected in the shape of $(x_{1},x_{2})$ trajectory,
see Fig. 5 (a). At the same time, one can recognize that the states remain in
vicinity to the origin, once reaching it and despite the further stepwise
changes of $\xi(t)$, see Fig. 5 (b) for the times $t>1.25$ sec.

##  VI CONCLUDING REMARKS

This paper described and analyzed the optimal terminal sliding mode control
aimed for the second-order motion systems with matched and bounded
perturbations. The control scheme uses the time-optimal trajectory of the
double-integrator dynamics, for which only one dedicated switching provides
the maximal possible acceleration and deceleration in presence of the inherent
control bounds. The proposed optimal terminal sliding mode has the single free
(and rather optional) design parameter $\alpha$, which allows for both, the
terminal and twisting control modes. The time-optimal twisting mode on the
boundary layer, i.e. for $\alpha=0.5$, is in particular interesting for those
motion control applications where a frequent sliding mode switching of the
control signal may be undesirable, like in the hydraulic control valves [21].
Here we however abstain from discussions on the pros and cons of the
discontinuous and continuous sliding mode controllers, see [22] for details.
We focus rather on the analysis and prove of existence and reachability
conditions of the proposed optimal terminal sliding mode control. The shown
developments and the determined optimal sliding surface are fully inline with
the Fuller’s problem. It was shown that for a motion dynamics with the given
inertial and actuator saturation parameters, there is no other time-optimal
surface than (13). All $\alpha>0$ provide a stable state convergence to the
origin, while the $\alpha$-values around boundary layer (i.e. $\alpha=0.5$)
are particularly interesting for the control application. These allow finding
a trade-off between a time-optimal convergence with minimal switching and
ability to capture multiplicative uncertainties of the inertial and saturation
parameters. Three illustrative examples of the bounded internal and external
perturbations demonstrated the control efficiency in a numerical setup.

## References

  * [1] V. T. Haimo, “Finite time controllers,” _SIAM Journal on Control and Optimization_ , vol. 24, no. 4, pp. 760–770, 1986.
  * [2] E. Moulay and W. Perruquetti, “Finite-time stability and stabilization: State of the art,” in _Advances in Variable Structure and Sliding Mode Control_. Springer, 2006, pp. 23–41.
  * [3] V. I. Utkin, _Sliding modes in contr. and optim._ Springer, 1992.
  * [4] W. Perruquetti and J.-P. Barbot, _Sliding mode control in engineering_. CRC press, 2002.
  * [5] Y. Shtessel, C. Edwards, L. Fridman, and A. Levant, _Sliding mode control and observation_. Springer, 2014.
  * [6] M. Zak, “Terminal attractors for addressable memory in neural networks,” _Physics Letters A_ , vol. 133, no. 1-2, pp. 18–22, 1988.
  * [7] S. T. Venkataraman and S. Gulati, “Control of nonlinear systems using terminal sliding modes,” in _American Control Conference_ , 1992, pp. 891–893.
  * [8] Z. Man, A. P. Paplinski, and H. R. Wu, “A robust mimo terminal sliding mode control scheme for rigid robotic manipulators,” _Transactions on automatic control_ , vol. 39, no. 12, pp. 2464–2469, 1994.
  * [9] Y. Feng, X. Yu, and Z. Man, “Non-singular terminal sliding mode control of rigid manipulators,” _Automatica_ , vol. 38, no. 12, pp. 2159–2167, 2002\.
  * [10] A. Fuller, “Relay control systems optimized for various performance criteria,” _IFAC Proceedings Volumes_ , vol. 1, no. 1, pp. 520–529, 1960\.
  * [11] F. Dinuzzo and A. Ferrara, “Higher order sliding mode controllers with optimal reaching,” _IEEE Transactions on Automatic Control_ , vol. 54, no. 9, pp. 2126–2136, 2009.
  * [12] A. Bryson, _Applied Optimal Control: Optimization, Estimation and Control_. CRC Press, 1975.
  * [13] A. Levant, “Sliding order and sliding accuracy in sliding mode control,” _International journal of control_ , vol. 58, no. 6, pp. 1247–1263, 1993\.
  * [14] L. Fridman, J. A. Moreno, B. Bandyopadhyay, S. Kamal, and A. Chalanga, “Continuous nested algorithms: The fifth generation of sliding mode controllers,” in _Recent advances in sliding modes: From control to intelligent mechatronics_. Springer, 2015, pp. 5–35.
  * [15] T. Sánchez and J. A. Moreno, “Lyapunov functions for twisting and terminal controllers,” in _13th International Workshop on Variable Structure Systems (VSS)_. IEEE, 2014, pp. 1–6.
  * [16] H. Geering, _Optimal control with eng. applic._ Springer, 2007.
  * [17] C. Edwards and S. Spurgeon, _Sliding mode control: theory and applications_. CRC Press, 1998.
  * [18] M. Ruderman, “On break-away forces in actuated motion systems with nonlinear friction,” _Mechatronics_ , vol. 44, pp. 1–5, 2017.
  * [19] M. Ruderman and M. Iwasaki, “Hybrid positioning control using time-optimal feed-forwarding and sliding mode feedback,” in _France-Japan and Europe-Asia Congress on Mecatronics_. IEEE, 2014, pp. 189–194.
  * [20] M. Ruderman, A. Ruderman, and T. Bertram, “Observer-based compensation of additive periodic torque disturbances in permanent magnet motors,” _IEEE Transactions on Industrial Informatics_ , vol. 9, no. 2, pp. 1130–1138, 2013.
  * [21] M. Ruderman, L. Fridman, and P. Pasolli, “Virtual sensing of load forces in hydraulic actuators using second-and higher-order sliding modes,” _Control Engineering Practice_ , vol. 92, p. 104151, 2019.
  * [22] U. Pérez-Ventura and L. Fridman, “When is it reasonable to implement the discontinuous sliding-mode controllers instead of the continuous ones? frequency domain criteria,” _International Journal of Robust and Nonlinear Control_ , vol. 29, no. 3, pp. 810–828, 2019.