11institutetext: Leibniz-Institut für Sonnenphysik, Schöneckstr. 6, 79104
Freiburg, Germany, 11email<EMAIL_ADDRESS>22institutetext: High
Altitude Observatory, NCAR, P.O. Box 3000, Boulder, CO 80307, USA

# On the (in)stability of sunspots††thanks: This paper is mainly based on Part
I of the Ph.D. thesis “On the decay of sunspots”, https://freidok.uni-
freiburg.de/data/165760.

H. Strecker 11 W. Schmidt 11 R. Schlichenmaier 11 M. Rempel 22

(Received: 22 December 2020 / Accepted: 04 March 2021)

###### Abstract

Context. The stability of sunspots is one of the long-standing unsolved
puzzles in the field of solar magnetism and the solar cycle. The thermal and
magnetic structure of the sunspot beneath the solar surface is not accessible
through observations, thus processes in these regions that contribute to the
decay of sunspots can only be studied through theoretical and numerical
studies.

Aims. We study the effects that destabilise and stabilise the flux tube of a
simulated sunspot in the upper convection zone. The depth-varying effects of
fluting instability, buoyancy forces, and timescales on the flux tube are
analysed.

Methods. We analysed a numerical simulation of a sunspot calculated with the
MURaM code. The simulation domain has a lateral extension of more than 98 Mm
$\times$ 98 Mm and extends almost 18 Mm below the solar surface. The analysed
data set of 30 hours shows a stable sunspot at the solar surface. We studied
the evolution of the flux tube at defined horizontal layers (1) by means of
the relative change in perimeter and area, that is, its compactness; and (2)
with a linear stability analysis.

Results. The simulation shows a corrugation along the perimeter of the flux
tube (sunspot) that proceeds fastest at a depth of about 8 Mm below the solar
surface. Towards the surface and towards deeper layers, the decrease in
compactness is damped. From the stability analysis, we find that above a depth
of 2 Mm, the sunspot is stabilised by buoyancy forces. The spot is least
stable at a depth of about 3 Mm because of the fluting instability. In deeper
layers, the flux tube is marginally unstable. The stability of the sunspot at
the surface affects the behaviour of the field lines in deeper layers by
magnetic tension. Therefore the fluting instability is damped at depths of
about 3 Mm, and the decrease in compactness is strongest at a depth of about 8
Mm. The more vertical orientation of the magnetic field and the longer
convective timescale lead to slower evolution of the corrugation process in
layers deeper than 10 Mm.

Conclusions. The formation of large intrusions of field-free plasma below the
surface destabilises the flux tube of the sunspot. This process is not visible
at the surface, where the sunspot is stabilised by buoyancy forces. The onset
of sunspot decay occurs in deeper layers, while the sunspot still appears
stable in the photosphere. The intrusions eventually lead to the disruption
and decay of the sunspot.

###### Key Words.:

Methods: simulations–Sun: sunspots – Sun: photosphere – convection zone

##  1 Introduction

The lifetime of sunspots covers a broad range of a few days to a few months
(Martínez Pillet, 2002; van Driel-Gesztelyi & Green, 2015). The observational
evolution and decay process of sunspots distinguishes a photometric and a
magnetic decay (Martínez Pillet, 2002). The origin of the differences in
lifetime and the flux removal is still not completely understood. In addition,
the onset of the decay process is not clear. In some cases, a sunspot already
starts to decay, that is, loose magnetic flux, while it might at the same time
gain magnetic flux by coalescence (Solanki, 2003). This means that the process
that leads to the decay of a sunspot operates already before the sunspot has
fully developed (McIntosh, 1981). The opaqueness of the photosphere inhibits
direct observations of processes below the solar surface. This limits our
ability of fully understanding sunspot evolution from direct observations.
Numerical models are used to obtain a more complete image of sunspots. A
variety of different models, especially static ones, has been presented over
the years (for an overview of different (static) sunspot models, we refer to
Jahn, 1997). The basic assumption for a magnetohydrostatic model is that the
sunspot is a monolithic, single magnetic flux tube in and below the
photosphere (Schlüter & Temesváry, 1958). For a better agreement of the
modelled sunspots with observations, for example, a sharp boundary between the
sunspot and the surroundings in intensity, Simon & Weiss (1970) implemented a
current sheet at the boundary of the flux tube. This current sheet, also
called magnetopause, balances the internal pressure with the pressure in the
surroundings. The discontinuity of the pressure across the current sheet is
correlated to a jump in the field strength (Jahn, 1992; Jahn & Schmidt, 1994).

To study the stability of sunspots in the convection zone, Meyer et al. (1977)
adapted a magnetohydrostatic model. They embedded a flux tube with an
axisymmetric meridional field $(B_{r},0,B_{z})$ in the non-magnetised plasma
of the upper convection zone. For the magnetic flux tube to be in equilibrium
with the surrounding plasma, they assumed the two regions to be separated by a
surface S (see Fig. 1). This surface is a current sheet that causes the abrupt
drop of the magnetic field and thereby horizontally balances the pressure of
the flux tube, $p_{\text{f}}$, and of the surrounding plasma, $p_{\text{p}}$.
Meyer et al. (1977) assumed that instabilities of the surroundings, for
instance, convection, do not affect the stability of this system. In addition,
the plasma in the tube was stably stratified and magnetohydrodynamically
stable. Then any perturbation leading to an instability depends on the
equilibrium balance at the surface, S. To describe the stability, they took
the decreasing pressure difference ($p_{\text{p}}-p_{\text{f}}$) with height
into account. In addition, the fanning-out of the field with height, which
causes a change in inclination angle, $\chi$, of the surface S to the vertical
(see Fig. 1) and a change in effect of the gravitational acceleration, ${\bf
g}$, onto the magnetised plasma. Finally, the curving of the field causes
magnetic tension forces, and the criterion for stability is written as

$(\rho_{\text{qs}}-\rho_{\text{s}})\cdot
g\cdot\text{sin}\chi-\frac{B^{2}}{4\pi R_{c}}>0,$ (1)

with $R_{c}$ the radius of curvature of the surface S. For a more detailed
derivation of the criterion, we refer to Meyer et al. (1977) and Strecker
(2020). Based on this inequality, a magnetic flux tube in the upper layers of
the convection zone has to balance two main effects to remain stable: (1) The
plasma of the magnetic flux tube is less dense than the surrounding plasma,
therefore the spot floats on the quiet Sun. (2) The concave curvature of the
magnetic field lines has a destabilising effect and causes the tube to become
vulnerable to interchange instability.

Figure 1: Geometry of a magnetic flux tube (V${}_{\text{f}}$) embedded in
field-free plasma (V${}_{\text{p}}$) (image adapted from Meyer et al., 1977).

Interchange instability is a well-known process in which a cylindrical
magnetic flux tube that is perturbed at its outer boundary is deformed.
Bundles of field lines are transferred into the non-magnetised neighbourhood.
On exchange, regions that were originally occupied by magnetic field lines are
occupied by non-magnetised plasma. The individual volumes remain unchanged in
this interchange process. The cylindrical flux tube obtains a rippled
interface, similar to a fluted column. Therefore the process is also called
fluting instability. The direction of the curving of the field lines towards
the non-magnetised plasma affects the stability. If the curvature is concave,
the rippled interface increases the cross-section and hence destabilises the
flux tube further. A detailed explanation of interchange instability can be
found in Priest (2014).

Schüssler (1984) developed the model of Meyer et al. (1977) further to study
the effect of velocity fields on the instability of flux tubes. He found that
converging convective flow cells surrounding the flux tube and flows along the
field lines within the flux tube do not have a stabilising effect. Instead,
their contribution can be destabilising. Only flows with an azimuthal
component around the magnetic flux tube have the ability to stabilise. Parker
(1979) made use of the stability criterion of Meyer et al. (1977). At a
sufficiently large depth below the surface, field lines are able to separate
from the spot through interchange instability. Parker suggested that a sunspot
is not a monolithic, single flux tube below the solar surface. Instead, at
about 1Mm below the photosphere, the monolithic structure divides into several
smaller flux bundles. Buoyancy and convective downdrafts allow these bundles
to bunch together towards the photosphere, and stabilise them against
interchange processes (Parker, 1979; Jahn, 1992). This model is often referred
to as the spaghetti or jelly fish model.

Studies performed by means of the photometric evolution of the sunspot take
advantage of the change in area of the sunspot in the photosphere in time. A
variety of studies have been performed, resulting in different decay processes
and scenarios. A linear decay law was found by Bumba (1963) for recurrent
sunspots and was found to fit the distribution in later studies. Martínez
Pillet et al. (1993) proposed a parabolic decay rate, with the area being a
quadratic function in time. However, they noted that it is difficult to
distinguish linear or quadratic decay rates based on the observations.
Petrovay & van Driel-Gesztelyi (1997) also proposed a parabolic decay law for
their observed sunspots. The determination of a decay law is important to
understand the decay process of the sunspots. A linear decay implies that the
loss of flux occurs everywhere within the sunspot through turbulent diffusion
of the magnetic field (e.g., Meyer et al., 1974; Solanki, 2003). Instead, a
parabolic decay law could be caused by erosion of the sunspot at its boundary,
for instance, by supergranular motions (Simon & Leighton, 1964) or turbulent
motions on a granular-size scale (Petrovay & Moreno-Insertis, 1997). These
models shift the question of an area change during sunspot decay to possible
scenarios regarding the loss of magnetic flux from the sunspots. A loss of
magnetic flux across the whole area of the sunspot is supported by the
appearance of moving magnetic features (MMFs). These appear in the non-
magnetised moat region surrounding a sunspot and were first reported by
Sheeley (1969). Harvey & Harvey (1973) were the first to propose that the
magnetic patches pull off magnetic flux from the sunspot while moving away.
However, later investigations found imbalances between the sunspot magnetic
flux loss and the magnetic flux transported by MMFs (e.g., Martínez Pillet,
2002; Kubo et al., 2007). Different types of MMFs were proposed to exist, and
were classified according to their polarity (Shine & Title, 2000). Unipolar
features of both polarities are observed as well as bipolar features. Martínez
Pillet (2002) proposed that only the unipolar features with the same polarity
as the sunspot are related to the loss of flux from the sunspot, that is,
sunspot decay. Unipolar features with opposite polarity as the sunspot and
bipolar features might be extensions of penumbral filaments. Such a scenario
was first proposed by Zhang et al. (1992). Thomas et al. (2002) described such
a scenario in further detail, with the bipolar features showing a sea-serpent
structure. Kubo et al. (2007), Kubo (2012), and Verma et al. (2012) also found
that the total magnetic flux of the MMFs is far higher than the flux loss of
the sunspot. Therefore the polarity of the MMFs has to be considered to
determine their relation to sunspot decay. In addition, Rempel (2015) found in
simulations that the emergence and submergence of horizontal field in the moat
region also contributes substantially to the flux budget of a sunspot.
Specifically, these contributions tend to offset the flux transport due to
MMFs. This might be an explanation for observations of MMFs in the
surroundings of stable sunspots that last for several months (van Driel-
Gesztelyi & Green, 2015). The relation of MMFs to the decay process of
sunspots is still unclear.

The knowledge that the flux tube of a sunspot evolves below the photosphere
might lead to insights into the decay process. However, the models described
above are static. The advent of new computers with more power and larger
computational domains enabled setting up realistic magnetohydrodynamic (MHD)
simulations in the past years. In this study we make use of a 3D MHD
simulation of a sunspot to study the evolution of a sunspot in the upper
convection zone and in the photosphere. In Sec. 2 we describe the setup and
characteristics of the simulation domain that hosts the sunspot. In addition,
we determine the flux tube boundary. We analyse the flux tube in Sec. 3 by
studying the evolution of its structure, and by performing a linear stability
analysis in Sec. 4. In both sections we describe the respective methods and
results separately. The results of both studies are combined and discussed in
Sec. 5. Finally, a picture of the stability of the flux tube in the upper
convection zone and photosphere and a decay scenario of sunspots is presented
in Sec. 6.

##  2 Data and methods

We analysed a sunspot in a 3D MHD simulation. It was calculated with the MURaM
code by Rempel (2015). He calculated two different simulations. We used the
simulation that is calculated with a modified numerical diffusivity. For a
more detailed description of the calculation and the specification, we refer
the reader to Section 2 of Rempel (2015).

