# Neutron star mass formula with nuclear saturation parameters for asymmetric
nuclear matter

Hajime Sotani<EMAIL_ADDRESS>Astrophysical Big Bang Laboratory,
RIKEN, Saitama 351-0198, Japan  Interdisciplinary Theoretical & Mathematical
Science Program (iTHEMS), RIKEN, Saitama 351-0198, Japan  Shinsuke Ota
Research Center for Nuclear Physics (RCNP), Osaka University, Ibaragi, Osaka
567-0047, Japan

###### Abstract

Low-mass neutron stars are directly associated with the nuclear saturation
parameters because their central density is definitely low. We have already
found a suitable combination of nuclear saturation parameters for expressing
the neutron star mass and gravitational redshift, i.e.,
$\eta\equiv(K_{0}L^{2})^{1/3}$ with the incompressibility for symmetric
nuclear matter, $K_{0}$, and the density-dependent nuclear symmetry energy,
$L$. In this study, we newly find another suitable combination given by
$\eta_{\tau}\equiv(-K_{\tau}L^{5})^{1/6}$ with the isospin dependence of
incompressibility for asymmetric nuclear matter, $K_{\tau}$, and derive the
empirical relations for the neutron star mass and gravitational redshift as a
function of $\eta_{\tau}$ and the normalized central number density. With
these empirical relations, one can evaluate the mass and gravitational
redshift of the neutron star, whose central number density is less than
threefold the saturation density, within $\sim 10\%$ accuracy, and the radius
within a few % accuracy. In addition, we discuss the neutron star mass and
radius constraints from the terrestrial experiments, using the empirical
relations, together with those from the astronomical observations.
Furthermore, we find a tight correlation between $\eta_{\tau}$ and $\eta$.
With this correlation, we derive the constraint on $K_{\tau}$ as $-348\leq
K_{\tau}\leq-237$ MeV, assuming that $L=60\pm 20$ and $K_{0}=240\pm 20$ MeV.

###### pacs:

04.40.Dg, 26.60.+c, 21.65.Ef

##  I Introduction

Neutron star is a massive remnant left after a supernova explosion, which
happens at the last moment of the massive star’s life. The density inside the
star becomes much higher than the standard nuclear density,
$\rho_{0}=2.7\times 10^{14}$ g/cm3, and the gravitational and magnetic fields
inside/around the star are significantly stronger than those observed in our
solar system [1]. The neutron star mass and radius strongly depend on the
equation of state (EOS) for neutron star matter under the $\beta$-equilibrium.
The mass of a neutron star model with a higher central density generally
becomes larger, even though the EOS is not fixed yet. So, the astronomical
observations of neutron stars or their activities tell us the information
about the EOS for a relatively higher density region, while the terrestrial
experiments tell us that for a lower density region (e.g., Fig. 2 in Ref.
[2]).

In practice, the discovery of $2M_{\odot}$ neutron stars [3, 4, 5, 6] has
excluded some soft EOSs, with which the expected maximum mass is less than the
observed mass. This discovery simultaneously reveals the problem that most of
the EOSs with hyperon are too soft to construct the $2M_{\odot}$ neutron star,
i.e., the so-called hyperon puzzle. Meanwhile, the gravitational wave event,
GW170817 [7], tells us the information on the tidal deformability of the
neutron star, which leads to the constraint that a $1.4M_{\odot}$ neutron star
radius should be less than 13.6 km [8]. We note that the constraint on neutron
star radius may become more stringent in view of the existing multi-messenger
observational data [9, 10]. The light bending due to the strong gravitational
field induced by the neutron star is also one of the important phenomena to
see the neutron star properties. That is, the pulsar light curve from the
rotating neutron star would be modified due to this relativistic effect and
one could get the neutron star properties by carefully observing it (e.g.,
[11, 12, 13, 14, 15, 16]). Based on this idea, the Neutron star Interior
Composition Explorer (NICER) is now operating on an International Space
Station (ISS) and it has already announced the constraint on two neutron stars
properties, i.e., PSR J0030+0451 [17, 18] and PSR J0740+6620 [19, 20].
Furthermore, the direct detection of the gravitational waves from the neutron
star in the future may enable us to extract the neutron star properties (e.g.,
[21, 22, 23, 24, 25, 26, 27, 28, 29]).

On the other hand, the EOS in a lower density region is also gradually
constrained through terrestrial nuclear experiments, but still, there are
large uncertainties in EOS parameters (or in neutron star properties)
constrained from terrestrial experiments. For instance, the fiducial value of
the density-dependent nuclear symmetry energy $L$ is $L\simeq 60\pm 20$ MeV
[30, 31], while the constraints of $L$ obtained recently seem to be
significantly larger than the fiducial value [32, 33]. This is because one has
to usually transform the experimental constraint to the EOS parameters, even
if the information determined via experiments is associated with some aspects
of nuclear EOS. Then, one can eventually discuss the neutron star mass and
radius as a solution of the Tolman-Oppenheimer-Volkoff (TOV) equation. Anyway,
the terrestrial experiments are definitely crucial for understanding the
neutron star EOS as well as the astronomical observation of neutron stars.

In order to discuss the neutron star mass and radius directly with the EOS
parameters constrained somehow from the experiments, it may be useful to
derive a kind of empirical formula expressing the neutron star properties as a
function of EOS parameters if it exists. In practice, we have already found
the empirical relations expressing the mass ($M$) and gravitational redshift
($z$) of a low-mass neutron star as a function of the suitable combination of
the nuclear saturation parameters, $\eta$, and the central density of neutron
star [34], with which the neutron star mass and radius expected from the
terrestrial experiments can be discussed as in Ref. [2]. On the other hand, in
this study, we consider deriving another type of empirical relations for the
low-mass neutron stars, focusing on the nuclear saturation parameters for
asymmetric nuclear matter, and then discuss the impact of the uncertainties in
the nuclear saturation parameters on the neutron star mass and radius.

This manuscript is organized as follows. In Sec. II, we briefly mention the
EOSs considered in this study and the saturation parameters in nuclear matter.
In Sec. III, we examine the neutron star models and derive the empirical
formulas for the neutron star mass and its gravitational redshift as a
function of a suitable combination of the nuclear saturation parameters we
derived in this study. Using the newly derived empirical formulas, in Sec. IV,
we discuss the neutron star mass and radius together with the constraints from
the astronomical and experimental observations. Finally, in Sec. V, we
conclude this study. Unless otherwise mentioned, we adopt geometric units in
the following, $c=G=1$, where $c$ and $G$ denote the speed of light and the
gravitational constant, respectively.

Table 1: EOS parameters adopted in this study, $K_{0}$, $n_{0}$, $L$, $Q$, $K_{sym}$, $Q_{sym}$, $K_{\tau}$, and $Q_{\tau}$ are listed, while $\eta$ and $\eta_{\tau}$ are specific combinations with them given by $\eta=\left(K_{0}L^{2}\right)^{1/3}$ and $\eta_{\tau}=\left(-{K_{\tau}}L^{5}\right)^{1/6}$.  EOS | $K_{0}$ | $n_{0}$ | $L$ | $Q$ | $K_{sym}$ | $Q_{sym}$ | $K_{\tau}$ | $Q_{\tau}$ | $\eta$ | $\eta_{\tau}$
---|---|---|---|---|---|---|---|---|---|---
| (MeV) | (fm-3) | (MeV) | (MeV) | (MeV) | (MeV) | (MeV) | (MeV) | (MeV) | (MeV)
OI-EOSs | 200 | 0.165 | 35.6 | -759 | -142 | 801 | -221 | 2017 | 63.3 | 48.2
|  | 0.165 | 67.8 | -761 | -27.6 | 589 | -176 | 2909 | 97.2 | 79.5
| 220 | 0.161 | 40.2 | -720 | -144 | 731 | -254 | 1915 | 70.9 | 54.7
|  | 0.161 | 77.6 | -722 | -9.83 | 486 | -221 | 2779 | 110 | 92.4
| 240 | 0.159 | 45.0 | -663 | -146 | 642 | -291 | 1760 | 78.6 | 61.4
|  | 0.158 | 88.2 | -664 | 10.5 | 363 | -274 | 2559 | 123 | 107
| 260 | 0.156 | 49.8 | -589 | -146 | 535 | -333 | 1551 | 86.4 | 68.4
|  | 0.155 | 99.2 | -590 | 32.6 | 219 | -338 | 2246 | 137 | 122
| 280 | 0.154 | 54.9 | -496 | -146 | 410 | -378 | 1285 | 94.5 | 75.7
|  | 0.153 | 111 | -498 | 57.4 | 54.4 | -412 | 1834 | 151 | 138
| 300 | 0.152 | 60.0 | -386 | -146 | 266 | -429 | 962 | 103 | 83.3
|  | 0.151 | 124 | -387 | 86.1 | -133 | -499 | 1310 | 167 | 157
KDE0v | 229 | 0.161 | 45.2 | -373 | -145 | 523 | -342 | 1187 | 77.6 | 63.4
KDE0v1 | 228 | 0.165 | 54.7 | -385 | -127 | 484 | -363 | 1317 | 88.0 | 75.0
SLy2 | 230 | 0.161 | 47.5 | -364 | -115 | 507 | -325 | 1183 | 80.3 | 65.4
SLy4 | 230 | 0.160 | 45.9 | -363 | -120 | 522 | -323 | 1175 | 78.7 | 63.7
SLy9 | 230 | 0.151 | 54.9 | -350 | -81.4 | 462 | -327 | 1215 | 88.4 | 73.9
SKa | 263 | 0.155 | 74.6 | -300 | -78.5 | 175 | -441 | 940 | 114 | 100
SkI3 | 258 | 0.158 | 101 | -304 | 73.0 | 212 | -412 | 1276 | 138 | 127
SkMp | 231 | 0.157 | 70.3 | -338 | -49.8 | 159 | -369 | 1086 | 105 | 92.7

