# Straintronics using the monolayer-Xene platform - a comparative study

Swastik Sahoo  Department of Electrical Engineering, Indian Institute of
Technology Bombay, Powai, Mumbai-400076, India  Namitha Anna Koshi  Indo-Korea
Science and Technology Center, Bengaluru-560064, India  Seung-Cheol Lee
Electronic Materials Research Center, KIST, Seoul 136-791, South Korea
Satadeep Bhattacharjee  Indo-Korea Science and Technology Center,
Bengaluru-560064, India  Bhaskaran Muralidharan<EMAIL_ADDRESS>Department
of Electrical Engineering, Indian Institute of Technology Bombay, Powai,
Mumbai-400076, India

###### Abstract

Monolayer silicene is a front runner in the 2D-Xene family, which also
comprises germanene, stanene, and phosphorene, to name a few, due to its
compatibility with current silicon fabrication technology. Here, we
investigate the utility of 2D-Xenes for straintronics using the ab-initio
density functional theory coupled with quantum transport based on the Landauer
formalism. With a rigorous bandstructure analysis, we show the effect of
strain on the K-point, and calculate the directional piezoresistances for the
buckled Xenes as per their critical strain limit. Further, we compare the
relevant gauge factors, and their sinusoidal dependences on the transport
angle akin to silicene and graphene. The strain-insensitive transport angles
corresponding to the zero gauge factors are $81^{\circ}$ and $34^{\circ}$ for
armchair and zigzag strains, respectively, for silicene and germanene. For
stanene as the strain limit is extended to $10\%$ and notable changes in the
fundamental parameters, the critical angle for stanene along armchair and
zigzag directions are $69^{\circ}$ and $34^{\circ}$ respectively. The small
values of gauge factors are attributed to their stable Dirac cones and strain-
independent valley degeneracies. We also explore conductance modulation, which
is quantized in nature and exhibits a similar pattern with other transport
parameters against a change in strain. Based on the obtained results, we
propose the buckled Xenes as an interconnect in flexible electronics and are
promising candidates for various applications in straintronics.

##  I Introduction

The prospect of exploiting the fascinating properties of graphene has
initiated a paradigm shift towards the expeditious development of its
analogous elemental monolayer sheets such as silicene [1], germanene [2],
stanene [3, 4], phosphorene [5, 6], arsenene [7] to name a few. The robust
electrical [8, 9], mechanical [10, 11, 12], and chemical [13, 14] properties
of graphene are touted to be the founding stone for investigation of Xenes as
mentioned above and for the materials among the group $III$ and $V$ elements
as well, which in due course will act as strong contenders for functional
nano-devices [15]. Group $IV$ materials form a buckled structure [16] whereas
group $V$ materials can exist both in buckled and puckered structures [15] as
compared to the planar structure of graphene.
Silicene, a monolayer of silicon (a group $IV$ element), is a prime contender
for various applications due of its consonance with the already certified
silicon industry [17]. Adding to this advantage, silicene also possesses
useful electro-mechanical properties [18, 19]. Germanene, the $2D$ counterpart
of germanium, is also endowed with some extra-ordinary properties like
silicene and they also possess a strong spin-orbit coupling [20, 21], giant
magneto-resistance [22], and a tunable bandgap due to applied electric fields
[23]. Like graphene, silicene, and germanene exhibit linear energy dispersion
relation near the Dirac points [24, 25], zero band-gap [26], and more
importantly, they are dynamically stable [25]. Among the $2D$ materials,
stanene excels with outstanding properties and potential applications, though
most of its electronic and structural properties are similar to other group
$IV$ counterparts. In group $V$ associates, black phosphorous is the most
stable crystal structure among the allotropes which include red, white, and
violet phosphorous [27, 28]. This forms two-dimensional layers which are held
together by weak van der Waals force [29] resulting in an orthorhombic
puckered structure with a honeycomb like lattice arrangement in each layer
[30, 31]. This material can be singularized from other competitors due to its
unique properties catering to various prospective listed as direct band gap
from bulk to monolayer [32] useful for optoelectronics [33], enhanced carrier
transport mobility [5, 6] for nanoelectronics [34], high on/off ratio for
transistor applications [5, 6], in-plane anisotropy [35] related to electrical
and optical properties, to name a few.
In spite of the above mentioned advantages, their respective applications in
straintronics remain largely unexplored. In this paper, we delve into the
strained lattice structure of 2D-Xenes monolayers in the quasi-ballistic
regime, which corresponds to a length-scale of around $100$ nm - $200$ nm)
[36]. The corresponding results are compared keeping in mind different nano
electro-mechanical systems (NEMS) applications and possible flexible
electronics applications, which are quite diverse, satisfying the fundamental
criterion of robust electronic and excellent mechanical response to strain
[37], followed by portability, and manufacturability [38].
The central goal of this paper is to compare silicene, germanene, stanene, and
phosphorene, based on their fundamental properties, band structures, and
transport properties in the presence of strain and explore different
applications of Xenes like piezoresistivity, conductance modulation, etc.,
using ab-initio density functional theory (DFT) and quantum transport theory.
Specifically, we calculate the directional piezoresistance for different
strain values along the armchair and zigzag directions. We typically obtaine a
smaller value of the directional piezoresistances and their sinusoidal
dependence. We show that the value of the gauge factor increases as we go up
for group $IV$ elements along their atomic numbers in the periodic table. The
strain-insensitive transport angles corresponding to the zero gauge factors
are $81^{\circ}$ and $34^{\circ}$ for armchair and zigzag strains,
respectively, for silicene and germanene. For stanene as the strain limit is
extended to $10\%$ and notable changes in the fundamental parameters, the
critical angle for stanene along armchair and zigzag directions are
$69^{\circ}$ and $34^{\circ}$ respectively. The small gauge factor of buckled
Xenes can be attributed to its robust Dirac cone and strain-independent valley
degeneracy. Based on the obtained results, we propose the buckled Xenes as an
interconnect in flexible electronics and are promising candidates for various
applications in straintronics.
The remainder of this paper is organized as follows. Section II represents a
briefing about the theory and fundamental properties of the Xenes. Also, the
tight binding parameters and band structures for the Xenes at different values
of strains are displayed. Section III illustrates the methodologies to obtain
the results. The depiction of these methods will go in parallel and
corresponding applications have been discussed in detail. Results will be
connoted forthwith, followed by a conclusion in section IV.

##  II Theory and properties

Figure 1: (a) Top view and side view of monolayer buckled Xene crystal lattice
applicable to silicene, germanene, and stanene with
$\overrightarrow{a_{1}^{s}}$ and $\overrightarrow{a_{2}^{s}}$ as the primitive
vectors, and $\delta$ as the buckling constant. A and B represent the non-co-
planar sub-lattices. ${t_{1}}$, ${t_{2}}$ and ${t_{3}}$ represent the nearest
neighbour tight binding parameters.(b) Schematic diagram depicting the high
symmetry path ${\Gamma MK\Gamma}$ in the $1^{st}$ Brillouin zone of buckled
strained Xene sheets. (c) Schematic presenting the setup and electron
transport in Xene sheets. (d) Side view, (e) Perspective view, (f) Top view of
phosphorene lattice. Blue- and -green-filled circles represent puckered-up and
puckered-down phosphorous atoms, respectively. $t_{1}$ to $t_{5}$ represent
the hopping parameters. (g) Sketch showing the high symmetry path ${Y\Gamma
XS}$ in the $1^{st}$ Brillouin zone of phosphorene.

In this section, we summarize the fundamental properties of Xenes in
comparative form. However, we highlight the commonalities in the dependent and
versatile properties of Xenes in the subsequent subsections.

The lattice structure and their corresponding reciprocal symmetry, i.e.,
Brillouin zone, is shown in Fig. 1. Fig. 1a, 1b and 1c are applicable for
buckled Xenes corresponding to silicene, germanene, and stanene, respectively.
The buckled honeycomb structure of Xenes consists of two triangular sub-
lattices, shown in Fig. 1a by red (denoted by A) and blue (denoted by B) dots,
respectively. These sub-lattices are non-co-planar, which gives the Xenes a
buckled honeycomb structure. The buckling constant $\delta$ differs for each
Xene and is enlisted in table 1 for reference. $t_{1}$, $t_{2}$, and $t_{3}$
represent nearest neighbor tight-binding (TB) parameters for buckled Xenes.
Figure 1b depicts the first-Brillouin zone of Xenes in reciprocal space. For
uniaxially strained Xenes, the high symmetry path is given by ${\Gamma
MK\Gamma}$. Using ab-initio calculations, we obtain the band structures of
strained Xene sheets along this path in the linear elastic regime (see Fig.
3). The simulation setup for the calculation of gauge factor for all the
monolayer buckled Xenes along different transport angles ($\theta$) is shown
in Fig. 1c. It consists of a buckled Xene sheet placed between two contacts, a
source, and a drain. A uniaxial tensile strain (s) is applied along the
armchair (AC) and zigzag (ZZ) direction varying from $0\%$ to $10\%$ depending
on the material. The resistance is measured along each transport direction‘
$\theta$’ ($0$∘ to $90$∘) for applied voltage ranging between $\pm 10$ mV. The
electron transport direction is shown as a consequence of this process. The
device dimension along the transport direction is nearly $100$ nm. Hence, the
electrons undergo quasi-ballistic transport.
The in-plane geometry of the black phosphorous crystal lattice is shown in
Fig. 1d, 1e and 1f. The side view and $3D$ view of the crystal lattice are
shown in Fig. 1d and 1e, respectively. Fig. 1f represents the alternate
stacking of phosphorous atoms along the AC direction, giving rise to the
puckered structure [39]. The five hopping parameters, $t_{1}$, $t_{2}$,
$t_{3}$, $t_{4}$, and $t_{5}$, are needed to describe phosphorene’s low energy
electronic structure within the TB model. $t_{1}$ and $t_{2}$ are the
prominent nearest neighbor hopping parameters and $t_{4}$ accounts for the
symmetricity between the conduction band and valence band. $t_{3}$ and $t_{5}$
are the minor parameters, used for dispersion correction [40]. The rectangle
on $k_{x}$ and $k_{y}$ plane represents the two-dimensional Brillouin zone of
the single puckered layer (see Fig. 1g). The high symmetry path is marked as
${Y\Gamma XS}$, and it is maintained throughout the calculation.
A uniaxial tensile strain is applied along both AC and ZZ directions of
buckled Xenes, which can be represented as [41]:

$s=s\begin{bmatrix}cos^{2}(\theta)-sin^{2}(\theta)&(1+\rho)cos(\theta)sin(\theta)\\\
(1+\rho)cos(\theta)sin(\theta)&sin^{2}(\theta)-cos^{2}(\theta)\end{bmatrix}$
(1)

Here, $\theta$ represents the transport angle, and it’s default value is
$0^{\circ}$ for ZZ strain and $90^{\circ}$ for AC strain. In this work, we
have varied the transport angle from $0^{\circ}$ to $90^{\circ}$, and
accordingly, the respective transport parameters are calculated.

###  II.1 Fundamental properties

We have listed the parameters related to the electronic and mechanical
properties of all the Xenes in tables 1 and 2 subsequently. More parameters
could have been documented, but we have kept only the relevant ones, keeping
in mind the aspect of the paper. Graphene is not integral to this work but is
listed as a reference material. Also, the dispersion relations for Xenes are
expressed below.

