# Isospin splitting of the Dirac mass probed by the relativistic Brueckner-
Hartree-Fock theory in the full Dirac space

Pianpian Qin  Department of Physics, Chongqing Key Laboratory for Strongly
Coupled Physics, Chongqing University, Chongqing 401331, China  Qiang Zhao
Center for Exotic Nuclear Studies, Institute for Basic Science,
Daejeon 34126, Korea  Hui Tong  Helmholtz-Institut f$\ddot{u}$r Strahlen- und
Kernphysik and Bethe Center for Theoretical Physics, Universit$\ddot{a}$t
Bonn, D-53115 Bonn, Germany  Chencan Wang  School of Physics and astronomy,
Sun Yat-Sen University, Zhuhai 519082, China  Sibo Wang<EMAIL_ADDRESS>Department of Physics, Chongqing Key Laboratory for Strongly Coupled Physics,
Chongqing University, Chongqing 401331, China

###### Abstract

The isospin splitting of the Dirac mass obtained with the relativistic
Brueckner-Hartree-Fock (RBHF) theory is thoroughly investigated. From the
perspective in the full Dirac space, the long-standing controversy between the
momentum-independence approximation (MIA) method and the projection method on
the isospin splitting of the Dirac mass in asymmetric nuclear matter (ANM) is
analyzed in detail. We find that, the assumption procedure of the MIA method,
which assumes that the single-particle potentials are momentum independent, is
not a sufficient condition that directly leads to the wrong sign of the
isospin splitting of the Dirac mass, while the extraction procedure of the MIA
method, which extracts the single-particle potentials from the single-particle
potential energy, leads to the wrong sign. By approximately solving the set of
equations involved in the extraction procedure, a formal expression of the
Dirac mass is obtained. The wrong isospin splitting of the Dirac mass is
mainly caused by that the extraction procedure forcely assumes the momentum
dependence of the single-particle potential energy to be a quadratic form
where the strength is solely determined by the constant scalar potential.

It is still a challenging task to understand the properties of nuclear many-
body system starting from the realistic nucleon-nucleon ($NN$) interaction.
The repulsive core of the realistic $NN$ interaction causes a strong
correlation for the many-body wave function and requires advanced many-body
methods that go beyond mean field Baldo and Maieron (2007); Barrett  _et al._
(2013); Hagen  _et al._ (2014); Carlson  _et al._ (2015); Hergert  _et al._
(2016); Shen  _et al._ (2019). The Brueckner-Hartree-Fock (BHF) theory
Brueckner  _et al._ (1954) is one of the representative nuclear many-body
methods, which is characteristic for its capacity to soften the realistic $NN$
interaction to an effective $G$ matrix in nuclear medium. The BHF theory can
be derived as the two-hole-line truncation to the general Bethe-Brueckner-
Goldstone expansion theory Day (1967); Baldo and Maieron (2007), where the
ground-state properties of the nuclear many-body system are calculated order
by order according to the number of independent hole lines contained in the
expansion diagrams Day (1967). Dating back to the 1960s, it is found that the
saturation points of symmetric nuclear matter (SNM) calculated by the BHF
theory with different two-body interactions are located on a so-called Coester
line Coester  _et al._ (1970), which deviates systematically from the
empirical values. Within the nonrelativistic framework, a way out of this
dilemma is to introduce three-body forces Grangé  _et al._ (1989); Baldo  _et
al._ (1997); Zuo  _et al._ (2002); Zhou  _et al._ (2004); Li  _et al._ (2006,
2008).

In 1980, major modification of the saturation properties of SNM is obtained by
including a relativistic description of the nucleon motion Anastasio  _et al._
(1980). With this pioneering work, much attention has been paid to developing
the relativistic Brueckner-Hartree-Fock (RBHF) theory Anastasio  _et al._
(1981); Brockmann and Machleidt (1984); Horowitz and Serot (1987); ter Haar
and Malfliet (1987). Starting from the Bonn potential Machleidt (1989), the
saturation points of SNM obtained by the RBHF theory have been shifted
remarkably close towards the empirical values, without introducing explicit
three-body forces. The success of the RBHF theory has been understood by the
fact that, relativistic effects contribute a particular part of the three-body
force Zuo  _et al._ (2002); Li  _et al._ (2006); Sammarruca  _et al._ (2012)
through virtual nucleon-antinucleon excitations in the intermediate states
(the so-called $Z$ diagrams) Jaroszewicz and Brodsky (1991); Brown (1987).

One of the essential points of the RBHF theory is using the Dirac equation to
describe the single-nucleon motion in the mean field, i.e., the so-called
single-particle potentials (SPPs). The SPP operator is generally divided into
a scalar and a vector component Serot and Walecka (1986). They can be obtained
in a self-consistent way from the effective $G$ matrix. In principle, the
calculations of SPPs and the $G$ matrix should be performed in the full Dirac
space, where both the positive-energy states (PESs) and negative-energy states
(NESs) are included. However, due to the complexity of these procedures, RBHF
calculations have been usually performed without NESs, by invoking approximate
methods to accomplish the Hartree-Fock calculation, such as the momentum-
independence approximation (MIA) method Brockmann and Machleidt (1990);
Sammarruca  _et al._ (2012); Li  _et al._ (1992); Alonso and Sammarruca
(2003); Tong  _et al._ (2018); Wang  _et al._ (2020) and the projection method
Gross-Boelting  _et al._ (1999); ter Haar and Malfliet (1987); van Dalen  _et
al._ (2004, 2005a, 2005b, 2007).

The MIA method assumes that the SPPs are independent of the momentum, and
extracts the SPPs from the single-particle potential energies calculated at
two selected momenta. In the projection method, the $G$ matrix elements are
projected onto a complete set of Lorentz invariant amplitudes Horowitz and
Serot (1987), from which the SPPs are calculated. Since the $G$ matrix coupled
to the NESs is not considered, the SPPs obtained with these two methods are
both suffered from ambiguities Nuppenau  _et al._ (1989); Gross-Boelting  _et
al._ (1999); Wang  _et al._ (2021). Moreover, two methods deviate
significantly for the isovector properties in asymmetric nuclear matter (ANM),
especially for the isospin splitting of the Dirac mass
$M^{*}_{D,n}-M^{*}_{D,p}$ Ulrych and Müther (1997); van Dalen  _et al._
(2005a); Van Giai  _et al._ (2010). The MIA method leads to
$M^{*}_{D,n}-M^{*}_{D,p}>0$, while the projection method obtains the opposite
sign $M^{*}_{D,n}-M^{*}_{D,p}<0$.

Recently, self-consistent RBHF calculations in the full Dirac space have been
achieved Wang  _et al._ (2021, 2022a, 2022b), which avoid the ambiguities
suffered from the RBHF calculations without NESs. In ANM, the full solution
predicts the sign $M^{*}_{D,n}-M^{*}_{D,p}<0$ Wang  _et al._ (2022b) and
clarifies the long-standing controversy between the MIA method and the
projection method on the isospin splitting of the Dirac mass. The RBHF theory
in the full Dirac space has also been applied to study the nonrelativistic
effective mass in nuclear matter Wang  _et al._ (2023), the properties of
208Pb with a liquid droplet model Tong  _et al._ (2023), and the neutron star
properties Tong  _et al._ (2022); Wang  _et al._ (2022c); Qu  _et al._ (2023).
In this work, we aim to thoroughly study the isospin splitting of the Dirac
mass obtained with the RBHF calculations, especially for the performance of
the MIA method, from the perspective in the full Dirac space.

In the RBHF theory, the single-particle motion for a nucleon inside the
infinite nuclear matter is described by the Dirac equation

$\left\\{\bm{\alpha}\cdot\bm{p}+\beta\left[M+\mathcal{U}_{\tau}(\bm{p})\right]\right\\}u_{\tau}(\bm{p},s)=E_{\bm{p},\tau}u_{\tau}(\bm{p},s).$
(1)

Here $\bm{\alpha}$ and $\beta$ are the Dirac matrices, $M$ is the nucleon
mass, $\bm{p}$ and $E_{\bm{p},\tau}$ are the momentum and single-particle
energy, $s$ denotes the spin, and $\tau=n,p$ denotes neutron $n$ and proton
$p$. The symbol $u_{\tau}$ in Eq. (1) represents the positive-energy Dirac
spinor, while the negative-energy Dirac spinor $v_{\tau}$ is obtained by
$v_{\tau}=\gamma_{5}u_{\tau}$. The SPP operator $\mathcal{U}_{\tau}(\bm{p})$
can be decomposed into the scalar potential $U_{S,\tau}(p)$, the timelike and
spacelike parts of the vector potential $U_{0,\tau}(p)$ and $U_{V,\tau}(p)$
Serot and Walecka (1986) as

$\mathcal{U}_{\tau}(\bm{p})=U_{S,\tau}(p)+\gamma^{0}U_{0,\tau}(p)+\bm{\gamma\cdot\hat{p}}U_{V,\tau}(p).$
(2)

Here $\hat{\bm{p}}=\bm{p}/|\bm{p}|=\bm{p}/p$ is the unit vector parallel to
the momentum $\bm{p}$. By calculating the matrix elements of
$\mathcal{U}_{\tau}(\bm{p})$ on the basis expanded by the PESs and NESs, i.e.,
$\Sigma^{++}_{\tau}(p)$, $\Sigma^{-+}_{\tau}(p)$, and $\Sigma^{--}_{\tau}(p)$,
the momentum dependent SPPs can be determined uniquely through Wang  _et al._
(2021); Tong  _et al._ (2022)

$\displaystyle U_{S,\tau}(p)=$ $\displaystyle\
\frac{\Sigma^{++}_{\tau}(p)-\Sigma^{--}_{\tau}(p)}{2},$ (3a) $\displaystyle
U_{0,\tau}(p)=$ $\displaystyle\
\frac{E^{*}_{\bm{p},\tau}}{M^{*}_{\bm{p},\tau}}\frac{\Sigma^{++}_{\tau}(p)+\Sigma^{--}_{\tau}(p)}{2}-\frac{p^{*}_{\tau}}{M^{*}_{\bm{p},\tau}}\Sigma^{-+}_{\tau}(p),$
(3b) $\displaystyle U_{V,\tau}(p)=$ $\displaystyle\
-\frac{p^{*}_{\tau}}{M^{*}_{\bm{p},\tau}}\frac{\Sigma^{++}_{\tau}(p)+\Sigma^{--}_{\tau}(p)}{2}+\frac{E^{*}_{\bm{p},\tau}}{M^{*}_{\bm{p},\tau}}\Sigma^{-+}_{\tau}(p).$
(3c)