##  II EOS for neutron star matter

To construct the neutron star models theoretically, one has to prepare the EOS
for neutron star matter. We note that, in order to discuss the neutron star
properties with the EOS parameters, one has to adopt the unified EOS. That is,
the EOS describing the neutron star core is constructed with the same nuclear
model as in the EOS for the neutron star crust. We note that, if one
constructs the neutron star model with the EOS (unlike the unified EOS), which
is assembled by connecting the EOS for the core region to the different EOS
for the crust region at an appropriate transition density, the radius of a
neutron star whose central density is around the transition density strongly
depends on the selection of the transition density (e.g., Ref. [35]). That is,
unless the unified EOSs are adopted, one can not discuss the empirical
relations, expressing the neutron star models, as we will derive in Sec. III.
In this study, to systematically see the EOS dependence of the neutron star
properties, we particularly adopt the phenomenological EOS proposed by
Oyamatsu and Iida [36, 37] (hereafter referred to as OI-EOS) and the EOSs with
the Skyrme-type interaction listed in Table 1.

OI-EOSs are based on the Padé-type potential energies and constructed in such
a way that the nucleus models with a simplified version of the extended
Thomas-Fermi theory should become consistent with the empirical masses and
radii of stable nuclei. We note that most of the OI-EOSs adopted here may not
fulfill either constraint of $2M_{\odot}$ observation or the $1.4M_{\odot}$
neutron star radius constrained from the GW170817, because the OI-EOSs are
constructed by focusing on the behavior around the saturation point.
Nevertheless, since we discuss the neutron star properties constructed with
the central density, which is not so high, and since the OI-EOSs are
definitely suitable for systematical study, we adopt the OI-EOSs in this
study. On the other hand, we also consider the EOSs with various Skyrme-type
interactions, such as KDE0v, KDE0v1 [38], SLy2, SLy4, SLy9 [39, 40], SKa [41],
SkI3 [42], and SkMp [43]. We note that the non-relativistic EOSs may break the
causality in a high density region, where the resultant neutron star models
are not realistic. But, in this study, we focus only on the density region
less than threefold the saturation density, where all the EOSs adopted in this
study do not break the causality.

Even though the EOS generally depends on the nuclear model, compositions, and
interaction, the bulk energy per nucleon for the uniform nuclear matter at
zero temperature can anyhow be expressed as a function of the baryon number
density, $n_{\rm b}$, and an asymmetry parameter, $\alpha$, such as

$\frac{E}{A}=w_{s}(n_{\rm b})+\alpha^{2}S(n_{\rm b})+{\cal O}(\alpha^{3}),$
(1)

where $n_{b}$ and $\alpha$ are given by $n_{\rm b}=n_{n}+n_{p}$ and
$\alpha=(n_{n}-n_{p})/n_{\rm b}$ with the neutron number density, $n_{n}$, and
the proton number density, $n_{p}$; $w_{s}$ corresponds to the energy per
nucleon of symmetric nuclear matter ($\alpha=0$); and $S$ denotes the density-
dependent symmetry energy. Additionally, $w_{s}$ and $S$ can be expanded
around the saturation density, $n_{0}$, of the symmetric nuclear matter as a
function of $u=(n_{\rm b}-n_{0})/(3n_{0})$;

$\displaystyle w_{s}(n_{\rm
b})=w_{0}+\frac{K_{0}}{2}u^{2}+\frac{Q}{6}u^{3}+{\cal O}(u^{4}),$ (2)
$\displaystyle S(n_{\rm
b})=S_{0}+Lu+\frac{K_{sym}}{2}u^{2}+\frac{Q_{sym}}{6}u^{3}+{\cal O}(u^{4}).$
(3)

The coefficients in this expansion correspond to the nuclear saturation
parameters, and each EOS has its own set of nuclear saturation parameters.
Among the nuclear saturation parameters, $n_{0}$, $w_{0}$, and $S_{0}$ are
especially well-constrained from the terrestrial experiments, i.e.,
$n_{0}\approx 0.15-0.16$ fm-3, $w_{0}\approx-15.8$ MeV [44], and $S_{0}\approx
31.6\pm 2.7$ MeV [31]. The constraint on $K_{0}$ and $L$ is relatively more
difficult. This is because one needs to obtain the nuclear data in a wide
density range to constrain $K_{0}$ and $L$, which are a kind of density
derivative. The current fiducial values of $K_{0}$ and $L$ are $K_{0}=240\pm
20$ MeV [45] and $L=60\pm 20$ MeV [30, 31]. Meanwhile, since the determination
of the values of $K_{0}$ and $L$ are strongly model-dependent, one may
consider $K_{0}=200-315$ MeV [46, 47] and $L=20-145$ MeV [2] as their
conservative values, which cover almost all predictions from various
experiments. The experimental constraints on the saturation parameters for
higher order terms, such as $Q$, $K_{sym}$, and $Q_{sym}$, are almost
disorganized, but they are theoretically evaluated as $-800\leq Q\leq 400$,
$-400\leq K_{sym}\leq 100$, and $-200\leq Q_{sym}\leq 800$ MeV [31].

On the other hand, considering $E/A$ as the energy per particle of infinite
asymmetric nuclear matter, one can also expand it around an isospin dependent
saturation density, $\tilde{n}_{0}(\alpha)\simeq
n_{0}\left[1-3(L/K_{0})\alpha^{2}\right]$ [48], such as

$\frac{E}{A}=\tilde{w}_{0}(\tilde{n}_{0})+\frac{\tilde{K}_{0}(\tilde{n}_{0})}{2}\tilde{u}^{2}+\frac{\tilde{Q}(\tilde{n}_{0})}{6}\tilde{u}^{3}+{\cal
O}(\tilde{u}^{4}),$ (4)

where $\tilde{u}$ is defined as $\tilde{u}=(n_{\rm
b}-\tilde{n}_{0})/(3\tilde{n}_{0})$ and the coefficients in this expansion are
related to the saturation parameters in Eqs. (2) and (3) through

$\displaystyle\tilde{w}_{0}(\tilde{n}_{0})=w_{0}+S_{0}\alpha^{2}+{\cal
O}(\alpha^{3}),$ (5)
$\displaystyle\tilde{K}_{0}(\tilde{n}_{0})=K_{0}+K_{\tau}\alpha^{2}+{\cal
O}(\alpha^{3}),$ (6)
$\displaystyle\tilde{Q}_{0}(\tilde{n}_{0})=Q+Q_{\tau}\alpha^{2}+{\cal
O}(\alpha^{3}).$ (7)

Here, $K_{\tau}$ and $Q_{\tau}$ are respectively the isospin dependence of
incompressibility and skewness coefficient at the saturation density,
$\tilde{n}_{0}(\alpha)$, given by

$\displaystyle K_{\tau}=K_{sym}-6L-\frac{Q}{K_{0}}L,$ (8) $\displaystyle
Q_{\tau}=Q_{sym}-\frac{9Q}{K_{0}}L.$ (9)

Compared to the constraints on $K_{0}$ and $L$, the constraints on $K_{\tau}$
are still very poor [49, 50]. The constraints obtained from the experiments
are $K_{\tau}=-550\pm 100$ MeV [51] from the analysis of the Isoscalar Giant
Monopole Resonance (ISGMR) in Sn isotopes and $K_{\tau}=-555\pm 75$ MeV from
the the analysis of the ISGMR in Cd isotopes [52], which are performed at the
Research Center for Nuclear Physics (RCNP), Osaka University, Japan. On the
other hand, using the experimental data in RCNP together with those in the
Texas A&M University (TAMU), the theoretical constraints are also obtained:
$K_{\tau}=-550\pm 30$ MeV from the analysis of the ISGMR in Pb and Sn isotopes
[53]; $-840<K_{\tau}<-350$ MeV from the reanalysis of ISGMR in Sn, Cd, and
data on $58\leq A\leq 208$ nuclei [54]; $K_{\tau}=-500_{-100}^{+120}$ MeV from
the analysis of the neutron-skin data for different anti-protonic atoms [55];
$K_{\tau}=-370\pm 120$ MeV from the analysis of the isotopic transport ratios
[56]; and $K_{\tau}=-500\pm 50$ MeV from the analysis of the ISGMR in Sn
isotopes [57]. In Fig. 1 we show the experimental and theoretical constraints
on $K_{\tau}$ obtained so far. In this study we particuraly adopt the result
obtained in Ref. [51], i.e., $K_{\tau}=-550\pm 100$ MeV, as a fiducial value
of $K_{\tau}$.

Figure 1: Experimental and theoretical constraints on $K_{\tau}$ obtained so
far. The constraint of $-348\leq K_{\tau}\leq-237$ MeV ($-718\leq
K_{\tau}\leq-245$ MeV) is obtained from the empirical relation given by Eq.
(12), assuming that $L=60\pm 20$ MeV [30, 31] and $K_{0}=240\pm 20$ MeV [45]
($L=20-145$ MeV [2] and $K_{0}=200-315$ MeV [46, 47]).