The simulation domain had a size of 98.308 $\times$ 98.308 $\times$ 18.432 Mm3
and was computed on a grid with a 48 $\times$ 48 $\times$ 24 km3 cell size.
The analysed data set had an increased grid cell size of 96 km in horizontal
direction and a grid spacing of 48 km in vertical direction. Thus, the
analysed 3D domain was composed of 1024 $\times$ 1024 grid cells horizontal
and 384 grid cells vertical. The data show the sunspot with an average
diameter of 30 Mm, the surrounding moat flow with an average extension of 10
Mm at the surface, and the surrounding quiet Sun. The vertical extension of
the simulation domain enabled us to cover sub-photospheric regions with their
own characteristics, such as convective motions. The simulation used the setup
of Rempel (2012), in which a more horizontal field is imposed through the top
boundary in order to create an extended penumbra.

The complete simulation covered a duration of 100 hours. In this analysis we
considered a section of 30 hours. This section covers the time range from
$t_{\text{sim}}=50$ h until $t_{\text{sim}}=79.75$ h, with a cadence of 15
minutes. For the analysis we used the 3D data of the temperature, $T$,
pressure, $p$, density, $\rho$, the three components of the magnetic field
($B_{x}$, $B_{y}$, $B_{z}$), and the velocity ($v_{x}$, $v_{y}$, $v_{z}$). In
addition, 2D maps of the magnetic field strength and intensity for selected
layers of constant optical depth, $\tau$, were used.

The vertical extension of the simulation domain of 18.432 Mm covers the
photosphere and uppermost part of the convection zone. The two regions are
attempted to be separated by a solar surface. To define a constant layer
within the gridded simulation domain and define it as the solar surface, we
determined an average $\tau$ = 1 layer. To do this, the optical depth along
vertical rays was calculated from the temperature and pressure of each grid
cell using the tabulated Rosseland mean opacity, which was also used to
advance the grey simulation in the first place. The Rosseland opacity was
computed from Kurucz (1993) opacity tables. The $z$-layer, in the following
termed solar surface, is defined as the mean $z$-position where $\tau$ = 1. A
restriction was made by only taking into account cells outside the sunspot
region, that is, the moat and quiet-Sun region. This layer was defined as $z$
= 0 Mm. Thus, the convection zone extends down to $z$ = -17.71 Mm with
negative $z$-values. Regions in the photosphere have positive $z$-values with
an extension up to 0.720 Mm above the solar surface.

In this analysis the evolution of the magnetic flux tube of the sunspot is
studied from the solar surface downward through the convective region. In this
region, we focus on defined depths in steps of $\Delta z$ = 0.625 Mm from $z$
= 0 Mm down to $z$ = -15 Mm.

###  2.1 Convective timescale

The evolution of the flux tube of the sunspot is affected by the surrounding
regions. We studied the time evolution of the flux tube at different depths in
time. Therefore the changing timescales of convective motion within the
convection zone had to be considered. In the following we determine the
convective timescale for different depths of the simulation domain using the
mixing length theory (see e.g. Stix, 2002). In addition, we compare the
obtained timescale with the corresponding values of the standard solar model.

Based on the mixing length theory, we can describe the evolution of convective
motions with the convective timescale, $t_{\text{conv}}$. It is calculated
from the mixing length, $l_{M}$, and the convective velocity, $v$, as

$t_{\text{conv}}=\frac{l_{m}}{v}=\frac{\alpha H_{p}}{v}.$ (2)

The mixing length can be approximated by $l_{M}=\alpha H_{p}$ with the mixing
length parameter, $\alpha$, which is of order unity, and the pressure scale
height, $H_{p}=p/(\rho g)$. Within the mixing length theory, the convective
velocity, $v$, is described as the average velocity of a rising bubble through
a defined layer. To determine the convective timescale within the simulation
domain, the thermodynamic quantities were determined in steps of $\Delta z$ =
0.625 Mm starting at $z$ = 0 Mm down to $z$ = -15 Mm, and additionally, the
depths $z$ = -16.25 Mm and $z$ = -17.5 Mm were considered. Only the region
outside the sunspot was used to determine the thermodynamic properties of the
simulated convection zone. There, in quiet-Sun convection, $\alpha=1$ is a
good approximation. For the velocity, the absolute value of the vertical
component $v_{z}$ was used. For each depth position, we first determined the
mean over nine layers (four layers above and four below the respective layer)
and finally the mean over all time steps. The resulting local timescale for
the different depths in the convective region of the simulation domain are
shown in Fig. 2 as red crosses. The increase in timescale with depth
represents a slower evolution of perturbations with deeper layers. Blue stars
in Fig. 2 show the convective timescale calculated with values based on the
standard solar model as provided by Stix (2002). A comparison of the values
from the simulation and the standard solar model shows that the thermodynamic
quantities of the two agree well.

Figure 2: Convective timescale $t_{\text{conv}}$ in the upper part of the
convection zone from the MHD simulation (red crosses) and a standard solar
model (blue stars) (values are taken from Stix, 2002).

###  2.2 Determination of the boundary of the sunspot flux tube

The sunspot region needs to be clearly defined when the structure of the
magnetic flux tube of the sunspot is to be analysed. To do this, we used the
magnetic field strength, $B=\sqrt{B_{x}^{2}+B_{y}^{2}+B_{z}^{2}}$. At each of
the 25 depth positions from $z$ = 0 Mm to $z$ = -15 Mm, an individual value of
the magnetic field strength was determined as a contour value to define the
boundary of the magnetic flux tube.

We applied two different methods to determine the contour values. The first
method is an estimate of the boundary by eye. The values $B_{\text{eye}}(z)$
were chosen to obtain the best match between the contours and the magnetic
flux tube at the individual depths. To this end, contours for different $B$
values were visualised on horizontal $B$ maps, clipped to $B\leq 6$ kG, for
arbitrary time steps. We then chose the contour that best represented the
visual contrast between the magnetic flux tube and the non-magnetic
surroundings. The obtained values, determined in steps of $\Delta z$ = 1.25
Mm, are shown in Fig. 3 as blue stars and are listed in Table 1 in Appendix A.
This manual method is neither consistent nor applicable to the large number of
contours.

As second method we applied an automatic procedure to determine the threshold
$B_{\rm c}(z)$. The maximum magnetic field strength $B^{\text{max}}(z,t)$ was
determined for each grid layer and each time step. $B^{\text{max}}(z,t)$ was
boxcar-smoothed with a width of 240 km in vertical direction and a width of 75
minutes in time around a defined depth position in steps of $\Delta z$ = 0.625
Mm. From the smoothed values, averages, $\langle
B^{\text{max}}_{\text{smooth}}\rangle_{t,z\pm 5}$, over all times and ten
adjacent depth layers were calculated. These values were divided by
$2\text{e}$ to yield the threshold $B_{\rm c}(z)$,

$B_{\text{c}}(z)=\frac{\langle B^{\text{max}}_{\text{smooth}}\rangle_{t,z\pm
5}}{2e}.$ (3)

We made tests for multiples of $e$ in visualising the different contours on
horizontal $B$ maps for three different times of 0 h, 15 h, and 27.75 h and
for depths in steps of $\Delta z=1.25$ Mm. The best representation of the
boundary of the magnetic flux tube throughout the different depth positions
was found for a fraction of $2e$ of $\langle
B^{\text{max}}_{\text{smooth}}\rangle_{t,z\pm 5}$. With this choice, the
contours from $B_{\text{c}}$ and $B_{\text{eye}}$ match fairly well, although
they are determined independent of each other, but differences exist in a few
cases.

The obtained $B_{\text{c}}(z)$-values are shown as red crosses in Fig. 3
overlaid on the temporal averaged values of the maximum magnetic field
strength for each vertical grid layer, $B^{\text{m,t}}(z)/2e$ (black dots).
The calculated contour values $B_{\text{c}}(z)$ are given in Table 1 in
Appendix A as well.

Figure 3: Contour values of the magnetic field strength ($B_{\text{eye}}(z)$
and $B_{\text{c}}(z)$) used to determine the sunspot boundary by eye (blue)
and from the maximum magnetic field strength (red).

We computed the area, $A$, and perimeter, that is, the length of the contour,
$C$, of the sunspot contours at each depth and time step as follows: We
smoothed each map of the magnetic field strength with a width of 5 in $x$\-
and $y$\- direction. We applied the values $B_{\text{eye}}(z)$ and
$B_{\text{c}}(z)$ to determine the contours.

##  3 Compactness of the sunspot

Because of the filamentary structure of the penumbra, the sunspot is not
perfectly circular, as shown in the intensity map in Fig. 4, which corresponds
to $t$ = 0 h. A contour, outlining the boundary of the sunspot, is drawn at an
intensity level $I_{\text{c}}=0.9I_{\text{n}}$, with the normalised intensity,
$I_{\text{n}}=I/\bar{I}_{\text{qs}}$, where $\bar{I}_{\text{qs}}$ denotes the
mean intensity of the quiet Sun. A similar structure can be found in maps of
the magnetic field strength at 1.25 Mm below the surface. This is shown for
the same time step in Fig. 5 a. The outer boundary of the flux tube, shown in
red, is determined with the automated procedure as described in Sec. 2.2. At
this depth, the flux tube is embedded in a mesh-like structure that covers the
moat. This structure terminates at the network boundary. This network boundary
is seen throughout all depth layers, as displayed in Fig. 5 b - f, meaning
that the sunspot is embedded in the centre of a moat cell (see also Rempel,
2015). We note that in Fig. 5 the field of view is reduced to focus on the
sunspot and the moat.

Field-free regions penetrate the flux tube. In deeper regions, the average
diameter of the flux tube decreases and the field strength increases. In
addition, the boundary of the tube becomes more ragged, and the indentations
at the boundary penetrate deeper into the inner part of the flux tube. This is
most prominent at depths of $z$ = -3.75 Mm and $z$ = -6.25 Mm, as shown in
Fig. 5 b and c, respectively. The flux tube is more compact in deeper regions,
that is, the outer boundary becomes smoother (see Figs. 5 d - f).

Figure 4: Intensity map at $\tau=1$ with a contour line at
$I_{\text{c}}=0.9I_{\text{n}}$ indicating the boundary of the sunspot. Figure
5: Maps of the magnetic field strength at $t=0$ h, shown at different depths
(indicated in red in each panel). The boundary of the sunspot (red line) is
defined by depth-dependent values of the magnetic field strength.

In time, the ragged structure of the flux tube becomes more pronounced. This
is shown in Fig. 6 for a depth of $z$ = -7.5 Mm with a cadence of 6 hours
between the snapshots, covering the 30 hours of the analysed data. An
animation with a cadence of 15 minutes for selected depths is provided as
additional material111The movies associated to Fig. 6 are available at
http://www.leibniz-
kis.de/fileadmin/user_upload/kis/media/promotionen/2020_strecker/BMaps.mp4.

Figure 6: Maps of the magnetic field strength at $z$ = -7.5 Mm at different
times (indicated in red in each panel). The boundary of the sunspot (red line)
is defined by the contour value $B_{\text{c}}=4781$ G.

At the beginning of the analysis, at $t$ = 0 h (see Fig. 6 a) the inner part
of the flux tube is mainly undisturbed. Within 6 hours, regions of weaker
field appear in the innermost part of the flux tube while the outer structure
becomes more ragged, as is shown in Fig. 6 b. The increasing raggedness causes
a degradation of the flux tube. This degradation process continues in time, as
Figs. 6 c - e show. At 29.75 h, the roundish structure of the flux tube has
completely vanished. The degradation process takes place in all regions deeper
than 1 Mm below the surface (see the animation provided as additional
material).

###  3.1 Method

For a quantitative description of the degradation process, we define the
compactness, $K$, as the ratio of the circumference, $2\sqrt{\pi A}$, that
corresponds to a circle with the area $A$, and the actual perimeter, $C$,
around the magnetised area $A$, that is, the length of the contour with area
$A$,

$K=2\sqrt{\pi}\frac{\sqrt{A}}{C}.$ (4)

This enables a comparison of different depths as the average diameter of the
flux tube decreases with depth. Thus, a circular flux tube has $K=1$. The
compactness $K(z,t)$ was calculated for all time steps and all selected
depths.

For an increasing raggedness, as discussed in context of Fig. 6, the length of
the contour increases. This causes a decrease in the compactness $K$ in time,
as described by Eq. (4).

###  3.2 Results