Table 1: Comparison of Lattice constant (a) in ($\AA$), buckling height ($\delta$), bond length (d) in ($\AA$), Spin-orbit coupling (SOC) in $meV$, Relaxation time near Dirac cone $(\tau_{d})$ in $ps$ , electron mobility ($\mu_{e}$) in $10^{5}cm^{2}V^{-1}s^{-1}$, hole mobility ($\mu_{h}$) in $10^{5}cm^{2}V^{-1}s^{-1}$ at $300$K for graphene, silicene, germanene, stanene, and phosphorene respectively.  Material | a | $\delta$ | d | SOC | $\tau_{d}$ | $\mu_{e}$ | $\mu_{h}$
---|---|---|---|---|---|---|---
Graphene | 2.46 [42] | 0 [9] | 1.42 [9] | 0.05 [43] | 2.24 [44] | 3.2 [44] | 3.51 [44]
Silicene | 3.83 [45] | 0.44 [46, 25] | 2.25 [25, 45] | 1.55 [20, 47] | 1.84 [44] | 2.57 [44] | 2.22 [44]
Germanene | 4.02 [20] | 0.64 [25, 45] | 2.38 [25, 45] | 23.9 [20, 47] | 5.26 [48] | 6.09 [48] | 6.39 [48]
Stanene | 4.67 [49] | 0.86 [49] | 2.69 [49] | 73 [47] | - | 27.7 [50] | 40.1 [50]
Phosphorene | 4.53 [6, 49] | 5.00 [6, 49] | 3.36 [6] | 20 [51] | - | 0.011 [52] | 0.007 [52]
Table 2: Comparison of elastic limit, isotropic limit, Poisson ratio $(\rho)$, Fermi velocity at Dirac point $(V_{f})$ in $(10^{5}ms^{-1})$ at $300$K for graphene, silicene, germanene, stanene, and phosphorene respectively. Material | Elastic limit | Isotropic region | $\rho$ | $V_{f}$
---|---|---|---|---
Graphene |  $20\%$ (AC) and $26\%$ (ZZ) [11, 53] |  $10\%$ (AC and ZZ) [41] | 0.14 [9] | 6.3-14.4 [54, 50]
Silicene |  $15\%$ (AC) and $16\%$ (ZZ) [55] |  $5\%$ (AC and ZZ) [56] | 0.31 [55] | 5.1 [54]
Germanene |  $16\%$ (AC) and $17\%$ (ZZ) [57] |  $5\%$ (AC and ZZ [58] | 0.33 [59] | 3.8 [54]
Stanene |  $18\%$ (AC and ZZ) [60] |  $10\%$ (AC and ZZ) [61] | 0.39 [59] | 8.3 [50]
Phosphorene |  $30\%$ (AC) and $27\%$ (ZZ) [62] |  $10\%$ (AC and ZZ) [63] | 0.20 (AC) and 0.70 (ZZ) [62] | 3.5 [64]

####  II.1.1 Buckled structured Xenes

In their most stable form, silicene, germanene, and stanene form buckled
hexagonal honeycomb structure [16] due to their larger bond length 1 compared
to planar-like graphene having smaller bandgap. This large bond length makes
the $\pi$ bond weaker leading to deviations from $sp^{2}$ hybridization. The
energy dispersion relation applicable to all buckled Xenes is depicted as
follows:

$E(k)=\pm|t_{2}+t_{1}e^{-(i\vec{k}\cdot\overrightarrow{a_{1}^{s}})}+t_{3}e^{-(i\vec{k}\cdot\overrightarrow{a_{2}^{s}})}|.$
(2)

Here, $\overrightarrow{a_{1}^{s}}$ and $\overrightarrow{a_{2}^{s}}$ represents
the primitive lattice vectors, shown in Fig. 1a along the nearest neighbour
tight binding parameters ($t_{1}$ to $t_{3}$).

####  II.1.2 Puckered structured Xene

Phosphorene, on the other hand, forms a puckered structure ascribed to the
covalent bond between each atom and three adjacent ones. This bond induces a
strong anisotropy in phosphorene resulting in modulation of electronic,
magnetic, optical, and transport properties [65, 66]. Black phosphorous exists
in four crystal structures: orthogonal, rhombic, simple cubic, and amorphous
[67]. At room temperature, phosphorene has an orthorhombic crystal structure,
as shown in Fig. 1.
The energy dispersion relation applicable to puckered Xenes [68] is
represented as:

$E(k)=f(k)\pm g(k).$ (3)

where, $f(k)$ given by $f(k)=4t_{4}cos({k_{x}a}/{2})cos({k_{y}b}/{2})$ and the
other component, $g(k)$ is depicted as:
$g(k)=2t_{1}e^{-k_{x}a_{1x}}cos({k_{y}b}/{2})+t_{2}e^{k_{x}a_{2x}}+2e^{k_{x}a_{3x}}cos({k_{y}b}/{2})+t_{5}e^{k_{x}a_{5x}}$
. Here, $+(-)$ represents conduction (valence) band.
The coefficients are given by $a_{1x}=1.41763$ $\AA$ and $a_{2x}=0.79732$
$\AA$. These represent the intraplanar and interplanar distance for nearest-
neighbour atoms along the x-direction, respectively. The other coefficients
are expressed as: $a_{3x}=a_{1x}+2a_{2x}$, $a_{4x}=a_{1x}+a_{2x}$ and
$a_{5x}=2a_{1x}+a_{2x}$. These parameters are obtained from the lattice
structure of phosphorene [68] and $t_{1}$ to $t_{5}$ represent TB hopping
parameters.

###  II.2 Computational details

The optimized band structure of silicene, germanene, stanene, and phosphorene
are calculated from DFT using Vienna ab-initio Simulation Package (VASP) [69,
70, 71]. The formalism of Perdew-Burke-Ernzerhof (PBE) [72] within generalized
gradient approximation (GGA) is used to treat exchange-correlation
interactions. The buckled silicene, germanene, and stanene monolayers belong
to the crystallographic space group $P-3m1(164)$. A vacuum separation of
$15(\AA)$ is employed for germanene and stanene. Black phosphorene
crystallizes in a puckered structure with space group $Pmna(53)$ different
from the hexagonal Xenes. In the case of phosphorene, a vacuum of thickness
$18(\AA)$ is used in our calculations. Plane-wave cut-off energies within VASP
are set to $500$ eV. Structures are optimized until each atom’s residual force
becomes less than $0.01$ eV/(Å). For geometry optimization of germanene and
stanene, the Brillouin zone is sampled using a $41\times 41\times 1$ and
$31\times 41\times 1$ k-mesh for 2-atom unit cell and 4-atom $1\times 1$
rectangular super cell respectively. Similarly, for a monolayer of black
phosphorus, a $31\times 41\times 1$ k-mesh is used for Brillouin zone
sampling. The van der Waals correction is considered by considering Grimme’s
semi-empirical DFT-D$2$ method [73]. For germanene and stanene, a $1\times 1$
rectangular supercell consisting of 4 atoms is constructed to study the effect
of uniaxial strain on their electronic structures. Further calculation details
are given in our previous work on silicene in [74].
The strain range for different Xenes is chosen differently depending on their
elastic limit and electronic properties, as listed in table 2. The critical
strain for each Xene can also be calculated quantitatively [75], based on the
continuum mechanics model. The formula is given as follows:

$S_{cr}=\frac{\pi^{2}t^{2}}{3(1-\rho^{2})W^{2}(\rho-(\frac{W}{2L})^{2})}$ (4)

Here, $t$ is the thickness of the material, $L$ and $W$ represent the two-
dimensional length and width of the material. $\rho$ represents Poisson ratio.
In case of silicene, by taking its thickness to be 0.313 nm [55] and the
values of length, width, and Poisson ratio from [74], critical strain obtained
is around $8\%$. In addition to this, there is a Dirac material-to-metal
transition for a uniaxial tensile strain of $7.5\%$ for silicene [55]. So, the
strain is restricted up to $7\%$ for silicene in our study. Similarly, for
germanene, the study is restricted to a strain of $5\%$, attributed to the
observed Dirac material-to-metal transition [76, 77] and for stanene, upto
$10\%$, considering the range of isotropic elastic limit [78].

Figure 2: Variation of the nearest neighbor tight-binding parameters as a
function of strain along the (a),(c),(e) armchair, and (b), (d), (f) zigzag
directions for silicene, germanene, and stanene, respectively. Figure 3: Band
structure of a) silicene, b) germanene, c) stanene along the path $\Gamma
MK\Gamma$ for AC strain adhering to the values of $0\%$, $2\%$ and $5\%$
respectively. The magnified view of the point (orange encircled) shows the
shifting of $K$ point away from the $K-\Gamma$ direction for each figure.

###  II.3 Portrayal of hopping parameters and band structures

This sub-section portrays the variation of nearest neighbor hopping parameters
against the change in the values of strain for individual Xenes and the strain
limit is dependent on their elastic limits and band structures.

Figure 4: Variation of the nearest neighbor tight-binding parameters as a
function of strain along (a)armchair and (b) zigzag directions for
phosphorene, respectively. (C) The band structure of phosphorene along the
path ${Y\Gamma XS}$ for AC strain adhering to the values of $0\%$, $2\%$, and
$5\%$.

We compute the Slater-Koster (SK) parameter, $V_{pp\pi}$, by fitting the DFT
band structure to a tight binding (TB) model using TB-Studio [79]. The
relation between the SK parameter of strained and pristine Xenes is provided
in [74]. The magnitude of hopping amplitude for pristine silicene, germanene,
and stanene are $1.03$ $0.99$ and $0.76$ eV, respectively. These are obtained
by benchmarking with our DFT results. In the case of black phosphorene, we use
the 5-parameter tight-binding model of Rudenko and Katsnelson [39] to describe
the energy dispersion near the Fermi level. The hopping parameters $t_{1}$ to
$t_{5}$ of pristine black phosphorene are obtained using Harrison’s method
[80]. The $\beta$ value is $2$ in the exponential relation linking the SK
parameter of strained and pristine black phosphorene [81]. The variation of
hopping amplitudes as a function of uniaxial strain is presented in Fig. 2 for
silicene, germanene, and stanene ($t_{1},t_{2},t_{3}$) and Fig. 4 for
phosphorene ($t_{1},t_{2},t_{3},t_{4},t_{5}$).
From Fig. 2, it can be illustrated that all the monolayer Xenes follow the
same pattern of hopping parameter variation along both AC and ZZ directions
despite having different magnitudes. We obtain the tight-binding (TB)
parameters for ten strain values irrespective of the elastic limit, and other
values are taken from the line obtained by fitting. For phosphorene, we have
shown the variation of hopping parameters with the applied uniaxial tensile
strain, which is the first of its kind and one of the major highlights of the
paper. From Fig. 4a and 4b, it can be illustrated that the hopping parameters
show minor variation compared to buckled Xene variation. Parameters $t_{3}$
and $t_{4}$ show slight variation with strain, whereas the rest of the
parameters are almost constant with the strain variation irrespective of the
direction.
The band structure of buckled and puckered Xenes computed from DFT
calculations are shown in Fig. 3 and Fig. 4c, respectively. The path shown in
the band structures is in good agreement with the high symmetry path shown in
the corresponding Brillouin zone. The region encircled in orange in each band
structure shows the shifting of the K-point away from $\Gamma$-point with an
increase in the uniaxial tensile strain value applied along the AC direction.
This shifting of K-point is also satisfied quantitatively (see Supplementary
material). The shifting of the K-point will be in the reverse direction when
the strain is in ZZ direction [74]. The band structure of black phosphorene is
presented in Fig. 4c. It is a direct band gap semiconductor with valence band
maximum $(VBM)$ and conduction band minimum $(CBM)$ at $\Gamma$. The valence
and conduction bands of black phosphorene in the $\Gamma-Y$ direction are
flatter than the $\Gamma-X$ direction, which could lead to anisotropic
effective mass in the AC and ZZ directions. With increasing tensile strain in
the AC direction, the band gap increases for strain up to $7\%$. Similarly,
when subjected to increasing tensile strain in the ZZ direction, the band gap
increases for strain up to $4\%$, after which it starts decreasing. This
increase in bandgap is consistent with a previous report by Peng et al. [62].
Beyond $8\%$ tensile strain, the CBM shifts slightly away from the
$\Gamma$-point.

##  III Results and discussion

In this section, we obtain the piezoresistance gauge factor of silicene,
germanene, and stanene along different transport angles using the theoretical
models discussed alongside and rationalize our findings in terms of the change
in transmission associated with the shifting and deformation of the Dirac
cones. This section demonstrates the variation of quantized conductance as a
function of strain. Further, we discuss the connotation of our results in
flexible electronic devices.

###  III.1 Piezoresistivity of Xenes