Since the neutron star properties are determined by solving the balance
equations, i.e., the TOV equation, together with the EOS depending on several
physical inputs, it may generally be difficult to predict the neutron star
properties as a function of the EOS parameters. Nevertheless, we could find a
suitable combination of the nuclear saturation parameters, $\eta$, for
expressing the low-mass neutron star [34]. Moreover, as mentioned in the
following section, we newly find another suitable combination of the
saturation parameters, $\eta_{\tau}$, for expressing the neutron stars. These
suitable combinations, $\eta$ and $\eta_{\tau}$, are given by

$\displaystyle\eta=\left(K_{0}L^{2}\right)^{1/3},$ (10)
$\displaystyle\eta_{\tau}=\left(-K_{\tau}L^{5}\right)^{1/6}.$ (11)

We note that the finding of $\eta$ and $\eta_{\tau}$ is just through trial and
error, and unfortunately we could not understand the physics behind these
quantities. Since one can express the neutron star properties with $\eta$ or
$\eta_{\tau}$, it is natural to expect the existence of a tight correlation
between $\eta$ and $\eta_{\tau}$. In fact, there is no correlation between
$K_{\tau}$ and $L$ as shown in the left panel of Fig. 2, while we find the
tight correlation between $\eta$ and $\eta_{\tau}$ as shown in the right panel
of Fig. 2, adopting the 304 models for OI-EOSs and 240 models for the Skyrme-
type EOSs listed in Ref. [58]. We note that, for plotting Fig. 2, we omit some
of the 240 models in Ref. [58] with the negative value of $L$. From this
figure, we can derive the empirical relation, such as

$\eta_{\tau}=61.439\eta_{100}+20.965\eta_{100}^{2},$ (12)

where $\eta_{100}$ is defined by $\eta_{100}\equiv\eta/(100\ {\rm MeV})$. From
this empirical relation together with the fiducial value of $\eta$, i.e.,
$70.6\leq\eta\leq 118.5$ MeV ($L=60\pm 20$ MeV [30, 31] and $K_{0}=240\pm 20$
MeV [45]), one can expect $53.8\leq\eta_{\tau}\leq 102.2$ MeV, which leads to
$-348\leq K_{\tau}\leq-237$ MeV. In a similar way, using
$30.4\leq\eta_{\tau}\leq 189.3$ MeV estimated with the conservative value of
$\eta$, i.e., $43.1\leq\eta\leq 187.8$ MeV ($L=20-145$ MeV [2] and
$K_{0}=200-315$ MeV [46, 47]), one can get the constraint on $K_{\tau}$ as
$-718\leq K_{\tau}\leq-245$ MeV. These constraints on $K_{\tau}$ are also
shown in Fig. 1. From Fig. 1, one can observe that the value of $K_{\tau}$
constrained with the conservative values of $K_{0}$ and $L$ through Eq. (12)
is more or less similar to the conservative value of $K_{\tau}$.

Figure 2: The left panel is the relation between $L$ and $K_{\tau}$ for the
OI-EOSs (open circles) and the Skyrme type EOS models (filled circles). The
right panel is the relation between $\eta$ and $\eta_{\tau}$. The thick-solid
line denote the fitting formula given by Eq. (12). With this empirical
relation together with the fiducial value of $\eta$, i.e., $70.6\leq\eta\leq
118.5$ MeV (shaded region), the expected value of $\eta_{\tau}$ is in the
range of $53.8\leq\eta_{\tau}\leq 102.2$ MeV (shaded region). In a similar
way, the value of $\eta_{\tau}$ is expected in the range of
$30.4\leq\eta_{\tau}\leq 189.3$ MeV (the region between two horizontal dashed
lines) with the conservative value of $\eta$, i.e., $43.1\leq\eta\leq 187.8$
MeV (the region between two vertical dashed lines).

##  III Neutron star mass formula

In the previous studies, we could derive the empirical formula expressing the
neutron star mass and gravitational redshift as a function of $\eta$ and
$u_{c}=\rho_{c}/\rho_{0}$ (or $n_{c}/n_{0}$), where $\rho_{c}$ and $\rho_{0}$
are respectively the central density of neutron star and the saturation
density, while $n_{c}$ is the baryon number density at the stellar center [34,
59]. With $\eta$, one can also discuss the rotational properties in low-mass
neutron stars [60] and the possible maximum mass of neutron stars [61, 62].

Now, we try to derive another type of empirical relation for the mass and
gravitational redshift of the neutron star in a similar way considered in Ref.
[34, 59]. The neutron star mass is determined by integrating the TOV equation,
assuming the central density, where the resultant mass also depends on the EOS
parameters. As a result of trial and error by hand, we find a suitable
combination of $K_{\tau}$ and $L$, i.e., $\eta_{\tau}$ given by Eq. (11), with
which the neutron star mass with a fixed central number density can be
characterized well, as shown in the top panel of Fig. 3. In the same way, we
find that the gravitational redshift, $z$, given by
$z=(1-2GM/Rc^{2})^{-1/2}-1$, for a fixed central number density is also well
expressed as a function of $\eta_{\tau}$, as shown in the bottom panel of Fig.
3. From this result, we can derive the fitting formulae for the neutron star
mass and gravitational redshift for each central number density as a function
of $\eta_{\tau}$, such as

$\displaystyle
M/M_{\odot}=a_{0}^{m}+a_{1}^{m}\eta_{\tau,100}+a_{2}^{m}\eta_{\tau,100}^{2}+a_{3}^{m}\eta_{\tau,100}^{3},$
(13) $\displaystyle
z=a_{0}^{z}+a_{1}^{z}\eta_{\tau,100}+a_{2}^{z}\eta_{\tau,100}^{2}+a_{3}^{z}\eta_{\tau,100}^{3},$
(14)

where $\eta_{\tau,100}$ is defined by $\eta_{\tau,100}\equiv\eta_{\tau}/(100\
{\rm MeV})$, while $a_{i}^{m}$ and $a_{i}^{z}$ for $i=0-3$ are coefficients
depending on $n_{c}/n_{0}$.

Figure 3: The mass (top panel) and gravitational redshift (bottom panel) for
the neutron star models with $n_{c}/n_{0}=1$, 2, and 3, constructed with
various EOSs, are plotted as a function of $\eta_{\tau}$. The thick-solid
lines are fitting lines given by Eqs. (13) and (14).

Moreover, the coefficients in Eqs. (13) and (14), i.e., $a_{i}^{m}$ and
$a_{i}^{z}$, are calculated by varying $n_{c}/n_{0}$, which are shown in Fig.
4. From this result, we can derive the fitting formulae for $a_{i}^{m}$ and
$a_{i}^{z}$ as a function of $n_{c}/n_{0}$, such as

$\displaystyle a_{i}^{m}=\sum_{j=0}^{4}a_{ij}^{m}(n_{c}/n_{0})^{j},$ (15)
$\displaystyle a_{i}^{z}=\sum_{j=0}^{4}a_{ij}^{z}(n_{c}/n_{0})^{j}.$ (16)

The concrete values of the coefficients in these equations, $a_{ij}^{m}$ and
$a_{ij}^{z}$, are listed in Table 2. At last, we can derive the empirical
relations for neutron star mass and gravitational redshift as a function of
$\eta_{\tau}$ and $n_{c}/n_{0}$, given by Eqs. (13) – (16).

Figure 4: The coefficients in the fitting formula (13) and (14) are shown as a function of $n_{c}/n_{0}$, where the top and bottom panels correspond to the results for $a_{i}^{m}$ and $a_{i}^{z}$, respectively. In the both panels, the solid lines denote the fitting of $a_{i}^{m}$ and $a_{i}^{z}$ given by Eqs. (15) and (16).  Table 2: Values of $a_{ij}^{m}$ and $a_{ij}^{z}$ in Eqs. (15) and (16). $j$ | 0 | 1 | 2 | 3 | 4
---|---|---|---|---|---
$a_{0j}^{m}$ | $-0.4498$ | $1.0714$ | $-0.5067$ | $-0.09871$ | $0.03648$
$a_{1j}^{m}$ | $2.2712$ | $-4.9277$ | $2.4996$ | $0.2726$ | $-0.1296$
$a_{2j}^{m}$ | $-2.4238$ | $4.6837$ | $-1.2923$ | $-0.7641$ | $0.1883$
$a_{3j}^{m}$ | $0.6795$ | $-1.1991$ | $0.1422$ | $0.3175$ | $-0.06662$
$a_{0j}^{z}$ | $-0.03506$ | $0.04173$ | $0.01848$ | $-0.03742$ | $0.006334$
$a_{1j}^{z}$ | $0.1561$ | $-0.3211$ | $0.09294$ | $0.09648$ | $-0.01888$
$a_{2j}^{z}$ | $-0.1884$ | $0.3743$ | $-0.03323$ | $-0.1200$ | $0.02129$
$a_{3j}^{z}$ | $0.05780$ | $-0.1110$ | $-0.001509$ | $0.04061$ | $-0.006882$

Here, we check how well the empirical relations derived in this study work. In
Fig. 5 we show the relative deviation of the neutron star mass (top panel),
gravitational redshift (middle panel), and radius (bottom panel) estimated
with the empirical relations from those for a TOV solution with various EOSs.
From this figure, one can observe that the neutron star mass and gravitational
redshift are estimated with the empirical relations within $\sim 10\%$
accuracy, while the radius for neutron star models with $n_{c}/n_{0}\geq 1.7$
is estimated within a few % accuracy.

Figure 5: The relative deviation of neutron star mass (top panel) and
gravitational redshift (middle panel) estimated with the empirical formulas
from those determined with the specific EOSs is shown as a function of
$n_{c}/n_{0}$. The bottom panel is the relative deviation of the neutron star
radius estimated with the empirical formulas for the stellar mass and
gravitational redshift from that of the TOV solution.