The compactness in time for eight selected depths is shown in Fig. 7. The
compactness evolution for all other analysed depths is shown in Fig. 15 and 16
in Appendix B. The determination of the sunspot boundary might affect the
results. Therefore we analysed a sunspot boundary with both contour values,
$B_{\text{c}}$ and $B_{\text{eye}}$ (see Sec. 2.2 and Fig. 3, red crosses and
blue stars, respectively). However, the analysis with $B_{\text{eye}}$ was
only performed for depths in steps of $\Delta z$ = 1.25 Mm. Therefore no
values are available at a depth of $z$ = -1.875 Mm, for instance. If not
explicitly noted otherwise, the following analysis describes the results
obtained with $B_{\text{c}}(z)$.

Figure 7: Compactness of the sunspot (normalised area-perimeter ratio) in
time, shown at different depths. The boundary of the sunspot is determined
with depth-dependent contour values $B_{\text{c}}(z)$ (red crosses) and
$B_{\text{eye}}(z)$ (blue stars). Degradation times are calculated from linear
fits (magenta and cyan lines).

At the surface, the compactness of the flux tube is almost constant (see Fig.
7 a). The scattering of the values is caused by the arbitrary change in
filamentary structure of the penumbra. At $z$ = -0.625 Mm the compactness is
much smaller than at other depths (see Fig. 15 a). Here, the area of the
sunspot is overestimated. The contour value $B_{\text{c}}$(-0.625 Mm) is not
able to reproduce the sunspot boundary. It includes small magnetic structures
of the surrounding region. Thus, a higher value for the area leads to a
smaller compactness (see Eq. (4)). In conclusion, these results and the visual
impression of the sunspot in maps of the magnetic field strength lead to the
picture of an overall stable sunspot above $z$ = -1 Mm. This agrees with the
evolution of the sunspots magnetic flux at the surface, which is almost
stationary, according to Rempel (2015). Below $z$ = -1 Mm, the compactness
decreases, as we show in Fig. 7 b - h and Fig. 15 b - h and 16 in Appendix B.

The initial compactness differs for different depths. A minimum compactness is
found 3.75 Mm below the surface (see Fig. 7 c). Towards higher and deeper
regions the compactness is larger. The maximum compactness is found at the
lowest analysed depth, that is, 15 Mm below the surface. This means that at
the beginning of the analysis, the sunspot structure is ragged most strongly
at 3.75 Mm below the surface. This agrees with the first visual impression
obtained from maps of the magnetic field strengths as described above and
shown in Fig. 5.

A similar compactness is found for an analysis with different contour values,
that is, $B_{\text{c}}(z)$ and $B_{\text{eye}}(z)$. The compactness evolutions
for the individual depths resemble each other. We can therefore assume that
the exact determination of the sunspot boundary does not substantially affect
the analysis of the sunspot compactness. The decrease in compactness in time
for all depths reproduces the increase in raggedness, as described above. The
evolution of the compactness can be approximated by a linear fit (magenta and
cyan lines in Fig. 7), despite the non-continuous decrease. The gradients of
the fits show that the compactness does not decrease at the same rate at
different depths.

##### Degradation time of the sunspot:

A decrease in compactness describes an increase in raggedness of the flux
tube. A higher raggedness causes an increase in the outer surface of the flux
tube. Thus, the contact area of the magnetic flux tube with the non-magnetic
surrounding plasma increases. This enables a stronger effect of the
surroundings, that is, plasma motions. We approximated the evolution of the
compactness by a linear fit (see Fig. 7) to calculate a degradation time of
the sunspot by the time at which the initial value of the fit has dropped by a
factor of $1/e$.

Above $z$ = -1 Mm the compactness is stable, as described and concluded above.
Therefore the following analysis focuses on layers deeper than 1 Mm below the
surface. The degradation time is calculated for the same depths as the
compactness. We also used the compactness obtained for both methods, that is,
a determination of the flux tube with $B_{\text{c}}$ and $B_{\text{eye}}$, to
calculate the degradation times. A similar distribution of the degradation
times with depth can be found for $B_{\text{c}}$ and $B_{\text{eye}}$ as shown
in Fig. 8 (red crosses and blue stars, respectively). Error bars show the
standard deviation obtained from the linear fits. Discrepancies between the
degradation times when applying different contour values for individual depths
can be ascribed to the evolution of the field strength within the flux tube.
For times $t\geq 20$ h at depths between 3 Mm and 5 Mm below the surface, the
contours of the flux tube obtained with the different methods are
significantly different. Fig. 9 shows a map of the magnetic field strength at
a depth of 3.75 Mm below the surface at $t$ = 20 h where the contour lines
obtained with the two methods differ significantly. This is one of the few
cases where the automatic and manual contours do not match: The higher value
of $B_{\text{c}}$ (red) excludes two regions that are included when
$B_{\text{eye}}$ (blue) is applied. This affects the determined compactness
and leads to the deviation of the degradation times for $B_{\text{c}}$ and
$B_{\text{eye}}$.

Figure 8: Degradation time of the sunspot at different depths from the
compactness of the flux tube. The analysis was made with different depth-
dependent contour values $B_{\text{c}}(z)$ and $B_{\text{eye}}(z)$ to
determine the sunspot boundary (see Sec. 2.2). Figure 9: Map of the magnetic
field strength at $z$ = -3.75 Mm, at $t$ = 20 h. The outlined contour of the
sunspot is determined with $B_{\text{c}}$ = 3500 G (red) and $B_{\text{eye}}$
= 2500 G (cyan).

The structure of the flux tube evolves fastest in regions between 6.25 Mm and
10 Mm below the surface. There, we measure the fastest degradation times of
less than three days, as Fig. 8 shows. The slowest decrease in compactness of
the flux tube is found at a depth of 15 Mm with a maximum degradation time of
5.4 days.

##  4 Stability of the sunspot

To understand the decay process of a sunspot, we took stabilising as well as
destabilising effects into account. To do this, we studied the stability of
the simulated sunspot in the convection zone.

###  4.1 Methods

The stability of a magnetic flux tube in the upper convection zone is affected
by two main effects: The stratification of the thermodynamic properties such
as density within the flux tube and in the direct surroundings, and the
curvature of the outer surface of the flux tube defining the sunspot. Meyer et
al. (1977) described the stability of a flux tube in the upper convection zone
and took both effects into account (see Sec. 1). We used Eq. (1) to study the
stability of the flux tube of the simulated sunspot in the convection zone. To
do this, we determined the following parameters: (1) The density difference in
the flux tube, $\rho_{\text{s}}$, and the surroundings, $\rho_{\text{qs}}$, as
a function of depth. (2) The inclination angle, $\chi$, of the outer surface
of the flux tube towards the vertical and its radius of curvature,
$R_{\text{c}}$. Both parameters are obtained from a surface function that
describes the radial extension of the sunspot with depth.

####  4.1.1 Density stratification in the upper convection zone

From observations we know that temperature, pressure, and density within a
sunspot are lower than in the surrounding quiet Sun. However, this only holds
for regions close to the surface. In deeper regions, these differences are
expected to decrease (e.g. Jahn & Schmidt, 1994; Schlichenmaier, 1997). We
determined spatial averages of the density for all defined depths $z$ and all
time steps for the areas covered by the flux tube of the sunspot and the quiet
Sun. The relative density difference between the two regions,
$(\rho_{\text{qs}}-\rho_{\text{s}})/\rho_{\text{qs}}$, decreases with
increasing depth, as Fig. 10 shows for time-averaged values. Because we
averaged the spot density over umbra and penumbra, the resulting density
contrast just beneath the surface is affected by the Wilson depression, which
leads to very low densities in the umbra. Together with a misinterpretation of
the flux tube contour, this leads to the discontinuities at $z>-0.625$ Mm.

Figure 10: Relative density difference between the quiet Sun and the magnetic
flux tube, time-averaged, at different depths.

####  4.1.2 Surface function of the flux tube

We define an equivalent radius of the flux tube as the radius, $r_{\rm
e}=\sqrt{A/\pi}$, of a circle with the same area (see Sec. 3.1). The
equivalent radius, $r_{\text{e}}(z,t)$, is calculated for all defined depths
$z$ and all 120 time steps. These individual values are shown in Fig. 11 as
rainbow-coloured dots. To approximate the depth dependence, we constructed an
asymptotic function,

$r_{\text{e}}(z,t)=\frac{a_{0}z^{2}+a_{1}}{a_{2}+a_{3}z^{3}}+a_{4}\,,$ (5)

with fitting parameters, $a_{i}$. The calculation was made for all 120 time
steps with $r_{\text{e}}(z,t)$ over a depth range from $z$ = -0.625 Mm to $z$
= -15 Mm. Fig. 11 shows the best-fit function (magenta line), obtained from
time-averaged equivalent radii (black stars). The obtained parameters for this
fit function are $a_{0}=-10.38$ Mm-1, $a_{1}=-2.74$ Mm, $a_{2}=-0.22$ Mm,
$a_{3}=1.85$ Mm-2, and $a_{4}=4.34$ Mm with $\chi^{2}_{\text{gof}}=0.002$ as
the goodness of the fit (gof).

Figure 11: Radial extension of a circular flux tube at different depths
obtained from the area $A$ of the simulated sunspot. The time-averaged values
(black stars) of all 120 time steps (rainbow-coloured dots, red to blue in
time) are fitted by an asymptotic function (magenta line).

Inclination angle $\chi$: The inclination angle $\chi$ is calculated from the
derivative of the surface function with

$\chi(z)=\tan^{\text{-1}}(r^{\prime}_{\text{e}}(z))=\tan^{\text{-1}}\left(\frac{z\cdot(2a_{0}a_{2}-a_{0}a_{3}z^{3}-3a_{1}a_{3}z)}{(a_{2}+a_{3}z^{3})^{2}}\right).$
(6)

At each depth position $z\leq$ -0.625 Mm and for each time step the spatial
derivative $r^{\prime}_{\text{e}}(z,t)$ is calculated. A minimum angle of less
than 1.7∘ is found at the lower-most depth, that is, 15 Mm below the surface.
The inclination angle slowly increases towards the surface and reaches 13∘ at
$z=-5$ Mm. Then, in the upper 5 Mm of the convection zone the increase
steepens and almost reaches 85∘ at $z$ = -0.625 Mm.

Radius of curvature $R_{\text{c}}$: To calculate the radius of curvature, the
first and second derivatives, $r^{\prime}_{\text{e}}(z,t)$ and
$r^{\prime\prime}_{\text{e}}(z,t)$, respectively, of the surface function have
to be calculated. It is then determined for each depth position $z\leq$ -0.625
Mm and time step individually with

$R_{\text{c}}(z)=\frac{(1+r^{\prime}_{\text{e}}(z)^{2})^{3/2}}{|r^{\prime\prime}_{\text{e}}(z)|}.$
(7)

The radius of curvature is largest at the lowest measurement point, that is,
15 Mm below the surface (see Fig. 12). The flux tube is almost vertical there,
thus, the curvature of the flux tube is smallest. The radius of curvature
decreases with decreasing depth and has a minimum at $z$ = -2.5 Mm. Upward of
the minimum, the radius of curvature increases again. Because of the
misinterpretation of the sunspot boundary at $z$ = -0.625 Mm (see Sec. 3.2),
the corresponding $R_{\text{c}}$ values show a large scatter, but the trend of
an upward-increasing curvature radius is obvious.

Figure 12: Radius of curvature at different depths, time-averaged (black
stars) over all 120 time steps (rainbow-coloured dots, blue to red in time).

Another parameter necessary to apply the stability criterion, that is, Eq.
(1), to the simulated sunspot is the magnetic field strength $B(z,t)$. It is
obtained as spatial average over the area of our flux tube contours for each
time step $t$ and each selected depth $z$.

###  4.2 Results

The stability of the sunspot at the different depth positions and for all time
steps was studied by evaluating Eq. (1) with $g$ = 274 m s-1. When the
difference on the left side of Eq. (1) is positive, then the flux tube is
stable. The difference obtained for the different depths and time steps is
visualised in Fig. 13 as dots.

Down to $z$ = -1.25 Mm the difference is positive (see Fig. 13), and this
means that the flux tube is stable. At $z$ = -1.875 Mm the time-averaged value
is still positive, although individual time steps are negative, indicating an
instability of the flux tube at these time steps. At greater depth, for $z$ ¡
-2.5 Mm, the flux tube is found to be unstable for all time steps. The change
from stable to unstable is almost always located at the same depth as the
smallest radius of curvature, $R_{\text{c}}$, of the sunspot flux tube. At
this depth the surface function has an inclination angle, $\chi$, between 42∘
and 62∘ with a temporal average of 46∘. The difference reaches a minimum at
$z$ = -3.125 Mm. There, the flux tube is most unstable. With increasing depth,
the negative difference approaches zero. However, the flux tube does not reach
a stable configuration down to the lowest measurement point.