The primary objective of this investigation is to calculate and compare the
piezoresistivity of $2D$ Xenes in the nanoscale regime obtained from ab initio
calculation using DFT and $V-I$ characteristics obtained from the quantum
transport model. The procedure to obtain the hopping parameters is discussed
in the previous section, and the transport model will be discussed here.
Phosphorene has a bandgap of $1.52$ eV [68] at zero strain. So, the band
counting method is not applicable here as it does not form a zero band gap or
linear energy dispersion relation at $CBM$ and $VBM$. Another process to
calculate the piezoresistance of these kinds of materials is the NEGF method,
which has already been calculated by Naurbaksh et’al [82]. However, we will
compare the data obtained using this method. Henceforth, our calculation has
been restricted to germanene and stanene. The results of silicene [74],
germanene, and stanene will be comprehensively compared in the subsequent
subsection.

Figure 5: The plots of transmission as a function of transport angle
$(\theta)$ for strain values $0\%$, $2\%$ and $5\%$ along the (a),(c)
armchair, and (b), (d) zigzag directions for silicene and germanene
respectively. Transmission plots for stanene akin to the strain values of
$0\%$, $5\%$, and $10\%$ along (e) armchair and (f) zigzag directions. The
green circle in all the plots represent the critical angle. Figure 6: The
plots of current as a function of transport angle $(\theta)$ for strain values
$0\%$, $2\%$ and $5\%$ along the (a),(c) armchair, and (b), (d) zigzag
directions for silicene and germanene respectively. Current plots for stanene
akin to the strain values of $0\%$, $5\%$, and $10\%$ along (e) armchair and
(f) zigzag directions. The green circle in all the plots represent the
critical angle.

After procurement of TB parameters, we use the Landauer formula [83] to
calculate the current-voltage characteristics of silicene in the quasi-
ballistic regime. This method is common to all the Xenes discussed here. The
Landauer formula is expressed as:

$I_{\theta}^{s}(V)=\frac{2q}{h}\int_{-\infty}^{\infty}{T^{s}}(E)[f_{1}(E-\mu_{1})-f_{2}(E-\mu_{2})]dE$
(5)

where $f_{1}$ and $f_{2}$ are the Fermi functions, and $\mu_{1}$ and $\mu_{2}$
are the electro-chemical potentials at the source and drain ends respectively,
$T^{s}$ is the transmission with strain in the Xenes, $q$ is the electronic
charge, $h$ is the Planck’s constant, $E$ is the energy well within the linear
regime, and $I_{\theta}^{s}(V)$ is the current along the transport direction
$(\theta)$ at an applied strain. Here, the voltage varies from $-10$ mV to
$+10$ mV. In general, the strained transmission $T^{s}(E)$ is the product of
transmission probability $T(E)$ and mode density $M_{\theta}^{s}(E)$. The
transmission probability depends on both the length $(L_{s})$ of the Xene
sheet and their mean free path $(\lambda)$ in a quasi-ballistic regime. The
mean free path in the linear regime depends on the relaxation time near the
Dirac cone [44], and as the calculation is in a quasi-ballistic regime, this
will be nearly equal to the length of the Xene sheet [36]. The fundamental
parameters in the ballistic regime are independent of this mean-free path.
Both regimes’ parameters are compared and tabulated in table 5. The mean free
path for silicene and germanene are calculated from [36, 84], and the same
method is applied for stanene to get the value. These values are satisfied in
the context of a quasi-ballistic regime. $L_{s}$ represents the length of the
strained Xene sheets, and with strain, its relation is depicted as
$(L_{s})=L(1+s\rho)$. The transmission probability is given as follows:

$T(E)=\frac{\lambda(E)}{\lambda(E)+L_{s}},$ (6)

Figure 7: The plots of resistance as a function of transport angle $(\theta)$
for strain values $0\%$, $2\%$ and $5\%$ along the (a),(c) armchair, and (b),
(d) zigzag directions for silicene and germanene respectively. Resistance
plots for stanene akin to the strain values of $0\%$, $5\%$, and $10\%$ along
(e) armchair and (f) zigzag directions. The red dotted line in the top row and
green circle in the bottom row for each figure represents strain insensitive
transport angle.

The transmission for silicene and germanene as a function of $(\theta)$ for
$0\%$, $2\%$ and $5\%$ AC strain is shown in Fig. 5a and 5c and for ZZ strain
in Fig. 5b and 5d respectively. The transmission of stanene is shown for the
strain values of $0\%$, $5\%$ and $10\%$ in Fig. 5e and 5f for AC and ZZ
strain, respectively. Here, $\theta$ varies from $0$∘ to $90$∘. From the
figures, we infer that all materials will exhibit opposite patterns concerning
AC and ZZ strain. Mode density at the Dirac cone is required to calculate the
effective number of transverse modes passing through the Dirac point, which is
accomplished through the band-counting method [85, 86]. The Dirac cone
degeneracy for silicene in the first Brillouin zone is two [87, 85]. So,
effectively the mode density is twice the number of transverse modes passing
through the Dirac cone. The pictorial representation of this calculation is
given in the supplementary material. Mathematically, mode density can be
expressed as $M_{\theta}^{s}(E)=2\times n_{\theta}^{s}(E)$. Here,
$n_{\theta}^{s}(E)$ is the number of transverse modes passing through the
constant energy surfaces of the Dirac cone at energy ‘$E$’. The superscript
’$s$’ indicates transverse modes are counted with applied strain on the Dirac
cone. These transverse modes have a separation of $2\pi/w_{s}$, where $w_{s}$
represents the width of strained Xene sheet and is equal to
$w_{s}=w(1+\frac{s}{100})$. Here $w$ represents the width of the unstrained
Xene sheet. Thus, transmission with strain is computed as:
$T^{s}(E)=T(E)*M_{\theta}^{s}(E)$ as explained before.
The transmission remains constant at zero strain for varying transport angles
due to circular constant energy surfaces in the Dirac cone of materials. Due
to AC strain, the value of $L_{x}$ (length of Dirac cone measured along the
x-axis) increases while $L_{y}$ (length of Dirac cone measured along the
y-axis) decreases. As a result, a circle like the Dirac cone is deformed into
an ellipsoid structure, and the number of transverse modes passing through
this will increase because of the large surface area, which increases the mode
density and eventually increases transmission. So, the transmission has a
higher value at $\theta=$ 0∘ than unstrained Xenes sheets, and it decreases
with an increase in the transport angle due to the change in the angle of the
transverse modes propagating through the Dirac cone. On the contrary, the
value of $L_{x}$ decreases while $L_{y}$ increases for ZZ strain which deforms
the Dirac cones into an oval shape, and mode density will reduce here. For the
above same reason, the transmission along $\theta=$ $0^{\circ}$ has a lower
value than the transmission of unstrained Xenes. However, due to a change in
the angle of the propagating modes, the transmission increases with an
increase in $\theta$.
The plots of current for silicene and germanene as a function of transport
angle for $0\%$, $2\%$ and $5\%$ AC strain is shown in Fig. 6a and 6c and for
ZZ strain in Fig. 6b and 6d respectively. The current plot for stanene is
shown for the strain values of $0\%$, $5\%$ and $10\%$ in Fig. 6e and 6f for
AC and ZZ strain, respectively. The plots will depict similar behavior to the
transmission plots as current is directly proportional to the transmission at
a particular energy [see Eq. (5)]. The green bubble in both the transmission
and current plot indicates an intersection point, where the strained and
unstrained lines of transmission and current coincide at a single value in
their respective plots, and the corresponding intersection angle is called the
critical angle. The critical angle is the angle at which the material is
insensitive to the applied strain. The transmission and value of current
increases as we go down group $IV$ in the periodic table, which can be
ascribed to the increase in radius of the Dirac cone (see Table. 5). The above
line can be a general statement applicable to group $IV$ elements.
The contacts are assumed to be ideal in our calculations. The energy range for
this Fermi-Dirac distribution is from $-0.2$ eV to $+0.2$ eV, which is well
within the linear regime. Once current in eq. (5) is obtained, the
piezoresistance is easily calculated from the dynamic change in voltage and
current values. Therefore, it is depicted as follows:

$R_{\theta}^{s}=\frac{dV}{dI_{\theta}^{s}(V)}$ (7)

The above equation is used to calculate the piezoresistance at different
transport angles. The plots of resistance versus $\theta$ as a function of AC
strain are shown in Figs. 7a, 7c and 7e for silicene, germanene, and stanene
respectively. Similarly, for the ZZ strain, the plots are shown in Fig. 7b, 7d
and 7f, respectively. It is obvious that resistance will show an opposite
pattern concerning current plots. However, the critical angle remains constant
for resistance as well, and the same can be depicted in the corresponding
diagram of AC strain through a significant red line drawn vertically. The
value of resistivity obtained in this work for $1~{}\mu m$ wide Xene sheets in
the quasi-ballistic regime are $2.42\text{\,}\mathrm{k\SIUnitSymbolOhm}$,
$1.092\text{\,}\mathrm{k\SIUnitSymbolOhm}$,
$0.815\text{\,}\mathrm{k\SIUnitSymbolOhm}$ for silicene, germanene and stanene
respectively (refer Table 5). It is now quite obvious that as the radius of
the Dirac cone increases, the corresponding resistance value will go down.
Using these values, we can calculate Angular Gauge Factor (AGF). AGF is
expressed as:

$(AGF)_{\theta}^{s}=\frac{(R_{\theta}^{s}-R_{\theta}^{0})}{s*R_{\theta}^{0}},$
(8)

where, ${R_{\theta}}^{0}$ is the resistance at zero strain and
${R_{\theta}}^{s}$ is the resistance at a particular value of strain $s$. As
the strain is applied throughout the transport angle $(\theta)$, the AGF is
the average of all GFs at different strains along a particular transport
angle. Thus, the average AGF is expressed as:

$(AGF)_{\theta}=\overline{(AGF)_{\theta}^{s}}$ (9)

The GFs of all the monolayer Xenes in the AC and ZZ strain are shown in Fig.
9. The result will be discussed in the subsequent subsection.

###  III.2 Conductance Modulation

Figure 8: Plots of conductance as a function of energy (eV) for the applied
strain having the values $0\%$, $2\%$ and $5\%$ for silicene, germanene and
values $0\%$, $5\%$ and $10\%$ for stanene along the (a),(c), (e) armchair
direction. Plots of conductance as a function of transport angle ($\theta$)
and applied strain for the values $0\%$, $2\%$ and $5\%$ for silicene,
germanene, and $0\%$, $5\%$ and $10\%$ for stanene along the (b), (d), (f)
zigzag direction. The green circle in bottom row plots represent the strain
insensitive transport angle.

Conductance modulation is the process of varying conductance concerning the
strained transmission as a function of energy. This section aims to show the
variation of conductance with respect to angular strain for buckled $2D$
Xenes. The transport properties of graphene taken from Landauer-Buttiker (LB)
formalism [88] can serve as a template for its equivalent $2D$ Xenes, which
are quasi-ballistic in nature. The corresponding ballistic conductance is
given as $G(E)=G_{0}*T^{s}(E)$, where $G_{0}$ is quantum conductance, a
fundamental parameter expressed as: $G_{0}=\frac{2e^{2}}{h}$. $T^{s}(E)$
represents strained transmission, defined in the previous section.
As mentioned earlier, conductance is related to energy through transmission,
so the plots for conductance concerning the energy window are shown in Fig.
8a, 8c and 8e with strain applied in AC direction. Similarly, plots for
conductance as a function of $\theta$ are shown in Fig. 8b, 8d and 8f with
strain applied in ZZ direction. These plots are applicable for monolayer
buckled Xenes; however, puckered structures are also expected to have
conductance modulation. Figures 8a, 8c and 8e depict that conductance
increases with an increase in energy and strain due to an increase in the
number of transverse modes passing through the Dirac cone in AC strain. The
increase in transverse modes is related to the deformation of the Dirac cone
from circle to ellipsoid [74]. In the case of ZZ strain, the deformation will
result in an oval shape which reduces the number of modes passing through the
Dirac cone [74]. Hence, the plots will exhibit opposite characteristics in the
ZZ strain. The increase in conductance for energy and strain is also observed
in graphene [89], which will validate our simulated plots. Variation of
conductance with $\theta$ is shown in Fig. 8b, 8d and 8f for silicene,
germanene and, stanene respectively. These plots resemble the plots of
transmission and current. Innately, the conductance below and above the
critical angle exhibits opposite behavior. Initially, conductance decreases
with increased strain, and as we go beyond the critical angle (shown by the
green circle), it starts increasing with strain. Zero strain conductance shows
an almost flat line, but the small changes in the values are discretized. As
the changes in plots are discrete, it is called quantized conductance.