##  IV Neutron star mass and radius relation

To simultaneously discuss the EOS parameters constrained from the terrestrial
experiments together with the EOSs for a higher density region constrained
from the astronomical observations, the neutron star mass and radius must be
reasonable. In this section, first, we briefly mention the constraints on the
neutron star mass and radius obtained from several astronomical observations.
Then, we discuss the impact of the EOS parameters constrained from the
terrestrial experiments on the neutron star mass and radius by using the
empirical formulae, Eqs. (13) – (16), derived in this study.

###  IV.1 Constraints from astronomical observations

The mass of a neutron star is one of the most important pieces of information
for constraining the EOS. In particular, the massive star is more important to
exclude the soft EOS, with which the expected maximum mass does not reach the
observed mass. Up to now, the maximum mass precisely determined is $M=2.08\pm
0.07M_{\odot}$ (68.3% credibility) for MSP J0740+6620 [5, 6]. Owing to the
NICER observations, the neutron star mass and radius for PSR J0030+0451 [17,
18] and MSP J0740+6620 [19, 20] are also constrained. The resultant $68\%$ and
$95\%$ confidential levels are shown in Fig. 6. In addition, the gravitational
wave observed from the binary neutron star merger, GW170817, told us the tidal
deformability of the neutron star, which leads to the constraint on the
$1.4M_{\odot}$ radius as $R_{1.4}\leq 13.6$ km [8]. This constraint must be
conservative, while more stringent constraints on the $1.4M_{\odot}$ radius
are also suggested from multimessenger observations and nuclear theory, i.e.,
$R_{1.4}=11.0^{+0.9}_{-0.6}$ km (90% confidence) [9] and
$R_{1.4}=11.75^{+0.86}_{-0.81}$ km (90% confidence) [10]. Furthermore, the
neutron star mass and radius are constrained from the x-ray bursts
observations [63], although this constraint may strongly depend on the
theoretical model. Meanwhile, as a theoretical constraint, the left-top region
on Fig. 6 should exclude due to the causality [64]. We collect these
constraints in Fig. 6. From this figure, one can see that the constraints from
the astronomical observations seem to be still uncertain, but they will
gradually become better. In the same figure, we also plot the neutron star
mass and radius relations constructed with some of the specific EOSs. We note
that the Shen EOS is ruled out from the $1.4M_{\odot}$ radius constraint, but
we show it because the Shen EOS is one of the standard EOSs adopted in many
astrophysical studies.

Figure 6: The expected region in the relation between neutron star mass and
radius is shown with the shaded region on the right-bottom side of the figure
by adopting the constraint on $K_{\tau}$ from experiments at RCNP, i.e.,
$K_{\tau}=-550\pm 100$ MeV [51], together with $L=60\pm 20$ MeV as a fiducial
value of $L$, which corresponds to $\eta_{\tau}=59.9-113$ MeV. We also show
the allowed region with the empirical relation as a function of $\eta$,
assuming the fiducial values of $L=60\pm 20$ MeV and $K_{0}=240\pm 20$ MeV.
For reference, we show the constraints obtained from the astronomical
observations, i.e., PSR J0030+0451 and MSP J0740+6620 with NICER, the shaded
region with the observation of x-ray bursts, and $1.4M_{\odot}$ neutron star
radius from GW170817 (see text for the details), together with the theoretical
mass-radius curves constructed with some EOSs. The shaded region on the left-
top side of the figure denotes the excluded region from the causality.

###  IV.2 Constraints from terrestrial experiments

Since the neutron star EOS cannot characterize only by the nuclear saturation
parameters in a higher density region, the neutron star mass and radius
discussed with the nuclear saturation parameters constrained from terrestrial
experiments are only for low-mass neutron stars, which are located at the
right-bottom side on the neutron star mass and radius relation. Anyway, to
connect the neutron star mass and radius with the nuclear saturation
parameters, the empirical relations for neutron star mass and radius (or
gravitational redshift) must be useful. With the empirical relation for the
neutron star mass and gravitational redshift as a function of $\eta$ [34], we
have already discussed the expected region in the neutron star mass and radius
relation [2], e.g., the allowed region with the fiducial values of $L$ and
$K_{0}$, i.e., $L=60\pm 20$ and $K_{0}=240\pm 20$ MeV, is shown as in Fig. 6,
where the central density is considered up to twice the saturation density. In
a similar way, with the empirical relations derived in this study as a
function of $\eta_{\tau}$, the expected region in the neutron star mass and
radius relation can be shown with the shaded region in Fig. 6, assuming that
$L=60\pm 20$ and $K_{\tau}=-550\pm 100$ MeV, where the central density is
considered up to threefold the saturation density and the stellar models with
twice the saturation density are also denoted with the dotted line in the
shaded region.

Figure 7: The expected region of the neutron star mass and radius, using the
empirical relations as a function of $\eta_{\tau}$. The shaded region
corresponds to the constraint on neutron star mass and radius with $L=60\pm
20$ and $K_{\tau}=-550\pm 100$ MeV, which is the same as the shaded region
shown in Fig. 6. In the left panel, the expected regions are additionally
shown for $L=60\pm 30$, $L=60\pm 10$, $L=60\pm 5$, and $L=60$ MeV together
with $K\tau=-550\pm 100$ MeV. In the right panel, the region enclosed with the
dashed, solid, and dotted lines correspond to the region expected with
$L=60\pm 20$ and $K_{\tau}=-550$ MeV, $L=60$ and $K_{\tau}=-550\pm 100$ MeV,
and $L=60$ and $-348\leq K_{\tau}\leq-237$ MeV respectively.

Finally, in Fig. 7, we show the expected region of the neutron star mass and
radius, adopting the specific value of $L$ and $K_{\tau}$, e.g., the shaded
region for $L=60\pm 20$ and $K_{\tau}=-550\pm 100$ MeV as in Fig. 6. In
addition, in the left panel, the expected regions are shown for $L=60\pm 30$,
$L=60\pm 10$, $L=60\pm 5$, and $L=60$ MeV with $K_{\tau}=-550\pm 100$ MeV,
while in the right panel, the region enclosed with the dashed line, the solid
line, and the dotted line correspond to the expected mass and radius for
$L=60\pm 20$ and $K_{\tau}=-550$ MeV; $L=60$ and $K_{\tau}=-550\pm 100$ MeV;
and $L=60$ and $-348\leq K_{\tau}\leq-237$ MeV, respectively. From this
figure, one can obviously see that the $L$ dependence is much stronger than
the $K_{\tau}$ dependence for determining the neutron star properties, but
still $K_{\tau}$ also plays an important role in neutron star mass and radius
relation.

##  V Conclusion

To discuss the neutron star mass and radius with the constraints on the
nuclear saturation parameters obtained from terrestrial experiments, the
empirical relations for the neutron star mass and radius (or gravitational
redshift) as a function of nuclear saturation parameters must be useful, if
they exist. In this study, we find a suitable combination, $\eta_{\tau}$, of
the nuclear saturation parameters for expressing neutron star mass and
gravitational redshift and derive their empirical relations, focusing on the
saturation parameters for asymmetric nuclear matter. With the resultant
empirical relations, which are applicable up to threefold the saturation
density of symmetric nuclear matter, one can evaluate the neutron star mass
and gravitational redshift within $\sim 10\%$ accuracy, and the radius within
a few % accuracies. In addition, with these empirical relations, we discuss
the expected region in the neutron star mass and radius relation, assuming the
experimental values of saturation parameters, together with the constraints
from the astronomical observations. Furthermore, we find a tight correlation
between $\eta_{\tau}$ and $\eta$, which is another suitable combination of
saturation parameters for symmetric nuclear matter. With this correlation, we
derive the constraint on the isospin dependence of incompressibility for
asymmetric nuclear matter, $K_{\tau}$, as $-348\leq K_{\tau}\leq-237$ MeV,
assuming the fiducial value of the density-dependent nuclear symmetry energy
and the incompressibility for symmetric nuclear matter. The value of
$K_{\tau}$ derived in this study shows the discrepancy from the values deduced
from experiments. This may be because the experiments have been done using
finite nuclei. As pointed out in Ref. [54], the surface term or the mass
number dependence especially in the expansion of the incompressibility in
finite nuclei should be evaluated in such a way as to appropriately deduce the
incompressibility and its isospin dependence of the nuclear matter. Systematic
measurement of the incompressibility in finite nuclei in a wide mass-number
range may provide an opportunity to evaluate such a finite system effect. In
addition, we found $\eta$ and $\eta_{\tau}$ through trial and error. Since
these are good quantities to express the low-mass neutron stars and we also
find a strong correlation between $\eta$ and $\eta_{\tau}$, there may exist a
physical meaning behind these quantities. If one will find such a physical
meaning, it must tell us a new aspect of nuclear physics, which has never been
known up to now.

###### Acknowledgements.

We are grateful to Hajime Togashi for preparing the EOS data adopted in this
study and also to David Blaschke and Tomoya Naito for valuable comments. This
work is supported in part by Japan Society for the Promotion of Science (JSPS)
KAKENHI Grant Numbers JP19KK0354, and JP21H01088, and by Pioneering Program of
RIKEN for Evolution of Matter in the Universe (r-EMU).