Here the plus and minus signs in the superscripts denote the PESs and NESs,
respectively. The effective quantities in Eq. (3) are defined as
$\bm{p}^{*}_{\tau}=\bm{p}+\hat{\bm{p}}U_{V,\tau}(p)$,
$M^{*}_{\bm{p},\tau}=M+U_{S,\tau}(p)$, and
$E^{*}_{\bm{p},\tau}=E_{\bm{p},\tau}-U_{0,\tau}(p)$. The effective mass
$M^{*}_{\bm{p},\tau}$ is also known as the Dirac mass $M^{*}_{D,\tau}$.

The quantities $\Sigma^{++}_{\tau}(p)$, $\Sigma^{-+}_{\tau}(p)$, and
$\Sigma^{--}_{\tau}(p)$ in Eq. (3) are calculated by summing up the effective
$G$ matrix elements with all the nucleons inside the Fermi sea

$\displaystyle\Sigma^{++}_{\tau}(p)=$ $\displaystyle\
\sum_{s^{\prime}\tau^{\prime}}\int^{k^{\tau^{\prime}}_{F}}_{0}\frac{d^{3}p^{\prime}}{(2\pi)^{3}}\frac{M^{*}_{\bm{p}^{\prime},\tau^{\prime}}}{E^{*}_{\bm{p}^{\prime},\tau^{\prime}}}\langle\bar{u}_{\tau}(\bm{p},1/2)\bar{u}_{\tau^{\prime}}(\bm{p}^{\prime},s^{\prime})|\bar{G}^{++++}(W)|u_{\tau}(\bm{p},1/2)u_{\tau^{\prime}}(\bm{p}^{\prime},s^{\prime})\rangle,$
(4a) $\displaystyle\Sigma^{-+}_{\tau}(p)=$ $\displaystyle\
\sum_{s^{\prime}\tau^{\prime}}\int^{k^{\tau^{\prime}}_{F}}_{0}\frac{d^{3}p^{\prime}}{(2\pi)^{3}}\frac{M^{*}_{\bm{p}^{\prime},\tau^{\prime}}}{E^{*}_{\bm{p}^{\prime},\tau^{\prime}}}\langle\bar{v}_{\tau}(\bm{p},1/2)\bar{u}_{\tau^{\prime}}(\bm{p}^{\prime},s^{\prime})|\bar{G}^{-+++}(W)|u_{\tau}(\bm{p},1/2)u_{\tau^{\prime}}(\bm{p}^{\prime},s^{\prime})\rangle,$
(4b) $\displaystyle\Sigma^{--}_{\tau}(p)=$ $\displaystyle\
\sum_{s^{\prime}\tau^{\prime}}\int^{k^{\tau^{\prime}}_{F}}_{0}\frac{d^{3}p^{\prime}}{(2\pi)^{3}}\frac{M^{*}_{\bm{p}^{\prime},\tau^{\prime}}}{E^{*}_{\bm{p}^{\prime},\tau^{\prime}}}\langle\bar{v}_{\tau}(\bm{p},1/2)\bar{u}_{\tau^{\prime}}(\bm{p}^{\prime},s^{\prime})|\bar{G}^{-+-+}(W)|v_{\tau}(\bm{p},1/2)u_{\tau^{\prime}}(\bm{p}^{\prime},s^{\prime})\rangle.$
(4c)

Here $\bar{G}$ is the antisymmetrized $G$ matrix, $k^{\tau}_{F}$ is the Fermi
momentum, and the starting energy is denoted by $W$. The additional factors
$M^{*}/E^{*}$ in Eq. (4) arise due to the normalization of the Dirac spinors,
i.e., $\bar{u}u=1,\bar{v}v=-1$. The effective interaction $G$ matrix is the
solution of the in-medium Thompson equation Brockmann and Machleidt (1990)
which describes the two-body scattering

$\begin{split}G_{\tau\tau^{\prime}}(\bm{q}^{\prime},\bm{q}|\bm{P},W)=&\
V_{\tau\tau^{\prime}}(\bm{q}^{\prime},\bm{q}|\bm{P})+\int\frac{d^{3}k}{(2\pi)^{3}}V_{\tau\tau^{\prime}}(\bm{q}^{\prime},\bm{k}|\bm{P})\\\
&\times\frac{M^{*}_{\bm{P}+\bm{k},\tau}M^{*}_{\bm{P}-\bm{k},\tau^{\prime}}}{E^{*}_{\bm{P}+\bm{k},\tau}E^{*}_{\bm{P}-\bm{k},\tau^{\prime}}}\frac{Q_{\tau\tau^{\prime}}(\bm{k},\bm{P})}{W-E_{\bm{P}+\bm{k},\tau}-E_{\bm{P}-\bm{k},\tau^{\prime}}}G_{\tau\tau^{\prime}}(\bm{k},\bm{q}|\bm{P},W).\end{split}$
(5)

The labels of PESs and NESs have been suppressed. In this work, the realistic
$NN$ interaction $V_{\tau\tau^{\prime}}$ is chosen as the Bonn potential
Machleidt (1989). Here $\bm{P}=\frac{1}{2}({\bm{k}}_{1}+{\bm{k}}_{2})$ and
$\bm{k}=\frac{1}{2}({\bm{k}}_{1}-{\bm{k}}_{2})$ are the center-of-mass and the
relative momenta of the two interacting nucleons with momenta ${\bm{k}}_{1}$
and ${\bm{k}}_{2}$. The initial, intermediate, and final relative momenta of
the two nucleons are $\bm{q},\bm{k}$, and $\bm{q}^{\prime}$, respectively. The
$NN$ scattering in the nuclear medium is restricted with the Pauli operator
$Q_{\tau\tau^{\prime}}(\bm{k},\bm{P})$.

Equations (1), (3), (4), and (5) constitute a coupled system that has to be
solved in a self-consistent way. After the convergence of SPPs, the single-
particle and bulk properties of nuclear matter can be calculated
straightforwardly van Dalen  _et al._ (2004); Katayama and Saito (2013); Tong
_et al._ (2018).

Figure 1: (Color online) The isospin splittings of Dirac masses
$(M^{*}_{D,n}-M^{*}_{D,p})/M$ as functions of the asymmetry parameter $\alpha$
obtained by the RBHF theory in the full Dirac space (red solid), the
projection method (olive dashed), and the MIA method (blue dotted). The
density is fixed at the normal nuclear saturation density $\rho=0.16\
\text{fm}^{-3}$. The $NN$ interaction Bonn A is used.

In SNM, the Dirac masses for nucleon calculated by the MIA method and the
projection method are both quantitatively close to the results obtained in the
full Dirac space Wang  _et al._ (2022b). However, the situation changes
dramatically in ANM. Figure 1 depicts the isospin splitting of the Dirac mass
$(M^{*}_{D,n}-M^{*}_{D,p})/M$ calculated with the RBHF theory as a function of
the asymmetry parameter $\alpha=(\rho_{n}-\rho_{p})/\rho$. The projection
method finds that in ANM there is $M^{*}_{D,n}-M^{*}_{D,p}<0$, while the MIA
method leads to the opposite sign $M^{*}_{D,n}-M^{*}_{D,p}>0$. This
contradiction, which is well known since 1997 Ulrych and Müther (1997), has
been clarified recently by considering the PESs and the NESs simultaneously in
Ref. Wang  _et al._ (2022b). Comparing to the RBHF results calculated in the
full Dirac space as also shown in Fig. 1, the projection method obtains a
qualitatively consistent isospin dependence of the Dirac mass with
overestimated amplitudes, while the MIA method gives a wrong sign.

It is still unclear why the MIA method succeeds in SNM but fails in ANM. To
reach the answer, we notice that this method has two essential procedures. The
first procedure can be called as the assumption procedure, which assumes that
the scalar potential and the timelike part of the vector potential are
momentum independent and the spacelike part of the vector potential is
negligible, i.e.,

$U_{S,\tau}(p)\approx U_{S,\tau},\quad U_{0,\tau}(p)\approx U_{0,\tau},\quad
U_{V,\tau}(p)\approx 0.$ (6)

The second procedure can be called as the extraction procedure, which extracts
the two constants $U_{S,\tau}$ and $U_{0,\tau}$ from the single-particle
potential energy $U_{\tau}(k)$ at two momenta $(k^{\tau}_{1},k^{\tau}_{2})$

$\displaystyle U_{\tau}(k^{\tau}_{1})=$ $\displaystyle\
\frac{M+U_{S,\tau}}{\sqrt{(k^{\tau}_{1})^{2}+(M+U_{S,\tau})^{2}}}U_{S,\tau}+U_{0,\tau},$
(7a) $\displaystyle U_{\tau}(k^{\tau}_{2})=$ $\displaystyle\
\frac{M+U_{S,\tau}}{\sqrt{(k^{\tau}_{2})^{2}+(M+U_{S,\tau})^{2}}}U_{S,\tau}+U_{0,\tau}.$
(7b)

The quantity $U_{\tau}(k)$ in Eq. (7) is calculated as
$(M^{*}_{\tau}/E^{*}_{\bm{k},\tau})\Sigma^{++}_{\tau}(k)$.

Figure 2: (Color online) Binding energies per nucleon $E/A$ of SNM and PNM as
functions of the nucleon density $\rho$ calculated in the full Dirac space
with the assumption procedure. For comparison, the self-consistent results
obtained by the RBHF theory in the full Dirac space are also shown.

The RBHF calculations in the full Dirac space provide us an opportunity to
analyze in detail the isospin splitting of the Dirac mass. In the following,
we try to study further by testing separately the two procedures of the MIA
method from the perspective in the full Dirac space.

Firstly, we apply the assumption procedure (6) in the full Dirac space, i.e.,
the quantities $U_{S,\tau}(p)$ in Eq. (3a) and $U_{0,\tau}(p)$ in Eq. (3b) are
assumed to be momentum independent and are calculated at the Fermi momentum
$k^{\tau}_{F}$, while $U_{V,\tau}(p)$ in Eq. (3c) is set as zero. The newly
obtained quantities $U_{S,\tau}$ and $U_{0,\tau}$ will be used to update the
Dirac spinors and $G$ matrix for next iteration. After the convergence of
SPPs, we calculate the binding energies per nucleon for SNM and pure neutron
matter (PNM). As can be seen in Fig. 2, the $E/A$ for SNM and PNM are very
close to those obtained by the RBHF theory in the full Dirac space. This
indicates that the assumption procedure is reasonable for describing the bulk
properties of nuclear matter.

Figure 3: (Color online) Dirac masses for neutron (solid lines) and proton
(dashed lines) as functions of the asymmetric parameter $\alpha$ at density
$\rho=0.16\ \text{fm}^{-3}$ calculated by the RBHF theory (a) in the full
Dirac space, (b) with the assumption procedure, and (c) with the extraction
procedure. Details can be found in the text.