###  III.3 Connotation of the results

Table 3 compares the piezoresistance GF of silicene, germanene, and stanene
with other prominent materials like phosphorene, graphene, and popular
semiconductors like silicon and germanium. From this table, we infer that the
GF of silicene and germanene are petite, whereas their 3D counterparts,
silicon, and germanium, have very high GF. Also, as we go along the periodic
table, the GF of the monolayer buckled Xenes will increase. This increase in
the value forms a pattern, which can be generalized to all the group $IV$
buckled Xenes. Phosphorene, a group $V$ element with puckered structure, has a
very high GF, which can act as a frontispiece for other group $V$ elements.

Table 3: Comparisons of the GFs of various materials. Material |  GF |  Reference
---|---|---
Silicene |  0.758 |  [74]
Germanene |  1.391 |  This work
Stanene |  2.349 |  This work
Phosphorene |  120 |  [82]
Graphene |  0.6 |  [86]
Silicon |  200 |  [90]
Germanium |  150 |  [90]
Figure 9: The plots of angular gauge factor along with its sinusoidal fit as a
function of transport angle along the (a),(c),(e) armchair and (b),(d),(f)
zigzag directions for silicene, germanene, and stanene respectively. The
highlighted point in each plot represents a critical angle.

Figures 9a, 9c and 9e and 9b, 9d and 9f depict AGFs as a function of the
transport angle for silicene, germanene and stanene for AC and ZZ directions
respectively. The AGFs have sinusoidal dependence on the angle ($\theta$) for
both the AC and ZZ strains, and it is evident by now that the variation of AGF
resembles a sinusoidal function similar to the one obtained by Sinha et al.
for graphene because of the similarity in deformation of Dirac points in the
reciprocal space. This variation in gauge factor goes from positive to
negative in ZZ strain (see Fig. 9b,d,f) as the strained resistance is more
than unstrained resistance(see Fig. 8b,d,f). The highlighted circle in each
plot represents the critical angle denoting zero AGF, which stands for the
angle insensitive to strain, and all the strained and unstrained parameters
will have only one single value approximately. After the critical angle, the
gauge factor becomes negative for ZZ strain.

Table 4: Values of constants P and Q for different Xenes to fit sinusoidal variation. | AC Strain | ZZ Strain
---|---|---
Material | P | Q | P | Q
Silicene | 0.249 | 0.237 | 1.079 | -0.239
Germanene | 0.710 | 0.268 | 0.820 | -0.275
Stanene | 1.376 | 1.009 | 1.7555 | -0.5075

In AC strain, the critical angle is the same ($81$∘) for both silicene and
germanene, pointing out that, despite the difference in quantifying cell
dimension and related aspects, still the two elements are very close to each
other. For stanene, the critical angle is $69$∘, attributed to a change in
applied strain and notable change in the parameters listed in table 5.
Interestingly, for ZZ strain, all the $3$ materials have the same critical
angle of $34$∘ (ref Fig. 9(b, d and f)) despite all the differences, which
draws an exciting conclusion for the rest of the group $IV$ elements. To fit
the sinusoidal plot, mathematically, AGF is in the form of

$\mathrm{AGF={-\\{P\cos(2\theta)+Q\\}}},$ (10)

Table 5: Comparison of various fundamental parameters of silicene, germanene, stanene, and phosphorene. | Silicene | Germanene | Stanene | Phosphorene
---|---|---|---|---
Parameters | $T(E)=1$ | T(E) $<$ 1 | $T(E)=1$ | T(E) $<$ 1 | $T(E)=1$ | T(E) $<$ 1 | $T(E)=1$
Dirac cone radius (1 meV) (Å-1) | $0.0002928$ | $0.0002928$ | $0.0002872$ | $0.0002872$ | $0.0003252$ | $0.0003252$ | $0.0002159$
Current ($\mu$A) | $24.76$ | $4.12$ | $24.30$ | $9.11$ | $27.598$ | $12.26$ |  $5.5$ [91]
Resistance (K$\Omega$) | $0.404$ | $2.428$ | $0.409$ | $1.092$ | $0.360$ | $0.815$ | -
Conductance($\mho$) | $0.0288$ | $0.00481$ | $0.03148$ | $0.0111$ | $0.03226$ | $0.01673$ |  $0.00061$ [91]

where P and Q are constants, and their values for buckled Xenes are listed
below in the table 4. Several combinations can be possible for the constants,
but the listed values suit the best for sinusoidal fit regarding the variation
of AGF concerning the transport angle.

Table 5 lists various fundamental parameters of buckled and puckered Xenes
like silicene, germanene, stanene, and phosphorene with and without
transmission probability at zero strain. Transport in a ballistic regime is
denoted as $T(E)=1$, stating it is independent of the mean free path, and
transmission for a quasi-ballistic regime is denoted as $T(E)<1$. The former
is called without transmission probability, and the latter is called with
transmission probability. Dirac cone radius is invariant to transmission
probability as it is deprived of any shifting and deformation at zero strain.
So, the radius will hold its value irrespective of the above-mentioned
regimes. However, the radius of the Dirac cones increases as the atomic number
of the constituent group $IV$ element increases, starting from graphene to
stanene at zero strain. Considering transmission to be equal to one, the
current value will increase. The corresponding effects can be seen for
resistance and conductance. We obtain the value of resistivity as
$0.404\text{\,}\mathrm{k\SIUnitSymbolOhm}$,
$0.409\text{\,}\mathrm{k\SIUnitSymbolOhm}$,
$0.360\text{\,}\mathrm{k\SIUnitSymbolOhm}$ for buckled Xenes in ballistic
regime. The resistivity value for stanene, germanene, and silicene is less
compared to the resistivity of graphene [85]. This small resistivity value is
because the number of TMs in other buckled Xenes is more, than in graphene
(for the sheet width) due to the larger size of its Dirac cones. For
phosphorene, we have considered only the ballistic regime. Its Dirac cone
radius is the smallest of all the listed values, resulting in the minimum
current value compared to other corresponding values taken only in the
ballistic regime. We have observed anisotropy in the shifting and deformation
of the Dirac cone for phosphorene in the presence of strain either in AC or ZZ
direction. The quantitative results are also listed (see Supplementary
information) for clarity. The proportion between the parameters, irrespective
of transmission probability, is maintained throughout, which persuades the
correctness of the result.
The discussed buckled Xenes are atomically-thin membranes as they are
monolayers and are expected to have a high value of adhesivity (like most of
the similar 2D materials) [92] and elasticity [55]. These significant values
lead to a future endeavor, where these Xenes are expected to have a very high-
pressure sensitivity like planar monolayer graphene [93] and
$\mathrm{PtSe_{2}}$ [94, 95] despite a low GF. These Xene sheets also show
high resistance to the change in strain. This resistant is realized when the
normalized resistance is plotted against strain for silicene and graphene
[74], which unveils their vital electronic characteristic to the change in
strain value. Several characteristics like levelheaded conductivity [96], high
mobility calculation [97], and high elastic limit [55] of buckled Xenes will
enable them to act as electrodes and interconnects in flexible electronic
devices. Also, they have the persuasion to act as potential topological
materials [98]. Based on results obtained in this voluminous model, we
propound that the Xenes can be strong contenders in flexible electronics and
show promising characteristics in any application related to straintronics.

##  IV Conclusion

In this paper, we investigated the piezoresistance of buckled Xenes like
germanene and stanene using ab-initio density function theory and quantum
transport theory. We compared the result with silicene and graphene. We
calculated the directional piezoresistance for different strain values along
the armchair and zigzag directions. We typically obtained a smaller value of
the directional piezoresistances and their sinusoidal dependence. The value of
the gauge factor increases as we go up for group $IV$ elements along their
atomic numbers in the periodic table. The strain-insensitive transport angles
corresponding to the zero gauge factors are $81^{\circ}$ and $34^{\circ}$ for
armchair and zigzag strains, respectively, for silicene and germanene. For
stanene as the strain limit is extended to $10\%$ and notable changes in the
fundamental parameters, the critical angle for stanene along armchair and
zigzag directions are $69^{\circ}$ and $34^{\circ}$ respectively. The small
gauge factor of buckled Xenes can be attributed to its robust Dirac cone and
strain-independent valley degeneracy. We have realized conductance
quantization in a quasi-ballistic regime for the buckled Xenes. We have also
shown the strained tight binding parameters of phosphorene, which is the first
of its kind in the case of puckered Xenes. Further, we are exploring the
spintronics part of $2D$ Xenes, where we will see the response of the
interfacial effects between the contacts and the Xenes against the applied
strain. Based on the obtained results, we propose the buckled Xenes as an
interconnect in flexible electronics and are promising candidates for various
applications in straintronics.

###### Acknowledgements.

We thank Dr.Abhinaba Sinha for his guidance and support towards preparing this
manuscript. The Research and Development work undertaken in the project under
the Visvesvaraya Ph.D. Scheme of Ministry of Electronics and Information
Technology, Government of India, is implemented by Digital India Corporation
(formerly Media Lab Asia). This work was also supported by the Science and
Engineering Research Board (SERB), Government of India, Grant No.
CRG/2021/003102.