## References

  * [1] S. L. Shapiro and S. A. Teukolsky, Black Holes, White Dwarfs, and Neutron Stars: The Physics of Compact Objects (Wiley-Interscience, New York, 1983).
  * [2] H. Sotani, N. Nishimura, and T. Naito, Prog. Theor. Exp. Phys. 2022, 041D01 (2022).
  * [3] P. Demorest, T. Pennucci, S. Ransom, M. Roberts, and J. Hessels, Nature 467, 1081 (2010).
  * [4] J. Antoniadis et al., Science 340, 6131 (2013).
  * [5] H. T. Cromartie et al., Nature Astronomy 4, 72 (2020).
  * [6] E. Fonseca et al., Astrophys. J. 915, L12 (2021).
  * [7] B. P. Abbott et al. (The LIGO Scientific Collaboration and the Virgo Collaboration), Phys. Rev. Lett. 119, 161101 (2017).
  * [8] E. Annala, T. Gorda, A. Kurkela, and A. Vuorinen, Phys. Rev. Lett. 120, 172703 (2018).
  * [9] C. D. Capano, I. Tews, S. M. Brown, B. Margalit, S. De, S. Kumar, D. A. Brown, B. Krishnan, and S. Reddy, Nature Astron. 4, 625 (2020).
  * [10] T. Dietrich, M. W. Coughlin, P. T. H. Pang, M. Bulla, J. Heinzel, L. Issa, I. Tews, and S. Antier, Science 370, 1450 (2020).
  * [11] K. R. Pechenick, C. Ftaclas, and J. M. Cohen, Astrophys. J. 274, 846 (1983).
  * [12] D. A. Leahy and L. Li, Mon. Not. R. Astron. Soc. 277, 1177 (1995).
  * [13] J. Poutanen and M. Gierlinski, Mon. Not. R. Astron. Soc. 343, 1301 (2003).
  * [14] D. Psaltis and F. Özel, Astrophys. J. 792, 87 (2014).
  * [15] H. Sotani and U. Miyamoto, Phys. Rev. D 98, 044017 (2018); 98, 103019 (2018).
  * [16] H. Sotani, Phys. Rev. D 101, 063013 (2020).
  * [17] T. E. Riley et al., Astrophys. J. 887, L21 (2019).
  * [18] M. C. Miller et al., Astrophys. J. 887, L24 (2019).
  * [19] T. E. Riley et al., Astrophys. J. 918, L27 (2021).
  * [20] M. C. Miller et al., Astrophys. J. 918, L28 (2021).
  * [21] N. Andersson and K. D. Kokkotas, Phys. Rev. Lett. 77, 4134 (1996).
  * [22] N. Andersson and K. D. Kokkotas, Mon. Not. R. Astron. Soc. 299, 1059 (1998).
  * [23] H. Sotani, K. Tominaga, and K. I. Maeda, Phys. Rev. D 65, 024010 (2001).
  * [24] H. Sotani and T. Harada, Phys. Rev. D 68, 024019 (2003); H. Sotani, K. Kohri, and T. Harada, ibid. 69, 084008 (2004).
  * [25] H. Sotani, N. Yasutake, T. Maruyama, and T. Tatsumi, Phys. Rev. D 83 024014 (2011).
  * [26] A. Passamonti and N. Andersson, Mon. Not. R. Astron. Soc. 419, 638 (2012).
  * [27] D. D. Doneva, E. Gaertig, K. D. Kokkotas, and C. Krüger, Phys. Rev. D 88, 044052 (2013).
  * [28] H. Sotani, Phys. Rev. D 102, 063023 (2020); 103021 (2020); 103, 123015 (2021).
  * [29] H. Sotani and A. Dohi, Phys. Rev. D 105, 023007 (2022).
  * [30] X. Viñas, M. Centelles, X. Roca-Maza, and M. Warda, Eur. Phys. J. A 50, 27 (2014).
  * [31] B.-A. Li, P. G. Krastev, D.-H. Wen, and N.-B. Zhang, Eur. Phys. J. A 55, 117 (2019).
  * [32] J. Estee et al. (S$\pi$RIT), Phys. Rev. Lett. 126, 162701 (2021).
  * [33] B. T. Reed, F. J. Fattoyev, C. J. Horowitz, and J. Piekarewicz, Phys. Rev. Lett. 126, 172503 (2021).
  * [34] H. Sotani, K. Iida, K. Oyamatsu, and A. Ohnishi, Prog. Theor. Exp. Phys. 2014, 051E01 (2014).
  * [35] M. Fortin, C. Providência, A. R. Raduta, F. Gulminelli, J. L. Zdunik, P. Haensel, and M. Bejger, Phys. Rev. C 94, 035804 (2016).
  * [36] K. Oyamatsu and K. Iida, Prog. Theor. Phys. 109, 631 (2003).
  * [37] K. Oyamatsu and K. Iida, Phys. Rev. C 75, 015801 (2007).
  * [38] B. K. Agrawal, S. Shlomo, and V. K. Au, Phys. Rev. C 72, 014310 (2005).
  * [39] F. Douchin and P. Haensel, Astron. Astrophys. 380, 151 (2001).
  * [40] E. Chabanat, Interactions effectives pour des conditions extremes d’isospin, Ph.D. thesis, University Claude Bernard Lyon-I (1995).
  * [41] H. S. Köhler, Nucl. Phys. A 258, 301 (1976).
  * [42] P. -G. Reinhard and H. Flocard, Nucl. Phys. A 584, 467 (1995).
  * [43] L. Bennour, P-H. Heenen, P. Bonche, J. Dobaczewski, and H. Flocard, Phys. Rev. C 40, 2834 (1989).
  * [44] M. Oertel, M. Hempel, T. Klähn, and S. Typel, Rev. Mod. Phys. 89, 015007 (2017).
  * [45] S. Shlomo, V. M. Kolomietz, and G. Colò, Eur. Phys. J. A 30, 23 (2006).
  * [46] J. R. Stone, N. J. Stone, and S. A. Moszkowski, Phys. Rev. C 89, 044316 (2014).
  * [47] Y. Wang, C. Guo, Q. Li, A. Le Fèvre, Y. Leifels, and W. Trautmann, Phys. Lett. B 778, 207 (2018).
  * [48] I. Vidaña, C. Providência, A. Polls, and A. Rios, Phys. Rev. C 80, 045806 (2009).
  * [49] U. Garg and G. Colò, Prog. Part. Nucl. Phys. 101, 55 (2018).
  * [50] X. Roca-Maza and N. Paar, Prog. Part. Nucl. Phys. 101, 96 (2018).
  * [51] T. Li et al. Phys. Rev. C 81, 034309 (2010).
  * [52] D. Patel et al. Phys. Lett. B 718, 447 (2012).
  * [53] P. Veselý, J. Toivanen, B. G. Carlsson, J. Dobaczewski, N. Michel, A. Pastore, Phys. Rev. C 86 024303 (2012).
  * [54] J. R. Stone, N. J. Stone, and S. A. Moszkowski, Phys. Rev. C 89, 044316 (2014).
  * [55] M. Centelles, X. Roca-Maza, X. Viñas, and M. Warda, Phys. Rev. Lett. 102, 122502 (2009).
  * [56] L.-W. Chen, B.-J. Cai, C. M. Ko, B.-A. Li, C. Shen, and J. Xu, Phys. Rev. C 80 014322 (2009).
  * [57] H. Sagawa, S. Yoshida, G.-M. Zeng, J.-Z. Gu, X.-Z. Zhang, Phys. Rev. C 76, 034327 (2007).
  * [58] M. Dutra, O. Lourenço, J. S. Sá Martins, A. Delfino, J. R. Stone, and P. D. Stevenson, Phys. Rev. C, 85, 035201 (2012).
  * [59] H. Sotani and H. Togashi, Phys. Rev. D 105, 063010 (2022).
  * [60] H. O. Silva, H. Sotani, and E. Berti, Mon. Not. R. Astron. Soc. 459, 4378 (2016).
  * [61] H. Sotani, Phys. Rev. C 95, 025802 (2017).
  * [62] H. Sotani and K. D. Kokkotas, Phys. Rev. D 95, 044032 (2017).
  * [63] A. W. Steiner, J. M. Lattimer, and E. F. Brown, Astrophys. J. Lett. 765, L5 (2013).
  * [64] J. M. Lattimer, Ann. Rev. Nucl. Part. Sci. 62, 485 (2012).


# Contingency Model Predictive Control for Automated Vehicles

John P. Alsterda1,2, Matthew Brown1 and J. Christian Gerdes1 1Department of
Mechanical Engineering, Stanford University, Stanford, CA 94305, USA
2Corresponding author<EMAIL_ADDRESS>
###### Abstract

We present Contingency Model Predictive Control (CMPC), a novel and
implementable control framework which tracks a desired path while
simultaneously maintaining a contingency plan – an alternate trajectory to
avert an identified potential emergency. In this way, CMPC anticipates events
that might take place, instead of reacting when emergencies occur. We
accomplish this by adding an additional prediction horizon in parallel to the
classical receding MPC horizon. The contingency horizon is constrained to
maintain a feasible avoidance solution; as such, CMPC is selectively robust to
this emergency while tracking the desired path as closely as possible. After
defining the framework mathematically, we demonstrate its effectiveness
experimentally by comparing its performance to a state-of-the-art
deterministic MPC. The controllers drive an automated research platform
through a left-hand turn which may be covered by ice. Contingency MPC prepares
for the potential loss of friction by purposefully and intuitively deviating
from the prescribed path to approach the turn more conservatively; this
deviation significantly mitigates the consequence of encountering ice.

##  I INTRODUCTION