Figure 13: Stability criterion of the sunspot flux tube at different depths
time-averaged (black stars) over all 120 time steps (rainbow-coloured dots,
blue to red in time). Positive values refer to a stable magnetic flux tube,
and negative values represent an unstable configuration.

Fig. 13 shows that the coloured dots for each depth do not follow the rainbow
pattern. This means that no time-dependent relation of the flux-tube stability
can be found.

##  5 Discussion

We applied the stability criterion that was derived with a simple magneto-
static model by Meyer et al. (1977) in a current numerical simulation. We find
that our simulated sunspot is stable in the upper 2 Mm below the photosphere.
Our analysis of the compactness of the flux tube shows that the compactness of
the sunspot is constant at the surface. For layers deeper than $z$=-1.25 Mm,
the compactness decreases in time. Meyer et al. (1977) have suspected that a
sunspot flux tube could be stable up to a depth of 2 Mm, although their
magneto-static model only extended down to 0.650 Mm below the photosphere. In
these layers, the inclined magnetic field is predominately kept together by
buoyancy forces. The density difference between the magnetic flux tube and the
surroundings stabilises the sunspot.

In deeper layers, the stability criterion yields an unstable configuration, as
shown in Fig. 13. This represents the increasing effect of fluting. The flux
tube is least stable at around $z$ = -3 Mm, whereas the fastest decrease in
compactness is found in deeper layers (see Fig. 8). We conjecture that the
interchange process is damped within this region. Above $z$ = -1.5 Mm, the
field lines are held together by buoyancy to a stable configuration. Magnetic
tension affects the behaviour of the field lines in deeper layers. These non-
local effects are not taken into account in the linear stability analysis. The
field lines can move more freely at larger depths. Therefore the decrease in
compactness is faster in layers below $z$ = -3 Mm. It proceeds most rapidly
between $z$ = -5 Mm and $z$ = -10 Mm (see Fig. 8).

In even deeper layers, the decrease in compactness is found to proceed more
slowly, and the instability of the flux tube approaches a marginally stable
state (see Fig. 13). This behaviour is caused by a combination of two effects:
(1) The convective timescale increases with depth, as described in Sec. 2.1
and shown in Fig. 2. This leads to the slower action of destabilising forces.
(2) The flux tube is more vertical and less strongly curved, which makes it
less susceptible to fluting, an effect that has been predicted by Meyer et al.
(1977).

The existence of the fluting stability alone is not sufficient for sunspot
decay. In addition, turbulent flows are required that remove flux from the
spot. Rempel (2015) found a stabilising effect of the moat cell that developed
around the simulated sunspot. Below the penumbra, the spot is surrounded by an
annulus of upflowing plasma that has a strongly reduced convective RMS
velocity and is almost devoid of downflows. This strongly inhibits turbulent
transport and in particular prevents the submergence of field lines below the
penumbra, which would be required for a decay of the penumbra.

A non-monolithic sunspot flux tube was proposed by Parker (1979). He described
the individual field lines of a sunspot that are to be kept together in the
photosphere. The observable monolithic structure would be divided into several
smaller flux bundles approximately 1 Mm below the visible photosphere.

The formation of intrusions at the sunspot boundary, which penetrate the flux
tube below the visible surface, was also found by Panja et al. (2021). They
calculated MHD simulations with different initial values for the vertical
magnetic field at the bottom boundary of their simulation domain. After about
$10$ hours, the sunspots simulated by Panja et al. (2021) showed a
fragmentation of the subsurface field that was substantially stronger than
what we found here. This might be due to their more shallow domain (about 5.5
Mm), which cuts out the deeper part of the spot that is less susceptible to
fluting. Despite this difference, they also find that the intrusions evolve
fastest in depths below $z$ = -5 Mm. Another critical difference is that their
setup did not impose a strongly inclined magnetic field at the top boundary.
They found that the fluting instability promotes penumbra formation, which is
not the case in our simulation. Switching to a potential field boundary
condition leads to a naked spot without penumbra in about half an hour, as
shown in Rempel (2015). The fluting instability does lead to a growth of the
penumbra over time at the expense of the umbra. This evolution was found in
Rempel (2015) and can be seen in intensity maps in the top panels of Fig. 14.
This points in the same direction as the findings of Panja et al. (2021): The
intrusions cause a conversion of magnetic flux from the umbra to the penumbra.

The observational study of the decay of a sunspot by Benko et al. (2018)
showed that the decrease in magnetic field strength within the umbra is one of
the first observable effects of sunspot decay. In time, regions of weaker
field are occupied by the penumbra or granulation. Benko et al. (2018)
concluded that the umbral area decreases during the decay phase of the
sunspot. In our simulation the intrusions of field free plasma that are most
prominent at a depth of a few Mm do lead to a visible imprint at the surface
(see Fig. 14). While we do not see fully developed light bridges, we find
penumbral intrusions and clusters of higher umbral dot density associated with
the subsurface fluting. This suggests that subsurface decay does not remain
hidden from the surface.

We conjecture that the degradation of the sunspot at its outer surface will in
time lead to the separation of the large magnetic flux tube into smaller ones.
Light bridges that are observed as the onset of the decay could be caused by
such intrusions. The weak magnetic field found within the sunspot in regions
above $z$ = -5 Mm, at $t$ ¿ 20 h is the first indication of such a scenario.
This agrees with the observations reported by McIntosh (1981) that the process
that leads to sunspot decay should already act before the sunspot has fully
developed. An interesting property of the simulated sunspot is that the total
flux remains constant in spite of the subsurface decay, which can be
attributed to the stabilising buoyancy of the penumbra and the stabilising
effect of the moat cell. We conjecture that this confinement will eventually
fail and at that time lead to a rapid disappearance of the already heavily
corrugated spot.

Figure 14: Comparison of the sunspot at the first ($t$ = 0h, left) and last
($t$ = 29.75h, right) time step in intensity maps at $\tau$ =1 (top row) and
the magnetic field strength at a depth of 1.25 Mm (bottom row). The outlined
contours are determined with $I_{\mathrm{c}}=0.9I_{\mathrm{n}}$ and
$B_{\mathrm{c}}=1669$ G.

##  6 Conclusion

We report the following main results from our study: (1) The convective
timescale (see Fig. 2), increases with depth and reaches a value of 8 hours at
$z=-15$ Mm, (2) the compactness of the sunspot decreases fastest at a depth
around $z=-8$ Mm (see Fig. 8), and (3) our stability analysis reveals that the
sunspot is stable down to approximately $-2$ Mm and most unstable around
$z=-3$ Mm (see Fig. 13). This leads to a consistent picture of sunspot
stability and decay.

Our analysis shows that sunspot (in)stability and decay is an interplay
between stabilising and destabilising forces as well as timescales. Buoyant
inclined field lines just beneath the photosphere stabilise the sunspot and
keep the field lines together. Non-local magnetic tension effects damp the
fluting instability in the adjacent layers. In even deeper layers, the fluting
instability steadily decreases the compactness of the flux tube. Thus, the
spot eventually loses its coherency in layers between $z$ = -5 Mm and $z$ =
-10 Mm. Deeper layers are less strongly affected because the destabilising
force is weaker, convective timescales (time of change) increase substantially
with depth, and the moat cell is stabilising.

We recall that this scenario of sunspot decay is based on simulated data. The
described process occurs within the upper convection zone, a region hidden
from real observations. For a more detailed clarification of the decay
process, we might therefore have to compare characteristics of observations
with those seen in the photospheric region of the simulated sunspot. Tracers
of the degradation process that are visible in the photosphere could then be
studied in greater detail in real observations. We propose to take light
bridges as tracers for the indentations at the boundary of the flux tube below
the solar surface into account. Panja et al. (2021) reported that upflows
within the intrusions become visible as light bridges at the surface. The role
of MMFs in the sunspot decay process has been discussed in many studies (see
Sec. 1). Their connection to the sunspot should be studied in simulations. In
simulations we can take advantage of tracing the connection of their field
lines below the surface. This might give insight into a possible role of these
features in the interchange and degradation process of a sunspot.

###### Acknowledgements.

We gratefully acknowledge fruitful discussions with Nazaret Bello González and
Markus Schmassmann, and we thank Rebecca Centeno Elliott for carefully reading
and commenting on the manuscript. We appreciate the useful comments and
suggestions of the anonymous referee. HS has been funded by the Deutsche
Forschungsgemeinschaft, under grant No. RE 3282 and acknowledges financial
support from the HAO visitor program. This material is based upon work
supported by the National Center for Atmospheric Research, which is a major
facility sponsored by the National Science Foundation under Cooperative
Agreement No. 1852977.

## References

  * Benko et al. (2018) Benko, M., González Manrique, S. J., Balthasar, H., et al. 2018, A&A, 620, A191
  * Bumba (1963) Bumba, V. 1963, Bull. astr. Inst. Czechosl., 14, 91
  * Harvey & Harvey (1973) Harvey, K. & Harvey, J. 1973, Sol. Phys., 28, 61
  * Jahn (1992) Jahn, K. 1992, in Sunspots: Theory and Observations, ed. J. H. Thomas & N. O. Weiss (Dordrecht: Springer Netherlands), 139–162
  * Jahn (1997) Jahn, K. 1997, in PASPC, Vol. 118, 1st Advances in Solar Physics Euroconference. Advances in Physics of Sunspots, ed. B. Schmieder, J. C. del Toro Iniesta, & M. Vazquez, 122
  * Jahn & Schmidt (1994) Jahn, K. & Schmidt, H. U. 1994, A&A, 290, 295
  * Kubo (2012) Kubo, M. 2012, in PASPC, Vol. 455, 4th Hinode Science Meeting: Unsolved Problems and Recent Insights, ed. L. Bellot Rubio, F. Reale, & M. Carlsson, 49
  * Kubo et al. (2007) Kubo, M., Shimizu, T., & Tsuneta, S. 2007, ApJ, 659, 812
  * Kurucz (1993) Kurucz, R. 1993, ATLAS9 Stellar Atmosphere Programs and 2 km/s grid. Kurucz CD-ROM No. 13. Cambridge, 13
  * Martínez Pillet (2002) Martínez Pillet, V. 2002, Astron. Nachr., 323, 342
  * Martínez Pillet et al. (1993) Martínez Pillet, V., Moreno-Insertis, F., & Vazquez, M. 1993, A&A, 274, 521
  * McIntosh (1981) McIntosh, P. S. 1981, in The Physics of Sunspots, ed. L. E. Cram & J. H. Thomas, 7–54
  * Meyer et al. (1977) Meyer, F., Schmidt, H. U., & Weiss, N. O. 1977, MNRAS, 179, 741
  * Meyer et al. (1974) Meyer, F., Schmidt, H. U., Weiss, N. O., & Wilson, P. R. 1974, MNRAS, 169, 35
  * Panja et al. (2021) Panja, M., Cameron, R. H., & Solanki, S. K. 2021, ApJ, 907, 102
  * Parker (1979) Parker, E. N. 1979, ApJ, 230, 905
  * Petrovay & Moreno-Insertis (1997) Petrovay, K. & Moreno-Insertis, F. 1997, ApJ, 485, 398
  * Petrovay & van Driel-Gesztelyi (1997) Petrovay, K. & van Driel-Gesztelyi, L. 1997, Sol. Phys., 176, 249
  * Priest (2014) Priest, E. 2014, Magnetohydrodynamics of the Sun (Cambridge University Press)
  * Rempel (2012) Rempel, M. 2012, ApJ, 750, 62
  * Rempel (2015) Rempel, M. 2015, ApJ, 814, 125
  * Schlichenmaier (1997) Schlichenmaier, R. 1997, PhD thesis, Ludwig-Maximilians-Universität München
  * Schlüter & Temesváry (1958) Schlüter, A. & Temesváry, S. 1958, Symposium - IAU, 6, 263
  * Schüssler (1984) Schüssler, M. 1984, A&A, 140, 453
  * Sheeley (1969) Sheeley, N. R. 1969, Sol. Phys., 9, 347
  * Shine & Title (2000) Shine, R. & Title, A. 2000, in Encyclopedia of Astronomy and Astrophysics, ed. P. Murdin (Bristol: Institute of Physics Publishing), 2038
  * Simon & Leighton (1964) Simon, G. W. & Leighton, R. B. 1964, ApJ, 140, 1120
  * Simon & Weiss (1970) Simon, G. W. & Weiss, N. O. 1970, Sol. Phys., 13, 85
  * Solanki (2003) Solanki, S. K. 2003, A&A Rev., 11, 153
  * Stix (2002) Stix, M. 2002, The Sun, A&A Library (Berlin Heidelberg: Springer)
  * Strecker (2020) Strecker, H. 2020, PhD thesis, Albert-Ludwigs-Universität Freiburg
  * Thomas et al. (2002) Thomas, J. H., Weiss, N. O., Tobias, S. M., & Brummell, N. H. 2002, Nature, 420, 390
  * van Driel-Gesztelyi & Green (2015) van Driel-Gesztelyi, L. & Green, L. M. 2015, Liv. Rev. Sol. Phys., 12, 1
  * Verma et al. (2012) Verma, M., Balthasar, H., Deng, N., et al. 2012, A&A, 538, A109
  * Zhang et al. (1992) Zhang, H., AI, G., Wang, H., Zirin, H., & Patterson, A. 1992, Sol. Phys., 140, 307