## References

  * Vogt  _et al._ [2012] P. Vogt, P. De Padova, C. Quaresima, J. Avila, E. Frantzeskakis, M. C. Asensio, A. Resta, B. Ealet, and G. Le Lay, Physical review letters 108, 155501 (2012).
  * Bianco  _et al._ [2013] E. Bianco, S. Butler, S. Jiang, O. D. Restrepo, W. Windl, and J. E. Goldberger, ACS nano 7, 4414 (2013).
  * Osaka  _et al._ [1994] T. Osaka, H. Omi, K. Yamamoto, and A. Ohtake, Physical Review B 50, 7567 (1994).
  * Zimmermann  _et al._ [1997] H. Zimmermann, R. C. Keller, P. Meisen, and M. Seelmann-Eggebert, Surface science 377, 904 (1997).
  * Li  _et al._ [2014] L. Li, Y. Yu, G. J. Ye, Q. Ge, X. Ou, H. Wu, D. Feng, X. H. Chen, and Y. Zhang, Nature nanotechnology 9, 372 (2014).
  * Liu  _et al._ [2014] H. Liu, A. T. Neal, Z. Zhu, Z. Luo, X. Xu, D. Tománek, and P. D. Ye, ACS nano 8, 4033 (2014).
  * Kamal and Ezawa [2015] C. Kamal and M. Ezawa, Physical Review B 91, 085423 (2015).
  * Novoselov  _et al._ [2004] K. S. Novoselov, A. K. Geim, S. V. Morozov, D.-e. Jiang, Y. Zhang, S. V. Dubonos, I. V. Grigorieva, and A. A. Firsov, science 306, 666 (2004).
  * Neto  _et al._ [2009] A. C. Neto, F. Guinea, N. M. Peres, K. S. Novoselov, and A. K. Geim, Reviews of modern physics 81, 109 (2009).
  * Lee  _et al._ [2008] C. Lee, X. Wei, J. W. Kysar, and J. Hone, science 321, 385 (2008).
  * Liu  _et al._ [2007] F. Liu, P. Ming, and J. Li, Physical Review B 76, 064120 (2007).
  * Kim  _et al._ [2009] K. S. Kim, Y. Zhao, H. Jang, S. Y. Lee, J. M. Kim, K. S. Kim, J.-H. Ahn, P. Kim, J.-Y. Choi, and B. H. Hong, nature 457, 706 (2009).
  * Geim [2009] A. K. Geim, science 324, 1530 (2009).
  * Rosas  _et al._ [2011] J. H. Rosas, R. R. Gutiérrez, A. Escobedo-Morales, and E. C. Anota, Journal of molecular modeling 17, 1133 (2011).
  * Yuan  _et al._ [2020] P. Yuan, T. Zhang, J. Sun, L. Liu, Y. Yao, and Y. Wang, Journal of Semiconductors 41, 081003 (2020).
  * Molle  _et al._ [2017] A. Molle, J. Goldberger, M. Houssa, Y. Xu, S.-C. Zhang, and D. Akinwande, Nature materials 16, 163 (2017).
  * Tao  _et al._ [2015] L. Tao, E. Cinquanta, D. Chiappe, C. Grazianetti, M. Fanciulli, M. Dubey, A. Molle, and D. Akinwande, Nature nanotechnology 10, 227 (2015).
  * Houssa  _et al._ [2010] M. Houssa, G. Pourtois, M. Heyns, V. Afanas’Ev, and A. Stesmans, Journal of The Electrochemical Society 158, H107 (2010).
  * Roman and Cranford [2014] R. E. Roman and S. W. Cranford, Computational materials science 82, 50 (2014).
  * Liu  _et al._ [2011] C.-C. Liu, W. Feng, and Y. Yao, Physical review letters 107, 076802 (2011).
  * Acun  _et al._ [2015] A. Acun, L. Zhang, P. Bampoulis, M. v. Farmanbar, A. van Houselt, A. Rudenko, M. Lingenfelder, G. Brocks, B. Poelsema, M. Katsnelson, _et al._ , Journal of physics: Condensed matter 27, 443002 (2015).
  * Rachel and Ezawa [2014] S. Rachel and M. Ezawa, Physical Review B 89, 195303 (2014).
  * Ni  _et al._ [2012] Z. Ni, Q. Liu, K. Tang, J. Zheng, J. Zhou, R. Qin, Z. Gao, D. Yu, and J. Lu, Nano letters 12, 113 (2012).
  * Guzmán-Verri and Voon [2007] G. G. Guzmán-Verri and L. L. Y. Voon, Physical Review B 76, 075131 (2007).
  * Cahangirov  _et al._ [2009] S. Cahangirov, M. Topsakal, E. Aktürk, H. Şahin, and S. Ciraci, Physical review letters 102, 236804 (2009).
  * Takeda and Shiraishi [1994] K. Takeda and K. Shiraishi, Physical Review B 50, 14916 (1994).
  * Hultgren  _et al._ [1935] R. Hultgren, N. Gingrich, and B. Warren, The Journal of Chemical Physics 3, 351 (1935).
  * Thurn and Kerbs [1966] H. Thurn and H. Kerbs, Angewandte Chemie International Edition in English 5, 1047 (1966).
  * Rodin  _et al._ [2014] A. Rodin, A. Carvalho, and A. C. Neto, Physical review letters 112, 176801 (2014).
  * Brown and Rundqvist [1965] A. Brown and S. Rundqvist, Acta Crystallographica 19, 684 (1965).
  * Island  _et al._ [2015] J. O. Island, G. A. Steele, H. S. van der Zant, and A. Castellanos-Gomez, 2D Materials 2, 011002 (2015).
  * Tran  _et al._ [2014] V. Tran, R. Soklaski, Y. Liang, and L. Yang, Physical Review B 89, 235319 (2014).
  * Low  _et al._ [2014a] T. Low, A. Rodin, A. Carvalho, Y. Jiang, H. Wang, F. Xia, and A. C. Neto, Physical Review B 90, 075434 (2014a).
  * Xia  _et al._ [2014] F. Xia, H. Wang, and Y. Jia, Nature communications 5, 1 (2014).
  * Low  _et al._ [2014b] T. Low, R. Roldán, H. Wang, F. Xia, P. Avouris, L. M. Moreno, and F. Guinea, Physical review letters 113, 106802 (2014b).
  * Abidin  _et al._ [2017] B. Abidin, K. Yeoh, D. S. Ong, and T. Yong, Journal of Physics D: Applied Physics 50, 425103 (2017).
  * Harris  _et al._ [2016] K. Harris, A. Elias, and H.-J. Chung, Journal of materials science 51, 2771 (2016).
  * Wong  _et al._ [2009] W. S. Wong, M. L. Chabinyc, T.-N. Ng, and A. Salleo, Flexible Electronics: Materials and Applications , 143 (2009).
  * Rudenko and Katsnelson [2014] A. N. Rudenko and M. I. Katsnelson, Physical Review B 89, 201408 (2014).
  * Yang  _et al._ [2016] M. Yang, H.-J. Duan, and R.-Q. Wang, Physica Scripta 91, 105801 (2016).
  * Pereira  _et al._ [2009] V. M. Pereira, A. C. Neto, and N. Peres, Physical Review B 80, 045401 (2009).
  * Reich  _et al._ [2002] S. Reich, J. Maultzsch, C. Thomsen, and P. Ordejon, Physical Review B 66, 035412 (2002).
  * Yao  _et al._ [2007] Y. Yao, F. Ye, X.-L. Qi, S.-C. Zhang, and Z. Fang, Physical Review B 75, 041401 (2007).
  * Shao  _et al._ [2013] Z.-G. Shao, X.-S. Ye, L. Yang, and C.-L. Wang, Journal of Applied Physics 114, 093712 (2013).
  * Qin  _et al._ [2012] R. Qin, C.-H. Wang, W. Zhu, and Y. Zhang, Aip Advances 2, 022159 (2012).
  * Lew Yan Voon  _et al._ [2010] L. Lew Yan Voon, E. Sandberg, R. Aga, and A. Farajian, Applied Physics Letters 97, 163114 (2010).
  * Yuhara and Le Lay [2020] J. Yuhara and G. Le Lay, Japanese Journal of Applied Physics 59, SN0801 (2020).
  * Ye  _et al._ [2014] X.-S. Ye, Z.-G. Shao, H. Zhao, L. Yang, and C.-L. Wang, Rsc Advances 4, 21216 (2014).
  * Balendhran  _et al._ [2015] S. Balendhran, S. Walia, H. Nili, S. Sriram, and M. Bhaskaran, small 11, 640 (2015).
  * Nakamura  _et al._ [2017] Y. Nakamura, T. Zhao, J. Xi, W. Shi, D. Wang, and Z. Shuai, Advanced Electronic Materials 3, 1700143 (2017).
  * Kurpas  _et al._ [2018] M. Kurpas, M. Gmitra, and J. Fabian, Journal of Physics D: Applied Physics 51, 174001 (2018).
  * Qiao  _et al._ [2014] J. Qiao, X. Kong, Z.-X. Hu, F. Yang, and W. Ji, Nature communications 5, 4475 (2014).
  * Choi  _et al._ [2010] S.-M. Choi, S.-H. Jhi, and Y.-W. Son, Physical Review B 81, 081407 (2010).
  * Xu  _et al._ [2013] M. Xu, T. Liang, M. Shi, and H. Chen, Chemical reviews 113, 3766 (2013).
  * Peng  _et al._ [2013] Q. Peng, X. Wen, and S. De, Rsc Advances 3, 13772 (2013).
  * Qin  _et al._ [2014] R. Qin, W. Zhu, Y. Zhang, and X. Deng, Nanoscale research letters 9, 1 (2014).
  * Ding  _et al._ [2018] N. Ding, H. Wang, L. Liu, W. Guo, X. Chen, and C.-M. L. Wu, Journal of Physics and Chemistry of Solids 113, 201 (2018).
  * Wang and Ding [2013] Y. Wang and Y. Ding, Solid State Communications 155, 6 (2013).
  * John and Merlin [2016] R. John and B. Merlin, Crystal Structure Theory and Applications 5, 43 (2016).
  * Mojumder  _et al._ [2015] S. Mojumder, A. A. Amin, and M. M. Islam, Journal of Applied Physics 118, 124305 (2015).
  * Van den Broek  _et al._ [2014] B. Van den Broek, M. Houssa, E. Scalise, G. Pourtois, V. Afanas‘ev, and A. Stesmans, 2D Materials 1, 021004 (2014).
  * Peng  _et al._ [2014] X. Peng, Q. Wei, and A. Copple, Physical Review B 90, 085402 (2014).
  * Sa  _et al._ [2014] B. Sa, Y.-L. Li, J. Qi, R. Ahuja, and Z. Sun, The Journal of Physical Chemistry C 118, 26560 (2014).
  * Fei  _et al._ [2015] R. Fei, V. Tran, and L. Yang, Physical Review B 91, 195319 (2015).
  * Fei and Yang [2014] R. Fei and L. Yang, Nano letters 14, 2884 (2014).
  * Wang  _et al._ [2015] L. Wang, A. Kutana, X. Zou, and B. I. Yakobson, Nanoscale 7, 9746 (2015).
  * Wang  _et al._ [2020] B. Wang, S. Zhong, Y. Ge, H. Wang, X. Luo, and H. Zhang, Nano Research 13, 891 (2020).
  * Yarmohammadi  _et al._ [2020] M. Yarmohammadi, M. Mortezaei, and K. Mirabbaszadeh, Physica E: Low-dimensional Systems and Nanostructures 124, 114323 (2020).
  * Kresse and Furthmüller [1996a] G. Kresse and J. Furthmüller, Computational materials science 6, 15 (1996a).
  * Kresse and Furthmüller [1996b] G. Kresse and J. Furthmüller, Physical review B 54, 11169 (1996b).
  * Kresse and Joubert [1999] G. Kresse and D. Joubert, Physical review b 59, 1758 (1999).
  * Perdew  _et al._ [1996] J. P. Perdew, K. Burke, and M. Ernzerhof, Physical review letters 77, 3865 (1996).
  * Grimme  _et al._ [2010] S. Grimme, J. Antony, S. Ehrlich, and H. Krieg, The Journal of chemical physics 132, 154104 (2010).
  * Sahoo  _et al._ [2022] S. Sahoo, A. Sinha, N. A. Koshi, S.-C. Lee, S. Bhattacharjee, and B. Muralidharan, Journal of Physics D: Applied Physics 55, 425301 (2022).
  * Zhang and Liu [2011] Y. Zhang and F. Liu, Applied Physics Letters 99, 241908 (2011).
  * Kaloni and Schwingenschlögl [2013] T. P. Kaloni and U. Schwingenschlögl, Chemical Physics Letters 583, 137 (2013).
  * Ni  _et al._ [2015] Z. Ni, E. Minamitani, Y. Ando, and S. Watanabe, Physical Chemistry Chemical Physics 17, 19039 (2015).
  * Shi and Singh [2017] Z. Shi and C. V. Singh, Nanoscale 9, 7055 (2017).
  * Nakhaee  _et al._ [2020] M. Nakhaee, S. A. Ketabi, and F. M. Peeters, Computer Physics Communications 254, 107379 (2020).
  * Harrison [2004] W. A. Harrison,  _Elementary electronic structure (revised edition)_ (World Scientific Publishing Company, 2004).
  * Midtvedt  _et al._ [2017] D. Midtvedt, C. H. Lewenkopf, and A. Croy, Journal of Physics: Condensed Matter 29, 185702 (2017).
  * Nourbakhsh and Asgari [2018] Z. Nourbakhsh and R. Asgari, Physical Review B 98, 125427 (2018).
  * Datta [2018] S. Datta,  _Lessons from Nanoelectronics: A New Perspective on Transport—Part B: Quantum Transport_ (World Scientific, 2018).
  * Gaddemane  _et al._ [2016] G. Gaddemane, W. G. Vandenberghe, and M. V. Fischetti, in  _2016 International Conference on Simulation of Semiconductor Processes and Devices (SISPAD)_ (IEEE, 2016) pp. 353–356.
  * Sinha  _et al._ [2019] A. Sinha, A. Sharma, A. Tulapurkar, V. R. Rao, and B. Muralidharan, Physical Review Materials 3, 124005 (2019).
  * Sinha  _et al._ [2020] A. Sinha, A. Sharma, P. Priyadarshi, A. Tulapurkar, and B. Muralidharan, Physical Review Research 2, 043041 (2020).
  * Liu  _et al._ [2013] Z. Liu, J. Wang, and J. Li, Physical Chemistry Chemical Physics 15, 18855 (2013).
  * Datta [1997] S. Datta,  _Electronic transport in mesoscopic systems_ (Cambridge university press, 1997).
  * Yan [2021] W. Yan, Results in Physics 20, 103723 (2021).
  * Smith [1954] C. S. Smith, Physical review 94, 42 (1954).
  * Rout  _et al._ [2019] C. S. Rout, D. J. Late, and H. Morgan,  _Fundamentals and sensing applications of 2D materials_ (Woodhead Publishing, 2019).
  * Megra and Suk [2019] Y. T. Megra and J. W. Suk, Journal of Physics D: Applied Physics 52, 364002 (2019).
  * Sinha  _et al._ [2022] A. Sinha, P. Priyadarshi, and B. Muralidharan, Journal of Applied Physics 132, 154501 (2022).
  * Smith  _et al._ [2013] A. Smith, F. Niklaus, A. Paussa, S. Vaziri, A. C. Fischer, M. Sterner, F. Forsberg, A. Delin, D. Esseni, P. Palestri, _et al._ , Nano letters 13, 3237 (2013).
  * Wagner  _et al._ [2018] S. Wagner, C. Yim, N. McEvoy, S. Kataria, V. Yokaribas, A. Kuc, S. Pindl, C.-P. Fritzen, T. Heine, G. S. Duesberg, _et al._ , Nano letters 18, 3738 (2018).
  * Zhang  _et al._ [2014] X. Zhang, H. Xie, M. Hu, H. Bao, S. Yue, G. Qin, and G. Su, Physical Review B 89, 054310 (2014).
  * Mir  _et al._ [2020] S. H. Mir, V. K. Yadav, and J. K. Singh, ACS omega 5, 14203 (2020).
  * Zhao and Wang [2020] A. Zhao and B. Wang, APL Materials 8, 030701 (2020).