In a future when automated vehicles are ubiquitous, we will expect them to
handle the surprises and challenges that humans tackle routinely. They must
avoid children who may run into the street, or maneuver safely when
encountering an icy surface. In control engineering, this concept is known as
robustness: a controller’s ability to safely operate when subject to
unexpected variation in system or environmental parameters. In the field of
Model Predictive Control (MPC), work to this end falls largely into two broad
categories: Robust MPC (RMPC) and Stochastic MPC (SMPC).

In two reviews, Mayne described some RMPC strategies, their limitations, and
direction for future development [1][2]. In these robust algorithms, uncertain
parameters’ values are first confined to a set between upper and lower bounds.
A control trajectory is then calculated which satisfies system constraints for
all possible combinations and realizations of these parameters. Thus RMPC is
sometimes dubbed ‘worst case’ MPC, prepared for the most extreme coincidences
of bad luck – however unlikely.

RMPC is often cast as a min-max problem, solving for commands that minimize
the maximum possible cost [3]. The problem is in general nonlinear and non-
convex, but Sartipizadeh et al. applied this technique on an approximate
convex hull of uncertainty to achieve computational feasibility [4]. D. He et
al. also ensured the tractability of their RMPC work by employing ‘quasi-min-
max’ MPC [5].

Figure 1: Volkswagen GTI Automated Test Platform

Another RMPC instantiation is Tube MPC, which reduces pessimism by solving for
control policies instead of open-loop control sequences [1]. Closed-loop
policies can react to future events and therefore pursue higher performance in
the near term. Tube MPC attempts to guarantee the set of possible disturbed
states is contained within problem constraints [6]. Kögel et al. presented a
formulation to further reduce conservatism by combining the effects of
estimation error and model uncertainty [7].

Robust MPC remains a principled strategy for systems that require guaranteed
stability, like a nuclear reactor [8]. However, it can be impracticably
cautious when uncertainties are numerous and large, as is the case driving on
legal roads. In this article, a selectively robust model predictive controller
is introduced which achieves responsible but practical conservatism in two
ways: First, contingency scenarios package uncertainties that are likely to
occur together. Whereas traditional RMPC assumes all uncertainties may agitate
simultaneously, engineers can leverage expert knowledge or parameter
covariance to focus on the most likely emergency events. Second, we separate
nominal and contingency planning into separate horizons, and define a unique
cost function for each. By separating nominal and emergency costs, the
contingency trajectory is unburdened by objectives typically found in MPC cost
functions. Whereas traditional RMPC would aim to maintain comfortable and
efficient emergency maneuvers, Contingency MPC preserves only a feasible
avoidance trajectory.

Stochastic MPC is another MPC technique to handle uncertainty. Control
engineers seeking robustness with less pessimism favor building parameter
ambiguity and stochastic uncertainty directly in their optimizations [9]. In
SMPC, uncertain quantities are characterized into statistical distributions,
and control inputs are calculated to minimize expected cost. In this way, the
stochastic algorithms do not place undue focus on highly improbable scenarios.
Farina et al. characterized the general problem statement, which requires a
disturbance model to be defined either by a probability density function (e.g.
the Gaussian) or some function which produces samples [10]. This assumption,
however, may not hold for some systems or applications.

For example, a well-defined distribution may not always exist for the
coefficient of friction along a roadway. We may not feel comfortable assigning
a probability density to predict a child’s behavior. On the contrary, if ice
is likely to exist somewhere on a winter day, we should drive in such a way
that a surprise icy patch can be handled safely. If children are seen playing
on the sidewalk, we should anticipate the possibility they run into the street
and approach with an avoidance strategy – a contingency plan.

One SMPC technique has been borrowed to formulate Contingency MPC. In Scenario
SMPC, a prediction horizon tree is built which branches from stage-nodes by
sampling from a parameter distribution. At each stage-node only one command is
calculated, an important constraint to prevent Scenario SMPC from improper
prognostication. Lucia et al. demonstrated this technique with a nonlinear
model, and Krishnamoorthy et al. worked to improve the distributed systems
over which Scenario SMPC can be decomposed [11] [12]. Contingency MPC’s
prediction horizon structure resembles a simple scenario tree, with a single
stage-node at $x_{0}$. At this stage, we also employ the single-command
constraint, shown in Fig. 2. This constraint is central to the function of
CMPC; the input $u^{0}$ is calculated to be valid for both nominal and
contingency horizons. Thus the vehicle will track the desired path as closely
as possible, while maintaining the ability to choose the contingency
trajectory.

This article offers a new option for control engineers who seek robustness
from foreseeable but uncertain events. In Section II, the general Contingency
MPC problem is defined. Section III narrows the general problem to an linear
convex program for use in automated vehicle (AV) control. In Section IV, we
present an experimental demonstration of CMPC. First the test platform turns
through a snow-covered corner; it then encounters ice performing the same
maneuver. Finally, Section V explores the potential for CMPC to tackle other
problems in control under uncertainty and outlines future steps to develop the
framework.

##  II CMPC FORMULATION

Model Predictive Control is a receding horizon optimal control technique whose
central strength is its ability to minimize future cost while explicitly
respecting model and environmental constraints [13]. At each time step, MPC
calculates a prediction horizon (open-loop trajectory) of discrete states $x$
and inputs $u$ to optimize a cost function, subject to constraints. Upon
computation, $u^{0}$ is deployed and a new horizon is optimized. Contingency
MPC’s receding horizons are illustrated in Fig. 2, where $x^{0}$ represents
the plant’s current state and two trajectories are calculated. Trajectory
$x_{nom}$ is analogous to the traditional deterministic MPC horizon,
containing the states we intend to guide the vehicle through by application of
control inputs $u_{nom}$. On the other hand, trajectory $x_{c}$ represents our
contingency plan of action, subject to a unique set of constraints and costs.

Figure 2: Contingency MPC prediction horizon. A nominal and emergency
trajectory are computed, which share a zeroth action.

###  II-A General Case

The following optimization problem calculates N control inputs
$\textbf{u}\in\mathbb{R}^{2\cdot m}$ to guide state vectors
$\textbf{x}\in\mathbb{R}^{2\cdot n}$ to minimize cost function $J$ while
abiding to equality constraints $H$ and inequality constraints $G$. N is the
number of prediction horizon stages, n the length of the nominal state vector,
and m the number of control inputs.

$J(x^{0})=\min_{\textbf{u}}\sum_{k=0}^{N}j^{k}(\textbf{x}^{k},\textbf{u}^{k})$
(1a)

where:

$\textbf{x}^{k}=\begin{bmatrix}[c]x_{nom}\\\ x_{c}\end{bmatrix}^{k}\hskip
14.22636pt;\hskip 14.22636pt\textbf{u}^{k}=\begin{bmatrix}[c]u_{nom}\\\
u_{c}\end{bmatrix}^{k}$ (1b)

such that:

$h(\textbf{x},\textbf{u})=0\hskip 11.38092pt\forall\hskip 11.38092pth\in H$
(1c)

$g(\textbf{x},\textbf{u})\leq 0\hskip 11.38092pt\forall\hskip 11.38092ptg\in
G$ (1d)

and for $k=0$: $\hskip 17.07182ptu^{0}_{nom}=u^{0}_{c}=u^{0}\hskip 68.28644pt$
(1e)

Equality constraints $H$ include the dynamics model:

$\textbf{x}^{k+1}=f(\textbf{x}^{k},\textbf{u}^{k})\hskip
14.22636pt\forall\hskip 14.22636ptk$ (1f)

As in typical MPC formulations, the cost function terms $j$ can be used to
encourage the nominal states $x_{nom}$ to follow a desired trajectory,
maintain comfort or smoothness, minimize fuel use, etc. However, these costs
need not be levied against the contingency trajectory. The separation of costs
allows the contingency trajectory to focus on maintaining a feasible emergency
maneuver, unburdened by costs unimportant during crisis. Constraint sets $H$
and $G$ encode the contingency scenario. For example, constraints on $x_{c}$
can be added to represent an obstacle which may enter the road. Alternatively,
to prepare for a dynamics-based contingency, the propagation model for
$x_{c}^{k+1}$ can be modified from the nominal case.

The resulting optimized horizon contains two input trajectories from $x^{0}$:
$u_{nom}$ and $u_{c}$. As in classical MPC, $u^{0}$ is deployed – the root
command for both $u_{nom}$ and $u_{c}$. Therefore the vehicle never chooses
between the nominal and contingency trajectories. $u_{0}$ tracks the desired
path to the greatest extent possible, while maintaining an avoidance maneuver
for the contingent event. If the contingency event occurs, CMPC is expected to
return safe trajectories in subsequent time-steps. Alternatively if the
contingency ceases to be relevant, the vehicle can proceed along the desired
path.

##  III Linear CMPC for AV Control

The general CMPC formulation above can be adapted to any MPC framework. For
the remainder of this article, the general form is narrowed to a convex linear
optimization extended from the deterministic MPC controller developed by Brown
et al., and modified to solve directly for steering angle $\delta$ as
presented by Zhang et al. [14] [15].

###  III-A Vehicle Dynamics Model

Figure 3: Bicycle model diagram

To adequately handle contingency scenarios which push the limits of vehicle
handling, dynamics must be modeled with an appropriate degree of fidelity. We
accomplish this using a planar bicycle model with successively linearized tire
forces to represent the vehicle [16]. Fig. 3 illustrates this model with three
position and velocity states, and forces at each axle. The position states are
local to a path with curvature $\kappa$, where $s$ represents the vehicle’s
longitudinal progress, $e$ its lateral error from the path, and $\Delta\psi$
is the heading angle error. Velocity states are composed of longitudinal speed
$U_{x}$, lateral speed $U_{y}$, and yaw rate $r$. The following differential
equations govern the states’ evolution:

$\dot{s}=U_{x}-U_{y}\Delta\psi\hskip 17.07182pt$ (2a)
$\dot{e}=U_{y}+U_{x}\Delta\psi\hskip 17.07182pt$ (2b)
$\dot{\Delta\psi}=r-\kappa U_{x}\hskip 43.53267pt$ (2c)
$\dot{U_{x}}=\frac{F_{xf}+F_{xr}}{m}+rU_{y}$ (2d)
$\dot{U_{y}}=\frac{F_{yf}+F_{yr}}{m}-rU_{x}$ (2e)
$\dot{r}=\frac{aF_{yf}-bF_{yr}}{I_{z}}\hskip 12.23447pt$ (2f)

Forces $F_{yf}$ and $F_{yr}$ are modeled by a Fiala brush tire curve which
relates lateral force to tire slip angles [17]:

$\delta+\alpha_{f}=\tan^{-1}\left(\frac{U_{y}+ar}{U_{x}}\right)\hskip
16.78702pt$ (3a)

$\alpha_{r}=\tan^{-1}\left(\frac{U_{y}-br}{U_{x}}\right)$ (3b)

In Fig. 4, tire models are plotted for asphalt, snow, and ice conditions,
which are the coefficients $\mu$ used in the experiments which follow.

Figure 4: Fiala brush tire models for asphalt, snow, and ice conditions

As a first demonstration of CMPC, this formulation calculates only the
steering angle $\delta$. Throttle and braking commands are calculated upstream
by a feedforward-feedback scheme, as developed by Funke et al. [18]. Including
longitudinal forces into CMPC is an opportunity for future development. Thus
the CMPC problem we formulate sees $F_{xf},F_{xr}$, $U_{x}$, and $s$ as a
known profile along the prediction horizon. The MPC state vector is then
$x=[U_{y}\enspace r\enspace\Delta\psi\enspace e]^{T}$.

To complete the MPC dynamics model, the equations of motion and tire model are
linearized with respect to operating points from the previous iteration’s
solution. The linear model was then discretized with respect to $\delta$: A
zero-order hold was used for the first ten 20 ms steps, and a first-order hold
was used for forty 300 ms steps to form a twelve second horizon. The number
and length of time steps were chosen to capture high frequency dynamics in the
short term, while extending the total horizon such that CMPC can see far
enough to prepare avoidance maneuvers. These choices were also influenced by
computation time, as reported in Section IV. After discretization, the
resulting affine dynamics equations take the form [15]:

$x^{k+1}=A^{k}x^{k}+B^{k}\delta^{k}+C^{k}$ (4a)

To support both nominal and contingency prediction horizons, (4a) is
duplicated for both $x_{nom}$ and $x_{c}$ state vectors:

$\displaystyle\textbf{x}^{k+1}=$ $\displaystyle\begin{bmatrix}[c]x_{nom}\\\
x_{c}\end{bmatrix}^{k+1}$ (4b) $\displaystyle=$
$\displaystyle\begin{bmatrix}[c]A_{nom}&0\;\\\
0&A_{c}\;\end{bmatrix}^{k}\\!\\!\textbf{x}^{k}+\begin{bmatrix}[c]B_{nom}&0\;\\\
0&B_{c}\;\end{bmatrix}^{k}\\!\\!\textbf{u}^{k}+\begin{bmatrix}[c]C_{nom}\\\
C_{c}\end{bmatrix}^{k}$

###  III-B Linear CMPC Problem Statement

The following optimization is extended from the deterministic MPC presented in
[14] to calculate a smooth trajectory which follows a desired path while
adhering to dynamics, stability, and environmental constraints.

$\displaystyle\min_{\textbf{u}}\;\;(\textbf{x}^{50})^{T}\begin{bmatrix}[c]\,0&0\\\
\,0&Q\end{bmatrix}\textbf{x}^{50}+\;\sum_{k=0}^{50}\,(\textbf{x}^{k})^{T}\begin{bmatrix}[c]Q&0\,\\\
0&0\,\end{bmatrix}$ $\displaystyle\textbf{x}^{k}$ (5a)
$\displaystyle+\;(\textbf{v}^{k})^{T}\begin{bmatrix}[c]R&0\,\\\
0&0\,\end{bmatrix}\textbf{v}^{k}+\,W\begin{bmatrix}[c]\,\sigma^{k}_{stab}\,\\\
\,\sigma^{k}_{env}\,\end{bmatrix}$

where:

$\displaystyle Q=\begin{bmatrix}[c]\,0&0&0&0\,\\\ \,0&0&0&0\,\\\
\,0&0&1&0\,\\\ \,0&0&0&1\,\\\ \end{bmatrix}$ $\displaystyle;\quad R=0.01$ (5b)
$\displaystyle\textbf{v}^{k}=\begin{bmatrix}[c]\,(u^{k}-u^{k-1})_{nom}\\\
\,(u^{k}-u^{k-1})_{c\hskip 14.22636pt}\\\ \end{bmatrix}$ $\displaystyle;\quad
W=\begin{bmatrix}[c]\;50&0\\\ \;0&500\\\ \end{bmatrix}$

$Q$ applies cost to heading and lateral state errors, $R$ penalizes change in
steering angle (slew rate), and $W$ discourages stability and environmental
constraint violation. Weights were tuned for experimental performance and
chosen to establish priority among the problem’s goals: Most important is
driving the constraints’ slack variables to zero to ensure stability and
obstacle avoidance. Next, the nominal trajectory should adhere to the desired
path. Lastly the nominal input trajectory should be smooth and comfortable for
passengers. This minimization is subject to (4b) and the following
constraints. First the steering angle and its slew rate are bounded:

$|\textbf{u}^{k}|\leq\begin{bmatrix}[c]\delta_{max}\\\
\delta_{max}\end{bmatrix}\quad\,;\,\quad|\textbf{v}^{k}|\leq\begin{bmatrix}[c]v_{max}\\\
v_{max}\end{bmatrix}\qquad\forall\quad k$ (5c)

Next, yaw rate and sideslip are confined to an invariant set, as developed by
Bobier et al. and Beal et al. [19] [20]. This stability envelope is
illustrated in Fig. 7 and enforces ${|r|\,\leq\,\mu g/U_{x}}$ and
$|\beta|\,\leq\,\alpha_{rear,peak}+br/U_{x}\,$:

$\begin{bmatrix}[c]H_{nom}&0\;\\\
0&H_{c}\;\end{bmatrix}_{stab}\\!\\!\\!\\!\\!\\!\textbf{x}^{k}\leq\begin{bmatrix}[c]G_{nom}\\\
G_{c}\end{bmatrix}_{stab}\\!\\!\\!\\!\\!\\!+\sigma_{stab}\quad\forall\quad k$
(5d)

To prevent collision with road edges, an environmental envelope encodes the
inequality $e_{min}\leq e\leq e_{max}$:

$\begin{bmatrix}[c]H_{nom}&0\;\\\
0&H_{c}\;\end{bmatrix}_{env}\\!\\!\\!\\!\\!\textbf{x}^{k}\leq\begin{bmatrix}[c]G_{nom}\\\
G_{c}\end{bmatrix}_{env}\\!\\!\\!\\!\\!+\sigma_{env}\quad\forall\quad k$ (5e)

Finally the nominal and contingency trajectories are coupled via their first
commands:

$u^{0}_{nom}=u^{0}_{c}$ (5f)

In the cost function and inequality constraints, block diagonal matrices are
employed to weight and constrain the nominal and contingency trajectories
independently. For example, the terminal cost on $\textbf{x}^{50}$ applies
penalty only to the contingency trajectory’s final angular and lateral
deviation, an important feature to mitigate the shortcomings of a finite
horizon. No path following or input-smoothing penalty is applied to the
contingency trajectory for $k<N$, however, because we only require it to obey
the imposed constraints.

##  IV Experimental Results and Analysis

Figure 5: Deterministic MPC and Contingency MPC experimentally navigated a
snow-covered road. Closed-loop behavior in green; nominal horizon in blue;
contingency horizon in red. CMPC took a more conservative route. For clarity,
Only ten of fifty horizon states are plotted from each trajectory.

Contingency MPC was demonstrated in experiment using the convex linear
formulation presented in Section III, deployed on an automated Volkswagen GTI
(Fig. 1). These results compare CMPC with a deterministic MPC based on the
nominal costs and constraints. In each comparison, CMPC purposefully deviated
from the desired path to maintain an avoidance strategy for the anticipated
contingency.

Figure 6: Deterministic MPC and Contingency MPC experimentally navigated an
ice-covered left turn. DMPC did not anticipate ice and lost control. CMPC
subtly deviated from the nominal path to increase robustness.

The first setting is shown in Fig. 5, a snow covered left-hand turn driven at
5 m/s – a feasible speed for the roughly homogeneous surface. Deterministic
MPC completed the maneuver while adhering closely to the desired path (the
lane’s centerline). The nominal dynamics modeled snow with a coefficient of
friction $\mu=0.25$. Next, Contingency MPC completed the same maneuver while
anticipating the potential for ice ($\mu=0.10$) to cover the road. The red
trajectory illustrates the contingency plan that CMPC computed for this time-
step. As the vehicle proceeded under CMPC, it deviated from the desired path
to keep the contingency trajectory within path boundaries. It shifted wide,
then steered early to make the turn more gentle and conservative. After
completing the turn, CMPC returned to the path. In this case the vehicle did
not encounter ice and the minor deviation was ultimately unnecessary.