##  Appendix A Determination of the flux tube boundary

Table 1: Threshold levels $B_{c}$ and $B_{\text{eye}}$ for determining the
boundary of the sunspot in maps of the magnetic field strength for defined
depths $z$.

Depth [Mm] | 0 | -0.625 | -1.25 | -1.875 | -2.5 | -3.125 | -3.75 | -4.375 | -5
---|---|---|---|---|---|---|---|---|---
$B_{c}$ [G] | 600 | 977 | 1669 | 2203 | 2711 | 3139 | 3500 | 3806 | 4130
$B_{\text{eye}}$ [G] | 1000 |  | 1700 |  | 2000 |  | 2500 |  | 3500
Depth [Mm] | -5.625 | -6.25 | -6.875 | -7.5 | -8.125 | -8.75 | -9.375 | -10 | -10.625
$B_{c}$ [G] | 4358 | 4547 | 4731 | 4782 | 4832 | 4928 | 5021 | 5036 | 5055
$B_{\text{eye}}$ [G] |  | 4000 |  | 4250 |  | 4500 |  | 4500 |
Depth [Mm] | -11.25 | -11.875 | -12.5 | -13.125 | -13.75 | -14.375 | -15 | -16.25 | -17.5
$B_{c}$ [G] | 5104 | 5156 | 5084 | 5014 | 5017 | 5086 | 5100 | 5309 | 5509
$B_{\text{eye}}$ [G] | 4500 |  | 4500 |  | 4500 |  | 4500 |  |

##  Appendix B Compactness of the flux tube

Figure 15: Compactness of the sunspot (normalised area-perimeter ratio) in
time, shown at different depths. The boundary of the sunspot is determined
with depth-dependent contour values $B_{c}(z)$ (red crosses) and
$B_{\text{eye}(z)}$ (blue stars). Same as Fig. 7. Figure 16: Same as Fig. 15,
but for different depths.


# Backup Control Barrier Functions: Formulation and Comparative Study

Yuxiao Chen1, Mrdjan Jankovic2, Mario Santillo2, and Aaron D. Ames1 1 Yuxiao
Chen and Aaron D. Ames are with Department of Mechanical and Civil
Engineering, California Institute of Technology, Pasadena, CA, USA
<EMAIL_ADDRESS>Mario Santillo and Mrdjan Jankovic are with Ford
Research and Advanced Engineering, Dearborn, MI, USA
<EMAIL_ADDRESS>
###### Abstract

The backup control barrier function (CBF) was recently proposed as a tractable
formulation that guarantees the feasibility of the CBF quadratic programming
(QP) via an implicitly defined control invariant set. The control invariant
set is based on a fixed backup policy and evaluated online by forward
integrating the dynamics under the backup policy. This paper is intended as a
tutorial of the backup CBF approach and a comparative study to some
benchmarks. First, the backup CBF approach is presented step by step with the
underlying math explained in detail. Second, we prove that the backup CBF
always has a relative degree 1 under mild assumptions. Third, the backup CBF
approach is compared with benchmarks such as Hamilton Jacobi PDE and Sum-of-
Squares on the computation of control invariant sets, which shows that one can
obtain a control invariant set close to the maximum control invariant set
under a good backup policy for many practical problems.

##  I Introduction

Control barrier functions (CBF) [1, 2] were proposed as a method that enforces
constraints on dynamic systems, which typically works as a supervisory
controller on top of a legacy controller. To guarantee the satisfaction of the
safety constraints, a CBF quadratic program (QP) is solved online. While CBF
is getting increasingly popular due to its simple implementation and strong
guarantee, the construction of CBF is sometimes overlooked.

To guarantee that the CBF QP is always feasible, a control invariant set is
needed, which is defined as a set in which any trajectory of the dynamic
system can stay indefinitely. The concept of control invariant sets has been
studied under various background and names, such as viability kernel [3],
infinite time reachable set [4], and various methods have been proposed to
compute the control invariant set depending on the system dynamics, see [5]
for an overview. Unfortunately, the computation of control invariant sets is
notoriously difficult. Even for simple cases such as linear or polynomial
dynamics, computation tools such as Minkowski operations [6], robust linear
program [7], and Sum-of-Squares [8] do not scale well. For general nonlinear
dynamic systems, the standard tool for computing invariant sets is Hamilton
Jacobi PDE [9], which typically cannot scale beyond systems with state
dimensions 4 due to the exponential complexity.

Due to the difficulty of synthesizing proper control barrier functions based
on control invariant sets, CBF QP has been implemented without control
invariant sets. One simple treatment is to assume infinite actuation power
[10], in which case the feasibility of the CBF QP can be guaranteed when a
simple sign condition is satisfied. To be specific, if one can show that when
the Lie derivative of the CBF w.r.t. the input dynamics is zero, the CBF
condition is satisfied by the intrinsic dynamics, the CBF QP is always
feasible under infinite actuation. Another simplification is to assume that
the velocity instead of the acceleration is under control [11], which can
guarantee the feasibility of the CBF QP for constraints on the position since
the CBF QP can always pick a velocity pointing away from the constraint. In
essence, assuming direct velocity control is similar to assuming infinite
actuation as an instantaneous change of velocity requires infinite force.
Obviously, these assumptions are not true in practice, and the safety of the
system is subject to parameter tuning. The simplified approaches might work
for low-speed cases as the change of velocity is not severe, but will not work
in general for highly dynamic applications.

Another issue of CBF without a control invariant set is the relative degree.
For example, a typical vehicle/robot model with acceleration input is a second
(or higher) order model with the position states and velocity states. If one
directly uses the constraint on the position as a control barrier function,
the relative degree of the CBF is 2, and its Lie derivative does not contain
the acceleration input. Several solutions have been proposed for the high
relative degree, such as Input-Output Linearization [12, 13] and backstepping
type formulations [14]. Again, the CBF performance is subject to parameter
tuning.

To the best of our knowledge, under limited input, there is not a generic
method that guarantees the feasibility of the CBF QP without a control
invariant set. All of the above-mentioned methods rely on heuristics and
tuning to work in practice. Even when the CBF works well in the test cases,
there is no guarantee of performance for cases not included in the test.

The idea of control barrier functions based on backup controllers (referred to
as backup CBF for the remainder of this paper) [15] was proposed based on a
simple observation that extends a small control invariant set, which is
typically easier to obtain, to a larger control invariant set by fixing a
backup controller. The backup CBF guarantees the feasibility of the CBF QP and
circumvents the difficult computation of a control invariant set by implicitly
representing the control invariant set. However, the implicit representation
calls for online integration of the dynamics, and the resulting control
invariant set is typically suboptimal in the sense that it is not the maximum
control invariant set. This paper is intended as a tutorial to the backup CBF
approach by demonstrating the method step by step on some simple examples. We
also provide explanations and visualizations on the advantage and
disadvantages of the backup CBF over some benchmark formulations.

##  II Preliminaries and a motivating example

###  II-A Control barrier functions

We begin by a brief review of control barrier functions. Consider the
following control affine system:

$\dot{x}=f(x)+g(x)u,\quad
x\in\mathbb{R}^{n},u\in\mathcal{U}\subseteq\mathbb{R}^{m},$ (1)

where $x$ is the system state and $u$ is the input. Suppose there exists a
function $h:\mathbb{R}^{n}\to\mathbb{R}$ that satisfies

$\displaystyle\forall~{}x\in{\mathcal{X}_{0}},$ $\displaystyle h(x)\geq 0$ (2)
$\displaystyle\forall~{}x\in{\mathcal{X}_{d}},$ $\displaystyle h(x)<0$
$\displaystyle\forall~{}x\in\left\\{{x\mid h(x)\geq 0}\right\\},$
$\displaystyle\exists~{}u\in\mathcal{U}\;\mathrm{s.t.}~{}\dot{h}+\alpha\left(h\right)\geq
0,$

where $\mathcal{X}_{0}$ is the set of initial states and $\mathcal{X}_{d}$ is
the danger set that we want to keep the state away from. $\alpha(\cdot)$ is a
class-$\mathcal{K}$ function, i.e., $\alpha(\cdot)$ is strictly increasing and
satisfies $\alpha(0)=0$. Then $h$ is called control barrier function, and for
any legacy controller, the CBF controller is a supervisory controller that
enforces the state to stay inside $\left\\{{x\mid h(x)\geq 0}\right\\}$ with
the following quadratic programming:

$\displaystyle u^{\star}=$
$\displaystyle\mathop{\arg\min}\limits_{u\in\mathcal{U}}\left\|{u-{u^{0}}}\right\|^{2}$
(3) $\displaystyle\mathrm{s.t.}{}$ $\displaystyle\nabla h\cdot
f\left({x,u}\right)+\alpha\left(h\right)\geq 0,$

where $u^{0}$ is the input of the legacy controller.

The third condition in (2) ensures that (3) is always feasible, yet it is
difficult to find an $h$ that satisfies it. For clarity, we refer to a
function that satisfies the first two conditions in (2) a CBF candidate, a
function that satisfies all three conditions a valid CBF.

The CBF condition is closely tied to the concept of a control invariant set,
which is defined as follows.

######  Definition 1.

A set $\mathcal{S}$ is a control invariant set if there exists a control law
$\pi:\mathbb{R}^{n}\to\mathcal{U}$ such that for all initial condition
$x(0)\in\mathcal{S}$, $\forall t\geq 0,x(t)\in\mathcal{S}$.

It is straightforward to see that for a valid CBF, $\\{x|h(x)\geq 0\\}$ is a
control invariant set. On the other direction, suppose for a CBF candidate
$h$, its 0-level set is a control invariant set. Immediately from Definition
1,

$h(x)=0\to\exists u\in\mathcal{U}~{}\mathrm{s.t.}\dot{h}(x,u)\geq 0,$

otherwise the state will exit $\mathcal{S}$, which contradicts Definition 1.
Then $h$ can be shown to be a valid control barrier function by picking an
$\alpha$ large enough given some continuity condition (Lipschitz continuity of
$h$ and $\dot{h}$). As a result, the third condition in (2) is also referred
to as the set invariance condition.

###  II-B Double integrator example

As a motivating example, consider a simple double integrator with

$\dot{x}=\begin{bmatrix}\dot{s}\\\ \dot{v}\end{bmatrix}=\begin{bmatrix}v\\\
u\end{bmatrix},\quad u\in[-u_{\max},u_{\max}],v\in[-v_{\max},v_{\max}].$

The safety constraint is $X\leq C$ with a constant $C$. Suppose one directly
takes the safety constraint as a CBF candidate: $h_{0}(x)=C-s$. There are two
issues. First, $h_{0}$ has relative degree 2, i.e., $\dot{h}_{0}$ is not a
function of $u$. Second, the set $\\{x|h_{0}(x)\geq 0\\}$ is not control
invariant. Consider the case where $s=C,v>0$, the limited input cannot stop
the state from crossing into the danger set due to inertia.

For this simple system, it is widely known that a simple valid CBF exists:

$h(x)=C-s-\mathds{1}_{v>0}\frac{v^{2}}{2u_{\max}},$ (4)

which is the safety constraint combined with the minimum stopping distance.
When $h(x)\geq 0$, one can always apply the maximum deceleration $u=-u_{\max}$
until $v=0$ and the safety constraint remains satisfied. Fig. 1 shows the
difference between $h$ and $h_{0}$ where $h_{0}$ includes a part of the state
space from which the safety constraint will eventually be violated.

Figure 1: Difference between $h$ and $h_{0}$ for the double integrator example

In fact, the above CBF is exactly the backup CBF with the backup policy being