# Long Horizon Risk-Averse Motion Planning: A Model-Predictive Approach

Chris van der Ploeg1,2, Robin Smit1, Arjan Teerhuis1, Emilia Silvas1,3

(Date: July 2022)

###### Abstract.

Developing safe automated vehicles that can be proactive, safe, and
comfortable in mixed traffic requires improved planning methods that are risk-
averse and that account for predictions of the motion of other road users. To
consider these criteria, in this article, we propose a non-linear model-
predictive trajectory generator scheme, which couples the longitudinal and
lateral motion of the vehicle to steer the vehicle with minimal risk, while
progressing towards the goal state. The proposed method takes into account the
infrastructure, surrounding objects, and predictions of the objects’ state
through artificial potential-based risk fields included in the cost function
of the model-predictive control (MPC) problem. This trajectory generator
enables anticipatory maneuvers, i.e., mitigating risk far before any safety-
critical intervention would be necessary. The method is proven in several case
studies representing both highways- and urban situations. The results show the
safe and efficient implementation of the MPC trajectory generator while
proving its real-time applicability.

1Netherlands Organisation for Applied Scientific Research, Integrated Vehicle
Safety Group, 5700 AT Helmond, The Netherlands.

2Eindhoven University of Technology, Dynamics and Control Group, Mechanical
Engineering Dept., P.O. Box 513, 5600 MB, Eindhoven, The Netherlands.(e-mail:
c.j.v.d.ploeg@tue.nl).

3Eindhoven University of Technology, Control Systems Technology Group,
Mechanical Engineering Dept., P.O. Box 513, 5600 MB, Eindhoven, The
Netherlands.

##  1\. Introduction

Increasing levels of vehicle automation are projected to reduce significantly
the number of fatalities caused by collisions and maximize driver comfort and
efficiency [1]. To achieve this, it is of great importance for automated
vehicles to negotiate and navigate in all traffic situations, dangerous or
not, with minimal risk. While minimal risk forms an essential argument in
terms of safety, progression towards the desired destination forms a large
additional factor for human acceptance of automated vehicles. For example, an
excessively slow-driving vehicle, which is arguably less susceptible to risks,
may severely reduce acceptance [2]. While planning the motion of a vehicle,
these factors of risk and progression need to be taken into account.

Although various automated driving architectures exist, there is yet no expert
consent on an optimal one. Nevertheless, in general, motion planning for an
automated vehicle can be divided into global (e.g., route planning) and local
planning (e.g., trajectory generation) [3]. The task of global planning, i.e.,
finding a route from some initial location to a destination, occurs on a large
road segment level and is generally done through path search algorithms such
as $\text{A}^{*}$ [4]. It is thereafter the task of a trajectory generator to
process these global points and create driveable, safe, and comfortable
trajectories to traverse the path towards the destination while satisfying the
requirements of minimal risk and progression.

Optimization-based methods, more specifically MPC algorithms, are often used
to solve the trajectory generation problem. These methods explicitly take into
account the vehicle model, as they can optimally balance multiple objectives
and various constraints. Several works deploy an MPC-based trajectory
generator, where the road boundaries and dynamical object avoidance are either
formulated through hard constraints [5] or through the use of artificial
potential fields (APFs) [6]. Note, that the notion of risk fits well with the
control-theoretic notion of artificial potential fields, first introduced in
[7]. These APFs allow the analytical formulation of repulsive fields, i.e.,
regions to stay away from, and attractive fields, i.e., regions to be
attracted towards. Several works in literature advocate the use of these APFs
as a replacement for hard anti-collision constraints, due to the potential
non-convexity of the constraints and, hence, a higher computational burden
[8]. Moreover, artificial potential fields circumvent the issue of the
potential infeasibility of the optimization problem caused by unavoidable
collisions, since they act in the problem cost function and, hence, a
collision would not lead to the infeasibility of the problem. In this case, a
minimal-risk form of damage mitigation would result in a solution, instead of
no solution at all.

Model-predictive APF-based trajectory generators have been proposed in several
different contexts in literature, ranging from emergency collision avoidance
[6] up to longer-term planning, e.g., highway maneuvering [9, 10, 11]. Many of
these solutions, however, only consider relatively short horizons (e.g., 1.6s
[8]) due to, e.g., the high fidelity of the used model. To make informed
decisions and actions, a trajectory generator is expected to look ahead for a
sufficient distance or time. For reference, a human driver would, depending on
the context, look ahead and anticipate on the environment $3-4$ seconds up to
tens of seconds [12]. Another shortcoming in the literature is related to
decision-making modules or pre-planned trajectories, which accompany the
model-predictive trajectory generator, and are deployed to decide and give
directions on, e.g., changing a lane [11, 8]. Involving separate components to
pre-determine actions or directions of the trajectory generator may result in
a loss of optimality (i.e., increased risk or reduced progression).

Our contributions: In summary, there exist various solutions for model-
predictive minimal risk trajectory planning. However, there does not yet exist
an approach that finds minimal risk decisions and actions in an optimal manner
(i.e., without the use of additional finite state machines or pre-defined
trajectories). As such, the contributions of our work are focused on the
following three points:

  1. (i)

Risk-averse trajectory planning without discrete decisions or pre-defined
trajectories: we propose a monolithic trajectory planner which can undertake
anticipatory complex risk-averse maneuvers towards a destination without the
need for discrete decisions or pre-defined trajectories to follow.

  2. (ii)

Multiple operational domains: the planner is simulated both in relevant
highway scenarios, as well as an urban VRU crossing scenario, incorporating
both timely and late object detection (e.g., bad weather effects on sensors).
This shows its applicability to a variety of operational domains and
conditions.

  3. (iii)

Real-time implementation: the algorithm is implemented in the CasADi [13]
optimization framework using an interior point solver (IPOPT [14]) and is
deployed in real-time, showing the feasibility of solving such a problem in
experimental settings.

The outline of our work is as follows. First, in Section 2, a mathematical
description of the problem statement is introduced. In Section 3, the
optimization problem is formulated through the notions of risk and
progression. In Section 4, a simulation study is presented in both highway and
urban scenarios, showing the expected minimal-risk results of our scenarios
and their robustness in the case of late object detections. Finally, Section 5
concludes the work.

##  2\. Preliminaries and Problem Statement

In various on-road situations (e.g., the lane change motion tasks shown in
Fig. 1(a)) the automated vehicle needs to plan a safe, comfortable and
drivable motion. To do so, the automated vehicle is expected to obey the
structured environment, (i.e., the road infrastructure, traffic rules, etc.),
interact with mixed-traffic (manually driven vehicles, cyclists, pedestrians),
and account for their future uncertain movements in various situations (e.g.
bad light, poor road markings, bad weather, dense traffic).

The vehicle, with coordinate system $(x,y)$, is expected to drive from an
initial state given by $(x_{0},y_{0},\theta_{0},v_{0})$, located on top of the
coordinate system (i.e., $x_{0},y_{0},\theta_{0}=0$) , towards the reference
state $r_{N}$ which is positioned relative to the initial condition of the
vehicle and contains a planar position and longitudinal velocity $v$ as a
reference, i.e., $(x_{N},y_{N},\theta_{N},v_{N})$. Here, $N$ denotes the
horizon of the trajectory generator to be planned towards the goal state,
i.e., the look-ahead time. Throughout this work, we assume that the reference
state, $\mathbf{r}_{N}$, is located in the center of the rightmost lane on the
road at a distance of $N\cdot t_{s}\cdot\bar{v}$, where $t_{s}$ represents the
sampling time between each trajectory element and $\bar{v}$ represents the
maximally allowed velocity, provided by the infrastructure. Furthermore, it is
assumed that this reference state is known a-priori, and could be provided by
a global planner (e.g., an $\text{A}^{*}$ algorithm [4]). The road
infrastructure is assumed in this work to be described by an integer set of
polynomials, representing the boundaries of the lanes, i.e.:

(1)
$\displaystyle\tilde{y}_{i}=c_{0}^{(i)}+c_{1}^{(i)}\tilde{x}_{i}+c_{2}^{(i)}\tilde{x}_{i}^{2}+c_{3}^{(i)}\tilde{x}_{i}^{3}$

where $i\in\mathbb{Z}$, the superscript (i) selects the polynomial constants
of the $i$-th line measurement, i.e., $c_{j}^{(i)},\>j\in[0,1,2,3]$, and
($\tilde{x},\tilde{y}$) denote the planar coordinates of the $i$-th line
relative to the vehicle initial condition. Note, that polynomial functions of
arbitrary order can be used following our methodology. Finally, the vehicle is
expected to safely drive along with the infrastructure in the presence of
objects (e.g., vulnerable road users (VRUs) or other vehicles). The state of
the objects $\mathbf{o}_{k}\forall k\in[0\ldots N]$ is denoted as
$(x^{o}_{k},y^{o}_{k},\theta^{o}_{k})$, which is measured and predicted with
respect to the initial condition of the vehicle. To include the predicted
motion of the object in motion planning, the state of the object is assumed to
be exactly measurable/predictable up to time $N$.

$x,v$
---

$-c_{0}^{2}$
---

$c_{0}^{1}$
---

$-c_{0}^{3}$
---

$r_{N}(x_{0},y_{0},\theta_{0})$
---
$o_{0}(x_{0},y_{0},\theta_{0})$
---
$o_{N}(x_{0},y_{0},\theta_{0})$
---
$\theta_{N}$
---

$y$
---

(a) A scenario in which the dynamical object at time $N$, i.e., $o_{N}$, does
not cover the ego vehicle reference position at time $N$.

$x,v_{x}$
---

$-c_{0}^{2}$
---

$c_{0}^{1}$
---

$y$
---

$r_{N}(x_{0},y_{0},\theta_{0})$
---
$o_{N}(x_{0},y_{0},\theta_{0})$
---
$o_{0}(x_{0},y_{0},\theta_{0})$
---

$-c_{0}^{3}$
---

(b) A scenario in which the dynamical object at time $N$, i.e., $o_{N}$,
covers the reference position at time instance $N$.

Figure 1. Two illustrations introduce the problem statement and show the
typical trade-off between driving towards a goal versus minimal-risk driving.