In Fig. 6, the vehicle turned left with the same geometry onto a polished ice
surface (cyan). Deterministic MPC’s performance is shown in the figure’s first
pane. Unaware of the ice ahead, it naively followed the centerline and could
not exert the lateral force necessary to track the path. It ran through the
boundary (grey), and recovered after sliding beyond the ice onto snow where it
regained traction. After encountering the ice, no avoidance maneuver could
have ensured safety. To escape, the vehicle must have behaved differently
before the emergency.

In panes 2 and 3 of Fig. 6, Contingency MPC demonstrated its effectiveness by
repeating the Fig. 5 approach and intuitively deviating from the desired path
to design a safer trajectory. The closed-loop behavior resembles the outside-
inside-outside cornering technique described by professional racecar driver
Taruffi, to maximize speed by minimizing curvature [21]. In this experiment
the goal was not to race; however the vehicle’s speed profile has been
designed to push the scenario’s limit. In order to cope with this speed (set a
priori), CMPC calculated steering commands to increase the radius of curvature
and decrease the tire-force necessary to stay inside the boundary. Thus, the
consequence of encountering the extremely slick zamboni-polished surface was
significantly mitigated by CMPC, compared to deterministic MPC. In closed-loop
(green), CMPC drifted only 20 cm beyond the road edge as it exited the turn.

Figure 7: CMPC prediction horizons and closed-loop trajectory on the $r$ –
$\beta$ phase plane. Stability envelopes for snow (blue) and ice (red) were
violated in closed-loop to prioritize staying within the road boundary.

In Fig. 7, velocity state data is presented from the CMPC test on ice,
corresponding to the same time-step as Fig. 6, pane 2. The parallelograms
illustrate the stability envelopes from (5d), inside which CMPC attempted to
confine the vehicle’s state. Measured velocity states are plotted in green,
which show the car’s actual closed-loop behavior differed significantly from
either plan. As it proceeded through the corner, CMPC evaluated the prospect
of breaching the stability and environmental envelopes. According to the cost
function weights defined in (5b), CMPC allowed a greater violation of the
stability guarantee as it prioritized avoidance of the road edge. Setting
$\sigma_{env}$ ten times greater than $\sigma_{stab}$ influenced CMPC to trade
the possibility of losing stability for a minimal collision.

The optimization for this experiment was solved using FORCES Pro on a single
i7 processor [22]. The fifty stage deterministic MPC problem converged in an
average of 4.5 ms, while Contingency MPC extended average solution times to 11
ms and up to 16 ms in some iterations. This experiment exhibited Contingency
MPC’s use in a dynamics-based scenario. The prospect of ice on the roadway was
encoded into the contingency horizon, which caused controller to act robustly
in anticipation of decreased friction.

##  V Discussion

The snow and ice experiments successfully demonstrated the concept and
efficacy of Contingency MPC. Furthermore, the extreme polished-ice
circumstance pushed this convex-linear implementation to its limit, providing
insight into how CMPC and similar controllers may be improved in the future.
One such limitation was the controller’s focus on steering. Indeed, the
simplest solution to the presented situation may have been to slow down.
Future work could incorporate longitudinal and lateral inputs to CMPC by
moving to a nonlinear optimization or by decoupling $F_{x}$ and $F_{y}$.

Furthermore, linearizing the low friction tire model presented difficulty and
contributed to the deviation illustrated in Fig. 7. As shown in Fig. 4, the
generous linear region present on asphalt nearly vanishes as friction
decreases. As a result, linearizations of the low friction model lose accuracy
quickly as one moves away from the operating point, allowing some CMPC
iterations to predict greater force than was physically available. This
problem may be abated by a strategy which maintains greater fidelity of the
nonlinear model.

##  VI CONCLUSION

In this article, Contingency MPC has been established as a credible strategy
to augment a deterministic model predictive controller with robustness. In
systems where potential emergency situations can be identified, CMPC maintains
an avoidance trajectory while following a desired path. Experimentally, the
controller successfully mitigated the effect of an abrupt decrease in road
surface friction by consciously sacrificing path tracking performance to
increase robustness. It leveraged the knowledge that ice may appear somewhere
in the environment, and thus planned trajectories which were prepared for its
occurrence.

Several avenues appear promising for future development: Applying the
framework to additional scenarios may shed valuable insight into the character
and range of CMPC. For example, pedestrians may cross unexpectedly or a
vehicle ahead may stop abruptly. To investigate its use within a larger
system, CMPC could be integrated with an online emergency recognition program
to identify contingencies in real-time. In more complex situations, several
contingency events may be simultaneously probable. For these cases, CMPC could
be augmented to maintain two or more contingency horizons.

## ACKNOWLEDGMENT

Thank you Volkswagen Group Research and Engineering Research Lab for
experimental support, and Embotech GmbH for academic licensing. Alsterda is
supported by Renault Group and the U.S. Dept. of Veterans Affairs G.I. Bill.
Brown is supported by Renault Group and the Graduate Research Fellowship from
the National Science Foundation.

## References

  * [1] D. Q. Mayne, “Model predictive control: recent developments and future promise,” Automatica, vol. 50, no. 12, pp. 2967–2986, 2014.
  * [2] D. Q. Mayne, “Robust and stochastic model predictive control: are we going in the right direction?,” Annual Reviews in Control, vol. 41, pp. 184–192, 2016.
  * [3] P. J. Campo and M. Morari, “Robust model predictive control,” in Proc. American Control Conference, 1987.
  * [4] H. Sartipizadeh and T. L. Vincent, “Uncertainty characterization for robust MPC using an approximate convex hull method,” in Proc. American Control Conference, 2016.
  * [5] D. F. He, H. Huang, and Q. X. Chen, “Quasi-min-max MPC for constrained nonlinear systems with guaranteed input-to-state stability,” Journal of the Franklin Institute, vol. 351, no. 6, pp. 3405–3423, 2014.
  * [6] F. Blanchini, “Control synthesis for discrete time systems with control and state bounds in the presence of disturbances,” Journal of Optimization Theory and Applications, vol. 65, no. 1, p. 29–40, 1990.
  * [7] M. Kögel and R. Findeisen, “Robust output feedback MPC for uncertain linear systems with reduced conservatism,” in Proc. 20th IFAC World Congress, 2017.
  * [8] H. Eliasi, M. B. Menhaj, and H. Davilu, “Robust nonlinear model predictive control for a PWR nuclear power plant,” Progress in Nuclear Energy, vol. 54, no. 1, pp. 177–185, 2012.
  * [9] M. B. Saltik, L. Özkan, J. H. Ludlage, S. Weiland, and P. M. Van den Hof, “An outlook on robust model predictive control algorithms: reflections on performance and computational aspects,” Journal of Process Control, vol. 61, pp. 77–102, 2018.
  * [10] M. Farina, L. Giulioni, and R. Scattolini, “Stochastic linear model predictive control with chance constraints - a review,” Journal of Process Control, vol. 44, pp. 53–67, 2016.
  * [11] S. Lucia, T. Finkler, and S. Engell, “Multi-stage nonlinear model predictive control applied to a semi-batch polymerization reactor under uncertainty,” Journal of Process Control, 2013.
  * [12] D. Krishnamoorthy, E. Suwartadi, B. Foss, S. Skogestad, and J. Jaschke, “Improving scenario decomposition for multistage MPC using a sensitivity-based path-following algorithm,” IEEE Control Systems Letters, vol. 2, no. 4, pp. 581 – 586, 2018.
  * [13] J. Richalet, A. Rault, J. L. Testud, and J. Papon, “Model predictive heuristic control: applications to industrial processes,” Automatica, vol. 14, no. 5, pp. 413–428, 1978.
  * [14] M. Brown, J. Funke, S. Erlien, and J. C. Gerdes, “Safe driving envelopes for path tracking in autonomous vehicles,” Control Engineering Practice, vol. 61, pp. 307–316, 2017.
  * [15] V. Zhang, S. M. Thornton, and J. C. Gerdes, “Tire modeling to enable model predictive control of automated vehicles from standstill to the limits of handling,” in Proc. 14th International Symposium on Advanced Vehicle Control, 2018.
  * [16] S. M. Erlien, J. Funke, and J. C. Gerdes, “Incorporating non-linear tire dynamics into a convex approach to shared steering control,” in Proc. American Control Conference, 2014.
  * [17] H. Pacejka, Tire and Vehicle Dynamics.  Butterworth-Heinemann, 3rd ed., 2012.
  * [18] J. Funke, M. Brown, S. M. Erlien, and J. C. Gerdes, “Collision avoidance and stabilization for autonomous vehicles in emergency scenarios,” IEEE Transactions on Control Systems Technology, vol. 25, no. 4, pp. 1204 – 1216, 2017.
  * [19] C. G. Bobier, S. Joe, and J. C. Gerdes, “Sliding surface envelope control: keeping the vehicle within a safe state-space boundary,” in Proc. Dynamic Systems and Control Conference, 2011.
  * [20] C. E. Beal and J. C. Gerdes, “Model predictive control for vehicle stabilization at the limits of handling,” IEEE Transactions on Control Systems Technology, vol. 21, no. 4, pp. 1258 – 1269, 2013.
  * [21] P. Taruffi, The Technique of Motor Racing.  Bentley Publishers, 2003.
  * [22] A. Domahidi and J. Jerez, “FORCES Professional,” 2014.