$\pi(x)=-\mathds{1}_{v>0}u_{\max}.$

We shall show later that the backup policy $\pi$ takes the system to a small
control invariant set $\mathcal{S}_{0}=\\{x|v=0\\}$, and $\\{x|h(x)\geq 0\\}$
with $h$ in (4) is the set of initial conditions that can be brought to
$\mathcal{S}_{0}$ while satisfying the safety constraint. In this view, $h$ is
induced from the backup policy $\pi$.

As mentioned in the Introduction, computing a control invariant set is very
difficult, especially for high-dimensional nonlinear systems. We shall show
that for complicated dynamic systems, a valid CBF can be fairly easily
obtained from a backup policy. Though the CBF induced from a backup policy
might not have the nice closed form as in the double integrator case, the CBF
QP can be solved online and is guaranteed to be feasible.

###  II-C Notation and preliminaries

Before going into the detail, some important notations are reviewed. Given the
dynamic system in (1), a control policy $\pi:\mathbb{R}^{n}\to\mathcal{U}$,
the closed-loop dynamics under $\pi$ is $\dot{x}=f_{\pi}(x)=f(x)+g(x)\pi(x)$.
We let
$\operatorname{\Phi_{f_{\pi}}}:\mathbb{R}^{n}\times(-\infty,\infty)\to\mathbb{R}^{n}$
denote the flow map, i.e., $\operatorname{\Phi_{f_{\pi}}}(x_{0},t)$ denotes
the solution $x(t_{0}+t)$ to the Initial Value Problem (IVP) at time $t_{0}+t$
with $x(t_{0})=x_{0}$ under the dynamics $f_{\pi}$. Note that since (1) is
time invariant, the initial time $t_{0}$ is irrelevant. Moreover, the flow map
is additive in $t$:

$\operatorname{\Phi_{f_{\pi}}}(\operatorname{\Phi_{f_{\pi}}}(x,t_{1}),t_{2})=\operatorname{\Phi_{f_{\pi}}}(x,t_{1}+t_{2}).$
(5)

Since the flow map can be written as

$\operatorname{\Phi_{f_{\pi}}}(x_{0},t)=x_{0}+\int_{0}^{t}f_{\pi}(x(\tau))d\tau,$
(6)

it is differentiable w.r.t. both $t$ and $x$. By inspection,
$\frac{\partial\operatorname{\Phi_{f_{\pi}}}(x_{0},t)}{\partial
t}=f_{\pi}(\operatorname{\Phi_{f_{\pi}}}(x_{0},t))$. To obtain the partial
derivative over $x_{0}$, take derivative on both sides:

$\frac{\partial\operatorname{\Phi_{f_{\pi}}}(x_{0},t)}{\partial
x_{0}}=\operatorname{\mathbf{I}}+\int_{0}^{t}\frac{df_{\pi}}{dx}\frac{\partial
x(\tau)}{\partial x_{0}}d\tau,$

from [16], we can define
$Q(t)=\frac{\partial\operatorname{\Phi_{f_{\pi}}}(x_{0},t)}{\partial x_{0}}$
as the sensitivity Jacobian, and taking derivative over $t$ on both sides
yields

$\dot{Q}(t)=\frac{df_{\pi}(x(t))}{dx(t)}Q(t),$ (7)

where $x(t)=\operatorname{\Phi_{f_{\pi}}}(x_{0},t)$. This fact will play a key
role on the derivation of the CBF condition.

######  Proposition 1.

For all $\tau\geq t$,

$\frac{d\Phi(x(t),\tau-t)}{dt}|_{u(t)=\pi(t)}=0.$ (8)

###### Proof.

Note that the flow map can be written as

$\operatorname{\Phi_{f_{\pi}}}(x(t),\tau-t)=x(t)+\int_{t}^{\tau}f_{\pi}(x(\overline{\tau}))d\overline{\tau},$

taking the derivative over $t$ gives

$\frac{d\Phi(x(t),\tau-t)}{dt}=\dot{x}(t)-f_{\pi}(x(t)),$

and the result follows as $\dot{x}|_{u=\pi(x)}=f_{\pi}(x)$.

∎

The physical interpretation of Proposition 1 is that the flow under the backup
strategy $\pi$ would not change with time if the current state follows $\pi$.

##  III Backup CBF

This section presents the backup CBF step by step.

###  III-A Definition and CBF QP

Consider a safe region $\mathcal{C}\doteq\\{x|h^{\mathcal{C}}\geq
0\\}\subseteq\mathbb{R}^{n}$, a control invariant set
$\mathcal{S}_{0}\doteq\\{x|h^{\mathcal{S}}\geq 0\\}\subseteq\mathcal{C}$ with
$h^{\mathcal{C}},h^{\mathcal{S}}:\mathbb{R}^{n}\to\mathbb{R}$ differentiable.
Given a backup controller $\pi$, define the $T$-time constrained reachable set
as

$\displaystyle\mathcal{S}$
$\displaystyle=\mathcal{R}(\mathcal{S}_{0},\mathcal{C},f_{\pi},T)$ (9)
$\displaystyle\doteq\\{x\in\mathbb{R}^{n}|\operatorname{\Phi_{f_{\pi}}}(x,T)\in\mathcal{S}_{0}\wedge\forall
t\in[0,T],\operatorname{\Phi_{f_{\pi}}}(x,t)\in\mathcal{C}\\},$

which contains the states from which the flow under $f_{\pi}$ will stay in
$\mathcal{C}$ and reach $\mathcal{S}_{0}$ by $T$. The theoretical root of the
backup CBF approach is the following simple observation.

######  Theorem 1.

[17] Given an initial control invariant set $\mathcal{S}_{0}$, a safe region
$\mathcal{C}\subseteq\mathbb{R}^{n}$, and a backup controller $\pi$, for all
$T>0$, the $T$-time constrained reachable set $\mathcal{S}$ defined in (9) is
a control invariant set and
$\mathcal{S}_{0}\subseteq\mathcal{S}\subseteq\mathcal{C}$.

The proof can be found in [17], we provide a brief proof for completeness.

###### Proof.