In Fig. 3(b), we show the Dirac masses of the neutron and proton obtained in
the full Dirac space with the assumption procedure. The density is fixed at
$\rho=0.16\ \text{fm}^{-3}$. The relation of $M^{*}_{D,n}-M^{*}_{D,p}<0$ is
found, which is consistent with the result obtained by the RBHF theory in the
full Dirac space, as shown in Fig. 3(a). This indicates that the assumption
procedure is not a sufficient condition that directly leads to the wrong sign
when one calculates the isospin splitting of the Dirac mass in ANM.

Secondly, we test the influence of the extraction procedure (7) from the
perspective in the full Dirac space. Starting from the converged $U_{\tau}(p)$
obtained by the RBHF theory in the full Dirac space, we extract the two
constants $U_{S,\tau}$ and $U_{0,\tau}$ by using Eq. (7) with two selected
momenta $(0.7k^{\tau}_{F},k^{\tau}_{F})$ and then calculate the Dirac mass.
The resulting isospin splitting of the Dirac mass is shown in Fig. 3(c). It is
seen that $M^{*}_{D,n}-M^{*}_{D,p}>0$ is obtained for the entire region of the
asymmetry parameter $\alpha$, which is opposite to results obtained in the
full Dirac space (Fig. 3(a)). This shows that it is the extraction procedure
(7) that leads to the wrong sign of the isospin splitting of the Dirac mass.

Figure 4: (Color online) (a) The single-particle potential energy
$U_{\tau}(k)$ obtained in the full Dirac space as a function of momentum $k$.
The cases for SNM and ANM with $\alpha=0.5$ and $\rho=0.16\ \text{fm}^{-3}$
are shown. The momenta $0.7k^{\tau}_{F}$ and $k^{\tau}_{F}$ are highlighted
with empty dots and solid squares respectively. (b) The Dirac mass
$M^{*}_{D}/M$ as a functional of dimensionless quantity $f$ in Eq. (12). (c)
The dimensionless function $f$ (11) obtained in the full Dirac space as a
function of momentum $k$. The notation for lines and symbols is the same as
panel (a).

To further investigate how the wrong isospin dependence of the Dirac mass
emerges from the extraction procedure, it is insightful to solve Eq. (7) to
get a formal expression of the Dirac mass. We start from the case of SNM and
suppress the isospin indexes for moment. For $M\simeq 1000$ MeV,
$U_{S}\simeq-400$ MeV, and $k_{F}=1.34\ \text{fm}^{-1}$, one notices that the
typical choice $k_{1}=0.7k_{F}$ and $k_{2}=k_{F}$ leads to
$\left[k_{1}/(M+U_{S})\right]^{2}\simeq 0.1\ll 1$,
$\left[k_{2}/(M+U_{S})\right]^{2}\simeq 0.2\ll 1$. This allows one to expand
the square root in Eq. (7) to the first order,

$U(k)=-\frac{1}{2}\frac{U_{S}}{(M+U_{S})^{2}}k^{2}+U_{0},$ (8)

which indicates a quadratic form of the single-particle potential energy
$U(k)$. With Eq. (8) it is straightforward to find the difference between
$U(k_{2})$ and $U(k_{1})$ as

$U(k_{2})-U(k_{1})=-\frac{1}{2}\frac{U_{S}}{(M+U_{S})^{2}}\left(k^{2}_{2}-k^{2}_{1}\right).$
(9)

For simplicity, we consider the limit that two momenta $k_{1}$ and $k_{2}$ are
very close and denote them as $k$. In this case Eq. (9) can be written as a
quadratic equation for $M^{*}_{D}/M$

$f(k)(M^{*}_{D}/M)^{2}+M^{*}_{D}/M-1=0,$ (10)

where the dimensionless function $f(k)$ is defined as

$f(k)\equiv\lim_{\begin{subarray}{c}k_{2}\rightarrow k\\\ k_{1}\rightarrow
k\end{subarray}}2M\frac{U(k_{2})-U(k_{1})}{k^{2}_{2}-k^{2}_{1}}=M\frac{U^{\prime}(k)}{k}.$
(11)

The derivative $U^{\prime}(k)$ describes the momentum dependence of the
single-particle potential energy $U(k)$.

In Fig. 4(a), the single-particle potential energy $U_{\tau}(k)$ obtained in
the full Dirac space for SNM and also ANM with $\alpha=0.5$ and $\rho=0.16\
\text{fm}^{-3}$ is shown as a function of momentum $k$. It is found that
$U_{\tau}(k)$ is a monotonic function of the $k$. Therefore, the function
$f(k)$ is positive definite and the solution of Eq. (10) is

$M^{*}_{D}/M=\frac{\sqrt{1+4f(k)}-1}{2f(k)}.$ (12)

From Fig. 4(b), the Dirac mass $M^{*}_{D}/M$ is decreasing with increasing
$f$. Therefore, the sign of $M^{*}_{D,n}-M^{*}_{D,p}$ is opposite to that of
$f_{n}(k_{n})-f_{p}(k_{p})$, i.e.,

$\begin{cases}M^{*}_{D,n}-M^{*}_{D,p}<0,\quad\text{if}\quad
f_{n}(k_{n})>f_{p}(k_{p}),\\\ M^{*}_{D,n}-M^{*}_{D,p}>0,\quad\text{if}\quad
f_{n}(k_{n})<f_{p}(k_{p}).\end{cases}$ (13)

In Fig. 4(c), the function $f_{\tau}(k)$ for nucleon $\tau$ obtained in the
full Dirac space is shown. It is found that $f_{\tau}(k)$ is decreasing with
the increasing $k$. Besides, generally there is $f_{n}(k)<f_{p}(k)$ for
$k<k^{p}_{F}$. This indicates that the momentum dependence of $U_{n}(k)$ is
weaker than that of $U_{p}(k)$.

In practice, one usually chooses $k^{\tau}_{1}=0.7k^{\tau}_{F}$ and
$k^{\tau}_{2}=k^{\tau}_{F}$ Tong  _et al._ (2018); Wang  _et al._ (2022b). For
ANM with $\alpha>0$, the Fermi momentum for neutron $k^{n}_{F}$ is larger than
the of proton, i.e., $k^{n}_{F}>k^{p}_{F}$. As shown by the symbols in panels
(a) and (c) in Fig. 4, this choice leads to $k_{n}>k_{p}$ and
$f_{n}(k_{n})<f_{p}(k_{p})$. Specifically, starting from the self-consistent
$U_{\tau}(k)$ obtained in the full Dirac space, the choices
$(0.7k^{n}_{F},k^{n}_{F})$ for neutron and $(0.7k^{p}_{F},k^{p}_{F})$ for
proton lead to $f_{n}(k_{n})\approx 0.6$ and $f_{p}(k_{p})\approx 0.8$.
According to the analysis in Eq. (13), the isospin splitting of the Dirac mass
$M^{*}_{D,n}-M^{*}_{D,p}>0$ emerges. This explains the reason why the MIA
method leads to the wrong isospin splitting of the Dirac mass in ANM.

Come back to Eq. (8), one can see that the extraction procedure forcely
assumes that the momentum dependence of the single-particle potential energy
to be a quadratic form, where the strength
$-\frac{1}{2}\frac{U_{S}}{(M+U_{S})^{2}}$ is solely determined by the scalar
potential $U_{S}$. In ANM, the momentum dependence of $U_{n}(k)$ is generally
weaker than that of $U_{p}(k)$, as shown in Fig. 4(c). It is exactly this fact
that leads to $U_{S,n}>U_{S,p}$ and the wrong sign of the isospin splitting of
the Dirac mass $M^{*}_{D,n}-M^{*}_{D,p}>0$.

In summary, the relativistic Brueckner-Hartree-Fock (RBHF) theory plays an
important role in deriving the nuclear many-body properties from the realistic
nucleon-nucleon interactions. In comparison to the results obtained self-
consistently in the full Dirac space, the momentum-independence approximation
(MIA) method leads to wrong isospin splitting of the Dirac mass in asymmetric
nuclear matter (ANM). The performance of this method is explored in detail
from the perspective in the full Dirac space. The assumption procedure of the
MIA method, which assumes that the single-particle potentials are momentum
independent, is not a sufficient condition that directly leads to the wrong
sign of the isospin splitting of the Dirac mass, while the extraction
procedure of the MIA method, which extracts the single-particle potentials
from the single-particle potential energy, is found to be responsible for the
wrong isospin splitting of the Dirac mass. By approximately solving the set of
equations involved in the extraction procedure, a formal expression of the
Dirac mass is obtained. With the typical choice of momenta adopted in
practical MIA calculations, the wrong isospin splitting of the Dirac mass is
found. We conclude that the wrong isospin splitting of the Dirac mass emerges
from the fact that the extraction procedure forcely assumes the momentum
dependence of the single-particle potential energy to be a quadratic form
where the strength is solely determined by the constant scalar potential. This
work improves greatly our understanding on the isospin splitting of the Dirac
mass with the RBHF theory.

###### Acknowledgements.

The authors thank the discussion with Professor Peter Ring. This work was
supported in part by the China Postdoctoral Science Foundation under Grant No.
2021M700610, the National Natural Science Foundation of China (NSFC) under
Grants No. 12205030, the Fundamental Research Funds for the Central
Universities under Grants No. 2020CDJQY-Z003 and No. 2021CDJZYJH-003, and the
Institute for Basic Science under Grant No. IBS-R031-D1. Part of this work was
achieved by using the supercomputer OCTOPUS at the Cybermedia Center, Osaka
University under the support of Research Center for Nuclear Physics of Osaka
University and the High Performance Computing Resources in the Research
Solution Center, Institute for Basic Science.