The trajectory of the vehicle is expected to minimize risks by abiding by the
traffic rules while respecting the presence of other road users. Furthermore,
the planner is expected to plan for a period that allows anticipating the
environment (i.e., tens of seconds down to 3-4 seconds [12]). Note that
minimizing risks while driving towards a reference $\mathbf{r}_{N}$ could
result in an optimal trajectory that does not reach the reference. To
illustrate this phenomenon, we reconsider Fig. 1(a) in a different scenario
depicted in Fig. 1(b). In this scenario, the vehicle is expected to avoid risk
by, e.g., steering the vehicle back into the center of the lane while avoiding
a dangerous situation with the object, meaning it will inherently not reach
its reference state $\mathbf{r}_{N}$.

Next to the planar movement of the vehicle (i.e., $x,y,\theta$), a notion of
time needs to be included to be able to incorporate the risk of moving
objects. As a result, the planar movement is augmented with the longitudinal
velocity of the vehicle, denoted by $v$. The main challenge of a trajectory
generator is, hence, to find a reference trajectory, $\mathcal{T}$, for the
vehicle to follow, which takes into account all the above aspects,

(2) $\mathcal{T}=(x_{k},y_{k},\theta_{k},v_{k})\>\forall k\in[0\ldots N]$

where the trajectory length is $N$ and satisfies kinematic vehicle behavior
and comfort-driven constraints while minimizing the risk to the passengers and
surrounding road users.

##  3\. Model-predictive planning for risk-averse driving

In this work, we approach the trajectory planning problem from a model-
predictive control (MPC) perspective. The characteristics of an MPC are in
line with our requirements, namely, minimizing risk (i.e., formulating the
cost as a function of risk), while satisfying kinematic vehicle dynamical
behavior and constraints (e.g., limited steering input or longitudinal
acceleration).

###  3.1. Formulation of the vehicle model

The model, used for MPC synthesis, considered in this work is the non-linear
kinematic bicycle model. This model is particularly suitable for larger
sampling times, i.e., larger time horizons while maintaining numerical
stability [9]. The model is depicted in Fig. 2 and is described in discrete
time, using forward Euler discretization with sampling time $t_{s}$. We denote
the state of the model as $\mathbf{x}=[x,y,\theta,v]$ and its inputs as
$\mathbf{u}=[\delta,a]$, and the state-update rules are denoted as follows

$x,v$
---

$y$
---

$L$
---
$L_{r}$
---
$L_{f}$
---
$\delta$
---
Figure 2. Kinematic bicycle model, with the planar reference frame in the
center of the rear axis.

(3) $\displaystyle\begin{aligned}
x_{k+1}=&x_{k}+t_{s}v_{k}\cos{\theta_{k}},\\\
y_{k+1}=&y_{k}+t_{s}v_{k}\sin{\theta_{k}},\\\
\theta_{k+1}=&\theta_{k}+t_{s}\frac{v_{k}}{L}\tan{\delta_{k}},\\\
v_{k+1}=&v_{k}+t_{s}a_{k}.\end{aligned}$

where the states $x,y,\theta$ denote the planar movement relative to the
initial condition, as denoted in Fig. 1(a), $v$ represents the velocity in the
direction of the $x$-axis, $\delta$ denotes the steering angle of the front
axle, $a$ denotes the acceleration in the direction of the $x$-axis and
finally $L$ denotes the wheelbase of the vehicle. The state vector
$\mathbf{x}$ and the input vector $\mathbf{u}$ belong to the sets

(4a) $\displaystyle\mathbf{X}=$
$\displaystyle\\{\mathbf{x}\in\mathbb{R}^{4}|\underline{\mathbf{x}}\leq\mathbf{x}\leq\overline{\mathbf{x}}\\},$
(4b) $\displaystyle\mathbf{U}=$
$\displaystyle\\{\mathbf{u}\in\mathbb{R}^{2}|\underline{\mathbf{u}}\leq\mathbf{u}\leq\overline{\mathbf{u}}\\},$

respectively. The bounds of the sets are governed by infrastructural, dynamic,
and comfort limitations.

###  3.2. Formulation of the risk

Similar to [15], we consider two sources of risk, namely, infrastructural risk
and object-imposed risk. The infrastructural risk is described as the safety
risk taken to travel across certain road boundaries [15]. We define the
infrastructural risk, at time instance $k$, by $\mathcal{U}^{I}_{k}$.
Furthermore, we define a risk-free situation, i.e.,
$\mathcal{U}^{I}_{k}\approx 0$ if and only if the vehicle is driving in the
center of a lane and assuming that all lanes are driveable. The risk-function
of the infrastructure, perceived by the vehicle at time instance $k$, can be
described as a set of structured Gaussian functions, using the polynomial road
description (1), as follows

(5) $\displaystyle\mathcal{U}^{I}_{k}=$
$\displaystyle\sum_{i=1}^{N_{d}}A_{I}e^{-\frac{\left(y_{k}-\tilde{y}_{k}^{(i)}(\tilde{x}^{(i)}_{k})\right)^{2}}{2\sigma^{2}}},\forall
i\in\mathbb{Z}^{+}$

where $N_{d}$ represents the number of road lines, $A_{I}$ represents the
amplitude of the risk and $\sigma$ represents the standard deviation of the
risk. The amplitude and standard deviation of the infrastructural risk are
assumed constant, though may depend on the type of lanes that are divided
(e.g., a lane with oncoming traffic may be riskier to drive on, hence the
amplitude may be higher). This risk-field description is a slight adaption of
the description in [16], where only a two-lane highway is considered. An
illustrational example of an infrastructural potential field is provided in
Fig. 3(a).

(a) Example infrastructural risk-field of a slightly curved two-lane highway.

(b) Example object risk-field.

Figure 3. Two illustrative examples of the proposed risk-field definitions.
$\Delta x$ denotes the $x$-distance from the vehicle’s initial position and
$\Delta y$ denotes the $y$-distance from the vehicle’s initial position.

The relative risk of dynamical objects is described using a repulsive
potential field formulated as a multivariate, potentially rotated, Gaussian
distribution as follows:

(6) $\displaystyle\mathcal{U}^{O}_{k}=$
$\displaystyle\sum_{n=1}^{N_{o}}A_{O}e^{-\frac{f_{o}^{(n)}(\mathbf{x}_{k},\mathbf{o}_{k})}{2}}$
$\displaystyle f_{o}^{(n)}=$
$\displaystyle\begin{bmatrix}x_{k}-\widehat{x}^{(n)}_{k}&y_{k}-\widehat{y}^{(n)}_{k}\end{bmatrix}R\Sigma^{-1}R^{\intercal}\begin{bmatrix}x_{k}-\widehat{x}^{(n)}_{k}\\\
y_{k}-\widehat{y}^{(n)}_{k}\end{bmatrix}$ $\displaystyle\Sigma=$
$\displaystyle\begin{bmatrix}{\sigma_{x}^{2}}&0\\\
0&{\sigma_{y}^{2}}\end{bmatrix},\quad
R=\begin{bmatrix}\cos\widehat{\theta}^{(n)}_{k}&-\sin\widehat{\theta}^{(n)}_{k}\\\
\sin\widehat{\theta}^{(n)}_{k}&\cos\widehat{\theta}^{(n)}_{k}\end{bmatrix}$

where
$(\widehat{x}_{k}^{(n)},\widehat{y}_{k}^{(n)},\widehat{\theta}_{k}^{(n)})$
represents the state of the $n$-th object, $N_{o}$ represents the number of
objects, $A_{O}\in\mathbb{R}^{N}_{d}$ represents for each object the maximal
risk, $\sigma_{x}$ represents the longitudinal standard deviation of the
perceived risk of the object and $\sigma_{y}$ represents the lateral standard
deviation of the perceived risk of the object. Note, again, that for
simplicity the amplitude and standard deviations are assumed constant, though
they can be varied based on the type of object (e.g., pedestrian versus a car)
or for example the velocity at which the ego vehicle is driving. An
illustrational example of such a risk field is provided in Fig. 3(b).

###  3.3. Formulation of the risk-averse trajectory generator

Using the definition of the vehicle dynamic model (3), its limitations (4) and
the definition of the risk field for the infrastructure (5) and objects (6),
the optimization problem for a risk-averse trajectory generator can be
composed. The cost and the constraints of a risk-averse trajectory generator
are formulated by the following model-predictive program:

$\displaystyle J^{\star}_{0\rightarrow N}(\mathbf{x}_{s})=$
$\displaystyle\min_{u,x}\sum_{k=1}^{N-1}\ell_{\mathbf{x}}\left(\mathbf{x}_{k}-\mathbf{r}_{k}\right)+\sum_{k=0}^{N-1}\ell_{\mathbf{u}}\left(\mathbf{u}_{k}\right)$
(7a) $\displaystyle+m$
$\displaystyle(\mathbf{x}_{N}-\mathbf{r}_{N})+\mathcal{U}_{k}^{O}+\mathcal{U}_{k}^{I}$
(7b) $\displaystyle\text{s.t.}\quad\mathbf{x}_{0}=$
$\displaystyle\mathbf{x}_{s}$ (7c) $\displaystyle\quad\mathbf{x}_{k+1}=$
$\displaystyle f(\mathbf{x}_{k},\mathbf{u}_{k}),\quad\forall k\geq 0$ (7d)
$\displaystyle\quad g(\mathbf{x}_{k},\mathbf{u}_{k})\leq$
$\displaystyle\>0,\quad\forall k\geq 0$

where, in (7a), the functions
$\ell_{\mathbf{x}}(\cdot),\>\ell_{\mathbf{u}}(\cdot)$ quadratically penalizes
the control error $\mathbf{x}_{k}-\mathbf{r}_{k}$, and the input
$\mathbf{u}_{k}$, respectively. Furthermore, the term $m(\cdot)$ quadratically
penalizes the terminal state at time $N$. The state $\mathbf{x}_{s}$ in (7b)
represents the state at present time. The function $f$ in (7c) represents the
non-linear dynamics as depicted in (3). Finally, the constraints in (7d)
enforce the boundedness of the state $\mathbf{x}$ and the input $\mathbf{u}$
through enforcement of the sets (4). Note, that the stage cost
$\ell_{\mathbf{x}}$ is not used in this work (i.e., equal to zero), due to the
fact that there is only the requirement to drive to a goal at time instance
$N$, hence, only the terminal cost $m(\cdot)$ is used to penalize the state.

##  4\. Simulation study

The non-linear program is implemented in a real-time operating optimization
framework CasADi [13] using the IPOPT Interior-point embedded solver [14]. The
trajectory generator is deployed in four different studies, where urban- and
highway driving in a risk-defined environment forms the basis. Throughout the
first three simulations, the same vehicle parameter set is used to show the
versatility of the algorithm without having to tune it for every separate use
case. For the urban case, slight adjustments are made. The parameter-set is
provided in Table 1, where it has to be noted, that any parameters which are
not denoted are either $0$ in the case of stage weights or not set in the case
of state bounds. In each of the test cases, the computational time is
denoted111The results are generated on an Intel Core i7-10850 2.7GHz
platform.. The optimization problem is initiated at each instance $t_{s}$,
hence, by real-time we mean that the MPC trajectory generator can generate a
new solution within a period of $t_{s}$ seconds. Finally, in the simulation,
we assume perfect tracking of the trajectory by the vehicle model.

Parameter |  Description |  Value |  Unit
---|---|---|---
$t_{s}$ |  Sampling time |  0.75 |  s
$N$ |  MPC horizon |  10 |  -
$\ell_{\mathbf{u}}$ |  Input weight |  $[1,100]$ |  -
$m$ |  Terminal weight |  $[1,0.01,0,0]$ |  -
$A_{I}$ |  Infrastructure risk amplitude |  100 |  -
$A_{O}$ |  Object risk amplitude |  1000 |  -
$\sigma$ |  Infrastructure risk deviation |  1.3 |  m
$\sigma_{x}$ |  Longitudinal risk deviation |  20 |  m
$\sigma_{y}$ |  Lateral risk deviation |  1.3 |  m
$\underline{a},\overline{a}$ |  Vehicle acceleration bounds |  $[-4,0.5]$ |  $\text{ms}^{-2}$
$\underline{\delta},\overline{\delta}$ |  Vehicle steering angle bounds |  $[-0.1,0.1]$ |  rad
$\underline{y},\overline{y}$ |  Lateral displacement bounds |  $[1,9.5]$ |  m
$\underline{v},\overline{v}$ |  Vehicle velocity bounds |  $[0,10]$ |  $\text{ms}^{-1}$
Table 1. Simulation parameters used for case I up to III and partially for
case IV.