By definition, $\mathcal{S}_{0}\subseteq\mathcal{S}$, and
$\mathcal{S}\in\mathcal{C}$ since
$\operatorname{\Phi_{f_{\pi}}}(x,0)\in\mathcal{C}$. By Definition 1, there
exists $\pi_{0}$ that keeps $\mathcal{S}_{0}$ invariant. Define
$\pi^{\prime}(x)=\left\\{{\begin{array}[]{*{20}{c}}{{\pi_{0}}(x),x\in{\mathcal{S}_{0}}}\\\
{\pi(x),x\notin{\mathcal{S}_{0}}}\end{array}}\right.$, it can be verified that
$\pi^{\prime}$ keeps $\mathcal{S}$ invariant. ∎

Figure 2: Enlarging a known control invariant set $\mathcal{S}_{0}$ with the
backup policy $\pi$

Fig. 2 shows theorem 1 pictorially. The yellow square is the state constraint
$\mathcal{C}$, given a small known control invariant set $\mathcal{S}_{0}$
(the blue ellipse), all states that can be driven to $\mathcal{S}_{0}$ under
$f_{\pi}$ (the closed-loop dynamics under the backup policy $\pi$) while
staying inside $\mathcal{C}$ in the meantime forms a larger control invariant
set $\mathcal{S}$ shown in brown. Note that the point on the top right is not
inside $\mathcal{S}$ because the trajectory under $f_{\pi}$ is not completely
contained $\mathcal{C}$.

######  Assumption 1.

For simplicity, we assume that the backup controller $\pi$ keeps
$\mathcal{S}_{0}$ invariant, and the class-$\mathcal{K}$ function $\alpha$ is
selected so that $\forall x\in\mathcal{S}_{0},\nabla
h^{\mathcal{S}}f_{\pi}(x)+\alpha(h(x))\geq 0$.

Regarding the horizon $T$, the following lemma is true.

######  Lemma 1.

For all $T_{1},T_{2}>0$, let
$\mathcal{S}=\mathcal{R}(\mathcal{S}_{0},\mathcal{C},f_{\pi},T_{1})$, then
$\mathcal{R}(\mathcal{S},\mathcal{C},f_{\pi},T_{2})=\mathcal{R}(\mathcal{S}_{0},\mathcal{C},f_{\pi},T_{1}+T_{2})$.

###### Proof.

By definition, $\forall
x\in\mathcal{R}(\mathcal{S},\mathcal{C},f_{\pi},T_{2}),$
$\operatorname{\Phi_{f_{\pi}}}(x,T_{2})\in\mathcal{S}=\mathcal{R}(\mathcal{S}_{0},\mathcal{C},f_{\pi},T_{2})$.
By (5),
$\operatorname{\Phi_{f_{\pi}}}(x,T_{1}+T_{2})=\operatorname{\Phi_{f_{\pi}}}(\operatorname{\Phi_{f_{\pi}}}(x,T_{2}),T_{1})\in\mathcal{S}_{0}$.
And by (9), the trajectory within $[0,T_{1}+T_{2}]$ is contained in
$\mathcal{C}$, and the conclusion follows. ∎

Lemma 1 states that $\mathcal{R}(\mathcal{S}_{0},\mathcal{C},f_{\pi},T)$
monotonically increases with $T$ in the set inclusion sense.

With a control invariant set $\mathcal{S}$ defined, a control barrier function
$h$ can be defined.

######  Lemma 2.

$\mathcal{S}$ is the 0-level set of the following function

$h(x)=\min\\{\mathop{\min}\limits_{t^{\prime}\in[0,T]}h^{\mathcal{C}}(\Phi_{f_{\pi}}(x,t^{\prime})),h^{\mathcal{S}}(\Phi_{f_{\pi}}(x,T))\\}.$
(10)

###### Proof.

First notice that by the continuity of the flow function $\Phi_{f_{\pi}}$ and
the $\min$ function, $h$ is continuous. For all $x\in\mathcal{S}$, by
definition, under the backup strategy $\pi$, the state evolution
$\Phi_{f_{\pi}}(x,t^{\prime})$ would satisfy the constraint and reach
$\mathcal{S}_{0}$ at time $T$, therefore $h(x)\geq 0$. On the other hand, for
all $x\notin\mathcal{S}$, under the backup strategy $\pi$, the state evolution
either violates the state constraint at some $t^{\prime}$, i.e., $\exists
t^{\prime}\in[0,T],h^{\mathcal{C}}(\Phi_{f_{\pi}}(x,t^{\prime}))<0$, or does
not reach $\mathcal{S}_{0}$ within the horizon $T$, i.e.,
$h^{\mathcal{S}}(\Phi_{f_{\pi}}(x,T))<0$, indicating that $h(x)<0$. Therefore,
$\mathcal{S}=\\{x|h(x)\geq 0\\}$. ∎

To implement the CBF $h$, $\dot{h}$ is needed. First $h$ needs to be written
as a function of time:

$h(t)=\min\\{\mathop{\min}\limits_{\tau\in[t,t+T]}h^{\mathcal{C}}(\Phi_{f_{\pi}}(x(t),\tau-t)),h^{\mathcal{S}}(\Phi_{f_{\pi}}(x(t),T))\\}.$
(11)

Since $h$ is defined as the minimum of multiple functions, it may not be
differentiable. Even in the differentiable case, the CBF condition may turn
out to be nonconvex in the control input $u$. Instead, the CBF condition is
enforced on every $\tau$, which is a sufficient condition for original CBF
condition.

First consider $h^{\mathcal{S}}(\operatorname{\Phi_{f_{\pi}}}(x(t),T))$. By
the chain rule,

$\frac{dh^{\mathcal{S}}(\operatorname{\Phi_{f_{\pi}}}(x(t),T))}{dt}=\frac{dh^{\mathcal{S}}}{dx}\frac{\partial\operatorname{\Phi_{f_{\pi}}}(x(t),T)}{\partial
x}(f(x(t))+g(x(t))u(t)),$

where the sensitivity Jacobian
$\frac{\partial\operatorname{\Phi_{f_{\pi}}}(x(t),T)}{\partial x}$ can be
calculated with (7). The computation of
$h^{\mathcal{C}}(\Phi_{f_{\pi}}(x(t)),\tau-t)$ is slightly different:

$\displaystyle\frac{dh^{\mathcal{C}}(\operatorname{\Phi_{f_{\pi}}}(x(t),\tau-t))}{dt}$
$\displaystyle=$
$\displaystyle\frac{dh^{\mathcal{S}}}{dx}(\frac{\partial\operatorname{\Phi_{f_{\pi}}}(x(t),\tau-t)}{\partial
x}\dot{x}-\frac{\partial\operatorname{\Phi_{f_{\pi}}}(x(t),\tau-t)}{\partial
t})$ $\displaystyle=$
$\displaystyle\frac{dh^{\mathcal{S}}}{dx}(\frac{\partial\operatorname{\Phi_{f_{\pi}}}(x(t),\tau-t)}{\partial
x}\dot{x}-f_{\pi}(\operatorname{\Phi_{f_{\pi}}}(x(t),\tau-t))),$

where $\dot{x}=(f(x(t))+g(x(t))u(t))$. By Proposition 1,
$\frac{d\operatorname{\Phi_{f_{\pi}}}(x(t),\tau-t)}{dt}|_{u=\pi(x)}=0$,
indicating that
$\frac{dh^{\mathcal{C}}(\operatorname{\Phi_{f_{\pi}}}(x(t),\tau-t))}{dt}=0$ if
$u(t)=\pi(x(t))$.

######  Remark 1.

The difference between the derivative of
$h^{\mathcal{S}}(\operatorname{\Phi_{f_{\pi}}}(x(t),T))$ and
$h^{\mathcal{C}}(\operatorname{\Phi_{f_{\pi}}}(x(t),\tau-t))$ is that the
former contains two parts, the change of the future state following the backup
strategy due to the change of current state, and the derivative due to $t+T$
increasing with $t$; whereas the latter only contains the first part.

The CBF condition is $\dot{h}+\alpha(h)\geq 0$. We shall impose this condition
on $h^{\mathcal{S}}(\operatorname{\Phi_{f_{\pi}}}_{f_{\pi}}(x,T))$ and
$h^{\mathcal{C}}(\operatorname{\Phi_{f_{\pi}}}_{f_{\pi}}(x,\tau-t))$ for every
$\tau\in[t,t+T]$ instead of only on the $\tau$ minimizing
$h^{\mathcal{C}}(\operatorname{\Phi_{f_{\pi}}}_{f_{\pi}}(x,\tau-t))$, and the
CBF QP is then

$\displaystyle\mathop{\min}\limits_{u\in\mathcal{U}}$
$\displaystyle||u-u^{0}||$ (12) $\displaystyle\mathrm{s.t.}$
$\displaystyle\forall\tau\in[t,t+T]$
$\displaystyle\frac{dh^{\mathcal{C}}(\operatorname{\Phi_{f_{\pi}}}(x,\tau-t))}{dt}(x,u)+\alpha(h^{\mathcal{C}}(\operatorname{\Phi_{f_{\pi}}}(x(t),\tau-t)))\geq
0$
$\displaystyle\frac{dh^{\mathcal{S}}(\operatorname{\Phi_{f_{\pi}}}(x(t),T))}{dt}(x,u)+\alpha(h^{\mathcal{S}}(\operatorname{\Phi_{f_{\pi}}}(x(t),T)))\geq
0.$

######  Proposition 2.

The constraint in (12) is a sufficient condition for $\dot{h}+\alpha(h)\geq
0$.

###### Proof.

For notational simplicity, let $\bar{\xi}^{\mathcal{C}}(x,\tau)\doteq
h^{\mathcal{C}}(\operatorname{\Phi_{f_{\pi}}}(x,\tau-t))$,
$\xi^{\mathcal{C}}(x)\doteq\mathop{\min}\limits_{\tau\in[t,t+T]}\bar{\xi}^{\mathcal{C}}(x,\tau)$,
and $\xi^{\mathcal{S}}\doteq
h^{\mathcal{S}}(\operatorname{\Phi_{f_{\pi}}}(x,T))$, then
$h(x)=\min\\{\xi^{\mathcal{C}}(x),\xi^{\mathcal{S}}(x)\\}$. First notice that

$\frac{d\xi^{\mathcal{C}}(x(t))}{dt}+\alpha(\xi^{\mathcal{C}}(x(t)))\geq\mathop{\min}\limits_{\tau\in[t,t+T]}\left[\dot{\bar{\xi}}^{\mathcal{C}}(x(t),\tau)+\alpha(\bar{\xi}^{\mathcal{C}}(x(t),\tau))\right]\geq
0,$

then

$\begin{aligned}
\dot{h}+\alpha(h)=&\frac{d\min\\{\xi^{\mathcal{C}}(x(t)),\xi^{\mathcal{S}}(x(t))\\}}{dt}+\alpha(\min\\{\xi^{\mathcal{C}}(x(t)),\xi^{\mathcal{S}}(x(t))\\})\\\
\geq&\min\\{\dot{\xi}^{\mathcal{C}}(x(t))+\alpha(\xi^{\mathcal{C}}(x(t))),\dot{\xi}^{\mathcal{S}}(x(t))+\alpha(\xi^{\mathcal{S}}(x(t)))\\}\geq
0.\end{aligned}$

∎

In practice,
$\mathop{\min}\limits_{\tau\in[t,t+T]}h^{\mathcal{C}}(\Phi_{f_{\pi}}(x(t)),\tau-t)$
may be difficult to evaluate, so a finite set of
$t=\tau_{0}<\tau_{1}...<\tau_{N}=t+T$ is used instead of $[t,t+T]$, and the
approximation error can be bounded given the Lipschitz constants of the flow
map.

######  Theorem 2.

Under Assumption 1, the backup CBF QP in (12) is always feasible for $h\geq
0$.

###### Proof.

Since $h(x)\geq 0$,
$\forall\tau\in[t,t+T],h^{\mathcal{C}}(\operatorname{\Phi_{f_{\pi}}}(x(t),\tau-t))\geq
0$, $h^{\mathcal{S}}(\operatorname{\Phi_{f_{\pi}}}(x(t),T))\geq 0$. By
proposition 1, if $u(t)=\pi(x(t))$, for all $\tau\in[t,t+T]$,
$\frac{dh^{\mathcal{C}}(\operatorname{\Phi_{f_{\pi}}}(x(t),\tau-t))}{dt}=0$,
$\frac{dh^{\mathcal{S}}(\operatorname{\Phi_{f_{\pi}}}(x(t),T))}{dt}=\frac{dh^{\mathcal{S}}}{dx}f_{\pi}(\operatorname{\Phi_{f_{\pi}}}(x(t),T))$.
Immediately, the first constraint in (12) is satisfied. By Assumption 1, since
$\operatorname{\Phi_{f_{\pi}}}(x(t),T)\in\mathcal{S}_{0}$, $\nabla
h^{\mathcal{S}}f_{\pi}(\operatorname{\Phi_{f_{\pi}}}(x(t),T))+\alpha(h^{\mathcal{S}}(\operatorname{\Phi_{f_{\pi}}}(x(t),T)))\geq
0$, indicating that the second constraint in (12) is satisfied. Since $\pi$ is
a feasible control policy, $\pi(x(t))\in\mathcal{U}$, thus a feasible solution
to (12). ∎

###  III-B Relative degree of backup CBFs

Another benefit of the backup CBF is that it is always relative degree one
given the system is weakly locally controllable. To show this, we first review
some definitions.

######  Definition 2.

A function $h:\mathbb{R}^{n}\to\mathbb{R}$ is said to have relative degree $r$
with respect to the dynamic system in (1) at a point $x_{0}$ if

$\begin{aligned} &\mathcal{L}_{g}\mathcal{L}_{f}^{k}h(x)=0,&\forall x\text{ in
a neighborhood of $x_{0}$, $\forall k\leq r-2$,}\\\
&\mathcal{L}_{g}\mathcal{L}_{f}^{r-1}h(x_{0})\neq 0&\end{aligned}$

Figure 3: $\\{x|h(x)\geq 0\\}$ under CBF computed with Hamilton Jacobi PDE
(green), backup CBF (red), and SOS (blue)

######  Definition 3.

[18] Given a dynamic system, $x_{1}$ is accessible from $x_{0}$ via
$\Omega\subseteq\mathbb{R}^{n}$ if there exists a control input signal
$\mathbf{u}:[t_{0},t_{1}]\to\mathcal{U}$ such that the trajectory starting at
$x_{0}$ reaches $x_{1}$ at $t_{1}$ and $\forall
t\in[t_{0},t_{1}],x(t)\in\Omega$. The accessible set of a point $x_{0}$ under
$\Omega$ is the set of all points accessible from $x_{0}$ under $\Omega$,
denoted as $A_{\Omega}(x_{0})$.

######  Definition 4.

[18] The dynamic system in (1) is locally weakly controllable if for every
$x_{0}\in\mathbb{R}^{n}$ and every neighborhood $\Omega$ of $x$,
$A_{\Omega}(x_{0})$ has a non-empty interior.

######  Theorem 3.

If the dynamic system (1) is locally weakly controllable, then for an
$h:\mathbb{R}^{n}\to\mathbb{R}$ that is not constant for any subset
$\Omega\subseteq\mathbb{R}^{n}$ with a non-empty interior,
$h\circ\operatorname{\Phi_{f_{\pi}}}(x,t)$ has relative degree 1 apart from
singular points for all $t>0$.

###### Proof.

We prove this by contradiction. Suppose there exists $x_{0}\in\mathbb{R}^{n}$
and a neighborhood $\Omega$ where $\forall
x\in\Omega,\mathcal{L}_{g}(h\circ\operatorname{\Phi_{f_{\pi}}})(x)=0$. Then we
can find a small $t_{0}>0$ such that the flow from $x_{0}$ stays inside
$\Omega$ within $t\in[0,t_{0}]$. Since the system is locally weakly
controllable, $A_{\Omega}(x_{0})$ is a compact set with a non-empty interior.
With a slight abuse of notation, let
$\operatorname{\Phi_{f_{\pi}}}(S,t)=\\{x|x=\operatorname{\Phi_{f_{\pi}}}(x_{0},t),x_{0}\in
S\\}$ be the image of the flow map acting on a subset
$S\subseteq\mathbb{R}^{n}$. For any $t>0$, since we only consider regular
points, $\operatorname{\Phi_{f_{\pi}}}(\cdot,t)$ has rank $n$, thus
$\operatorname{\Phi_{f_{\pi}}}(A_{\Omega}(x_{0}),t)$ has non-empty interior.
However, since $\mathcal{L}_{g}h(x)=0$ for all $x\in\Omega$, $h(x)$ is
constant inside $\operatorname{\Phi_{f_{\pi}}}(A_{\Omega}(x_{0}),t)$, which
contradicts the assumption that $h$ is not constant in any set with non-empty
interior. Thus, $h\circ\operatorname{\Phi_{f_{\pi}}}(x,t+t_{0})$ has relative
degree 1. Since $t_{0}$ can be chosen arbitrarily small,
$h\circ\operatorname{\Phi_{f_{\pi}}}(x,t)$ has relative degree 1 for all
$t>0$. ∎

We cannot show that the CBF defined in (10) has relative degree 1 due to the
$\min$ function over $t\in[0,T]$, however, note that the constraints in the
CBF QP in (12) are on all
$h^{\mathcal{C}}\circ\operatorname{\Phi_{f_{\pi}}}(x,\tau-t)$ with
$\tau\in[t,t+T]$ and $h^{\mathcal{S}}\circ\operatorname{\Phi_{f_{\pi}}}(x,T)$.
Therefore, as long as $h^{\mathcal{S}}$ and $h^{\mathcal{C}}$ are not constant
for any subset with a non-empty interior, the CBF QP can be solved as if $h$
has relative degree 1.

Theorem 3 shows that the backup CBF is sufficient without any high order
extension or backstepping, it naturally bridges the potentially high relative
degree functions $h^{\mathcal{C}}$ and $h^{\mathcal{S}}$ with the input
dynamics via the flow map.

Figure 4: $\\{x|h(x)\geq 0\\}$ with an aggressive backup policy

##  IV Comparative study

This section presents the comparative study of the backup CBF against some
benchmark methods.

One major concern of the backup CBF is that since the control invariant set is
induced by a fixed backup policy, how conservative is the resulting CBF? The
conservatism of the CBF can be measured by the size of the set $\\{x|h(x)\geq
0\\}$, which is the induced control invariant set. We compare the backup CBF
with the Hamilton Jacobi (HJ) PDE result, which is a close approximation of
the maximum control invariant set, and a Sum-of-Squares (SOS) result. The
Hamilton Jacobi formulation follows [9] and is computed with the level-set
toolbox. In particular, the following HJ PDE is solved:

$-\frac{\partial h}{\partial
t}=H(x,p)=\max\limits_{u\in\mathcal{U}}p^{\intercal}f(x,u),$

where $p=\frac{\partial h}{\partial x}$ is the momentum vector.
$\\{x|h(x,t)\geq 0\\}$ is then the backward reachable set at time $t$. For a
sufficiently large $t$, $\lim_{t\to\infty}h(x,t)$ is an inner approximation of
the maximum control invariant set.

The SOS approach follows [19, 20], which computes a valid CBF under a fixed
controller. For fairness, we chose the same backup controller $\pi$ for the
backup CBF and the SOS program.

###  IV-A Lane keeping with Dubin’s car model

To demonstrate the result, a simple lane keeping problem under Dubin’s car
model is considered:

$\dot{x}=\begin{bmatrix}\dot{Y}&\dot{v}&\dot{\psi}\end{bmatrix}^{\intercal}=f(x,u)=\begin{bmatrix}v\sin(\psi)&a&r\end{bmatrix}^{\intercal}$
(13)

where the state $x$ consists of lateral position $Y$, velocity $v$, and
heading angle $\psi$, and the input consists of acceleration $a$ and yaw rate
$r$, The state constraint considered is $\mathcal{C}=\\{|Y|\leq
Y_{\max}\wedge|\psi|\leq\psi_{\max}\\}$ with
$Y_{\max}=1.8m,\psi_{\max}=\pi/3$, and $\mathcal{U}=\\{a,r||a|\leq
a_{\max}\wedge|r|\leq r_{\max}\\}$.

The backup policy is chosen as

$\mathbf{u}(x)=\begin{bmatrix}a\\\
r\end{bmatrix}=\begin{bmatrix}\operatorname{\mathbf{Sat}}_{a_{\max}}(k_{v}(v_{des}-v))\\\
\operatorname{\mathbf{Sat}}_{r_{\max}}(k_{y}[Y;\psi]),\end{bmatrix}$ (14)

where $\operatorname{\mathbf{Sat}}$ is the saturation function, $k_{v}>0$ is a
constant and $k_{y}$ is calculated via LQR. Under the same controller $\pi$, a
CBF is synthesized with the following SOS program:

$\displaystyle\mathop{\max}{}$ $\displaystyle R$ (15)
$\displaystyle\mathrm{s.t.}{}$ $\displaystyle
h(x)-s_{1}(x)(Q(x)-R^{2})\in\Sigma[x]$
$\displaystyle-h(x)-s_{2}(x)(W^{2}-Y^{2})\in\Sigma[x]$ $\displaystyle\nabla
h\cdot f(x,\mathbf{u}(x))+\alpha(h)-s_{3}(x)(Q(x)-R_{0}^{2})\in\Sigma[x]$
$\displaystyle s_{1}(x),s_{2}(x),s_{3}(x)\in\Sigma[x],$

with $Q(x)\doteq\hat{x}^{\intercal}Q\hat{x}$, and
$\hat{x}=[Y,v-v_{des},\psi]^{\intercal}$, where $Q$ is a PSD matrix.
$\Sigma[x]$ is the set of sum-of-squares polynomial of $x$,
$s_{1},s_{2},s_{3}$ are SOS multipliers for the Positivstellensatz procedure
that enforces positive definiteness of a polynomial on a semialgebraic set.
Due to the cross product term in line 2, (15) cannot be directly solved via
SOS, a line search is used to find the largest $R$ that renders the SOS
program feasible. $R^{0}$ is a radius that contains the region of state space
we are interested in, (15) essentially tries to fit the largest ball of $Q(x)$
inside $\\{h(x)\geq 0\\}$ where $h$ is the CBF that satisfies (2). For
simplicity, $\sin\psi\approx\psi-\psi^{3}/3$, and the saturation is lifted,
which should lead to a larger control invariant set.

| Dimension | HJ Offline | HJ online | SOS offline | SOS online | Backup online integration | Backup CBF QP
---|---|---|---|---|---|---|---
Dubin’s car | 3 | 37.3s | 0.5ms | 3.6s | 0.4ms | 1.4ms | 3.6ms
Quadrotor | 16 | NA | NA | NA | NA | 4.92ms | 28.33ms
TABLE I: Computation time of backup CBF, CBF based on HJ, and CBF based on SOS

Fig. 3 shows the comparison of the CBF computed with HJ, backup CBF, and SOS
where the upper plot shows the 3D surface and the lower plot shows the slicing
of the set at $v=5m/s$. As expected, Hamilton Jacobi PDE generates the largest
control invariant set, the one from backup CBF is smaller, and the SOS one is
the smallest since $h$ is restricted to be a polynomial (4-th order in this
case).

Obviously, one can choose a better backup policy and greatly improve the
result. For example, if we change $v_{des}=0$, and switch to a more aggressive
LQR design, the control invariant set computed from backup CBF is almost
identical to the one from Hamilton Jacobi, as shown in Fig. 4.

###  IV-B Aeroplane avoidance

One classic example in safety-critical control is the aeroplane avoidance with
the following dynamics:

$\begin{bmatrix}\Delta\dot{X}\\\ \Delta\dot{Y}\\\
\Delta\dot{\psi}\end{bmatrix}=\begin{bmatrix}-v_{a}+v_{b}\cos(\Delta\Psi)+u\Delta
Y\\\ v_{b}\sin(\Delta\psi)-u\Delta X\\\ -u\end{bmatrix},$ (16)

where $\Delta X$, $\Delta Y$, and $\Delta\psi$ are the difference of $X,Y$
coordinates and heading angles plane $a$ and $b$, $u$ is the turning rate of
plane $a$, $v_{a}$ and $v_{b}$ are their velocities, assumed to be constant.

The backup CBF approach is not good at handling malicious disturbance,
therefore we assume that plane $b$ maintains a fixed orientation. The safety
constraint is defined as $\Delta X^{2}+\Delta Y^{2}\geq R^{2}$, where $R$ is
the minimum distance to maintain between the two planes. Fig. 5 shows the
danger set ($\\{x|h(x)\leq 0\\}$) computed with the HJI PDE and the one
computed with a simple backup policy:

$u=-u_{\max}\mathbf{sign}(\Delta Y).$ (17)

The two danger sets are almost identical.

Figure 5: Danger set of the aeroplane collision avoidance

###  IV-C Computation

Technically, it is not fair to compare the computation time of the backup CBF
approach with the benchmarks as the former does not require any offline
computation while the online computation is more complicated than normal CBF
QP, while CBF based on explicit control invariant set requires computing the
control invariant set offline.

Table I shows the computation time of the CBF QP with CBF constructed by the
HJ PDE, SOS, and the backup CBF approach. The computation time for backup CBF
consists of two parts, the integration of the backup policy, and the CBF QP.
As shown by the result, the majority of solver time is spent on the CBF QP,
which is 1 magnitude larger than the time used to solve the CBF QP with
explicit form. This is mainly because (12) enforces the CBF condition on every
$\tau\in[t,t+T]$ (replaced with a sequence of $\tau_{i}$ in practice), which
significantly increases the number of constraints. The integration time is
relatively small compared to the QP time and can be further reduced. For
instance, [21] uses a Koopman operator approach to further simplify the online
integration. We also apply the backup CBF on a 16-dimensional quadrotor model,
which is way beyond the limit of HJ (maximum dimension 4-5) and SOS (maximum
dimension around 8-10), and the backup CBF can still be implemented with a
reasonable loop rate.

##  V Conclusion

We give a tutorial of the backup CBF approach and compare it to the HJ PDE
approach and the SOS approach as benchmarks. The result shows that the
implicit control invariant set induced by the backup policy is close to the
maximum control invariant set under a properly chosen backup policy in many
practical problems. The backup CBF is much more scalable than the benchmark
methods with explicitly computed control invariant sets and is applicable to
general nonlinear dynamics. Furthermore, since the backup CBF has relative
degree 1 under mild assumptions, it is a better choice than the high-order
CBFs that bear no feasibility guarantee.

## References

  * [1] A. D. Ames, J. W. Grizzle, and P. Tabuada, “Control barrier function based quadratic programs with application to adaptive cruise control,” in _53rd IEEE Conference on Decision and Control_. IEEE, 2014, pp. 6271–6278.
  * [2] P. Wieland and F. Allgöwer, “Constructive safety using control barrier functions,” _IFAC Proceedings Volumes_ , vol. 40, no. 12, pp. 462–467, 2007\.
  * [3] J.-P. Aubin, A. M. Bayen, and P. Saint-Pierre, _Viability theory: new directions_. Springer Science & Business Media, 2011.
  * [4] D. Bertsekas, “Infinite time reachability of state-space regions by using feedback control,” _IEEE Transactions on Automatic Control_ , vol. 17, no. 5, pp. 604–613, 1972.
  * [5] F. Blanchini, “Set invariance in control,” _Automatica_ , vol. 35, no. 11, pp. 1747–1767, 1999.
  * [6] S. Rakovic, P. Grieder, M. Kvasnica, D. Mayne, and M. Morari, “Computation of invariant sets for piecewise affine discrete time systems subject to bounded disturbances,” in _2004 43rd IEEE Conference on Decision and Control (CDC)(IEEE Cat. No. 04CH37601)_ , vol. 2. IEEE, 2004, pp. 1418–1423.
  * [7] Y. Chen, H. Peng, J. Grizzle, and N. Ozay, “Data-driven computation of minimal robust control invariant set,” in _2018 IEEE Conference on Decision and Control (CDC)_. IEEE, 2018, pp. 4052–4058.
  * [8] M. Korda, D. Henrion, and C. N. Jones, “Convex computation of the maximum controlled invariant set for polynomial control systems,” _SIAM Journal on Control and Optimization_ , vol. 52, no. 5, pp. 2944–2969, 2014.
  * [9] I. M. Mitchell, A. M. Bayen, and C. J. Tomlin, “A time-dependent hamilton-jacobi formulation of reachable sets for continuous dynamic games,” _IEEE Transactions on automatic control_ , vol. 50, no. 7, pp. 947–957, 2005\.
  * [10] A. Taylor, A. Singletary, Y. Yue, and A. Ames, “Learning for safety-critical control with control barrier functions,” in _Learning for Dynamics and Control_. PMLR, 2020, pp. 708–717.
  * [11] L. Lindemann and D. V. Dimarogonas, “Control barrier functions for signal temporal logic tasks,” _IEEE control systems letters_ , vol. 3, no. 1, pp. 96–101, 2018.
  * [12] Q. Nguyen and K. Sreenath, “Exponential control barrier functions for enforcing high relative-degree safety-critical constraints,” in _2016 American Control Conference (ACC)_. IEEE, 2016, pp. 322–328.
  * [13] X. Xu, “Constrained control of input–output linearizable systems using control sharing barrier functions,” _Automatica_ , vol. 87, pp. 195–201, 2018.
  * [14] W. Xiao and C. Belta, “Control barrier functions for systems with high relative degree,” in _2019 IEEE 58th Conference on Decision and Control (CDC)_. IEEE, 2019, pp. 474–479.
  * [15] T. Gurriet, M. Mote, A. Singletary, P. Nilsson, E. Feron, and A. D. Ames, “A scalable safety critical control framework for nonlinear systems,” _IEEE Access_ , vol. 8, pp. 187 249–187 275, 2020.
  * [16] H. Seywald and R. Kumar, “Desensitized optimal trajectories,” _Analytical Mechanics Associates Rept_ , pp. 03–16, 2003.
  * [17] T. Gurriet, M. Mote, A. D. Ames, and E. Feron, “An online approach to active set invariance,” in _2018 IEEE Conference on Decision and Control (CDC)_. IEEE, 2018, pp. 3592–3599.
  * [18] R. Hermann and A. Krener, “Nonlinear controllability and observability,” _IEEE Transactions on automatic control_ , vol. 22, no. 5, pp. 728–740, 1977\.
  * [19] X. Xu, J. W. Grizzle, P. Tabuada, and A. D. Ames, “Correctness guarantees for the composition of lane keeping and adaptive cruise control,” _IEEE Transactions on Automation Science and Engineering_ , vol. 15, no. 3, pp. 1216–1229, 2017.
  * [20] Y. Chen, A. Hereid, H. Peng, and J. Grizzle, “Enhancing the performance of a safe controller via supervised learning for truck lateral control,” _Journal of Dynamic Systems, Measurement, and Control_ , vol. 141, no. 10, 2019.
  * [21] C. Folkestad, Y. Chen, A. D. Ames, and J. W. Burdick, “Data-driven safety-critical control: Synthesizing control barrier functions with koopman operators,” _IEEE Control Systems Letters_ , 2020.