## References

  * Baldo and Maieron (2007) M. Baldo and C. Maieron, Journal of Physics G: Nuclear and Particle Physics 34, R243 (2007).
  * Barrett  _et al._ (2013) B. R. Barrett, P. Navratil, and J. P. Vary, Prog. Part. Nucl. Phys. 69, 131 (2013).
  * Hagen  _et al._ (2014) G. Hagen, T. Papenbrock, M. Hjorth-Jensen, and D. J. Dean, Rep. Prog. Phys. 77, 096302 (2014).
  * Carlson  _et al._ (2015) J. Carlson, S. Gandolfi, F. Pederiva, S. C. Pieper, R. Schiavilla, K. E. Schmidt, and R. B. Wiringa, Rev. Mod. Phys. 87, 1067 (2015).
  * Hergert  _et al._ (2016) H. Hergert, S. Bogner, T. Morris, A. Schwenk, and K. Tsukiyama, Phys. Rep. 621, 165 (2016).
  * Shen  _et al._ (2019) S. Shen, H. Liang, W. Long, J. Meng, and P. Ring, Prog. Part. Nucl. Phys. 109, 103713 (2019).
  * Brueckner  _et al._ (1954) K. A. Brueckner, C. A. Levinson, and H. M. Mahmoud, Phys. Rev. 95, 217 (1954).
  * Day (1967) B. D. Day, Reviews of Modern Physics 39, 719 (1967).
  * Coester  _et al._ (1970) F. Coester, S. Cohen, B. Day, and C. M. Vincent, Physical Review C 1, 769 (1970).
  * Grangé  _et al._ (1989) P. Grangé, A. Lejeune, M. Martzolff, and J.-F. Mathiot, Phys. Rev. C 40, 1040 (1989).
  * Baldo  _et al._ (1997) M. Baldo, I. Bombaci, and G. F. Burgio, Astron. Astrophys. 328, 274 (1997).
  * Zuo  _et al._ (2002) W. Zuo, A. Lejeune, U. Lombardo, and J. Mathiot, Nuclear Physics A 706, 418 (2002).
  * Zhou  _et al._ (2004) X. R. Zhou, G. F. Burgio, U. Lombardo, H.-J. Schulze, and W. Zuo, Phys. Rev. C 69, 018801 (2004).
  * Li  _et al._ (2006) Z. H. Li, U. Lombardo, H.-J. Schulze, W. Zuo, L. W. Chen, and H. R. Ma, Phys. Rev. C 74, 047304 (2006).
  * Li  _et al._ (2008) Z. H. Li, U. Lombardo, H.-J. Schulze, and W. Zuo, Phys. Rev. C 77, 034316 (2008).
  * Anastasio  _et al._ (1980) M. R. Anastasio, L. S. Celenza, and C. M. Shakin, Phys. Rev. Lett. 45, 2096 (1980).
  * Anastasio  _et al._ (1981) M. R. Anastasio, L. S. Celenza, and C. M. Shakin, Phys. Rev. C 23, 2273 (1981).
  * Brockmann and Machleidt (1984) R. Brockmann and R. Machleidt, Physics Letters B 149, 283 (1984).
  * Horowitz and Serot (1987) C. J. Horowitz and B. D. Serot, Nucl. Phys. A 464, 613 (1987).
  * ter Haar and Malfliet (1987) B. ter Haar and R. Malfliet, Physical Review Letters 59, 1652 (1987).
  * Machleidt (1989) R. Machleidt, Adv. Nucl. Phys. 19, 189 (1989).
  * Sammarruca  _et al._ (2012) F. Sammarruca, B. Chen, L. Coraggio, N. Itaco, and R. Machleidt, Phys. Rev. C 86, 054317 (2012).
  * Jaroszewicz and Brodsky (1991) T. Jaroszewicz and S. J. Brodsky, Phys. Rev. C 43, 1946 (1991).
  * Brown (1987) G. E. Brown, Physica Scripta 36, 209 (1987).
  * Serot and Walecka (1986) B. D. Serot and J. D. Walecka, Adv. Nucl. Phys. 16, 1 (1986).
  * Brockmann and Machleidt (1990) R. Brockmann and R. Machleidt, Phys. Rev. C 42, 1965 (1990).
  * Li  _et al._ (1992) G. Q. Li, R. Machleidt, and R. Brockmann, Phys. Rev. C 45, 2782 (1992).
  * Alonso and Sammarruca (2003) D. Alonso and F. Sammarruca, Physical Review C 67, 054301 (2003).
  * Tong  _et al._ (2018) H. Tong, X.-L. Ren, P. Ring, S.-H. Shen, S.-B. Wang, and J. Meng, Phys. Rev. C 98, 054302 (2018).
  * Wang  _et al._ (2020) C. Wang, J. Hu, Y. Zhang, and H. Shen, Journal of Physics G: Nuclear and Particle Physics 47, 105108 (2020).
  * Gross-Boelting  _et al._ (1999) T. Gross-Boelting, C. Fuchs, and A. Faessler, Nuclear Physics A 648, 105 (1999).
  * van Dalen  _et al._ (2004) E. N. E. van Dalen, C. Fuchs, and A. Faessler, Nuclear Physics A 744, 227 (2004).
  * van Dalen  _et al._ (2005a) E. N. E. van Dalen, C. Fuchs, and A. Faessler, Physical Review Letters 95, 022302 (2005a).
  * van Dalen  _et al._ (2005b) E. N. E. van Dalen, C. Fuchs, and A. Faessler, Phys. Rev. C 72, 065803 (2005b).
  * van Dalen  _et al._ (2007) E. N. E. van Dalen, C. Fuchs, and A. Faessler, The European Physical Journal A 31, 29 (2007).
  * Nuppenau  _et al._ (1989) C. Nuppenau, Y. Lee, and A. MacKellar, Nuclear Physics A 504, 839 (1989).
  * Wang  _et al._ (2021) S. Wang, Q. Zhao, P. Ring, and J. Meng, Phys. Rev. C 103, 054319 (2021).
  * Ulrych and Müther (1997) S. Ulrych and H. Müther, Phys. Rev. C 56, 1788 (1997).
  * Van Giai  _et al._ (2010) N. Van Giai, B. V. Carlson, Z. Ma, and H. Wolter, J. Phys. G 37, 064043 (2010).
  * Wang  _et al._ (2022a) S. Wang, H. Tong, and C. Wang, Phys. Rev. C 105, 054309 (2022a).
  * Wang  _et al._ (2022b) S. Wang, H. Tong, Q. Zhao, C. Wang, P. Ring, and J. Meng, Phys. Rev. C 106, L021305 (2022b).
  * Wang  _et al._ (2023) S. Wang, H. Tong, Q. Zhao, C. Wang, P. Ring, and J. Meng, Phys. Rev. C 108, L031303 (2023).
  * Tong  _et al._ (2023) H. Tong, J. Gao, C. Wang, and S. Wang, Phys. Rev. C 107, 034302 (2023).
  * Tong  _et al._ (2022) H. Tong, C. Wang, and S. Wang, The Astrophysical Journal 930, 137 (2022).
  * Wang  _et al._ (2022c) S. Wang, C. Wang, and H. Tong, Phys. Rev. C 106, 045804 (2022c).
  * Qu  _et al._ (2023) X. Qu, H. Tong, C. Wang, and S. Wang, Sci. China-Phys. Mech. Astron. 66, 242011 (2023).
  * Katayama and Saito (2013) T. Katayama and K. Saito, Phys. Rev. C 88, 035805 (2013), 1307.2067 .


# Spatio-temporal Motion Planning for Autonomous Vehicles with Trapezoidal
Prism Corridors and Bézier Curves

Srujan Deolasee1, Qin Lin2, Jialun Li3, and John M. Dolan4 *This work was
supported by the National Science Foundation.1Srujan Deolasee is with
Department of Computer Science & Information Systems, Birla Institute of
Technology and Science, Pilani, Rajasthan, India<EMAIL_ADDRESS>pilani.ac.in2Qin Lin is with the Department of Electrical Engineering and
Computer Science, Cleveland State University, Cleveland, OH, the USA
<EMAIL_ADDRESS>(corresponding author)3Jialun Li is with Dajiang
Innovation Technology Co., Ltd, Shenzhen, China<EMAIL_ADDRESS>M.
Dolan is with the Robotics Institute, Carnegie Mellon University, Pittsburgh,
PA, the USA<EMAIL_ADDRESS>
###### Abstract

Safety-guaranteed motion planning is critical for self-driving cars to
generate collision-free trajectories. A layered motion planning approach with
decoupled path and speed planning is widely used for this purpose. This
approach is prone to be suboptimal in the presence of dynamic obstacles.
Spatial-temporal approaches deal with path planning and speed planning
simultaneously; however, the existing methods only support simple-shaped
corridors like cuboids, which restrict the search space for optimization in
complex scenarios. We propose to use trapezoidal prism-shaped corridors for
optimization, which significantly enlarges the solution space compared to the
existing cuboidal corridors-based method. Finally, a piecewise Bézier curve
optimization is conducted in our proposed corridors. This formulation
theoretically guarantees the safety of the continuous-time trajectory. We
validate the efficiency and effectiveness of the proposed approach in
numerical and CommonRoad simulations.

###### Index Terms:

Autonomous vehicle, motion planning, trajectory optimization

##  I INTRODUCTION

Motion planning is one of the key modules in autonomous driving systems. The
task of motion planning in a dynamic traffic environment is to generate
trajectories for a low-level controller to follow considering collision-free
safety constraints, dynamic feasibility, and comfort. A Frenet frame [1] is
commonly used for motion planning due to the significant advantage of its
independence from complex road geometry. The lateral motion (in the $L$
direction) and longitudinal motion (in the $S$ direction) can be projected
onto the reference, which is usually the centerline of the road with an
arbitrary shape. Including the time dimension $T$, a 3D $S-L-T$ coordinate
system can be established for insightful and convenient planning.

Path-speed (or Layered planning) is a practical real-time solution to
decompose a planning problem into two stages: path planning and speed planning
[2, 3, 4]. A path ($S-L$) is generated in the first stage in a static or low-
speed environment. The generation of the speed profile ($S-T$ or $L-T$) in the
speed planning stage allows an AV to respond to dynamic obstacles. The
significant limitation of layered planning is that it is prone to be
suboptimal in the presence of dynamic obstacles in complicated scenarios.

Spatio-temporal planning considers spatial and temporal maneuvers
simultaneously [5, 6, 7]. This method of direct optimization in the 3D $S-L-T$
space is generally superior to the layered planning approach due to the larger
search space. See the motivating example illustrated in Fig. 1 (discussed in
detail in Section III.A): in driving scenarios involving even small deviations
along the lateral direction, coupled longitudinal and lateral planning helps
guarantee global optimality.

(a)

(b)

(c)

(d)

Figure 1: Yielding example. (a) bird’s-eye view; (b) our spatio-temporal
planning in 3D $S-L-T$ graph; (c) another angle of view for car C; (d) another
angle of view for car A.

Ensuring safety is the vital objective of motion planning. Many existing speed
planning methods use discrete time instants to impose safety constraints.
However, a provable safety guarantee independent of sampling time in
continuous time space is preferable. To address this problem, the spatial
corridor is widely applied in trajectory generation. We are motivated by these
efforts to further extend the spatial corridor to the spatio-temporal domain
to cope better with dynamic obstacles. The convex hull property of Bézier
polynomials is leveraged to enforce that the continuous trajectory always
falls into a safe spatio-temporal region. In addition, such an optimization
problem’s solution space is enlarged via our proposed trapezoidal-prism-shaped
corridors.