###  4.1. Case I: single overtake

In the first case, we consider a scenario where the vehicle drives in the
rightmost lane and approaches several slow-driving vehicles. The objects are
driving at a velocity of $5\text{ms}^{-1}$ (bottom lane) and $2\text{ms}^{-1}$
(middle lane), respectively as typical in highway traffic jam situations. The
results are shown in Fig. 4 at several interesting points in time. First, the
vehicle approaches the slow-driving vehicle in the bottom lane and overtakes
the vehicle. Then, the vehicle has three choices, it can either pass the
object on the top or bottom lane, or it can slow down and drive behind the
vehicle. In this case, it is of minimal cost (i.e., a combination of risk and
benefit) to return to the original lane, such that the vehicle can, once
again, drive towards the goal point. The average computational time of the
planner in this scenario is $30.4\text{ms}$ with a standard deviation of
$5.1\text{ms}$.

Figure 4. Simulation results of case I with two vehicles driving slowly on the
two rightmost lanes, motivating a single lane change.

###  4.2. Case II: double overtake

In the second case, we consider the scenario where the vehicle starts on the
rightmost lane, overtakes a slow vehicle, and then decides whether to slow
down, or to once again overtake due to a static vehicle in the middle lane.
The objects are driving at a velocity of $5\text{ms}^{-1}$ (bottom lane) and
$0\text{ms}^{-1}$ (middle lane), respectively. The results are shown in Fig.
5. The first overtake maneuver is identical to case I. Once the vehicle
approaches the static vehicle, it can be observed that a lane change to the
right would compromise the taken risk due to the proximity of the right-lane
object prediction at horizon $N$ ahead. As a result, the MPC takes the
anticipatory decision to take another lane change to the left, at the cost of
a higher lateral distance error to the desired goal point. The average
computational time of the planner in this scenario is $30.2\text{ms}$ with a
standard deviation of $4.1\text{ms}$.

Figure 5. Simulation results of case II with two vehicles driving slowly on
the two rightmost lanes, motivating a double lane change.

###  4.3. Case III: blocked lanes, slowly driving participants

In the third case, we consider the scenario where the three-lane road is
blocked by slowly driving objects. The objects are all driving at a velocity
of $8\text{ms}^{-1}$ and their longitudinal position initial condition
deviates to challenge the algorithm. The results are shown in Fig. 6. It can
be observed that the MPC notices that there is only one real sensible
solution, which is to stay inside the lane and brake to the same velocity as
the succeeding vehicles. Note, that the steady-state euclidean distance
between the vehicles can be tuned using the longitudinal parameter
$\sigma_{x}$. The average computational time of the planner in this scenario
is $28.2\text{ms}$ with a standard deviation of $1.6\text{ms}$.

Figure 6. Simulation results of case III with three vehicles slowly driving
and blocking the lanes.

###  4.4. Case IV: VRU crossing with potential late object detection

Whereas case I to III presented highway maneuvers, the fourth and final case
will cover an urban driving setting, showing the applicability of the MPC
trajectory generator for design domains other than highway driving. In this
case, we consider the vehicle driving along a two-lane road. On this two-lane
road, a VRU is about to cross the road as depicted in Fig. 7. This case is
deemed to be highly relevant and largely responsible for VRU incidents (up to
$51\%$ [17]). The proposed simulations are twofold. First, we study how the
MPC trajectory generator acts when it sees the VRU and can anticipate its
actions based on these observations. Subsequently, we study how the trajectory
generator acts when the VRU observation is obstructed and the detection of the
VRU is late (advocated in [17] as a relevant root cause in this particular
scenario). Transitioning from a highway scenario towards an urban scenario
requires an adaptation of MPC parameters. The revised parameters, different
from Table 1, are proposed in Table 2. In an urban scenario, the maneuvers
likely require a higher tracking fidelity, for this reason, the sampling time
is reduced and, as a result, the horizon is increased to achieve a similar
look-ahead time concerning the previous cases. Moreover, the reduced
longitudinal velocity results in a change of bound and a change in the
longitudinal standard deviation of the object risk-field. Finally, an urban
scenario with two-lane roads usually allows bi-directional traffic. As a
result, the infrastructural potential field amplitude is increased (as shown
in Table 1 and also in Fig. 7 by the infrastructural risk field color). The
VRU crosses the road at a velocity of $1\text{ms}^{-1}$ and the ego vehicle
starts at an initial velocity of $8.33\text{ms}^{-1}$. First, we consider the
results without late VRU detection, as depicted in Fig. 8. It is clear that,
from $23\text{s}$ onward the vehicle starts to anticipate and brake, down to a
velocity of $3\text{ms}^{-1}$, to wait for the VRU to pass. Once the VRU has
passed, the vehicle accelerates back to its original velocity, which is a
natural and intuitive expected behavior of the vehicle. Now, consider a late
detection of the VRU (coming, e.g., from bad weather conditions) of which the
results are depicted in Fig. 9. The detection occurs at around $30\text{s}$
(i.e., the top graph), a time at which there is only $13\text{m}$ of distance
left concerning the object. At this time, the trajectory generator path
indicates that braking is no longer an option and an evasive maneuver towards
the left lane is necessary. Once the evasive maneuver is finished, the vehicle
returns to its lane while a safe distance has been kept with the VRU. The
average computational time of the planner for the early VRU detection is
$43.5\text{ms}$ with a standard deviation of $4.2\text{ms}$. The average
computational time of the planner in the late detection case is
$48.2\text{ms}$ with a standard deviation of $8.1\text{ms}$. The increased
mean computational time can be explained largely by the fact that, due to the
higher required fidelity of the path, the horizon is doubled.

Figure 7. Initial setup for the VRU simulation case IV. Parameter |  Description |  Value |  Unit
---|---|---|---
$t_{s}$ |  Sampling time |  0.38 |  s
$N$ |  MPC horizon |  20 |  -
$A_{I}$ |  Infrastructure risk amplitude |  200 |  -
$\sigma_{x}$ |  Longitudinal risk deviation |  10 |  m
$\underline{y},\overline{y}$ |  Lateral displacement bounds |  $[1,6]$ |  m
$\underline{v},\overline{v}$ |  Vehicle velocity bounds |  $[0,8.33]$ |  $\text{ms}^{-1}$
Table 2. Changed parameters w.r.t. Table 1 for case IV.

Figure 8. Simulation results of case IV with timely object detection.

Figure 9. Simulation results of case IV with late object detection.

##  5\. Conclusion

In this work, we have proposed a monolithic model-predictive trajectory
generator that, as opposed to many works in literature, does not use a set of
pre-defined decisions or trajectories to decide on actions to progress towards
a goal point while minimizing risk. This offers a simpler solution to
trajectory generation, which can be used for longer horizons while being
computationally feasible. The proposed method explicitly considers predicted
trajectories of other road users and it’s shown effective in various use
cases. The method is shown to be effective by demonstrating it in four
relevant scenarios in which the approach shows the ability to make logical and
safe decisions in both highway and urban scenarios. Furthermore, the planner
is challenged in a more difficult urban scenario, where the late detection of
a pedestrian (for example happening due to bad weather) is compared with the
timely detection. The results show that this approach can function safely
under both circumstances. Future work includes the adaptation of the MPC
problem formulation and the potential fields to incorporate even more complex
situations and be robust for measurement uncertainty, induced by predictions
or faults (for example, the lack of fault isolability in object fault
detection [18]). Moreover, adaptive potential fields will be studied to make
the trajectory generator suitable for miscellaneous environments.

##  6\. Acknowledgments

This work is supported by the EU Horizon 2020 R&D program under grant
agreement No. 861570, project SAFE-UP (proactive SAFEty systems and tools for
a constantly UPgrading road environment).

## References

  * [1] A. Papadoulis, M. Quddus and M. Imprialou  “Evaluating the safety impact of connected and autonomous vehicles on motorways”  In _Accident Analysis and Prevention_ 124, 2019, pp. 12–22
  * [2] J. De Freitas et al.  “From driverless dilemmas to more practical commonsense tests for automated vehicles”  In _Proceedings of the National Academy of Sciences of the United States of America_ 118.11, 2021
  * [3] D. González, J. Pérez, V. Milanés and F. Nashashibi  “A Review of Motion Planning Techniques for Automated Vehicles”  In _IEEE Transactions on Intelligent Transportation Systems_ 17.4, 2016, pp. 1135–1145
  * [4] M. Likhachev and D. Ferguson  “Planning long dynamically feasible maneuvers for autonomous vehicles”  In _International Journal of Robotics Research_ 28.8, 2009, pp. 933–945
  * [5] B. Gutjahr, L. Gröll and M. Werling  “Lateral Vehicle Trajectory Optimization Using Constrained Linear Time-Varying MPC”  In _IEEE Transactions on Intelligent Transportation Systems_ 18.6, 2017, pp. 1586–1595
  * [6] J. Ji, A. Khajepour, Sr. Melek and Y. Huang  “Path planning and tracking for vehicle collision avoidance based on model predictive control with multiconstraints”  In _IEEE Transactions on Vehicular Technology_ 66.2, 2017, pp. 952–964
  * [7] O. Khatib  “Real-time obstacle avoidance for manipulators and mobile robots”  In _IEEE International Conference on Robotics and Automation_ , 1985, pp. 500–505
  * [8] S. Dixit et al.  “Trajectory Planning for Autonomous High-Speed Overtaking using MPC with Terminal Set Constraints”  In _IEEE Conference on Intelligent Transportation Systems_ , 2018, pp. 1061–1068
  * [9] Vincent A. Laurense and J. Gerdes  “Long-Horizon Vehicle Motion Planning and Control Through Serially Cascaded Model Complexity”  In _IEEE Transactions on Control Systems Technology_ 30.1, 2022, pp. 166–179
  * [10] Z. Huang, Q. Wu, J. Ma and S. Fan  “An APF and MPC combined collaborative driving controller using vehicular communication technologies”  In _Chaos, Solitons and Fractals_ 89, 2016, pp. 232–242
  * [11] P. Hang et al.  “An Integrated Framework of Decision Making and Motion Planning for Autonomous Vehicles Considering Social Behaviors”  In _IEEE Transactions on Vehicular Technology_ 69.12, 2020, pp. 14458–14469
  * [12] S.C. Calvert, W.J. Schakel and J.W.C. Lint  “A generic multi-scale framework for microscopic traffic simulation part II – Anticipation Reliance as compensation mechanism for potential task overload”  In _Transportation Research Part B: Methodological_ 140, 2020, pp. 42–63
  * [13] Joel A E Andersson et al.  “CasADi – A software framework for nonlinear optimization and optimal control”  In _Mathematical Programming Computation_ , 2018
  * [14] A. Wächter and L.T. Biegler  “On the implementation of an interior-point filter line-search algorithm for large-scale nonlinear programming”  In _Mathematical Programming_ 106.1, 2006, pp. 25–57
  * [15] L. Li et al.  “Dynamic Driving Risk Potential Field Model under the Connected and Automated Vehicles Environment and Its Application in Car-Following Modeling”  In _IEEE Transactions Intelligent Transportation Systems_ 23.1, 2022, pp. 122–141
  * [16] F. Siebenrock, M. Gunther and S. Hohmann  “LTV-MPC based trajectory planning considering uncertain object prediction through adaptive potential fields”  In _IEEE Conference on Control Technology and Applications_ , 2020, pp. 666–672
  * [17] L. Yue et al.  “In-depth approach for identifying crash causation patterns and its implications for pedestrian crash prevention”  In _Journal of Safety Research_ 73, 2020, pp. 119–132
  * [18] Chris Ploeg et al.  “Anomaly Detection from Cyber Threats via Infrastructure to Automated Vehicle”  In _European Control Conference_ , 2021, pp. 1788–1794