The main contributions of our work can be briefly described as follows:

  1. 1.

We propose an efficient convexification algorithm to construct 3D convex-
feasible regions consisting of trapezoidal-prism-shaped corridors.

  2. 2.

We provide a sufficient condition on coefficients of the Bézier polynomials to
theoretically guarantee the trajectory’s safety in trapezoidal-prism
corridors. Compared with existing cuboidal corridors [5], the condition is
relaxed and the solution space is significantly enlarged, which leads to a
higher chance of finding an optimal solution.

The remainder of this paper is structured as follows. We review related works
in Sec. II. We introduce necessary notations and background materials in Sec.
III. The 3D convex safe region construction is presented in Sec. IV. In Sec.
V, we present our optimization formulation. The simulation results and
analysis can be found in Sec. VI. We make concluding remarks in Sec. VII.

##  II RELATED WORKS

###  II-A Speed Planning

Speed planning techniques can be classified into three categories: 1) search
and optimization; 2) sampling lattices and selecting the minimum-cost
trajectory; 3) approximated optimization. The Search and optimization method
searches for the best candidate speed profile and optimizes the curve for
smoothness; see the post-optimization method[3], the Baidu EM motion planner
[8], and the Piecewise-Jerk Speed Optimization [9]. The Sampling approach
samples different speed lattices combined with path lattices. The generated
local spatial-temporal trajectories are evaluated and the one with minimum
cost is selected. Related works can be found in [2, 3, 4]. Most works in the
first and the second categories conduct search and optimization directly in
the $S-T$ graph. Approximated optimization considers a vehicle dynamic model
in a sequential optimization problem; see convex feasible set algorithm [10]
and optimal control methods, such as model predictive control (MPC)[11] and
constrained iterative linear quadratic regulator (CiLQR)[12]. The advantage of
these approaches is that they mitigate the planning and control inconsistency
problem since the dynamic model has already been considered in the planning
layer. However, the disadvantage is the high computation cost.

###  II-B Corridor generation for Autonomous Vehicles

The spatial corridor is widely used in trajectory generation. Some previous
works generate the corridors in a static environment and cannot deal with
dynamic obstacles [13], [14]. Liu et al. [15] find a convex feasible set
around the reference trajectory, but the computation complexity restricts the
method for real-time applications. Our previous work proposes the use of
trapezoidal corridors for convexifying 2D space in the $S-T$ graph [16]. Zhang
et al. present a general convex spatio-temporal corridors-based approach [17].
Xu et al. propose using a modified vertical cell decomposition approach for
speed planning in [18]. All these methods suffer from the limitations of the
layered planning approach discussed in the previous section. Ding et al. use
the spatio-temporal semantic corridor (SSC) method to uniformly express
obstacles and traffic rules in the 3D $S-L-T$ space [5]. However, restricting
the shape of the corridors to simple cuboids drastically limits the search
space for optimization in complex scenarios. Our proposed method of extending
trapezoid-shaped 2D corridors in $S-T$ to 3D $S-L-T$ space significantly
enlarges the solution space for trajectory optimization. This enables us to
extend the spatial corridor to the spatio-temporal domain to cope with dynamic
obstacles while meeting the real-time requirement.

###  II-C Bézier Polynomials-Based Planning

Previously, monomial basis polynomials have been used to generate trajectories
[1], [8]. However, these methods often fail to represent highly constrained
maneuvers in the presence of dynamic obstacles. They also fail to give safety
guarantees between sample points, as the constraints are only enforced/checked
on a finite set of sampled points. In [19], a smooth and continuous speed
profile is computed by proper curve concatenation without optimization and
dynamic obstacles. Bézier polynomials combined with rectangular corridors was
originated in the area of unmanned aerial vehicles (UAVs) [20], [21]. Ding et
al. extended it for motion planning of unmanned ground vehicles (UGVs) [5].
The significant limitation is that the proposed cuboidal corridor
representation fails to make the most of free space for optimization. In our
work, we propose to use time-dependent trapezoidal prism-shaped corridors and
give sufficient conditions to enforce Bézier curves in these time-dependent
corridors for safety. It is theoretically proved that the trapezoidal prism-
shaped corridors can enlarge the solution space for improved optimization.

##  III $S-L-T$ GRAPH AND TRAJECTORY REPRESENTATION

In this section, we briefly introduce background materials on the $S-L-T$
graph, Bézier polynomials, and trajectory representations using piecewise
Bézier polynomials.

###  III-A Representing Dynamic Agents in $S-L-T$ graph

The $S-L-T$ graph represents all traffic participants’ positions at each time
step including the past, current, and prediction. As an example, in Fig. 1a,
we take the case of two cars moving at constant speeds for simplicity. The
scenario is described as follows: car A and the ego vehicle are driving in a
lane which has a static obstacle (e.g., a construction site). A lane change
maneuver is enforced for both vehicles. We list the following typical
entities:

1) Static obstacle (6-zero-slope-faces type): As the most simple entity, all
the six faces have zero slopes, see the bird’s-eye view of the yellow block in
Fig. 1a and the cuboid in the $S-L-T$ plot in Fig. 1b.

2) Moving car with only longitudinal motion (4-zero-slope-faces type): Car C
moving straight forward is an example shown in Fig. 1a. The top and bottom
faces of Car C in Fig. 1b have zero slopes. The two faces (see another angle
of view for left and right faces in Fig. 1c) are perpendicular to the $S-L$
plane without slopes. The side length of the parallelogram along the $L$ axis
is the width of the vehicle plus the safety region. The side length of the
parallelogram along the $S$ axis is the length of the vehicle plus half the
length of the ego vehicle as a safety region.

3) Moving car with longitudinal and lateral motions (2-zero-slope-faces type):
Car A moving left and forward is an example shown in Fig. 1a. As we can see in
Fig.1b, only the top and down faces are zero-slope.

The 3D free space in the $S-L-T$ graph is non-convex in general. We propose an
over-approximation of the _2-zero-slope-faces type_ parallelepiped, e.g., car
A in Fig. 1b, into a _4-zero-slope-faces type_ parallelepiped, see the pink
inflated space in Fig. 1b and Fig. 1d. There are two significant benefits of
doing so: 1) we will show that in the presence of such parallelepipeds, we can
extend our 2D corridor construction algorithm [16] to construct 3D trapezoidal
prism-shaped corridors efficiently; 2) though we pay the cost of losing some
space due to over-approximation, the safety corridor is still significantly
larger than the state-of-the-art cuboidal corridors. Note that the
transformation between 2D and 3D in our method is without the loss of search
space. Thus, in summary, we make a good trade-off between exactness and
efficiency. The faces chosen for the over-approximation step are decided by a
simple minimization of the volume of search space compromised in the process.
Intuitively, if the lateral velocity is less than the longitudinal velocity of
the vehicle, the corresponding faces are chosen for over-approximation.

###  III-B Bézier Polynomials and Properties

A Bézier polynomial is a polynomial function represented by linear
combinations of Bernstein bases. The $n$th-order Bézier polynomial is written
as

$B(t)=c_{0}b_{n}^{0}(t)+c_{1}b_{n}^{1}(t)+\cdots+c_{n}b_{n}^{n}(t)=\sum_{i=0}^{n}c_{i}b_{n}^{i}(t)$

where the Bernstein bases satisfy $b_{n}^{i}(t)=C_{i}^{n}\cdot t^{i}\cdot(1-$
$t)^{n-i},t\in[0,1]$. The coefficients of the polynomial $c_{i}(i=$
$0,1,\ldots,n$) are also called control points. Compared to monomial
polynomials, Bézier curves have the following properties:

  * •

The time interval is defined on $t\in[0,1]$.

  * •

The Bézier polynomial starts at control point $B(0)=c_{0}$ and ends at
$B(1)=c_{n}$.

  * •

Convex hull property: The Bézier curve $B(t)$ is confined within the convex
hull of control points.

  * •

Hodograph property: By the hodograph property, the derivative of
$B(t),\dot{B}(t)$, can also be written as a Bézier polynomial with control
points $c_{i}^{1}=n\cdot\left(c_{i+1}-c_{i}\right),i=$ $0,1,\ldots,n-1$. By
applying the convex hull property to the derivative Bézier curve, the entire
dynamical profile of the original curve $B(t)$ can be confined within a given
dynamical range.

###  III-C Trajectory Representation using Bézier Polynomials

To mitigate the numerical instability issue, piecewise Bézier polynomials with
lower orders are used instead of using a high-order Bézier polynomial for the
whole planning horizon. Each piece of the trajectory is associated with one
trapezoidal-prism corridor. Note that $B(t)$ is defined on a fixed time
interval $[0,1]$. For a whole trajectory with $m+1$ pieces, in each piece
$\left[T_{k},T_{k+1}\right](k=0,1,\ldots,m)$, we use a scaling transformation
and translation transformation in the time domain to map it into the interval
$[0,1]$ [13]. Then, the whole piece-wise trajectory in one dimension
$\sigma\in\\{s,l\\}$ is:

$f^{\sigma}(t)=\left\\{\begin{array}[]{c}h_{0}B_{0}\left(\frac{t-T_{0}}{h_{0}}\right),t\in\left[0,T_{1}\right]\\\
h_{1}B_{1}\left(\frac{t-T_{1}}{h_{1}}\right),t\in\left[T_{1},T_{2}\right]\\\
\vdots\\\
h_{m}B_{m}\left(\frac{t-T_{m}}{h_{m}}\right),t\in\left[T_{m},T_{m+1}\right].\end{array}\right.$

where $h_{i}$ is the scaling transformation factor and $T_{i}$ is the
translation transformation factor for $i=0,1,\ldots,m$ with setting $T_{0}=0$.

##  IV CORRIDOR GENERATION

In this section, a convexification algorithm is introduced to construct convex
corridors from the original non-convex optimization problem for real-time
solving. A reference trajectory is often used to provide a warm start to the
optimization process. In this work, we use simple piecewise functions for
generating valid reference waypoints in the configuration space of the ego
vehicle.

###  IV-A Piecewise Convex Safe Regions Representations

Suppose the whole safe region is divided into $m+1$ pieces with time intervals
$\left[T_{0},T_{1}\right],\ldots,\left[T_{m},T\right]$ and $T=T_{m+1}$, with
each interval corresponding to a convex safe region. The details of such a
convexification algorithm will be introduced in the next section. The $k$-th
convex safe region in $S-L-T$ space can be represented as

$\mathcal{S}_{k}=\\{\ (t_{i},s_{i},l_{i})\ |$
$\underline{p^{k}_{0}}+h_{k}\underline{p^{k}_{1}}\frac{t_{i}-T_{k}}{h_{k}}\leq
s_{i}\leq\overline{p^{k}_{0}}+h_{k}\overline{p^{k}_{1}}\frac{t_{i}-T_{k}}{h_{k}},$
$l_{beg}\leq l_{i}\leq l_{end},\ t_{i}\in[T_{k},T_{k+1}]\ \\}$

where $s_{i}$ and $l_{i}$ are the longitudinal and lateral coordinates of the
$i^{th}$ control point respectively,
$\underline{p_{0}^{k}},\underline{p_{1}^{k}}$ are bias and skew of the lower
bound and $\overline{p_{0}^{k}},\overline{p_{1}^{k}}$ are those of the upper
bound. $h_{k}$ denotes the length of the $k$-th time interval and satisfies
$h_{k}=T_{k+1}-T_{k},k=0,1,\ldots,m$.

Then, the whole safe region is the union of a set of piecewise-safe sub-
regions: $\mathcal{S}=\mathcal{S}_{0}\cup\cdots\cup\mathcal{S}_{m}$. The speed
planning is safe if $\forall
t_{0}\in[0,T],s\left(t_{0}\right)\in\mathcal{S},l\left(t_{0}\right)\in\mathcal{S}$,
which is equivalent to for
$t_{0}\in\left[T_{k},T_{k+1}\right],s\left(t_{0}\right)\in\mathcal{S}_{k},l\left(t_{0}\right)\in\mathcal{S}_{k},k=0,1,\ldots,m$,
i.e.,

$\underline{p_{0}^{k}}+h_{k}\underline{p_{1}^{k}}\frac{t_{0}-T_{k}}{h_{k}}\leq
s\left(t_{0}\right)\leq\overline{p_{0}^{k}}+h_{k}\overline{p_{1}^{k}}\frac{t_{0}-T_{k}}{h_{k}}$
$l_{beg}\leq l\left(t_{0}\right)\leq l_{end}$

###  IV-B Construction of Piecewise-Convex Safe Regions

(a) $S-L-T$ graph

(b) over-approximation (red)

(c) $S-T$ graph for $l\in[1,3)$

(d) $S-T$ graph for $l\in[3,6.7)$

Figure 2: Scenario in Fig. 1b after over-approximating car A

Algo. 1 outlines the 3D trapezoidal corridor generation process. The original
non-convex space is sliced along the $L$ axis at the starting or ending $L$
coordinates of any obstacles in the $S-L-T$ graph. This generates 3D chunks of
the non-convex space which can be projected in a 2D $S-T$ graph without the
loss of any search space. As an example, in Fig. 2a, any slices at $L$
coordinates in the range $[1,3)$ will give us the 2D $S-T$ cross-section as
seen in Fig. 2c. Similarly, any slice at an $L$ coordinate between $[3,6.7)$
will give us the $S-T$ graph as seen in Fig. 2d. The inputs to Algo. 1 are the
upper and lower bounds in the $S$ and $L$ direction w.r.t. the ego vehicle and
the road. The bounds are measured over a time horizon using a discrete time
interval $\Delta$. For each slice, we construct 2D convex corridors in the
corresponding $S-T$ graph. In this work, we extend the 2D convex trapezoidal
corridor generation algorithm in our previous work [16]. The new algorithm is
presented as Algo. 2.

Algo. 2 outlines the construction of 2D piecewise-convex safe regions in any
given $S-T$ cross-section along the $L$ axis. The lower and upper bounds in
the $S$ direction serve as inputs to this algorithm. We refer readers to [16]
for further details about the working of Algo. 2. A key modification in our
work is that in the subroutine $SingleRegionCaculate()$, we also initialize
the upper and lower boundaries of the regions in the $L$ direction (Algo. 3.
Lines 6,7). Our over-approximation step and the design of Algo. 1 guarantee
that these boundaries are the same for all 2D convex regions generated by
Algo. 2. Thus, we essentially get 2D trapezoidal-shaped corridors dragged
along the $L$ axis to form 3D trapezoidal prism-shaped convex corridors.
Finally, $RegionSplit()$ is used to check the length of each 2D convex region.
If it is above a user-defined threshold (e.g., $1\ s$ in our experimental
setting), it will be split into multiple sub-regions, for which the time
intervals are all below the threshold. This refinement operation aims to avoid
underfitting.

This 2D corridor generation process is repeated for all distinct obstacle
boundaries in our $S-L-T$ graph (Line 1, Algo. 1). The initialization of
bounds along the $L$ axis ensures that we get 3D trapezoidal-shaped convex
corridors. Since the length of the corridors in the $L$ direction is given by
the starting or ending of the obstacles in the $S-L-T$ space, we can guarantee
the safety of all the corridors generated using Algo. 1. Note that for the
space divided by obstacles, we select the unique space enclosing the reference
trajectory using the $SelectCorridors()$ method. For yielding to car A, the
corridors lying in the green region of Fig. 2c and 2d are chosen according to
the 3D reference waypoints. Note that once the corridors are selected, the
whole optimization is solved as a single problem and not decomposed into
individual optimizations for separate corridors.

Input: $obs,lb_{s}[obs],up_{s}[obs],lb_{l}[obs],up_{l}[obs],nums,\Delta$

1

Output: 3D_Region

2

3Initialize: for  _$i\leftarrow 0\textbf{ to }obs$_ do

4 $corridor=Convexify2D(lb_{s}[i],ub_{s}[i],lb_{l}[i],ub_{l}[i],nums,\Delta)$;

5 $new\\_corridor.append(corridor)$;

6 end for

7$final\\_corridor=SelectCorridors(new\\_corridor)$;

$\textbf{Return }final\\_corridor$

Algorithm 1 Piecewise 3D Convex Regions Generation

Input: $lb_{s},ub_{s},lb_{l},ub_{l},nums,\Delta$

1

Output: $regions$

2

$\textbf{Initialize: }regions[0],i=0,j=1$ /* i and j are counters for meta-
pieces and resulting convex regions, respectively */

3 $SingleRegionCalculate(region,0,lb_{s}[0],ub_{s}[0],lb_{s}[1],$

4 $ub_{s}[1],lb_{l},ub_{l})$

5 $regions.append(region)$

6 for  _$i\leftarrow 2\textbf{ to }nums-1$_ do

$lskew=(lb_{s}[i]-lb_{s}[i-1])/\Delta$ /* lower bound’s skew for two
consecutive meta-pieces */

$uskew=(ub_{s}[i]-ub_{s}[i-1])/\Delta$ /* upper bound’s skew */

7 if  _$||lskew-regions[j-1].lskew|| >\epsilon\textbf{ or }||uskew-
regions[j-1].uskew||>\epsilon$_ then

8 $regions[j-1].t_{end}=i$

9 $regions[j-1].t=(regions[j-1].t_{end}-regions[j-1].t_{beg})*\Delta$

10 $SingleRegionCalculate(region,j,lb_{s}[i],ub_{s}[i],$

11 $lb_{s}[i+1],ub_{s}[i+1],lb_{l},ub_{l})$

12 $regions.append(region)$

13 $j\leftarrow j+1$

14 end if

15

16 end for

17$regions[j-1].t_{end}=nums-1$

18 $regions[j-1].t=(regions[j-1].t_{end}-regions[j-1].t_{beg})*\Delta$

19 $RegionSplit(regions)$

Return $regions$

Algorithm 2 Convexify2D

Input: $region,j,lb_{s}[i],ub_{s}[i],lb_{s}[i+1],ub_{s}[i+1],lb_{l},ub_{l}$

Output: $region$

1 $region.t_{beg}=j$

2 $region.lskew=(lb_{s}[i+1]-lb_{s}[i])/\Delta$

3 $region.lbias=lb_{s}[i]$

4 $region.uskew=(ub_{s}[i+1]-ub_{s}[i])/\Delta$

5 $region.ubias=ub_{s}[i]$

6 $region.l_{beg}=lb_{l}$

7 $region.l_{end}=ub_{l}$

Algorithm 3 Single Region Caculate

##  V PIECEWISE BÉZIER POLYNOMIAL OPTIMIZATION

In this section, we discuss more about the limitations of using the cuboidal
corridors. We then discuss the safety enforcement in our trapezoidal prism-
shaped corridors. The formulation of quadratic optimization using the newly
designed convex solution space is introduced thereafter.

###  V-A Limitations of Safety Enforcement in Cuboidal Corridors

As discussed in Sec. III C, the convex hull property of the Bézier curves is
used to enforce that the trajectory in the $S-L-T$ graph stays in the safe
region $\mathcal{S}$. We first formally define a corridor for our trajectory
generation:

######  Definition 1.

Let the coefficients of the Bézier Polynomial be
$c_{i}\in\Omega,i=0,1,\ldots,n$. Each control point has two dimensions -
$\\{S,L\\}$. These control points lying in the safe region $\mathcal{S}$ form
a subset $\mathcal{S}^{cor}\subseteq\mathcal{S}$, which is called a corridor.

Ding et al. presented the construction of cuboidal corridors in the $S-L-T$
graph [5]. Constraints of the control points of cuboidal corridors are given
by the following proposition:

######  Proposition 1.

If a trajectory has control points in each time interval satisfying
$c^{k}_{i}\in\Omega^{k}_{cub}=\\{c^{k}|\underline{p_{0}^{k}}+h_{k}\underline{p_{1}^{k}}\leq
c^{k,s}\leq\overline{p_{0}^{k}},l^{k}_{beg}\leq c^{k,l}\leq
l^{k}_{end},i=0,1,\ldots,n,k=0,1,\ldots,m\\},f(t)$ is guaranteed to be safe,
and the upper bounds and lower bounds form cuboidal corridors
$\mathcal{S}^{cub}$.

The proof of safety enforcement in rectangular corridors can be found in [16],
and can be straightforwardly extended to the third dimension $L$ for cuboidal
corridors.

The optimization fails if any lower bound
($\underline{p}_{0}^{k}+h_{k}\underline{p_{1}^{k}}$) is greater than the upper
bound ($\overline{p_{0}^{k}}$). In order to avoid this, the time interval of
the $k$-th corridor must satisfy
$h_{k}\leq\frac{\overline{p_{0}^{k}}-p_{0}^{k}}{\underline{p_{1}^{k}}}$.

In [5], Ding et al. propose a seed generation and cube inflation method to
adjust time intervals. However, this method generates a significant number of
corridors and optimized variables, which leads to a high computation cost. In
common driving scenarios (Fig. 1), we have $\underline{p_{1}^{k}}>0$ or
$\overline{p_{1}^{k}}>0$, due to which the cuboidal corridors fail to cover
all the safe regions. As a result, the search space is sub-optimal and the
constraints on control points to enforce the trajectory are overtightened.

###  V-B Safety Enforcement in Trapezoidal-Prism Corridors

The sufficient conditions of control points $c_{i}$ to keep the longitudinal
and the lateral trajectory safe and in our proposed trapezoidal-prism
corridors are built upon the following lemma.

######  Lemma 1.

Let $M\in\mathbb{R}^{(n+1)\times(n+1)}$ denote the transition matrix from the
Bernstein basis
$\left\\{b_{n}^{0}(t),b_{n}^{1}(t),\ldots,b_{n}^{n}(t)\right\\}$ to the
monomial basis $\left\\{1,t,t^{2},\ldots,t^{n}\right\\}$. We have
$M_{i,0}=1,0\leq$ $M_{i,j}\leq 1,i=0,1,\ldots,n,j=0,1,\ldots,n$.

The proof can be found in [16]. We leverage the following theorem meant for 2D
trapezoidal corridors to construct 3D trapezoidal prism-shaped corridors.

######  Theorem 1.

For a trajectory, if it has control points in each time interval satisfying
$c_{i}^{k}\in\Omega^{k}$, where
$\Omega^{k}=\\{c^{k}|\underline{p_{0}^{k}}+h_{k}\underline{p_{1}^{k}}M_{i,1}\underline{p_{1}^{k}}\leq
c_{i}^{k,s}\leq\overline{p_{0}^{k}}+h_{k}\overline{p_{1}^{k}}M_{i,1},l^{k}_{beg}\leq
c_{i}^{k,l}\leq l^{k}_{end},i=0,1,\ldots,n,k=0,1,\ldots,m\\}$, $f(t)$ is
guaranteed to be safe. The upper and lower bounds in the $S$ and $L$
directions help form a trapezoidal prism-shaped corridor $\mathcal{S}^{trp}$.

The proof for the 2D case of the above theorem can be found in [16], and can
be easily extended to another dimension $L$.

In Theorem 1, conditions on $c_{i}^{s}$ are
$\underline{p_{0}^{k}}+h_{k}\underline{p_{1}^{k}}M_{i,1}\underline{p_{1}^{k}}\leq
c_{i}^{k,s}\leq\overline{p_{0}^{k}}+h_{k}\overline{p_{1}^{k}}M_{i,1}$.
Compared to the safety enforcement in cuboidal corridors in Proposition 1, we
have
$\underline{p_{0}^{k}}+h_{k}\underline{p_{1}^{k}}M_{i,1}\leq\underline{p_{0}^{k}}+h_{k}\underline{p_{1}^{k}}$
and
$\overline{p_{0}^{k}}+h_{k}\overline{p_{1}^{k}}M_{i,1}\geq\overline{p_{0}^{k}}$.
The advantage of having trapezoidal corridors is twofold: i) By the proof of
$\underline{p_{0}^{k}}+h_{k}\underline{p_{1}^{k}}M_{i,1}\underline{p_{1}^{k}}\leq
c_{i}^{k,s}\leq\overline{p_{0}^{k}}+h_{k}\overline{p_{1}^{k}}M_{i,1}$, the
lower boundaries are guaranteed to be smaller than the upper boundaries all
the time. Recall that for the rectangular corridors, we need to always check
$h_{k}\leq\frac{\overline{p_{0}^{k}}-p_{0}^{k}}{\underline{p_{1}^{k}}}$; ii)
The constraints are relaxed, therefore the solution space is enlarged compared
with the rectangular corridors (see the illustration for the comparison in
Fig. 3).

###  V-C Trajectory Optimization Formulation

The objective function is established as

$\displaystyle J=J_{s}+J_{l}$ (1) $\displaystyle
J_{s}=w_{1}\int_{0}^{T}\left(s(t)-s^{r}(t)\right)^{2}\mathrm{~{}d}t+w_{2}\int_{0}^{T}\left(\dot{s}(t)-v_{s}^{r}\right)^{2}\mathrm{~{}d}t$
$\displaystyle+w_{3}\int_{0}^{T}\ddot{s}(t)^{2}\mathrm{~{}d}t+w_{4}\int_{0}^{T}\dddot{s}(t)^{2}\mathrm{~{}d}t+w_{5}\left(s(T)-s^{r}(T)\right)^{2}$
$\displaystyle
J_{l}=w_{6}\int_{0}^{T}\left(l(t)-l^{r}(t)\right)^{2}\mathrm{~{}d}t+w_{7}\int_{0}^{T}\left(\dot{l}(t)-v_{l}^{r}\right)^{2}\mathrm{~{}d}t$
$\displaystyle+w_{8}\int_{0}^{T}\ddot{l}(t)^{2}\mathrm{~{}d}t+w_{9}\int_{0}^{T}\dddot{l}(t)^{2}\mathrm{~{}d}t+w_{10}\left(l(T)-l^{r}(T)\right)^{2}$

where $s^{r}(t)$ and $l^{r}(t)$ are the reference longitudinal and lateral
trajectories, and $v^{r}_{s}$ and $v^{r}_{l}$ are the reference velocities in
the two directions. For $J_{s}$ and $J_{l}$, their first terms penalize the
deviation from the reference; the second ones penalize the deviation between
the actual and reference speed; the third and the fourth terms penalize
acceleration and jerk, respectively. The last terms penalize the deviation of
the ending position from the reference. We used Optuna[22] for tuning all the
10 parameters.

The optimization considers the following constraints:

  * •

Boundary Constraints: The piecewise curve starts from fixed position, speed,
and acceleration, i.e.,

$c_{i}^{0,l}h_{k}^{(1-l)}=\left.\frac{\mathrm{d}^{l}f(t)}{\mathrm{d}t^{l}}\right|_{t=0},\quad
l=0,1,2,$

where $c_{i}^{k,l}$ is the control point for the $l$th-order derivative of the
$k$-th Bézier curve. Note that $c_{i}^{k,l}$ has two dimensions: $\\{S,L\\}$.

  * •

Continuity Constraints: The piecewise curve must be continuous at the
connected time points for position, speed, and acceleration.

$c_{n}^{k,l}h_{k}^{(1-l)}=c_{0}^{k+1,l}h_{k+1}^{(1-l)},l=0,1,2,k=0,1,\ldots,m-1$

  * •

Safety Constraints: With our proposed trapezoidal-prism corridors, safety
constraints for the longitudinal dimension of the control point can be given
as

$\underline{p_{0}^{k}}+h_{k}\underline{p_{1}^{k}}M_{i,1}\leq
c_{i}^{k,0}\leq\overline{p_{0}^{k}}+h_{k}\overline{p_{1}^{k}}M_{i,1},k=0,1,\ldots,m$

and those for the lateral dimension of the control point can be given as

$l_{beg}\leq c_{i}^{k,0}\leq l_{end}$

  * •

Physical Constraints: The physical constraints under consideration include the
limit of a vehicle’s velocity, acceleration, and jerk. We can use the
hodograph property of a Bézier curve to calculate velocity, acceleration, and
jerk. The constraints are given by

$\displaystyle\beta^{k,1}$ $\displaystyle\leq
c_{i}^{k,l}\leq\overline{\beta^{k,1}}$ $\displaystyle\beta^{l}$
$\displaystyle\leq c_{i}^{k,l}\leq\overline{\beta^{l}},l=2,3$

where $k=0,1,\ldots,m$ and it follows that $c_{i}^{k,l+1}=(n-$
$l)\left(c_{i+1}^{k,l}-c_{i}^{k,l}\right)$. The upper bounds
$\overline{\beta^{k,1}}$ are determined by speed limits on road and
centripetal acceleration constraints. Let $a_{cm}$ be the maximum acceleration
permitted and $\kappa_{k}$ the maximum curvature of the path for
$t\in\left[T_{k},T_{k+1}\right]$ (see [23] for details). The lateral
acceleration constraints are given by

$c_{i}^{k,l}\leq\overline{\beta^{k,1}}=\sqrt{\frac{a_{cm}}{\kappa_{k}}}.$

The bounds on longitudinal and lateral accelerations and jerks are constant
for different pieces of speed profiles.

Then, the trajectory optimization process can be formulated as a quadratic
programming (QP) problem as

$\displaystyle\mathbf{P}:\quad\min_{\mathbf{c}}$
$\displaystyle\frac{1}{2}\mathbf{c}^{T}\mathbf{Q}_{\mathbf{c}}\mathbf{c}+\mathbf{q}_{\mathbf{c}}^{T}\mathbf{c}+\mathrm{const}$
s.t. $\displaystyle\mathbf{A}_{eq}\mathbf{c}=\mathbf{b}_{eq}$
$\displaystyle\mathbf{A}_{ie}\mathbf{c}\leq\mathbf{b}_{ie}.$

We refer readers to the appendix of our previous work [24] for the detailed
formulation process. This problem can be solved in real-time by a modern
solver such as OSQP[25].

##  VI SIMULATIONS AND RESULTS ANALYSIS

Our framework has been implemented using C++11. All simulations are carried
out on a personal computer with a 2.60GHz Intel i10-10750H processor.

###  VI-A Numerical Simulations

We conduct numerical simulations to compare the proposed approach with
cuboidal corridors [5]. The planning horizon is 7 s. Different road scenarios
are as follows:

####  VI-A1 Merging into another lane due to road construction

Consider the scenario in Fig. 1. We project different stations of the vehicles
onto the $S-L-T$ graph. The initial velocity and the acceleration of the ego
vehicle are $v_{s}(0)=7.0\ \textnormal{m/s},v_{l}(0)=0\ \textnormal{m/s}$ and
$a_{s}(0)=0\ \textnormal{m/s\textsuperscript{2}},a_{l}(0)=0\
\textnormal{m/s}\textsuperscript{2}$, respectively.

(a)

(b)  (c) Longitudinal Acceleration profiles (d) Lateral Acceleration profiles
(f) Lateral Velocity profiles Figure 3: Piecewise Bézier polynomial and its
dynamic profile

Generated trajectory using both methods (green: trapezoidal, red: cuboidal)

(e) Longitudinal Velocity profiles

Fig. 3a and Fig. 3b show Bézier curves generated by cuboidal (red, [5]) and
trapezoidal-prism (green, ours) corridors for the scenario presented in Fig.
1. From Fig. 3c, we observe that the maximum acceleration required for our
method is less than that needed by the cuboidal corridors approach. The
superiority of using trapezoidal corridors is more clear from Fig. 3d, which
records the lateral acceleration of both the methods. We observe that our
approach yields a smoother acceleration plot with minimal jerk and the lower
maximum acceleration. We also test the maximum initial conditions of both the
methods for the same scenario to show the effect of the enlarged search space.
While using trapezoidal corridors, we can generate a trajectory for $a_{s}=2$
m/s2, $v_{s}=10.5$ m/s, $a_{l}=1.2$ m/s2, $v_{l}=2$ m/s where the bounds on
longitudinal acceleration are $[-3,2]$ m/s2 and those on lateral acceleration
are $[-2,2]$ m/s2. Using the cuboidal corridors fails to generate a trajectory
for these initial conditions and is only successful when the initial velocity
in the longitudinal direction is reduced to 9 m/s.

####  VI-A2 Overtaking a low-speed vehicle in front

We test our planner on overtaking a slowly moving car in front by lane
changing twice (second time to merge back into the original lane of the ego
vehicle). In layered planning techniques, these kinds of scenarios are
typically tackled by considering the obstacle to be static for a few seconds.
Hence, this approach proves to be conservative. The differences in the
longitudinal acceleration graphs between the two corridor generation
techniques can be seen in Fig. 4. Clearly, using the trapezoidal corridors
generates a trajectory with much lower acceleration. Here, the vehicle in
front is assumed to be moving with $v_{s}=5$ m/s and the ego vehicle’s initial
condition is $v_{s}=7$ m/s.

(a) Longitudinal acceleration

(b) matplotlib animation

Figure 4: Overtaking scenario

####  VI-A3 Unprotected left turn

As shown in Fig. 5a, there are two cars coming from the front which obstruct
the ego vehicle from making a left turn without yielding to them. As seen in
Fig. 5b, our planner can successfully find a trajectory while meeting all the
safety and dynamic feasibility constraints. Since the ego vehicle needs to
yield to the cars in front, we also test the maximum initial velocity
($v_{s}=1$ m/s) and acceleration ($a_{s}=0.5$ m/s2) in the longitudinal
direction for this case. If the distance between Car A and Car B is sufficient
for the ego vehicle to go in between them, our planner finds the corresponding
trajectory (Fig. 5c). In this scenario, the additional search space obtained
by trapezoidal corridors is not used at all, as the trajectory passing through
the enlarged search space can only result in a lane change, which is not
desired. Hence, both the trajectories obtained are the same and overlap each
other, as seen in Fig. 5.

(a) ego vehicle is shown in blue with the corresponding reference trajectory

(b)

(c)

Figure 5: Unprotected Left Turn at an Intersection- trajectories obtained from
both corridor construction methods (green: trapezoidal, red: cuboidal) are
identical

###  VI-B CommonRoad Simulations

The simulations in this part are conducted on the CommonRoad platform [26],
which provides an interactive simulated and non-interactive real traffic
environment for validating motion planning algorithms. A given scenario is
considered “solved” when the ego vehicle reaches the desired goal region while
satisfying all the constraints. We visualize the bird’s-eye view simulation of
the lane change scenario in Fig. 1. The results obtained using our approach
can be seen in Fig. 6. The cuboidal corridor approach did not yield a
collision-free trajectory as it failed to replan owing to the lack of search
space. It also had significantly high acceleration as can be observed in Fig.
3c.

(a) $t=0\ s$

(b) $t=2\ s$

(c) $t=4\ s$

(d) $t=7\ s$

(e) $t=9\ s$

(f) $t=12\ s$

Figure 6: CommonRoad Simulation for “Merging into another lane due to road
construction” scenario. The ego vehicle (initial location is represented by
the green arrow) is shown in green, and the other cars are shown in blue. The
construction site (static obstacle) is shown in red.

##  VII CONCLUSION

In this paper, we propose a novel convexification algorithm for generating
safety corridors in the $S-L-T$ space. We show that our method of trapezoidal
prism-shaped corridors enlarges the solution space as compared to the existing
cuboidal corridors-based method. We provide the sufficient conditions of
control points in the trapezoidal corridors to provably guarantee the safety
of trajectories represented by Bézier polynomials. Finally, we formulate the
trajectory optimization as a QP problem. The numerical and CommonRoad
simulations show that the proposed approach is superior in terms of optimality
and low failure rates. Future work includes using a dynamic programming-based
approach to generate a comfort-optimal reference trajectory in the $S-L-T$
space.

## References

  * [1] M. Werling, S. Kammel, J. Ziegler, and L. Gröll, “Optimal trajectories for time-critical street scenarios using discretized terminal manifolds,” _The International Journal of Robotics Research_ , vol. 31, no. 3, pp. 346–359, 2012.
  * [2] T. Gu, J. M. Dolan, and J.-W. Lee, “Runtime-bounded tunable motion planning for autonomous driving,” in _2016 IEEE Intelligent Vehicles Symposium (IV)_. IEEE, 2016, pp. 1301–1306.
  * [3] W. Xu, J. Wei, J. M. Dolan, H. Zhao, and H. Zha, “A real-time motion planner with trajectory optimization for autonomous vehicles,” in _2012 IEEE International Conference on Robotics and Automation_. IEEE, 2012, pp. 2061–2067.
  * [4] X. Li, Z. Sun, D. Cao, Z. He, and Q. Zhu, “Real-time trajectory planning for autonomous urban driving: Framework, algorithms, and verifications,” _IEEE/ASME Transactions on mechatronics_ , vol. 21, no. 2, pp. 740–753, 2015\.
  * [5] W. Ding, L. Zhang, J. Chen, and S. Shen, “Safe trajectory generation for complex urban environments using spatio-temporal semantic corridor,” pp. 2997–3004, 2019.
  * [6] T. Mercy, W. Van Loock, and G. Pipeleers, “Real-time motion planning in the presence of moving obstacles,” in _2016 European Control Conference (ECC)_. IEEE, 2016, pp. 1586–1591.
  * [7] J. Ziegler, P. Bender, T. Dang, and C. Stiller, “Trajectory planning for Bertha—A local, continuous method,” in _2014 IEEE intelligent vehicles symposium proceedings_. IEEE, 2014, pp. 450–457.
  * [8] H. Fan, F. Zhu, C. Liu, L. Zhang, L. Zhuang, D. Li, W. Zhu, J. Hu, H. Li, and Q. Kong, “Baidu apollo EM motion planner,” _arXiv preprint arXiv:1807.08048_ , 2018.
  * [9] J. Zhou, R. He, Y. Wang, S. Jiang, Z. Zhu, J. Hu, J. Miao, and Q. Luo, “Autonomous driving trajectory optimization with dual-loop iterative anchoring path smoothing and piecewise-jerk speed optimization,” _IEEE Robotics and Automation Letters_ , vol. 6, no. 2, pp. 439–446, 2020.
  * [10] C. Liu, W. Zhan, and M. Tomizuka, “Speed profile planning in dynamic environments via temporal optimization,” in _2017 IEEE Intelligent Vehicles Symposium (IV)_. IEEE, 2017, pp. 154–159.
  * [11] X. Qian, I. Navarro, A. de La Fortelle, and F. Moutarde, “Motion planning for urban autonomous driving using bézier curves and mpc,” in _2016 IEEE 19th international conference on intelligent transportation systems (ITSC)_. Ieee, 2016, pp. 826–833.
  * [12] Y. Pan, Q. Lin, H. Shah, and J. M. Dolan, “Safe planning for self-driving via adaptive constrained ilqr,” in _2020 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS)_. IEEE, 2020, pp. 2377–2383.
  * [13] Z. Zhu, E. Schmerling, and M. Pavone, “A convex optimization approach to smooth trajectories for motion planning with car-like robots,” in _2015 54th IEEE conference on decision and control (CDC)_. IEEE, 2015, pp. 835–842.
  * [14] S. M. Erlien, S. Fujita, and J. C. Gerdes, “Safe driving envelopes for shared control of ground vehicles,” _IFAC Proceedings Volumes_ , vol. 46, no. 21, pp. 831–836, 2013.
  * [15] C. Liu, C.-Y. Lin, and M. Tomizuka, “The convex feasible set algorithm for real time optimization in motion planning,” _SIAM Journal on Control and optimization_ , vol. 56, no. 4, pp. 2712–2733, 2018.
  * [16] J. Li, X. Xie, J. He, Q. Lin, and J. Dolan, “Motion planning by search in derivative space and convex optimization with enlarged solution space,” in _2022 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS)_ , 07 2022.
  * [17] W. Zhang, P. Yadmellat, and Z. Gao, “A sufficient condition for convex hull property in general convex spatio-temporal corridors,” in _2022 IEEE Intelligent Vehicles Symposium (IV)_. IEEE, 2022, pp. 1033–1039.
  * [18] W. Xu and J. M. Dolan, “Speed planning in dynamic environments over a fixed path for autonomous vehicles,” in _2022 International Conference on Robotics and Automation (ICRA)_. IEEE, 2022, pp. 3321–3327.
  * [19] D. González, V. Milanés, J. Pérez, and F. Nashashibi, “Speed profile generation based on quintic bézier curves for enhanced passenger comfort,” in _2016 IEEE 19th international conference on intelligent transportation systems (ITSC)_. IEEE, 2016, pp. 814–819.
  * [20] F. Gao, W. Wu, Y. Lin, and S. Shen, “Online safe trajectory generation for quadrotors using fast marching method and bernstein basis polynomial,” in _2018 IEEE International Conference on Robotics and Automation (ICRA)_ , pp. 344–351.
  * [21] B. Zhou, F. Gao, L. Wang, C. Liu, and S. Shen, “Robust and efficient quadrotor trajectory generation for fast autonomous flight,” _IEEE Robotics and Automation Letters_ , vol. 4, no. 4, pp. 3529–3536, 2019.
  * [22] T. Akiba, S. Sano, T. Yanase, T. Ohta, and M. Koyama, “Optuna: A next-generation hyperparameter optimization framework,” in _Proceedings of the 25rd ACM SIGKDD International Conference on Knowledge Discovery and Data Mining_ , 2019.
  * [23] Y. Zhang, H. Sun, J. Zhou, J. Hu, and J. Miao, “Optimal trajectory generation for autonomous vehicles under centripetal acceleration constraints for in-lane driving scenarios,” in _2019 IEEE Intelligent Transportation Systems Conference (ITSC)_ , 2019, pp. 3619–3626.
  * [24] J. L. et al., “Motion Planning by Search in Derivative Space and Convex Optimization with Enlarged Solution Space,” https://github.com/JialunLi/Derivate-Space-Search-Supplement-Material-, 2022, [extended version].
  * [25] B. Stellato, G. Banjac, P. Goulart, A. Bemporad, and S. Boyd, “OSQP: an operator splitting solver for quadratic programs,” _Mathematical Programming Computation_ , vol. 12, no. 4, pp. 637–672, 2020.
  * [26] M. Althoff, M. Koschi, and S. Manzinger, “Commonroad: Composable benchmarks for motion planning on roads,” in _2017 IEEE Intelligent Vehicles Symposium (IV)_ , 2017, pp. 719–726.