# Surface plasmonics of Weyl semimetals

Xin Lu  Laboratoire de Physique des Solides, Univ. Paris-Sud, Université Paris
Saclay, CNRS, UMR 8502, F-91405 Orsay Cedex, France  Dibya Kanti Mukherjee
Laboratoire de Physique des Solides, Univ. Paris-Sud, Université Paris Saclay,
CNRS, UMR 8502, F-91405 Orsay Cedex, France  Mark O. Goerbig  Laboratoire de
Physique des Solides, Univ. Paris-Sud, Université Paris Saclay, CNRS, UMR
8502, F-91405 Orsay Cedex, France

###### Abstract

Smooth interfaces of topological systems are known to host massive surface
states along with the topologically protected chiral one. We show that in Weyl
semimetals these massive states, along with the chiral Fermi arc, strongly
alter the form of the Fermi-arc plasmon, Most saliently, they yield further
collective plasmonic modes that are absent in a conventional interfaces. The
plasmon modes are completely anisotropic as a consequence of the underlying
anisotropy in the surface model and expected to have a clear-cut experimental
signature, e.g. in electron-energy loss spectroscopy.

_Introduction._ —Weyl semimetals (WSMs) are often considered as a three-
dimensional version of graphene since their low-energy
$\mathbf{k}\cdot\mathbf{p}$ Hamiltonian is described by the massless Weyl
equation with linear energy dispersion known as Weyl cones [1]. WSMs must have
an even number of Weyl cones: at least four for inversion symmetry broken WSM
and two for time-reversal symmetry broken ones. Previous studies have shown
that the bulk dielectric properties such as Friedel oscillations and
(magneto-)plasmon in WSMs [2, 3, 4, 5, 6] are different from those of graphene
[7, 8] due to the increased dimensionality. For example, the bulk plasmon’s
dispersion in WSMs is gapped and parabolic in momentum while it follows a
gapless square-root dispersion in graphene. The chiral anomaly in WSMs [9] may
be probed by the bulk plasmon [10] whose dispersion depends on the chirality-
resolved chemical potential. As metals, WSMs can also host surface plasmons
[11] and surface plasmon polaritons [12, 13] by using Maxwell’s equations in
the bulk with a topological Chern-Simons $\theta$-term [14]. In particular, in
ferromagnetic WSMs, due to its gapless spectrum and large Berry curvature [15,
16], electromagnetic waves propagate non-reciprocally, i.e., one direction is
preferred.

Localized surface states can also give birth to original surface plasmons. In
WSMs, topologically protected Fermi-arc (FA) states connecting two Weyl cones
emerge on the surface due to the bulk-edge correspondence: the presence of
topologically protected edge states is dictated by the topological invariant
of the twisted bulk band structure. The FA states have been shown to induce a
chiral linear FA surface plasmon with total non-reciprocity [17, 18, 19, 20,
21, 22, 23], i.e. it propagates only in one direction determined by the the
chirality of the FA dispersion. Linear dispersion and total non-reciprocity
are highly desirable for further plasmonic applications. Furthermore, a smooth
surface of a topological material is known to host massive states called
Volkov-Pankratov (VP) states [24, 25, 26, 27, 28, 29], along with the
protected topological chiral states. These gapped bands can be visualized as
pseudo-Landau levels of the system where the smoothness of the interface is
modeled as a pseudo-magnetic field [30]. Though not protected topologically,
their presence may heavily modify the transport [31] and magneto-optical
properties [32, 33] of surfaces of topological materials.

Figure 1: Various optical excitations involving the surface bands for a fixed
transverse momentum. The chiral FA is denoted by the red line whereas the
massive VP bands are shown in green. As discussed in the main text, for
$q_{z}=0$, only $n\to\pm n$ excitations are allowed. Inset: global view on
surface bands connecting two Weyl nodes. Along the transverse momentum, the
quasi-1D surface bands do not disperse.

Here, we show within a simple single-boundary model how the FA and the VP
states conspire to give rise to new plasmon modes on a smooth surface of a WSM
applying random phase approximations (RPA). We confirm that the FA plasmon is
chiral and exhibits strong anisotropy and a singularity at zero momentum [17,
18, 20, 21, 22, 23] because the two-dimensional (2D) dispersion of the FA band
evolves into an effectively one-dimensional (1D) one: the energy disperses
linearly perpendicular to the $z$-direction connecting two Weyl nodes and
remains almost constant along $z$. A spectacular consequence of this
anisotropy is the finite gap which the FA plasmon acquires at $q_{z}=0$ and
that vanishes when the longitudinal wavevector $q_{z}\neq 0$. Moreover, a VP
intraband plasmon appears when the chemical potential is above the minimum of
the first VP band. Somewhat surprisingly, this plasmon is also non-reciprocal
in spite of the $k_{y}\leftrightarrow-k_{y}$ symmetric dispersion of the VP
bands. As we show below, this is due to coupling to the chiral FA state. We
also find a gapped plasmon mode that stems from excitations between two VP
bands of same band index and that we call VP interband plasmon.

_Smooth WSM surface_ —We consider a smooth interface in the $x$-direction
between a time-reversal breaking WSM and a trivial insulator modeled by the
Hamiltonian [34, 35, 33]

$\displaystyle H$
$\displaystyle=v(k_{x}\sigma_{x}+k_{y}\sigma_{y})+\left(\frac{k_{z}^{2}}{2m}-\Delta+2\frac{\Delta}{\ell}x\right)\sigma_{z},$
(1)

where all the material-related parameters are positive and henceforth we use
$\hbar=1$ for notational simplicity. Without the $x$-dependent term, this is
the simplest model for a time-reversal breaking two-node WSM with Weyl nodes
at $\mathbf{k}=\eta\sqrt{2m\Delta}\hat{z}$ with opposite chirality $\eta=\pm
1$. The spatially variant gap parameter describes explicitly how the inverted
band gap at the center of Brillouin zone is closed and reopened across the
interface of width $\ell$ from WSM ($x<0$) to trivial insulator ($x>\ell$).
The smoothness of the surface can be viewed as an effective chiral pseudo-
magnetic field $\textbf{B}_{P}=-\eta 2\Delta/ev\ell\hat{y}$ that couples to
the two Weyl nodes of opposite chirality with respective signs. Hamiltonian
(1) can be diagonalized by introducing creation and annihilation operators
constructed from linear combinations of the $k_{x}$ and the $x$ dependent
terms. Thus, the effective surface bands are reminiscent of Landau levels
following the dispersion [36]

$\displaystyle E_{n}^{\lambda}(k_{y})=\lambda
v\sqrt{k_{y}^{2}+\frac{2n}{\ell_{S}^{2}}}=\lambda\sqrt{v^{2}k_{y}^{2}+ne_{0}^{2}},$
(2)

for $n\geq 1$, where $\lambda=\pm$ is the band index and the smoothness of the
surface has been encoded in an effective magnetic length
$\ell_{S}=1/\sqrt{eB_{P}}$. The VP band gap $e_{0}=\sqrt{2}v/\ell_{S}$, which
is the separation between the $n=0$ and $n=1$ bands at $\textbf{k}=0$, sets
the characteristic energy scale of this surface model. The FA is naturally
described here by the $n=0$ band with

$\displaystyle E_{0}(k_{y})=vk_{y}$ (3)

and breaks the symmetry $k_{y}\to-k_{y}$, its counterpart with opposite sign
of the dispersion being localized at the other surface of the WSM that we do
not consider here. The FA state is independent of the surface details such as
its smoothness, i.e. the band dispersion does not depend on $\ell$, indicating
its topological nature. However, the $n\geq 1$ VP bands depend strongly on the
surface modeling. In the sharp-surface limit ($\ell\to 0$), the VP bands rise
up in energy and eventually merge with the bulk states when
$v\sqrt{2n}/\ell_{S}\sim\Delta$, while only the FA state survives.

From Eq. (2), we can see that the VP bands are completely flat in the
$k_{z}$-direction until they hybridize with the bulk Weyl cones, as shown in
the inset of Fig. 1. In spite of being embedded in a 2D $(k_{y},k_{z})$
manifold, the VP bands are effectively 1D and thus exhibit van Hove
singularities in the density of states at the band extrema. The underlying 2D
nature and the $k_{z}$-dependence is encoded in the location of the surface
states: along the interface, they have a Gaussian profile of a characteristic
width $\ell_{S}$ centered at $\langle x\rangle=B_{P}(\Delta-k_{z}^{2}/2m)$.
This, as shown below, results in non-diagonal overlap matrix elements for
excitations in the $k_{z}$-direction.

_Quasi-two-dimensional RPA._ – In order to analyze the behavior of surface
electrons, consider the non-interacting dynamical polarization

$\displaystyle\chi^{(0)}_{i,j}(\mathbf{q},\omega)=$
$\displaystyle\sum_{i,j}\frac{1}{V}\sum_{\mathbf{k}}\frac{f_{D}(E_{i}(\mathbf{k}))-f_{D}(E_{j}(\mathbf{k+q}))}{\omega+E_{i}(\mathbf{k})-E_{j}(\mathbf{k+q})+i\delta}$
$\displaystyle\times|F_{i,j}(\mathbf{k},\mathbf{k+q})|^{2},$ (4)

where $\delta=0^{+}$ and the $i,j$ indices are shorthand notations for both
band labels $n$ and $\lambda$. In general, the overlap matrix element
$F_{i,j}$ is not diagonal because of the abovementioned $k_{z}$-dependence of
the eigenstates so that $\chi^{(0)}_{i,j}$ is generally a tensor. However, for
$q_{z}=0$, the particle-hole excitations are also 1D and $F_{i,j}$ becomes
diagonal, meaning that only excitations from $n$ to $\pm n$ are possible. The
RPA dielectric function then retrieves its usual form

$\displaystyle\epsilon^{\text{RPA}}(q_{y},\omega)=1-V_{\text{2D}}(q_{y})\chi^{(0)}(q_{y},\omega),$
(5)

where $V_{\text{2D}}(q_{y})=e^{2}/2\epsilon_{0}\epsilon_{r}|q_{y}|$, and
$\chi^{(0)}$ is the non-interacting charge susceptibility, in terms of the
environmental dielectric constant $\epsilon_{r}$ [36].

When $q_{z}\neq 0$, the overlap matrix element
$F_{i,j}(\mathbf{k},\mathbf{k}+\mathbf{q})$ is more involved. However, in the
long-wavelength limit, the off-diagonal term $F_{i,j}$ is proportional to
$q_{z}^{|n_{i}-n_{j}|}$ [36] so that the $n\to\pm n$ excitations still remain
the leading contributions to the charge susceptibility. Nevertheless, due to
the complicated form of the off-diagonal terms, we cannot factorize the
Coulomb interaction operator [36] even in the long-wavelength limit when
several VP bands are present. We therefore consider only the chiral FA and the
$n=\pm 1$ VP bands (threeband model), where [36]

$\displaystyle\chi^{(0)}(\mathbf{q},\omega)=\sum_{i,j}\chi^{(0)}_{i,j}(\mathbf{q},\omega),$
(6)

and $\chi^{(0)}_{i,j}$ are the contributions by the excitation
$(n_{i},\lambda_{i})\to(n_{j},\lambda_{j})$. Accordingly, we can generalize
$(q_{y},\omega)$ in Eq. (5) to $(\mathbf{q},\omega)$.

_Plasmons_ —We summarize our results for $q_{z}=0$ in Fig. 2 and for
$q_{z}\ell_{S}=0.2$ in Fig. 3 where we numerically calculate the profile of
$-\text{Im}(\chi^{(0)})$ in the $(q_{y}>0,\omega>0)$-plane, for different
values of the chemical potential $\mu$ and a given disorder amplitude
$\delta=0.01$. Indeed, the imaginary part of $\chi^{(0)}$ bears important
information about possible electronic excitations and therefore damping of the
plasmon modes, indicated by the red dashed lines (zeros of the real part of
$\epsilon^{RPA}$). The plasmon modes are only long-lived and undamped in the
black regions where $\text{Im}(\chi^{(0)})=0$. Within the abovementioned
threeband model, one obtains three particle-hole continua, with
$-\text{Im}(\chi^{(0)})\neq 0$, shown in Fig. 1. The FA particle-hole
continuum extends linearly from in the $(q_{y}>0,\omega>0)$-plane. A second
particle-hole spectrum is delimited from below by
$\omega>\sqrt{4e_{0}^{2}+v^{2}q_{y}^{2}}$ due to interband excitations
involving the VP conduction and valence bands $n=\pm 1$, respectively.
However, it vanishes at small momenta because the eigenstates associated with
the VP conduction and valence bands are orthogonal at $\mathbf{q}=0$.

Figure 2: Profile of the imaginary part of the non-interacting dynamical
polarization $-\text{Im}(\chi^{(0)})$ in the $(q_{y},\omega)$ phase space at
$\mu=0.01e_{0},0.90e_{0},1.01e_{0}$ and $1.30e_{0}$ for $q_{z}=0$. The zeros
of the real part of $\epsilon^{\text{RPA}}$ (red dashed lines) indicate the
plasmon modes.

As we increase $\mu$ above the VP conduction band [see Figs. 2(c), (d) and
3(c), (d)], the poles of the FA excitations remain unchanged because of the
linear FA dispersion whereas that of the interband particle-hole continua gets
heavily modified due to Pauli blocking at the conduction band minima. At low
frequencies, intraband excitations of the $n=+1$ VP band induce a third
particle-hole continuum. With $\mu>e_{0}$ just above the conduction band
minimum, the VP band is approximately parabolic, and its quasi-1D character is
apparent in the form of the particle-hole intraband spectrum with its typical
exclusion dome for $q\ell_{S}<2k_{F}$ with $k_{F}$ defined as
$\mu=\sqrt{v^{2}k_{F}^{2}+e_{0}^{2}}$.

Figure 3: Same as Fig. 2 for $q_{z}\ell_{S}=0.2$. The zeros of the real part
of $\epsilon^{\text{RPA}}$ (red dashed lines) indicate the plasmon modes.

In the $q_{z}=0$ limit, two plasmon modes are present for $\mu<e_{0}$ as we
can see in Figs. 2(a), (b). The first one is the linearly dispersing FA
plasmon with a gap at $q_{y}=0$, in agreement with theoretical approaches
using classical electrodynamics [17], hydrodynamic description [22], or
quantum-mechanical calculations [18, 20, 21, 23]. From the zeros of the real
part of the equation $\epsilon^{\text{RPA}}(q_{y},\omega)=0$, we find the FA-
plasmon disersion

$\displaystyle\omega\approx\text{sgn}(q_{y})\frac{k_{0}e^{2}}{4\pi^{2}\epsilon_{0}\epsilon_{r}}+\left(1+\frac{2k_{F}}{\sqrt{k_{F}^{2}+\frac{2}{\ell_{S}^{2}}}}\delta_{n_{F},1}\right)vq_{y},$
(7)

where $2k_{0}=2\sqrt{2m\Delta}$ is the separation between two Weyl nodes in
the bulk, $\epsilon_{0}$ is the vacuum and $\epsilon_{r}$ the relative
permittivity, while $n_{F}$ is the integer part of the ratio between $\mu$ and
$e_{0}$. Let us first focus on the case where $n_{F}=0$. For positive
$\omega$, the FA plasmon is allowed to propagate only in the direction of
positive $q_{y}$, due to the chirality of the FA. For the usual Coulomb
potential, recall that the 1D and 2D plasmon dispersions are linear and
square-root, respectively. In spite of the the quasi-1D nature of the FA, the
Coulomb potential remains 2D here, and one might naively expect a square-root
plasmon dispersion. Surprisingly, this is not the case, and one finds a linear
gapped plasmon mode thanks to its chiral nature. As simple it is, Eq. (7)
accurately describes the mode found numerically in Fig. 2(a), (b), even when
the $n=\pm 1$ VP states are retained in the calculation. We emphasize that the
experimentally measurable FA plasmon gap,

$\displaystyle\gamma=\frac{k_{0}e^{2}}{4\pi^{2}\epsilon_{0}\epsilon_{r}},$ (8)

yields direct information about the separation $2k_{0}$ of the bulk Weyl
nodes.

VP intraband excitations significantly modify the dispersion of the FA plasmon
when $n_{F}=1$. The dispersion of FA plasmon at large $q_{y}$ remains almost
linear with the same velocity $v$. At small $q_{y}$, although the gap sticks
to the same value $\gamma$ in Eq. (8), the FA plasmon acquires an enhanced
velocity that can be further boosted by increasing the chemical potential as
shown in Fig. 2(c),(d). However, the change of velocity for finite $q_{y}$, as
seen in Fig. 2, is not captured by Eq. (7).

Since the bands are effectively 1D, the particle-hole continua are independent
of $q_{z}$ if the coupling between FA and VP bands is omitted, as shown in
Fig. 3. It is indeed legitimate to neglect this coupling in the long-
wavelength limit where it scales as $\sim(q_{z}\ell_{S})^{2}$. However, the
plasmon dispersion gets strongly modified for $q_{y}<q_{z}$ due the
$q_{z}$-dependence of the Coulomb interaction. Neglecting a small
hybridization between the VP bands and the FA [36], the FA-plasmon gap at
small momenta becomes

$\displaystyle\gamma^{\prime}\approx\gamma\frac{q_{y}}{\sqrt{q_{y}^{2}+q_{z}^{2}}}.$
(9)

When $q_{z}=0$, Eqs. (9) and (8) coincide, and the FA plasmon is gapped as
shown earlier. However when $q_{z}\neq 0$, the gap vanishes at $q_{y}=0$ as a
consequence of the strong anisotropy of the FA state, which is only quasi-1D
but embedded in a 2D manifold. When $q_{y}\gg q_{z}$, the FA plasmon disperses
again linearly with slope $v$. This is further validated by our numerical
calculations [see Fig. 3]: the FA plasmon gap vanishes when
$\textbf{q}=q_{z}\hat{z}$. This singular behaviour of the gap at
$\mathbf{q}=0$ is also reported in Refs. [17, 18, 20, 21, 22, 23].

The second plasmon mode in Fig.2(a), (b) is the VP interband plasmon, which
stems mostly from the $n=\pm 1$ interband excitations. It is also gapped and
starts at a finite momentum for the same orthogonality reason that makes the
spectral weight of the particle-hole continuum vanishingly small at $q_{y}\sim
0$, which makes sustained plasmonic osciallations impossible. The interband VP
plasmon mode lies in the VP interband particle-hole region and is thus Landau-
damped. However, since the amplitude of $-\text{Im}(\chi^{(0)})$ drops at high
energy, this plasmon may be visible as an additional bump in Electron Energy
Loss Spectroscopy (EELS), as we show later.

It is interesting to point out that the spectrum in Fig. 2(b), where
$\mu=0.9e_{0}$, is exactly the same as that of Fig. 2(a), where the chemical
potential is very close to the charge neutrality point. Indeed the interband
excitations between the VP bands are unchanged as long as the chemical
potential remains between the two VP bands.

A third plasmon mode emerges when $\mu>e_{0}$ [see Figs. 2(c), (d) and 3(c),
(d)]. Interestingly, this mode exists in a region delimited by the particle-
hole continua of the FA and the VP conduction bands. It starts at small but
finite momentum and its energy disperses along with the upper boundary of the
intraband continuum and eventually gets merged in it at larger momentum. One
may naively think that this VP intraband plasmon originates only from
intraband band excitations and has a square-root dispersion at small momenta
[36]. However, our numerical calculations invalidate this picture, and one
needs to take into account the other particle-hole continua, namely the linear
one associated with the FA, which prohibits such a square-root dependence of
an undamped plasmon. Moreover, remote VP interband excitations do not only
modify $\epsilon_{r}$ in the low-energy modes because of the diverging density
of states when the chemical potential crosses a VP conduction band. This
modifies significantly dynamical screening and, as shown in Fig. 2, the VP
intraband plasmon acquires positive energy only at non-zero finite momentum
and disperses linearly with a velocity smaller than $v$. Increasing $\mu$ from
$1.01$ to $1.3e_{0}$, the exclusion dome at low frequencies becomes wider. The
available phase space for the VP intraband plasmon between the FA and the VP
intraband continua has reduced even further so that this plasmon might be less
visible at larger values of $\mu$.

_Non-reciprocity._ — To show what one can see in experiments, we plot in Fig.
4 electron loss function $-\Im[1/\epsilon^{\text{RPA}}]$, measurable by EELS,
in $(q_{y},\omega>0)$-plane with intensity indicated by colorbar. One of the
intriguing properties of FA is the non-reciprocity of FA plasmon, reflecting
the chiral nature of FA state. Therefore, we should also study
$-\Im[1/\epsilon^{\text{RPA}}]$ for $(q_{y}<0,\omega>0)$. The result for
$\omega<0$ can be found easily by reversing simultaneously the sign of
$\omega$ and $\mathbf{q}$ in known results. As shown in Fig. 4 where
$\mu=1.01e_{0}$, the FA plasmon is completely absent when $q_{y}<0$ as well as
the corresponding particle-hole continuum. Being non-reciprocal, FA plasmon
only propagates in one direction with fixed velocity, highly desirable for
applications. Strikingly, also the VP intraband plasmon is non-reciprocal even
if it involves the $k_{y}\leftrightarrow-k_{y}$ symmetry of the $n=\pm$ VP
bands [see Eq. (2)]: it has a different dispersion for $q_{y}<0$, which can be
calculated analytically there

$\displaystyle\omega\approx\text{sgn}(-q_{y})v\left(|q_{y}|+\frac{2k_{F}}{\sqrt{k_{F}^{2}+\frac{2n}{\ell_{S}^{2}}}}\sqrt{q_{y}^{2}+q_{z}^{2}}\right).$
(10)

Contrary to $q_{y}>0$, it starts from the origin of $(q_{y},\omega)$ and
disperses with a velocity larger than $v$. This non-reciprocity is a
consequence of the hybridization with the FA mode and particle-hole continuum,
which is in close vicinity of the intraband VP plasmon for $q_{y}>0$ but
further well separated in energy for $q_{y}<0$ [36]. The chirality of the FA
modes thus induces a non-reciprocity in the other excitations due to their
mutual coupling. This can also be seen in the VP interband plasmon, where the
starting point moves to higher energies and larger momenta. As anticipated
above, the VP interband plasmon is submerged amid particle-hole continuum but
nevertheless visible on EELS.

Figure 4: Electron loss function at $\mu=1.01e_{0}$ for
$q_{y}\ell_{S}\in[-2.0,2.0]$: left for $q_{z}=0$ and right for
$q_{z}\ell_{S}=0.2$. Yellow dashed lines show the symmetric VP intraband
plasmon in the absence of the FA state. Three plasmon modes are all non-
reciprocal and visible.

_Discussions._ — We have investigated the effect of surface smoothness on the
charge oscillation spectrum of a WSM surface. Within RPA calculations, we
observe the emergence of two collective modes stabilized by the inter and
intra VP band excitations, in addition to the FA plasmon. The plasmons exhibit
anisotropy and non-reciprocity inherited from the underlying surface model.
Our findings could be verified experimentally, e.g. in EELS, which in addition
to a proof of these plasmons could probe the chirality of the FA.
Furthermore,the plasmon gap in Eq. (8) gives us a direct experimental measure
of the separation between the Weyl nodes.

We acknowledge financial support from Agence Nationale de la Recherche (ANR
project “Dirac3D”) under Grant No. ANR-17-CE30-0023.

## References

  * Armitage  _et al._ [2018] N. P. Armitage, E. J. Mele, and A. Vishwanath, Rev. Mod. Phys. 90, 015001 (2018).
  * Das Sarma and Hwang [2009] S. Das Sarma and E. H. Hwang, Phys. Rev. Lett. 102, 206412 (2009).
  * Lv and Zhang [2013] M. Lv and S.-C. Zhang, International Journal of Modern Physics B 27, 1350177 (2013), https://doi.org/10.1142/S0217979213501774 .
  * Panfilov  _et al._ [2014] I. Panfilov, A. A. Burkov, and D. A. Pesin, Phys. Rev. B 89, 245103 (2014).
  * Hofmann and Das Sarma [2015] J. Hofmann and S. Das Sarma, Phys. Rev. B 91, 241108 (2015).
  * Giri  _et al._ [2020] D. Giri, D. K. Mukherjee, S. Verma, H. A. Fertig, and A. Kundu, Dipolar optical plasmon in thin-film Weyl semimetals (2020), arXiv:2011.06862 [cond-mat.mes-hall] .
  * Hwang and Das Sarma [2007] E. H. Hwang and S. Das Sarma, Phys. Rev. B 75, 205418 (2007).
  * Wunsch  _et al._ [2006] B. Wunsch, T. Stauber, F. Sols, and F. Guinea, New Journal of Physics 8, 318 (2006).
  * Son and Spivak [2013] D. T. Son and B. Z. Spivak, Phys. Rev. B 88, 104412 (2013).
  * Zhou  _et al._ [2015] J. Zhou, H.-R. Chang, and D. Xiao, Phys. Rev. B 91, 035114 (2015).
  * Ritchie [1957] R. H. Ritchie, Phys. Rev. 106, 874 (1957).
  * Hofmann and Das Sarma [2016] J. Hofmann and S. Das Sarma, Phys. Rev. B 93, 241402 (2016).
  * Tamaya  _et al._ [2019] T. Tamaya, T. Kato, K. Tsuchikawa, S. Konabe, and S. Kawabata, Journal of Physics: Condensed Matter 31, 305001 (2019).
  * Zyuzin and Burkov [2012] A. A. Zyuzin and A. A. Burkov, Phys. Rev. B 86, 115133 (2012).
  * Kotov and Lozovik [2018] O. V. Kotov and Y. E. Lozovik, Phys. Rev. B 98, 195446 (2018).
  * Pellegrino  _et al._ [2015] F. M. D. Pellegrino, M. I. Katsnelson, and M. Polini, Phys. Rev. B 92, 201407 (2015).
  * Song and Rudner [2017] J. C. W. Song and M. S. Rudner, Phys. Rev. B 96, 205443 (2017).
  * Andolina  _et al._ [2018] G. M. Andolina, F. M. D. Pellegrino, F. H. L. Koppens, and M. Polini, Phys. Rev. B 97, 125431 (2018).
  * Lošić [2018] Ž. B. Lošić, Journal of Physics: Condensed Matter 30, 365003 (2018).
  * Adinehvand  _et al._ [2019] F. Adinehvand, Z. Faraei, T. Farajollahpour, and S. A. Jafari, Phys. Rev. B 100, 195408 (2019).
  * Chen  _et al._ [2019] Q. Chen, A. R. Kutayiah, I. Oladyshkin, M. Tokman, and A. Belyanin, Phys. Rev. B 99, 075137 (2019).
  * Gorbar  _et al._ [2019] E. V. Gorbar, V. A. Miransky, I. A. Shovkovy, and P. O. Sukhachov, Phys. Rev. B 99, 155120 (2019).
  * Ghosh and Timm [2020] S. Ghosh and C. Timm, Phys. Rev. B 101, 165402 (2020).
  * Volkov and Pankratov [1985] B. Volkov and O. Pankratov, JETP Lett. 42, 178 (1985).
  * Pankratov  _et al._ [1987] O. Pankratov, S. Pakhomov, and B. Volkov, Solid state communications 61, 93 (1987).
  * Tchoumakov  _et al._ [2017a] S. Tchoumakov, V. Jouffrey, A. Inhofer, E. Bocquillon, B. Plaçais, D. Carpentier, and M. O. Goerbig, Phys. Rev. B 96, 201302 (2017a).
  * Inhofer  _et al._ [2017] A. Inhofer, S. Tchoumakov, B. A. Assaf, G. Fève, J. M. Berroir, V. Jouffrey, D. Carpentier, M. O. Goerbig, B. Plaçais, K. Bendias, D. M. Mahler, E. Bocquillon, R. Schlereth, C. Brüne, H. Buhmann, and L. W. Molenkamp, Phys. Rev. B 96, 195104 (2017).
  * Alspaugh  _et al._ [2020] D. J. Alspaugh, D. E. Sheehy, M. O. Goerbig, and P. Simon, Phys. Rev. Research 2, 023146 (2020).
  * van den Berg  _et al._ [2020a] T. L. van den Berg, A. De Martino, M. R. Calvo, and D. Bercioux, Phys. Rev. Research 2, 023373 (2020a).
  * Lu and Goerbig [2020] X. Lu and M. O. Goerbig, Phys. Rev. B 102, 155311 (2020).
  * van den Berg  _et al._ [2020b] T. L. van den Berg, M. R. Calvo, and D. Bercioux, Phys. Rev. Research 2, 013171 (2020b).
  * Lu and Goerbig [2019] X. Lu and M. O. Goerbig, EPL (Europhysics Letters) 126, 67004 (2019).
  * Mukherjee  _et al._ [2019] D. K. Mukherjee, D. Carpentier, and M. O. Goerbig, Physical Review B 100, 195412 (2019).
  * Okugawa and Murakami [2014] R. Okugawa and S. Murakami, Phys. Rev. B 89, 235315 (2014).
  * Tchoumakov  _et al._ [2017b] S. Tchoumakov, M. Civelli, and M. O. Goerbig, Phys. Rev. B 95, 125306 (2017b).
  * [36] See supplementary material for more details.

Supplemental Materials for “Surface plasmonics of Weyl semimetals”

##  I S1. Solution of the Hamiltonian for the smooth surface of a time-
reversal-broken Weyl semimetal

In all calculations below, we use $\hbar=1$ for notational simplicity. The
basic Hamiltonian for an interface between a Weyl semimetal (WSM) and a
(trivial) band insulator (“vacuum”)

$\displaystyle
H=v(k_{x}\sigma_{x}+k_{y}\sigma_{y})+\left(\frac{k_{z}^{2}}{2m}-\Delta+\frac{2\Delta}{\ell}x\right)\sigma_{z}$
(S1)

reads, after a unitary transformation,

$\displaystyle
H_{T}=v\begin{bmatrix}-k_{y}&\frac{\sqrt{2}}{\ell_{S}}\hat{a}\\\
\frac{\sqrt{2}}{\ell_{S}}\hat{a}^{\dagger}&k_{y}\end{bmatrix}$ (S2)

where

$\displaystyle\ell_{S}$ $\displaystyle=\sqrt{\ell\frac{v}{2\Delta}}$ (S3)
$\displaystyle\hat{a}^{\dagger}$
$\displaystyle=\frac{\ell_{S}}{\sqrt{2}}\left(k_{x}+i\frac{x-\langle
x\rangle}{\ell_{S}^{2}}\right),$ (S4)

and $\langle x\rangle=(\Delta-k_{z}^{2}/2m)\ell/2\Delta$ determines the
average position of the surface state. It is important to note that $\langle
x\rangle$ depends on the surface momenta and we explicitly mark this
dependence on the spinor components $|n,k_{z}\rangle$. The eigenstates are
thus of the form as

$\displaystyle|\psi^{\lambda}_{n}\rangle$
$\displaystyle=\frac{1}{\sqrt{2}}\begin{pmatrix}u_{n,\lambda}(k_{y})|n-1,k_{z}\rangle\\\
\lambda
v_{n,\lambda}(k_{y})|n,k_{z}\rangle\end{pmatrix}=\frac{1}{\sqrt{2}}\begin{pmatrix}\sqrt{1-\lambda\frac{k_{y}}{\sqrt{k_{y}^{2}+\frac{2n}{\ell_{S}^{2}}}}}|n-1,k_{z}\rangle\\\
\lambda\sqrt{1+\lambda\frac{k_{y}}{\sqrt{k_{y}^{2}+\frac{2n}{\ell_{S}^{2}}}}}|n,k_{z}\rangle\end{pmatrix}$
if $n\geq 1$ (S5) $\displaystyle|\psi_{0}\rangle$
$\displaystyle=\begin{pmatrix}0\\\ |0,k_{z}\rangle\end{pmatrix}$ if $n=0$ (S6)

where $\lambda=\pm 1$ is band index and $|n,k_{z}\rangle$ is the eigenstates
of quantum harmonic oscillator defined by $\hat{a}^{\dagger}(k_{z})$ and
$\hat{a}(k_{z})$. The energy is thus

$\displaystyle E_{n}^{\lambda}$ $\displaystyle=\lambda
v\sqrt{k_{y}^{2}+\frac{2n}{\ell_{S}^{2}}}$ if $n\geq 1$ (S7) $\displaystyle
E_{0}$ $\displaystyle=vk_{y}$ if $n=0$ (S8)

Note that the FA state ($n=0$) has linear dispersion on $k_{y}$ breaking the
parity symmetry $k_{y}\to-k_{y}$, which justifies the name of chiral state,
and the VP states preserve this symmetry. Although the eigenstates live on a
two-dimensional manifold $(k_{y},k_{z})$, their energies disperse only in the
$k_{y}$-direction, the direction perpendicular to the interface and the line
connecting the Weyl nodes at $k_{z}^{W}=\pm k_{0}=\pm\sqrt{2\Delta m}$ in
reciprocal space. The bands are thus quasi-1D and the density of states of the
VP states diverges at the band extremum. The dependence of $k_{z}$ is only
encoded in the eigenstates, more precisely the cyclotron center of
$|n\rangle$.

##  II S2. Non-interacting dynamical polarization

In this section, we give the analytical expressions of the non-interacting
charge susceptibility. In the multi-band model, the charge susceptibility
$\chi^{(0)}_{1,2}$ is in general a tensor

$\displaystyle\chi^{(0)}_{1,2}(\mathbf{q},\omega)=\sum_{1,2}\frac{1}{V}\sum_{\mathbf{k}}\frac{f_{D}(E_{1}(\mathbf{k}))-f_{D}(E_{2}(\mathbf{k+q}))}{\omega+E_{1}(\mathbf{k})-E_{2}(\mathbf{k+q})+i\eta}\times|F_{1,2}(\mathbf{k},\mathbf{k+q})|^{2}$
(S9)

where the number indices are shorthand notation of all degrees of freedom,
$j=(n_{j},\lambda_{j})$ except the momentum $\mathbf{k}$, i.e. $E_{1}$ is the
energy of the band with band order $n_{1}$ and band index $\lambda_{1}$, and
$F_{1,2}(\mathbf{k},\mathbf{k+q})$ is the overlap function

$\displaystyle|F_{1,2}(\mathbf{k},\mathbf{k+q})|^{2}$ $\displaystyle=|\langle
1,\mathbf{k}|2,\mathbf{k+q}\rangle|^{2}$ (S10)
$\displaystyle=\frac{1}{4}|u_{1}^{*}(k_{y})u_{2}(k_{y}+q_{y})R_{n_{1}-1,n_{2}-1}(k_{z},k_{z}+q_{z})+\lambda_{1}\lambda_{2}v_{1}^{*}(k_{y})v_{2}(k_{y}+q_{y})R_{n_{1},n_{2}}(k_{z},k_{z}+q_{z})|^{2}$
(S11)

where

$\displaystyle R_{n_{1},n_{2}}(k_{z},k_{z}+q_{z})$ $\displaystyle=\langle
n_{1},k_{z}|n_{2},k_{z}+q_{z}\rangle$ (S12)
$\displaystyle=\sqrt{\frac{n_{2}!}{n_{1}!}}\alpha^{n_{1}-n_{2}}L_{n_{2}}^{(n_{1}-n_{2})}(|\alpha|^{2})e^{-\frac{|\alpha|^{2}}{2}}\quad\text{if
$n_{1}\geq n_{2}$}$ (S13)
$\displaystyle=\sqrt{\frac{n_{1}!}{n_{2}!}}\alpha^{n_{2}-n_{1}}L_{n_{1}}^{(n_{2}-n_{1})}(|\alpha|^{2})e^{-\frac{|\alpha|^{2}}{2}}\quad\text{if
$n_{2}>n_{1}$},$ (S14)

and

$\displaystyle\alpha=\alpha(k_{z},k_{z}+q_{z})=i\sqrt{\frac{l}{8m^{2}v(\Delta+\Delta^{\prime})}}(2k_{z}+q_{z})q_{z}=i\sqrt{\frac{\Delta}{2(\Delta+\Delta^{\prime})}\frac{\Delta
l}{v}}\frac{(2k_{z}+q_{z})q_{z}}{k_{0}^{2}}.$ (S15)

This is precisely a manifestation of the strong anisotropy of the edge states
in WSM, and we can already anticipate a highly anisotropic dynamical
polarization.

Even if the overlap function is complicated, the expressions can be simplified
significantly in the long-wavelength limit. Let us consider $q_{z}/k_{0}\to
0$, $k_{0}$ appears to be a cut-off in the long-wavelength limit. When
$q_{z}=0$, $R_{n_{1},n_{2}}$ is simplified to $\delta_{n_{1},n_{2}}$ due to
the orthogonality of the wave functions. Therefore,

$\displaystyle|F_{1,2}(\mathbf{k},\mathbf{k+q})|^{2}=\frac{1}{4}\delta_{n_{1},n_{2}}|u_{1}^{*}u_{2}+\lambda_{1}\lambda_{2}v_{1}^{*}v_{2}|^{2}$
(S16)

and one notices that only excitations $n\to n$ are allowed. Suppose that
$q_{z}\neq 0$ and for example $n_{1}\geq n_{2}$

$\displaystyle
R_{n_{1},n_{2}}(k_{z},k_{z}+q_{z})\sim\sqrt{\frac{n_{2}!}{n_{1}!}}\binom{n_{2}}{n_{1}}\alpha^{n_{1}-n_{2}}\propto
q_{z}^{n_{1}-n_{2}}.$ (S17)

The $n\to n$ excitations are therfore still the leading contribution to the
charge susceptibility, and the coupling strength of other excitations decays
as $q_{z}^{|n_{1}-n_{2}|}$.

Based on the arguments exposed above, we consider in the rest of this section
$q_{z}=0$. Accordingly, the tensor $\chi^{(0)}_{1,2}$ is diagonal and can be
treated as a scalar. We can thus divide it into contributions from different
excitation invoking different bands,

$\displaystyle\chi^{(0)}=\sum_{n\geq
1}\chi_{n}(q_{y},\omega)+\chi_{0}(q_{y},\omega).$ (S18)

For $n=0$,

$\displaystyle\chi_{0}(q_{y},\omega)=\frac{1}{V}\sum_{\mathbf{k}}\frac{f_{D}(E_{0}(\mathbf{k}))-f_{D}(E_{0}(\mathbf{k+q}))}{\omega+E_{0}(\mathbf{k})-E_{0}(\mathbf{k+q})+i\eta},$
(S19)

while for $n\geq 1$

$\displaystyle\chi_{n}(q_{y},\omega)=\sum_{\lambda_{1},\lambda_{2}}\frac{1}{V}\sum_{\mathbf{k}}\frac{f_{D}(E_{n}^{\lambda_{1}}(\mathbf{k}))-f_{D}(E_{n}^{\lambda_{2}}(\mathbf{k+q}))}{\omega+E_{n}^{\lambda_{1}}(\mathbf{k})-E_{n}^{\lambda_{2}}(\mathbf{k+q})+i\eta}\times|F_{\lambda_{1},\lambda_{2}}^{n}(\mathbf{k},\mathbf{k+q})|^{2},$
(S20)

where

$\displaystyle|F_{\lambda_{1},\lambda_{2}}^{n}(\mathbf{k},\mathbf{k+q})|^{2}=\frac{1}{2}\left(1+\lambda_{1}\lambda_{2}\frac{k_{y}(k_{y}+q_{y})+\frac{2n}{\ell_{S}^{2}}}{\sqrt{\left(k_{y}^{2}+\frac{2n}{\ell_{S}^{2}}\right)\left((k_{y}+q_{y})^{2}+\frac{2n}{\ell_{S}^{2}}\right)}}\right).$
(S21)

To obtain the plasmon mode, we need first calculate the imaginary part of each
susceptibility for $n=0$ and $n\geq 1$ from which their real part can then be
deduced by the Kramers-Kronig relations. The single-pair excitations (SPE)
resides in the domain in the $(q_{y},\omega)$-space where the imaginary part
is not vanishing. We assume $T=0$ in the following calculations and we can
consider only $\omega>0$.

###  II.1 Fermi arc $n=0$

Because of the chirality of the FA state, particle-hole excitations with
$\omega>0$ require $q_{y}>0$,

$\displaystyle\chi_{0}(q_{y}>0,\omega>0)$
$\displaystyle=\frac{k_{0}}{2\pi^{2}}\frac{q_{y}}{\omega-vq_{y}+i\eta}$ (S22)
$\displaystyle=\frac{k_{0}q_{y}}{2\pi^{2}}\left(\frac{\omega-vq_{y}}{(\omega-
vq_{y})^{2}+\eta^{2}}-i\frac{\eta}{(\omega-vq_{y})^{2}+\eta^{2}}\right).$
(S23)

Note that this expression is valid for $\omega$ and $q_{y}$ of arbitrary sign
and that $\chi_{0}$ is independent of the chemical potential due to the linear
dispersion of the FA state and a well-defined SPE $\omega=vq_{y}$. Using RPA
and $V(q)=e^{2}/2\epsilon_{0}q$ which we will justify in the following
section, we can easily find the FA plasmon

$\displaystyle\omega=vq_{y}+\frac{k_{0}e^{2}}{4\pi^{2}\epsilon_{0}}\text{sgn}(q_{y}),$
(S24)

which is gapped and disperses linearly when $q_{z}=0$.

###  II.2 Massive states $n\geq 1$

The VP states are not chiral such that the sign of $q_{y}$ is irrelevant.

####  II.2.1 Charge neutral point

When the chemical potential is at the charge neutral point, only the interband
excitations are possible. Since the VP states are gapped, the following
results are also valid for finite chemical potential if only it resides
between two band extrema of the VP states. The imaginary part of $\chi_{n}$
reads

$\displaystyle\Im[\chi_{n}^{-}](q_{y},\omega>0)=-\frac{k_{0}}{\pi
v}\frac{\frac{4n}{\ell_{S}^{2}}q_{y}^{2}}{\left(\frac{\omega^{2}}{v^{2}}-q_{y}^{2}\right)^{2}\sqrt{1-\frac{\frac{8n}{\ell_{S}^{2}}}{\frac{\omega^{2}}{v^{2}}-q_{y}^{2}}}}\Theta\left(\frac{\omega^{2}}{v^{2}}-q_{y}^{2}-\frac{8n}{\ell_{S}^{2}}\right)$
(S25)

where the minus sign on $\chi_{n}$ refers to the contribution from the
$\lambda=-$ valence band. The particle-hole spectrum
$\omega>v\sqrt{q_{y}^{2}+\frac{8n}{\ell_{S}^{2}}}$ is above the FA plasmon for
realistic value of $k_{0}$. Note that the imaginary part is proportional to
$q_{y}^{2}$ and vanishing at $q_{y}=0$ due to the vanishing overlap function
when $q_{y}=0$. This is simply because the eigenstates of Bloch Hamiltonian
for a given $k$ are orthogonal. The corresponding real part reads

$\displaystyle\Re[\chi_{n}^{-}]$
$\displaystyle=-\frac{4k_{0}\frac{2n}{\ell_{S}^{2}}q_{y}^{2}}{\pi^{2}v}\left(\frac{1}{\frac{4n}{\ell_{S}^{2}}\left(q_{y}^{2}-\frac{\omega^{2}}{v^{2}}\right)}+\frac{\pi-2\arctan
A}{\left|q_{y}^{2}-\frac{\omega^{2}}{v^{2}}\right|^{\frac{3}{2}}\left|\frac{8n}{\ell_{S}^{2}}+q_{y}^{2}-\frac{\omega^{2}}{v^{2}}\right|^{\frac{1}{2}}}\right)\quad\text{if
$q_{y}^{2}\leq\frac{\omega^{2}}{v^{2}}\leq
q_{y}^{2}+\frac{4n}{\ell_{S}^{2}}$}$ (S26)
$\displaystyle=-\frac{4k_{0}\frac{2n}{\ell_{S}^{2}}q_{y}^{2}}{\pi^{2}v}\left(\frac{1}{\frac{4n}{\ell_{S}^{2}}\left(q_{y}^{2}-\frac{\omega^{2}}{v^{2}}\right)}-\frac{\log|\frac{A+1}{A-1}|}{\left|q_{y}^{2}-\frac{\omega^{2}}{v^{2}}\right|^{\frac{3}{2}}\left|\frac{8n}{\ell_{S}^{2}}+q_{y}^{2}-\frac{\omega^{2}}{v^{2}}\right|^{\frac{1}{2}}}\right)\quad\text{otherwise}$
(S27)

where

$\displaystyle
A=\sqrt{\left|\frac{\frac{8n}{\ell_{S}^{2}}+q_{y}^{2}-\frac{\omega^{2}}{v^{2}}}{q_{y}^{2}-\frac{\omega^{2}}{v^{2}}}\right|}.$
(S28)

Using RPA and $V(q)=e^{2}/2\epsilon_{0}q$, we can prove the existence of a
plasmon mode damped by the interband SPE. Let’s check two limits. When
$\omega^{2}/v^{2}-q_{y}^{2}\gg 8n/\ell_{S}^{2}$,

$\displaystyle 1-V(q_{y})\Re[\chi_{n}^{-}](\omega,q_{y})\to 1^{-}.$ (S29)

When $\omega^{2}/v^{2}-q_{y}^{2}\to^{+}8n/\ell_{S}^{2}$,

$\displaystyle 1-V(q_{y})\Re[\chi_{n}^{-}](\omega,q_{y})$
$\displaystyle=1-\frac{e^{2}}{2\pi^{2}\epsilon_{0}v}\frac{k_{0}|q_{y}|}{\frac{2n}{\ell_{S}^{2}}},$
(S30)

which can only be satisfied for non-zero values of $|q_{y}|$. This plasmon
emerges only at finite $|q_{y}|$ for the same reason as the vanishing
imaginary part.

####  II.2.2 Band-crossing chemical potential

Now we consider the case of a chemical potential crossing one of the VP bands.
Since VP states are particle-hole symmetric, we consider a chemical potential
in the conduction band ($n=+1$),

$\displaystyle\mu=v\sqrt{k_{F}^{2}+\frac{2n}{\ell_{S}^{2}}}>0,$ (S31)

where we define a Fermi wavevector $k_{F}$. Now intraband excitations in the
$\lambda=+$ conduction band are possible and we have to suppress a part of
interband contributions in the susceptibility due to Pauli blocking. Since the
latter does not change the boundary of the particle-hole spectrum (checked by
numerical calculations), we give here only the expressions for the intraband
excitations. The imaginary part for intraband excitations reads

$\displaystyle\Im[\chi_{n,\text{intra}}^{+}](q_{y},\omega>0)$
$\displaystyle=-\frac{k_{0}}{4\pi
v}\frac{q_{y}^{2}\frac{8n}{\ell_{S}^{2}}}{\left|q_{y}^{2}-\frac{\omega^{2}}{v^{2}}\right|^{\frac{3}{2}}\left|\frac{8n}{\ell_{S}^{2}}+q_{y}^{2}-\frac{\omega^{2}}{v^{2}}\right|^{\frac{1}{2}}}$
(S32)
$\displaystyle\times\left(\Theta\left(2k_{F}-\left|q_{y}-\frac{\omega}{v}\sqrt{1+\frac{\frac{8n}{\ell_{S}^{2}}}{q_{y}^{2}-\frac{\omega^{2}}{v^{2}}}}\right|\right)-\Theta\left(2k_{F}-\left|q_{y}+\frac{\omega}{v}\sqrt{1+\frac{\frac{8n}{\ell_{S}^{2}}}{q_{y}^{2}-\frac{\omega^{2}}{v^{2}}}}\right|\right)\right).$
(S33)

For the intraband excitations in the $n=+1$ VP band, we find a particle-hole
region between two parallel parabolic-like boundaries excluding a dome from 0
to $2k_{F}$, as it should be for the Lindhard function for a one-dimensional
parabolic band, as shown in Figure S1 where we gather also the other
excitations (intra-FA and inter-VP band excitations from $-1$ to $+1$). This
is due to the one-dimensional dispersion of the massive VP states. We can
anticipate that when the chemical potential is high the Dirac linear
dispersion of the VP states is probed. So the two parallel parabolic-like
boundaries becomes linear and parallel to the line defined by $\omega=vq_{y}$.

Figure S1: Particle-hole region

The Kramers-Konig relations yield a real part of the intraband contributions
that reads

$\displaystyle\Re[\chi_{n,\text{intra}}^{+}]$
$\displaystyle=\frac{k_{0}\frac{2n}{\ell_{S}^{2}}q_{y}^{2}}{\pi^{2}v}\left(g_{1}\left(q_{y}^{2}-\left|\sqrt{k_{F}^{2}+\frac{2n}{\ell_{S}^{2}}}-\sqrt{|k_{F}-|q_{y}||^{2}+\frac{2n}{\ell_{S}^{2}}}\right|\right)-g_{1}\left(q_{y}^{2}-\left|\sqrt{k_{F}^{2}+\frac{2n}{\ell_{S}^{2}}}-\sqrt{|k_{F}+|q_{y}||^{2}+\frac{2n}{\ell_{S}^{2}}}\right|\right)\right)$
if $q_{y}^{2}\leq\frac{\omega^{2}}{v^{2}}\leq
q_{y}^{2}+\frac{4n}{\ell_{S}^{2}}$ (S34)
$\displaystyle=\frac{k_{0}\frac{2n}{\ell_{S}^{2}}q_{y}^{2}}{\pi^{2}v}\left(g_{2}\left(q_{y}^{2}-\left|\sqrt{k_{F}^{2}+\frac{2n}{\ell_{S}^{2}}}-\sqrt{|k_{F}-|q_{y}||^{2}+\frac{2n}{\ell_{S}^{2}}}\right|\right)-g_{2}\left(q_{y}^{2}-\left|\sqrt{k_{F}^{2}+\frac{2n}{\ell_{S}^{2}}}-\sqrt{|k_{F}+|q_{y}||^{2}+\frac{2n}{\ell_{S}^{2}}}\right|\right)\right)$
otherwise (S35)

where

$\displaystyle g_{1}(x)$
$\displaystyle=\frac{1}{\frac{4n}{\ell_{S}^{2}}\left(q_{y}^{2}-\frac{\omega^{2}}{v^{2}}\right)}\sqrt{\frac{x+\frac{8n}{\ell_{S}^{2}}}{x}}-\frac{2\arctan
f(x)}{\left|q_{y}^{2}-\frac{\omega^{2}}{v^{2}}\right|^{\frac{3}{2}}\left|\frac{8n}{\ell_{S}^{2}}+q_{y}^{2}-\frac{\omega^{2}}{v^{2}}\right|^{\frac{1}{2}}}$
(S36) $\displaystyle g_{2}(x)$
$\displaystyle=\frac{1}{\frac{4n}{\ell_{S}^{2}}\left(q_{y}^{2}-\frac{\omega^{2}}{v^{2}}\right)}\sqrt{\frac{x+\frac{8n}{\ell_{S}^{2}}}{x}}-\frac{\log|\frac{f(x)+1}{f(x)-1}|}{\left|q_{y}^{2}-\frac{\omega^{2}}{v^{2}}\right|^{\frac{3}{2}}\left|\frac{8n}{\ell_{S}^{2}}+q_{y}^{2}-\frac{\omega^{2}}{v^{2}}\right|^{\frac{1}{2}}}$
(S37) $\displaystyle f(x)$
$\displaystyle=\sqrt{\left|\frac{\frac{8n}{\ell_{S}^{2}}+q_{y}^{2}-\frac{\omega^{2}}{v^{2}}}{q_{y}^{2}-\frac{\omega^{2}}{v^{2}}}\right|\cdot\frac{x+\frac{8n}{\ell_{S}^{2}}}{x}}.$
(S38)

In the long wavelength limit ($q_{y}\to 0$),

$\displaystyle\Re[\chi_{n,\text{intra}}^{+}]=\frac{k_{0}}{\pi^{2}v}\frac{k_{F}}{\sqrt{k_{F}^{2}+\frac{2n}{\ell_{S}^{2}}}}\frac{v^{2}q_{y}^{2}}{\omega^{2}-v^{2}q_{y}^{2}}.$
(S39)

Recall that for a spinless 1D parabolic band, these expressions may be
simplified, and one obtains

$\displaystyle\Re[\chi_{\text{para}}]=\frac{k_{F}q_{y}^{2}}{\pi m\omega^{2}}.$
(S40)

At the band minimum where $k_{F}\to 0$, $\Re[\chi_{n,\text{intra}}^{+}]$
becomes indeed

$\displaystyle\Re[\chi_{n,\text{intra}}^{+}]=\frac{2k_{0}}{2\pi}\frac{k_{F}q_{y}^{2}}{\pi
m_{D}\omega^{2}},$ (S41)

where $m_{D}$ is the Dirac mass $\sqrt{2}/v\ell_{S}$. The strong anisotropy is
explicitly unveiled in this form: dispersed in the transverse direction and
flat in the longitudinal direction, effectively 1D band. However, the
prefactor $2k_{0}/2\pi$ reminds us that the underlying manifold is 2D. Indeed,
$2k_{0}$ is the extension of the VP bands in the $k_{z}$-direction delimited
by the two Weyl nodes at $k_{z}=\pm k_{0}$, and $2k_{0}/2\pi$ is thus its
contribution to the density of states. It is also interesting to calculate the
dynamical polarization for a 1D linear band $E=\pm vk_{y}$ for a positive
chemical potential and the excitation momentum $|q_{y}|<k_{F}$, that is

$\displaystyle\Re[\chi_{\text{lin}}]=\frac{k_{0}}{\pi^{2}v}\frac{v^{2}q_{y}^{2}}{\omega^{2}-v^{2}q_{y}^{2}},$
(S42)

which coincides with $\Re[\chi_{n,\text{intra}}^{+}]$ when $k_{F}\gg
2n/\ell_{S}^{2}$.

Using RPA and $V(q)=e^{2}/2\epsilon_{0}q$, a square-root intraband plasmon is
obtained:

$\displaystyle\omega=\sqrt{\frac{e^{2}vk_{0}}{2\pi^{2}\epsilon_{0}}\frac{k_{F}}{\sqrt{k_{F}^{2}+\frac{2n}{\ell_{S}^{2}}}}q_{y}+v^{2}q_{y}^{2}}.$
(S43)

The intraband plasmon is pushed up with increasing chemical potential and
$k_{F}$.

##  III S3. Quasi-two-dimensional random phase approximation

In this section we show how to use the random-phase approximation (RPA) in our
quasi-2D Hamiltonian. We show $\ell_{S}$ as another length scale cut-off to
consider our Hamiltonian as a quasi-2D system. We have already shown the
complicated overlap function when $q_{z}\neq 0$. So the matrix element of the
Coulomb interaction must be at least as messy as the overlap function. We
first show the matrix element of the Coulomb interaction in the most general
form without any approximation. Then, we explicitly show how the long-
wavelength approximation simplifies the calculations considerably.

###  III.1 Coulomb interaction matrix element

Consider the wavefunction of the eigenstates

$\displaystyle\langle\mathbf{r}|\psi_{n,\lambda}\rangle$
$\displaystyle=\frac{1}{\sqrt{S}}e^{i\mathbf{k_{\parallel}}\cdot\mathbf{r_{\parallel}}}\chi_{n,\lambda}(x,k_{\parallel})$
(S44)
$\displaystyle=\frac{1}{\sqrt{2S}}e^{i\mathbf{k_{\parallel}}\cdot\mathbf{r_{\parallel}}}\begin{pmatrix}u_{n,\lambda}(k_{y})\langle
x|n-1\rangle(k_{z})\\\ \lambda v_{n,\lambda}(k_{y})\langle
x|n\rangle(k_{z})\end{pmatrix},$ (S45)

where $S$ is the size of the surface and in terms of which the Coulomb
interaction operator reads

$\displaystyle\hat{V}=\frac{1}{2}\frac{e^{2}}{4\pi\epsilon_{0}S}\sum_{1,2,3,4}\sum_{k_{1},k_{2},q}\iint
dxdx^{\prime}\frac{2\pi
e^{-|q(x-x^{\prime})|}}{|q|}\chi_{1}^{\dagger}(x,k_{1})\chi_{3}(x,k_{1}+q)\chi_{2}^{\dagger}(x^{\prime},k_{2})\chi_{4}(x^{\prime},k_{2}-q)\hat{c}^{\dagger}_{1,k_{1}}\hat{c}^{\dagger}_{2,k_{2}}\hat{c}_{4,k_{2}-q}\hat{c}_{2,k_{1}+q}.$
(S46)

For smooth interfaces with $\ell\gg\ell_{S}$, one can show that the dominant
contribution to the Coulomb operator comes from the region where $q(\langle
x\rangle-\langle x^{\prime}\rangle)\ll 1$, such that the effective Coulomb
potential reads

$\displaystyle\hat{V}\approx\frac{1}{2}\frac{e^{2}}{2\epsilon_{0}S|q|}\sum_{1,2,3,4}\sum_{k_{1},k_{2},q}F_{1,3}(k_{1},k_{1}+q)F_{2,4}(k_{2},k_{2}-q)\hat{c}^{\dagger}_{1,k_{1}}\hat{c}^{\dagger}_{2,k_{2}}\hat{c}_{4,k_{2}-q}\hat{c}_{2,k_{1}+q}.$
(S47)

The complicated form of the overlap function renders further simplification
impossible unless we suppose $q_{z}=0$. Then we can define

$\displaystyle\hat{V}\approx\frac{1}{2}\sum_{1,2,3,4}\sum_{q}V_{13,24}(q)\frac{1}{S}\sum_{k_{1},k_{2}}F^{n_{1}}_{\lambda_{1},\lambda_{3}}(k_{1},k_{1}+q)F^{n_{2}}_{\lambda_{2},\lambda_{4}}(k_{2},k_{2}-q)\hat{c}^{\dagger}_{1,k_{1}}\hat{c}^{\dagger}_{2,k_{2}}\hat{c}_{4,k_{2}-q}\hat{c}_{2,k_{1}+q},$
(S48)

where the diagonal two-dimensional Coulomb interaction matrix element is
defined as

$\displaystyle
V_{13,24}(q)=\frac{e^{2}}{2\epsilon_{0}|q|}\delta_{n_{1},n_{3}}\delta_{n_{2},n_{4}}.$
(S49)

Thus the interaction tensor becomes a scalar, and the RPA polarizability reads

$\displaystyle\chi(q_{y},\omega)=\frac{\chi^{(0)}(q_{y},\omega)}{1-V_{\text{2D}}(|q_{y}|)\chi^{(0)}(q_{y},\omega)},$
(S50)

where

$\displaystyle\chi^{(0)}$ $\displaystyle=\sum_{n\geq
1}\chi_{n}(q_{y},\omega)+\chi_{0}(q_{y},\omega)$ (S51)
$\displaystyle\text{and}\quad V_{\text{2D}}(|q_{y}|)$
$\displaystyle=\frac{e^{2}}{2\epsilon_{0}|q_{y}|}$ (S52)

When $q_{z}\neq 0$, we can in general not simplify the Coulomb interaction
operator to the previous form. Luckily, when we consider only a few bands, the
long-wavelength limit renders the expressions tractable.

###  III.2 Three-band model

Let us consider a three-band model with the chiral FA state and the two $n=\pm
1$ VP states. As simple as it is, this toy model gives already many promising
results on plasmon. As shown before, it is difficult to deal with overlap
functions when $q_{z}\neq 0$ because of $R_{n_{1},n_{2}}\neq 0$ for $n_{1}\neq
n_{2}$. In the long wavelength limit, i.e. $q/k_{0}\ll 1$, we can write them
explicitly for $n_{1},n_{2}\in\\{0,1\\}$ as

$\displaystyle R_{0,0}(k_{z},k_{z}+q_{z})$
$\displaystyle=R_{1,1}(k_{z},k_{z}+q_{z})=1+O(q^{2})\approx 1$ (S53)
$\displaystyle R_{0,1}(k_{z},k_{z}+q_{z})$
$\displaystyle=R_{1,0}(k_{z},k_{z}+q_{z})=\alpha(k_{z},k_{z}+q_{z})+O(q^{3})\approx|\alpha|,$
(S54)

where $\alpha\sim q$. Then we can rewrite overlap functions as

$\displaystyle|F_{1,2}(\mathbf{k},\mathbf{k+q})|^{2}=\frac{1}{4}|u_{1}^{*}u_{2}+\lambda_{1}\lambda_{2}v_{1}^{*}v_{2}|^{2}|R_{n_{1},n_{2}}(k_{z},k_{z}+q_{z})|^{2}=|G_{1,2}(k_{y},k_{y}+q_{y})|^{2}|R_{n_{1},n_{2}}(k_{z},k_{z}+q_{z})|^{2},$
(S55)

where we decouple the $q_{y}$ and $q_{z}$ directions. On the RPA level, the
final expression for $\chi$ reads

$\displaystyle\chi(q,\omega)=\frac{\chi^{(0)}(q,\omega)}{1-V_{\text{2D}}(|q|)\chi^{(0)}(q,\omega)},$
(S56)

where

$\displaystyle\chi^{(0)}$ $\displaystyle=\sum_{1,2}\chi^{0}_{1,2}(q,\omega),$
(S57) $\displaystyle\chi^{(0)}_{1,2}(q,\omega)$
$\displaystyle=\sum_{1,2}\frac{1}{V}\sum_{\mathbf{k}}\frac{f_{D}(E_{1}(k_{y}))-f_{D}(E_{2}(k_{y}+q_{y}))}{\omega+E_{1}(k_{y})-E_{2}(k_{y}+q_{y})+i\eta}\times|G_{1,2}(k_{y},k_{y}+q_{y})|^{2}|R_{n_{1},n_{2}}(k_{z},k_{z}+q_{z})|^{2},$
(S58) $\displaystyle\text{and}\quad V_{\text{2D}}(|q|)$
$\displaystyle=\frac{e^{2}}{2\epsilon_{0}|q|}.$ (S59)

Before we calculate the full dynamical polarization including all the
contributions in the three-band model, let us write the dispersion for FA
plasmon and VP intraband plasmon when $q_{z}\neq 0$. In the long-wavelength
limit, the coupling between states of different $n$, for example $n=0$ and
$n=1$, is negligible because it is at least proportional to $q_{z}^{2}$. Thus,
we simply need to replace $q_{y}$ by $q\cos\theta$ with
$\cos\theta=q_{y}/\sqrt{q_{y}^{2}+q_{z}^{2}}$ in the previous expressions of
dynamical polarization. For the FA state,

$\displaystyle\Re[\chi_{0}]$
$\displaystyle=\frac{k_{0}}{2\pi^{2}}\frac{q\cos\theta}{\omega-vq\cos\theta},$
(S60)

and the FA plasmon dispersion reads

$\displaystyle\omega=vq\cos\theta+\frac{k_{0}e^{2}}{4\pi^{2}\epsilon_{0}}\cos\theta.$
(S61)

For the $n=1$ VP state,

$\displaystyle\Re[\chi_{n,\text{intra}}^{+}]=\frac{k_{0}}{\pi^{2}v}\frac{k_{F}}{\sqrt{k_{F}^{2}+\frac{2n}{\ell_{S}^{2}}}}\frac{v^{2}q^{2}\cos^{2}\theta}{\omega^{2}},$
(S62)

and the VP intraband plasmon dispersion reads

$\displaystyle\omega=\cos\theta\sqrt{\frac{e^{2}vk_{0}}{2\pi^{2}\epsilon_{0}}\frac{k_{F}}{\sqrt{k_{F}^{2}+\frac{2n}{\ell_{S}^{2}}}}q_{y}}.$
(S63)

As shown in the main text of the paper, this expression cannot explain the
linear dispersion of the VP intraband plasmon, and we need to take into
account the influence of the FA state.

In the long-wavelength limit, assuming $q_{z}=0$, we have

$\displaystyle\Re[\chi_{\text{tot}}]=\frac{k_{0}}{\pi^{2}v}\left(\frac{k_{F}}{\sqrt{k_{F}^{2}+\frac{2n}{\ell_{S}^{2}}}}\cdot\frac{q_{y}^{2}}{\frac{\omega^{2}}{v^{2}}-q_{y}^{2}}+\frac{1}{2}\cdot\frac{q_{y}}{\frac{\omega}{v}-q_{y}}\right).$
(S64)

The first term in the parenthesis comes from the VP intraband excitations and
the second one from the intra-FA excitations. Here we do not include
$\Re[\chi_{n}^{-}]$ since it describes an insulator playing the role of
substrate, which only change $\epsilon_{0}$ to $\epsilon_{0}\epsilon_{r}$. For
now, we do not take into the interband part in $\Re[\chi_{n}^{+}]$ either. As
we will show later, this is essential for VP intraband plasmon. Using RPA and
$V(q_{y})=e^{2}/2\epsilon_{0}|q_{y}|$, the dielectric function reads

$\displaystyle\epsilon=1-\frac{e^{2}k_{0}}{2\epsilon_{0}\pi^{2}v}\left(\frac{k_{F}}{\sqrt{k_{F}^{2}+\frac{2n}{\ell_{S}^{2}}}}\cdot\frac{|q_{y}|}{\frac{\omega^{2}}{v^{2}}-q_{y}^{2}}+\frac{1}{2}\cdot\frac{\text{sgn}(q_{y})}{\frac{\omega}{v}-q_{y}}\right).$
(S65)

To find the plasmon mode, we need to solve the following equation

$\displaystyle\omega^{2}-\text{sgn}(q_{y})\gamma\omega-\gamma
v|q_{y}|\left(1+\frac{2k_{F}}{\sqrt{k_{F}^{2}+\frac{2n}{\ell_{S}^{2}}}}\right)-v^{2}|q_{y}|^{2}=0$
(S66)

where $\gamma=e^{2}k_{0}/4\epsilon_{0}\pi^{2}$ is the gap of FA plasmon. This
equation is readily solved, and one finds

$\displaystyle\omega=\frac{\text{sgn}(q_{y})\gamma\pm\sqrt{\gamma^{2}+4\gamma
v|q_{y}|\left(1+\frac{2k_{F}}{\sqrt{k_{F}^{2}+\frac{2n}{\ell_{S}^{2}}}}\right)+4v^{2}|q_{y}|^{2}}}{2}.$
(S67)

As $q_{y}$ is small, we have

$\displaystyle\omega=\text{sgn}(q_{y})\gamma+\left(1+\frac{2k_{F}}{\sqrt{k_{F}^{2}+\frac{2n}{\ell_{S}^{2}}}}\right)vq_{y}\quad\text{or}\quad-\left(1+\frac{2k_{F}}{\sqrt{k_{F}^{2}+\frac{2n}{\ell_{S}^{2}}}}\right)vq_{y}.$
(S68)

Most saliently, the VP intraband plasmon has linear behavior compared to a
square-root-like dispersion expected naively in our previous calculations,
where we consider only the VP intraband excitations without taking into
account the FA plasmon. For $q_{z}\neq 0$, we can find the plasmon

$\displaystyle\omega=\frac{\gamma\cos\theta\pm\gamma|\cos\theta|}{2}\pm
vq\left(|\cos\theta|+\frac{2k_{F}}{\sqrt{k_{F}^{2}+\frac{2n}{\ell_{S}^{2}}}}\right)$
(S69)

where we omit the coupling between FA and VP band which is only a correction
of $q_{z}^{2}$.

This result can explain the enhanced slope in the dispersion of the FA plasmon
at small $|q_{y}|$ and the VP intraband plasmon when $q_{y}<0$. However, we do
not find an intraband plasmon for $q_{y}>0$, which we observe in our numerical
calculations. Indeed, as we have anticipated, the interband part in
$\Re[\chi_{n,inter}^{+}]$ accounts for this plasmon,

$\displaystyle\Re[\chi_{n,inter}^{+}]$
$\displaystyle=\frac{1}{V}\sum_{\mathbf{k}}\frac{f_{D}(E_{n}^{+}(\mathbf{k}))}{\omega+E_{n}^{+}(\mathbf{k})-E_{n}^{-}(\mathbf{k+q})+i\eta}\times|F_{+,-}^{n}(\mathbf{k},\mathbf{k+q})|^{2}$
(S70)
$\displaystyle+\frac{1}{V}\sum_{\mathbf{k}}\frac{-f_{D}(E_{n}^{+}(\mathbf{k+q}))}{\omega+E_{n}^{-}(\mathbf{k})-E_{n}^{+}(\mathbf{k+q})+i\eta}\times|F_{-,+}^{n}(\mathbf{k},\mathbf{k+q})|^{2}.$
(S71)

To simplify the discussion, let us assume that $k_{F}\gg\sqrt{2n}/\ell_{S}$ so
that the energy dispersion is $\omega=\pm v|q_{y}|$. Also,
$k_{F}/\sqrt{k_{F}^{2}+2n/\ell_{S}^{2}}\approx 1$. Neglecting the overlap
function which is only a correction of $q_{y}^{2}$, we have

$\displaystyle\Re[\chi_{n,inter}^{+}]$
$\displaystyle=\frac{k_{0}}{2\pi^{2}v}\int_{-k_{F}}^{k_{F}}dk_{y}\left(\frac{1}{\frac{\omega}{v}+|k_{y}|+|k_{y}+q_{y}|}-\frac{1}{\frac{\omega}{v}-|k_{y}|-|k_{y}-q_{y}|}\right)$
(S72)
$\displaystyle=\frac{k_{0}}{2\pi^{2}v}\left(\log\left(\frac{4v^{2}k_{F}^{2}}{v^{2}q_{y}^{2}-\omega^{2}}\right)-\frac{2q_{y}^{2}}{\frac{\omega^{2}}{v^{2}}-q_{y}^{2}}\right).$
(S73)

where $q_{y}>0$. Then,

$\displaystyle\Re[\chi_{n}^{+}]=\Re[\chi_{n,inter}^{+}]+\Re[\chi_{n,intra}^{+}]=\frac{k_{0}}{2\pi^{2}v}\log\left(\frac{4v^{2}k_{F}^{2}}{v^{2}q_{y}^{2}-\omega^{2}}\right)$
(S74)

where we set $0<\omega<vq_{y}$ as we have observed numerically. Therefore the
RPA dielectric function reads

$\displaystyle\epsilon=1-\frac{e^{2}k_{0}}{4\epsilon_{0}\pi^{2}v}\left(\frac{1}{|q_{y}|}\log\left(\frac{4v^{2}k_{F}^{2}}{v^{2}q_{y}^{2}-\omega^{2}}\right)+\frac{\text{sgn}(q_{y})}{\frac{\omega}{v}-q_{y}}\right).$
(S75)

This is a transcendental equation that yields a plasmon but cannot be solved
analytically. If we omit the second term in the parenthesis, we find

$\displaystyle\omega=v\sqrt{\left(1-\frac{2v^{2}k_{F}^{2}}{\gamma^{2}}\right)q^{2}+4k_{F}^{2}\frac{v}{\gamma}q-4k_{F}^{2}},$
(S76)

where we have a velocity smaller than $v$. Although this formula cannot
quantitatively match our numerical calculations shown in the main text, it
demonstrates qualitatively the importance of $\Re[\chi_{n,inter}^{+}]$.


# Finding Needles in Haystack: Formal Generative Models for Efficient Massive
Parallel Simulations
††thanks: This work is part of the project “KImaDiZ”, supported by the German
Aerospace Center (DLR) with funds of the German Federal Ministry of Economics
and Technology (BMWi), support code 50 RA 1934.

Osama Maqbool  Institute of Man-Machine-Interaction
RWTH Aachen University
Aachen, Germany
<EMAIL_ADDRESS>Jürgen Roßmann  Institute for Man-Machine
Interaction
RWTH Aachen University
Aachen, Germany
<EMAIL_ADDRESS>
###### Abstract

The increase in complexity of autonomous systems is accompanied by a need of
data-driven development and validation strategies. Advances in computer
graphics and cloud clusters have opened the way to massive parallel high
fidelity simulations to qualitatively address the large number of operational
scenarios. However, exploration of all possible scenarios is still
prohibitively expensive and outcomes of scenarios are generally unknown
apriori. To this end, the authors propose a method based on bayesian
optimization to efficiently learn generative models on scenarios that would
deliver desired outcomes (e.g. collisions) with high probability. The
methodology is integrated in an end-to-end framework, which uses the
OpenSCENARIO standard to describe scenarios, and deploys highly configurable
digital twins of the scenario participants on a Virtual Test Bed cluster.

###### Index Terms:

massive parallel simulations, bayesian optimization, virtual test beds,
experimentable digital twins

##  I Introduction

The advent of intelligent vehicles has brought with it increasing levels of
system complexity. Vehicles with machine learning components additionally
incorporate black-boxes and uncertainty within the system, transforming an
already difficult problem to a non-deterministic one. These systems therefore
require data-driven strategies in addition to classical approaches to deliver
statistical metrics on safety and reliability. Simulations offer a natural
supplement to real-world tests, allowing reproduction of expensive or
dangerous scenarios virtually and being scaled as needed. Advancements in
computer graphics has made high fidelity simulations possible, which are
especially beneficial for generating large volumes of realistic sensor data
required for machine learning based perception systems. More recently, the
availability of cloud computing resources, e.g. Microsoft Azure, has offset
the procurement effort of on-premise compute clusters enabling large-scale
parallel simulations for a wider community.

Figure 1: Visualizations of cut-in-from-left simulation. Transparent “ghosts”
illustrate parallel variants at the same time-point.

Nevertheless, the space of all possible operational scenarios remains
prohibitively large, and usage of on-premise or cloud-based compute clusters
is additionally a cost-incurring process. This warrants a motivation to only
simulate scenarios that would be meaningful to the actual development process.
For instance, the training process of a machine learning based accident-
prevention system would require simulations either violating a safety metric
or close to the point of violation. Assuming knowledge of the system, this can
be achieved by constraints on the type of scenarios, i.e. inputs to the
simulation [1] [2]. Intelligent vehicular systems, however are black-box
systems and outcomes of scenarios are rarely known apriori. Although extensive
literature exists on adversarial validation of systems [3], i.e. iterative
search within the scenario space to yield failure outcomes, these are
generally falsification approaches and do not suffice to generate the large
quantity of scenarios required for comprehensive coverage of system behaviors.

This paper proposes an efficient and flexible methodology to learn generative
models over the scenario space with respect to given outcomes. Specifying
outcome metrics as cost functions on simulation traces, we employ bayesian
optimization [4] to learn a surrogate function as belief over the cost
function with reasonable confidence around minima. The surrogate is then used
for fitting generative models. The methodology is formalized in an end-to-end
framework which uses the OpenSCENARIO standard [5] to specify scenarios and
brings them to life by constructing digital twins of scenario participants in
a virtual test bed (VTB) cluster. The setup is tested on a highway scenario,
some examples of which are illustrated in Fig. 1.

The rest of the paper is structured as follows: Sec. II provides a survey of
related work followed by a formal problem definition in Sec. III and
introduction to theoretical methods in Sec. IV. Sec. V explains the full
architecture that employs the proposed methodology, and Sec. VI applies the
methodology to an example application, finally followed by the conclusion.

##  II Related Work

The scenario-based architecture used by the authors is based on the PEGASUS
project [6], which aimed towards standardized safety qualification of
autonomous vehicles and has spawned various works for formal scenario
description [7][8] and derivations of critical scenarios via data- [9] or
knowledge-driven [10] methods. Scenic [1] provides another formal description
language for scenarios, which is similarly used for formal safety
qualifications [11].

Stepping outside the formalized scenario domain, many approaches tackle the
validation problem of systems with machine learning (ML) components via
falsification - typically optimization-based search policies to find counter-
examples for systems. Among them are image generation methods, such as domain
randomization in [12] and [13], or similarly executed point clouds generation
approaches as in [14]. Such approaches attack the object-detection or planning
ML directly, but the exploited parameters lack semantic information about the
scenario and do not offer a holistic interpretation with respect to the
system. To this end, VERIFAI [15] introduced a comprehensive validation tool
to find falsifying examples via semantic scenario parameters that is usable
with external simulators. [16] proposed a compositional falsification
framework that isolates feature spaces for the non-ML and ML sub-systems to
generate counter-examples separately. [17] uses signal temporal logic for
formalized safety specification and automated falsification of ML-based
systems with graph-based bayesian optimization.

The above approaches are effective for efficient falsification but do not
offer generative models for massive parallel simulations. Within this area,
data-driven approaches are typically employed [18], which use traffic
databases to directly sample [19] or learn targeted probability distribution
models [20]. [20] explicitly allows the incorporation of desired outcomes that
dictate the created distributions. [21] uses a model prior in addition to
traffic databases to synthesize both knowledge and data in subsequent models.
Data-driven methods, while critical for bridging real and virtual domains,
have inherently limited potential of exploring novel scenarios via intelligent
search algorithms, unlike their rule-based counterparts introduced above. Our
proposed methodology integrates an efficient optimization routine in a
formalized rule-based scenario framework, which yields a generative model
rather than singular counter-examples, which can be used to generate a large
number of counter-examples and allows convenient storage and reuse of insights
achieved by the optimization run.

##  III Problem Statement

A scenario $\mathbf{S}$ is composed of participating agents $\mathbb{A}$ and
actions $\mathbb{U}$, where each action $u\in\mathbb{U}$ is defined for
certain agents $\mathbf{a}\in\mathbb{A}$ and controlling parameters
$p\in\mathbb{P}$:

$\mathbf{S}:=\\{\mathbb{A},\mathbb{U}\\},\>\>\>\>\>\>u:=f(\mathbf{a},p).$ (1)

For instance, Sec. VI specifies a cut-in scenario with two vehicles as agents
and cut-in maneuver as action. The parameters such as action trigger times,
relative initial positions and velocities control the action exection and form
the basis of scenario variation. A simulation model $M$ maps the scenario
agents to sub-models and derives initial states $\underline{s}_{0}$ and input
trajectories $U(t)$ as dictated by scenario actions and parameters, such that

$\displaystyle\underline{s}(t)$ $\displaystyle=M(\underline{s}_{0},U(t),t),$
(2) $\displaystyle\mathbf{T}$
$\displaystyle=\\{\underline{s}(t_{0}),...,\underline{s}(t_{N})\\},\>\>t_{0}\leq
t\leq t_{N}$ (3)

where $\underline{s}(t)$ are system states, and $\mathbf{T}$ is a trace
containing all state trajectories. An outcome specification is defined as a
formula $\phi$ of predicates over $\mathbf{T}$, such that ($\vDash$ indicates
is true for)

$\phi\vDash\mathbf{T}\Leftrightarrow\xi(\mathbf{T})=0.$ (4)

$\xi$ is a positive real-valued cost function over the trace. The objective is
to find a distribution $\mathbb{P}_{\theta}$ over the scenario parameter space
$\mathbb{P}$ such that the resulting simulation and its subsequent trace has a
high probability to minimize $\xi$,

$\underset{\theta}{\mathrm{argmax}}\>\>p_{\theta}(\xi=0).$ (5)

##  IV Background

This section presents a theoretical overview of the methods used in Sec. V to
tackle the objective in (5).

###  IV-A Signal Temporal Logic

Signal temporal logic (STL) [22] is used within the formal verification
community for monitoring continuous temporal signals of complex systems and is
used to express the outcome specifications in (4). A specification (formula)
is expressed as a combination of multiple predicates:

$\phi:=true\>|\>\pi\>|\>\lnot\>\phi\>|\>\phi\land\psi\>|\>\phi\>\mathcal{U}_{[a,b]}\psi,$
(6)

where $\psi$ is another formula and $\mathcal{U}_{[a,b]}$ an until operator
over the time interval $[a,b]$. A predicate $\pi$ is defined as a real
function at a single time-point of the trace. STL allows formal and intuitive
outcome specifications with complex boolean elements such as implication, if-
and-only-if, as well as timed qualifiers such as eventually ($\mathbf{F}$) or
always ($\mathbf{G}$) true:

$\displaystyle\mathbf{F}_{[a,b]}\phi$
$\displaystyle=true\>\>\mathcal{U}_{[a,b]}\>\>\phi$ (7) $\displaystyle=\exists
t^{\prime}\in[t+a,t+b]\>\>such\>that\>\>\phi\vDash t^{\prime},$ (8)
$\displaystyle\mathbf{G}_{[a,b]}\phi$
$\displaystyle=\lnot(\mathbf{F}_{[a,b]}\lnot\phi)$ (9) $\displaystyle=\forall
t^{\prime}\in[t+a,t+b]\>\>such\>that\>\>\phi\vDash t^{\prime}.$ (10)

E.g. Sec. VI uses the time-to-collision based predicate $(1-ttc>0)$ to define
a specification over the whole trace: $\mathbf{F}_{[0,T]}(1-ttc>0)$. The real
value of predicates $\pi$ affords a robustness metric $\rho$ to the
specification in addition to truth-false inferences. E.g. the robustness of
the $\mathbf{F}_{[0,T]}(1-ttc>0)$ is the maximum value of $(1-ttc)$ throughout
the simulation. Robustness for compound specifications can be calculated
automatically:

$\displaystyle\rho^{\phi\land\psi}(t)$
$\displaystyle:=\mathrm{min}(\rho^{\phi}(t),\rho^{\psi}(t)),$ (11)
$\displaystyle\rho^{\phi\>\mathcal{U}_{[a,b]}\>\psi}(t)$
$\displaystyle:=\underset{\tau\in
t+[a,b]}{\mathrm{max}}\left(\mathrm{min}\left(\rho^{\psi}(\tau),\underset{s\in[t,\tau]}{\mathrm{min}}\rho^{\phi}(s)\right)\right),$
(12)

where $\rho^{\phi}$ indicates robustness of specification $\phi$. Robustness
metric impart both qualitative and quantitative semantics to specifications
and open the door for mathematical algorithms to automatically search for
falsifying traces.

###  IV-B Bayesian Optimization

Bayesian optimization (BO) is known for its sample-efficiency and flexibility
in optimizing black-box functions. BO maintains a surrogate function as belief
over the black-box function and updates it via adaptive measurements. Typical
choice for the surrogate function is a Gaussian Process (GP), but other
methods such as random forests are also used [23].

####  IV-B1 Gaussian Processes

Following the formalism in [24], GPs model a prior function
$f(x)\sim\mathcal{GP}(m(x),k(x,x^{\prime}))$ as a joint gaussian distribution
over the continuous function values. $m(x)$ is the prior mean and
$k(x,x^{\prime})$ is a kernel that defines covariance between any two points
and encodes prior assumptions on the target function. We use the commonly used
Matern kernel which provides explicit parameters to control the smoothness of
the fitted function. The GP prior is conditioned on measurements
$y(\mathbf{x})$ to create a posterior belief, also a GP. For a zero prior mean
and non-noisy measurements (since simulations deliver ground-truth data), the
posterior GP is

$\displaystyle f^{*}|x^{*},\mathbf{X},y(\mathbf{x})$
$\displaystyle\sim\mathcal{N}(\mu(x^{*}),\Sigma(x^{*},\mathbf{x}))$ (13)
$\displaystyle\mu(x^{*})$
$\displaystyle=k(x^{*},\mathbf{x})K_{\mathbf{x}\mathbf{x}}^{-1}y(\mathbf{x}),$
(14) $\displaystyle\Sigma(x^{*},\mathbf{x})$
$\displaystyle=k(x^{*},x^{*})-k(x^{*},\mathbf{x})K_{\mathbf{x}\mathbf{x}}^{-1}k(\mathbf{x},x^{*}),$
(15)

where the covariance $\Sigma(x^{*},\mathbf{x})$ is between the new points
$x^{*}$ and previously measured $\mathbf{x}$ and $K_{\mathbf{x}\mathbf{x}}$ is
the kernel matrix whose entries $K_{\mathbf{x}\mathbf{x}}(i,j)=k(x_{i},x_{j})$
are evaluated between each of the previously measured points $\mathbf{x}$.

####  IV-B2 Acquisition Function

Given a surrogate prior GP, BO optimizes an acquisition function iteratively
to choose the next candidate $x^{*}$ to compute a surrogate posterior. Among
the variety of acquisition functions available, we use the Thompson Sampling
(TS) [25] method which chooses the candidate $x^{*}$ by

$x^{*}=\underset{x\in D}{\mathrm{argmin}}\>\>f^{*},$ (16)

where $f^{*}$ is sampled from the surrogate prior. As it is a distribution,
samples $f^{*}$ from the same prior vary from one-another, thus implicitly
imparting both exploration and exploitation to (16). The posterior calculated
with the candidate $x^{*}$ and its measurement $y(x^{*})$ is used as prior for
the next iteration, and the process continues until $f^{*}_{min}$ converges.

###  IV-C Gaussian Mixture Models

Gaussian mixture model are a weighted sum of $N$ gaussian distributions

$p(x|\theta)=\sum_{i=1}^{N}w_{i}\mathcal{N}(x;\mu_{i},\Sigma_{i}),$ (17)

commonly used as an unsupervised learning technique. The algorithm learns a
number of normal distributions based on clusters of given data. Typically, an
expectation maximization algorithm is employed to learn the maximum likelihood
estimates of hyper-parameters $\theta$ (means $\mu_{i}$, covariance
$\sigma_{i}$) [26]. We employ bayesian gaussian mixture models from the python
scikit-learn library [27] which also learns the optimal number of
distributions to fit on the data.

Figure 2: Framework for scenario-based cluster simulations

##  V Framework

This section presents the framework to solve (5), illustrated in Fig. 2. The
overall workflow can be reiterated as such: given a scenario template and its
parameter space, the desired outcomes of scenario simulations are specified as
formal specifications. An optimizer learns a distribution over the parameter
space via successive simulations so that sampling from said distribution has
high likelihood of conforming to given specifications, i.e. delivering desired
simulation outcomes. The overall formalism is roughly derived from [28] with
emphasis on the scenario variation aspect.

####  V-1 Scenarios

Scenarios specify the actors and events in a human-readable manner independent
of the simulation framework. We specify the static scenario (road layout etc.)
with ASAM OpenDRIVE [29] format and the dynamic scenario (actor behaviour)
with OpenSCENARIO [5]. Our framework additionally differentiates between the
abstract, logical and concrete scenarios. The abstract scenario, referring
both the static and dynamic scenario, is complete with respect to actors and
events but defines no values for the scenario parameters. The logical scenario
allows human designers to specify a scenario parameter space using a
formalized meta-model illustrated in Fig. 3 from [2]. Given parameter ranges
and inter-parameter relations specified per the Constraint and probability
distributions per the Distribution template, the Variation Engine use Markov-
Chain Monte-Carlo methods to efficiently sample the resulting parameter space.
This paper focuses on the case where the designer desires certain simulation
outcomes but does not know the corresponding parameter distributions or
constraints beforehand. These outcomes are specified within the logical
scenario as-are, the subsequent adaptive sampling algorithm derives the
optimal parameter distribution and saves it as an instance of the Distribution
template (see the yellow highlighted boxes in Fig. 3). The Variation Engine
can then sample the parameter space to generate concrete scenarios. Concrete
scenarios are complete in all aspects, and can be converted to simulation
models for the virtual test beds (VTB).

Figure 3: Meta-model for logical scenario based on [2]

####  V-2 Outcome Specification

An outcome specification is defined with the STL formalism introduced in Sec.
IV-A over one or more predicates - real-valued functions that can be evaluated
over a simulation trace. For each concrete scenario and its subsequent
simulation trace, STL formalism allows the evaluation of a robustness value,
whose sign indicates conformity and value indicates the extent of conformity.
We additionally use a cost metric $\xi$ that affords control over the
subsequent optimization algorithms:

$\xi(\mathbf{T})=\begin{cases}-\rho^{\phi_{\mathbf{T}}},&\text{if
}\rho^{\phi_{\mathbf{T}}}<0\\\ 0,&\text{otherwise}.\end{cases}$ (18)

The cost metric offers mere convenience and is irrelevant for the optimization
routine. (18) uses it to invert robustness for minimization and to assert that
all scenarios conforming to the specification are equally relevant. Since
scenario distributions rather than singular examples are the goal, this
encourage the adaptive sampling algorithm to fit on a diverse range of
conforming scenarios rather than the most robust scenario.

####  V-3 Adaptive Sampling

Illustrated in Fig. 2, bayesian optimization (BO) initializes with rasterized
parameter samples and their evaluated costs $(p_{0},\xi_{0})$ to fit a
Gaussian Process (GP) $f$ as a surrogate function as in (13). Since $\xi$ must
be evaluated on VTB, it is an expensive process, and the idea is not just to
find the minimum $\xi$, but to estimate an $f$ that forms a reasonable belief
on $\xi$ near minima. This is further encouraged by a batch variant of
Thompson Sampling (TS), wherein $N$ surrogate $f$ are optimized
simultaneously. In each successive iteration, BO optimizes the TS acquisition
function in (16) to chose the next best candidate parameters. The candidates
and their costs $(p_{t}^{*},\xi_{t}^{*})$ thus create a new posterior GP as in
(13).

The converged surrogate forms the basis to infer scenarios likely to yield the
desired outcome, as parameters with low cost on a sufficiently-fit surrogate
are likely to have lower costs on the VTB as well. A large number of
parameters corresponding to minimum cost range of the surrogate are therefore
fit on a gaussian mixture model in (17). The GMM is presented as
$\mathbb{P}_{\theta}$ the solution to (5), and serves as an intuitive and
lightweight description for the scenario distribution. It is saved within the
logical scenario as an instance of the Distribution template via a set of
hyper-parameter, which can also be tuned to control the variance of concrete
scenarios.

#### Example

An example for the adaptive sampling method is illustrated in Fig. 4 for a
two-dimensional Griewank function [30]. Fig. 4(a) shows that the surrogate
function regresses well on the original function around the minima even with
existence of multiple minima. This is relevant since all parameter regions
corresponding to an outcome are typically desired. The experiment was carried
out with an initial draw of 11 samples, and converged in 8 iterations, each
with a batch size of 5. The GMM was fitted to parameter values corresponding
to the surrogate cost range $[0-0.25]$. Samples from the GMM and their actual
function evaluation is illustrated in Fig. 4(b).

(a) BO surrogate estimation (b) Samples from fitted GMM

Figure 4: Adaptive sampling example on the Griewank function Figure 5: (Left)
Bayesian optimization iterations for cut-in-from-left scenario, initialized
with a rasterized grid (not shown). (Right) Samples from fitted GMM vs uniform
samples evaluated on the VTB

####  V-4 Simulation on VTB Clusters

Concrete scenarios, sampled from GMM and any other constraints, are imported
in VTB instances running on a compute cluster. VTB serves as a cross-domain
platform offering modular simulation algorithms (e.g. rigid body dynamics,
realistic sensor simulation). For each concrete scenario, a VTB instance
constructs a simulation model by mapping actors to digital twins from a
central repository, and actions to input trajectories and initial states of
the digital twins. We use the framework developed by [28] within the VEROSIM
platform [31]. Each concrete scenario and subsequent VTB instance are fully
independent and self-contained, constituting a freely scalable “embarassingly
parallel” problem. The simulation traces are logged in an SQLite data-bank
[32] accessible throughout the framework cluster.

##  VI Transfer in Application

The presented application is from a highway Traffic Jam Chauffeur (TJC)
research project, where a machine-learning based TJC module must be
comprehensively validated against its operational scenarios. We present one of
the scenarios, cut-in-from-left, wherein a test vehicle Vehicle_1 behind the
ego vehicle cuts in its lane from left, and the goal is to find a distribution
of critically dangerous scenarios that the ego vehicle may encounter.

Scenario 1 Abstract Scenario: Cut-In-From-Left

0: $dS,dV,T$

0:

1: Actor ego pos $\leftarrow 1000m$, vel $\leftarrow 16.667ms^{-1}$, lane
$\leftarrow 0$

2: Actor: vehicle_1 rel_pos $\leftarrow dS$, rel_vel $\leftarrow dV$, rel_lane
$\leftarrow-1$

2:

2: Cut-In

3: Actors $\leftarrow$ vehicle_1

4: Trigger_Time $\leftarrow T$

The abstract scenario is defined with OpenSCENARIO formalism and summarized in
Scenario 1. The initial relative position and velocity of vehicle_1, and
trigger time for the cut-in maneuver are undefined parameters, whose
corresponding ranges are described in the Logical Scenario 2, which also
states an appropriate outcome specification for critical scenarios: the time-
to-collision between two vehicles must be below $1s$ at-least once during the
simulation. The parameter distribution to achieve high conformity to the
specification is specified with a gaussian mixture model (GMM) with unknown
hyper-parameters $\theta$.

Figure 6: Vehicle tracks for 200 samples from GMM for $Spec1$. Scenario 2
Logical Scenario: Cut-In-From-Left

0:

1: $dS\leftarrow$ range $[-30,0]$

2: $dV\leftarrow$ range $[0.5,2.0]$

3: $T\leftarrow$ range $[0.5,3.0]$

3:

4: $Spec1:\mathbf{F}(1-ttc)>0$

5: Distribution : GMM $\leftarrow\theta$

The adaptive sampling methodology creates a an initial set of concrete
scenarios with 64 rasterized samples of the parameter space. For each concrete
scenario, a VTB instance delivers a time-to-collision trace calculated at each
simulation time-step via projections of the position and velocity vectors of
both vehicles. An STL parser evaluates the outcome specification $Spec1$ to
calculate the cost of each concrete scenario. The bayesian optimization (BO)
routine then iteratively suggests new samples to evaluate on the VTB while
maintaining a surrogate Guassian Process (GP). Fig. 5 (left) illustrates BO
for 6 iterations each with a batch size of 10, suggesting samples in optimal
yet diverse regions. The converged surrogate GP is used to infer parameter
candidates with cost (time-to-collision) less than 1, which are then fitted to
a GMM to derive $\theta$. Fig. 5 (center) depicts 200 samples from the learned
$\theta$ that were evaluated on the VTB. The high percentage of evaluations
conforming to $Spec1$ (zero cost) is evident as compared to 200 uniform
samples evaluated on the VTB (Fig. 5 right). The VTB evaluations on GMM
samples are also visualized interactively using methodologies from [33] in
Fig. 1 and 6. Fig. 6 visualizes the vehicle tracks for all 200 simulations,
whereas Fig. 1 shows selected examples in detail.

##  VII Conclusion

This paper proposed generative models to efficiently test relevant use-cases
for on-cluster massive parallel simulations. The authors proposed a formal
framework based on formal specifications, human-readable scenarios and modular
virtual test beds to set up the optimization and simulation tool chain, and
presented a bayesian optimization based adaptive sampling algorithm to learn
the optimal generative models. The framework was tested on an example scenario
and compared with uniform simulations, where the presented method showed
substantially better performance.

## References

  * [1] D. J. Fremont, T. Dreossi, S. Ghosh, X. Yue, A. L. Sangiovanni-Vincentelli, and S. A. Seshia, “Scenic: a language for scenario specification and scene generation,” in _Proceedings of the 40th ACM SIGPLAN Conference on Programming Language Design and Implementation_ , 2019, pp. 63–78.
  * [2] O. Maqbool and J. Roßmann, “Formal scenario-driven logical spaces for randomized synthetic data generation,” in _Proceedings of the 10th International Conference on Model-Driven Engineering and Software Development, MODELSWARD 2022, Online Streaming, February 6-8, 2022_ , L. F. Pires, S. Hammoudi, and E. Seidewitz, Eds. SCITEPRESS, 2022, pp. 203–210. [Online]. Available: https://doi.org/10.5220/0010816400003119
  * [3] S. Riedmaier, T. Ponn, D. Ludwig, B. Schick, and F. Diermeyer, “Survey on scenario-based safety assessment of automated vehicles,” _IEEE Access_ , vol. 8, pp. 87 456–87 477, 2020.
  * [4] D. R. Jones, M. Schonlau, and W. J. Welch, “Efficient global optimization of expensive black-box functions,” _Journal of Global optimization_ , vol. 13, no. 4, pp. 455–492, 1998.
  * [5] “Asam openscenario,” https://www.asam.net/standards/detail/openscenario/, accessed: 2022-12-05.
  * [6] H. Winner, K. Lemmer, T. Form, and J. Mazzega, “Pegasus—first steps for the safe introduction of automated driving,” in _Road Vehicle Automation 5_. Springer, 2019, pp. 185–195.
  * [7] B. Kramer, C. Neurohr, M. Büker, E. Böde, M. Fränzle, and W. Damm, “Identification and quantification of hazardous scenarios for automated driving,” in _International Symposium on Model-Based Safety and Assessment_. Springer, 2020, pp. 163–178.
  * [8] T. Menzel, G. Bagschik, and M. Maurer, “Scenarios for development, test and validation of automated vehicles,” in _2018 IEEE Intelligent Vehicles Symposium (IV)_. IEEE, 2018, pp. 1821–1827.
  * [9] A. Pütz, A. Zlocki, J. Bock, and L. Eckstein, “System validation of highly automated vehicles with a database of relevant traffic scenarios,” _situations_ , vol. 1, p. E5, 2017.
  * [10] G. Bagschik, T. Menzel, C. Körner, and M. Maurer, “Wissensbasierte szenariengenerierung für betriebsszenarien auf deutschen autobahnen,” in _Workshop Fahrerassistenzsysteme und automatisiertes Fahren. Bd_ , vol. 12, 2018, p. 12.
  * [11] D. J. Fremont, E. Kim, Y. V. Pant, S. A. Seshia, A. Acharya, X. Bruso, P. Wells, S. Lemke, Q. Lu, and S. Mehta, “Formal scenario-based testing of autonomous vehicles: From simulation to the real world,” in _2020 IEEE 23rd International Conference on Intelligent Transportation Systems (ITSC)_. IEEE, 2020, pp. 1–8.
  * [12] A. Prakash, S. Boochoon, M. Brophy, D. Acuna, E. Cameracci, G. State, O. Shapira, and S. Birchfield, “Structured domain randomization: Bridging the reality gap by context-aware synthetic data,” in _2019 International Conference on Robotics and Automation (ICRA)_ , 2019, pp. 7249–7255.
  * [13] R. Khirodkar, D. Yoo, and K. M. Kitani, “Vadra: Visual adversarial domain randomization and augmentation,” 2018.
  * [14] M. Abdelfattah, K. Yuan, Z. J. Wang, and R. Ward, “Towards universal physical attacks on cascaded camera-lidar 3d object detection models,” in _2021 IEEE International Conference on Image Processing (ICIP)_. IEEE, 2021, pp. 3592–3596.
  * [15] T. Dreossi, D. J. Fremont, S. Ghosh, E. Kim, H. Ravanbakhsh, M. Vazquez-Chanlatte, and S. A. Seshia, “Verifai: A toolkit for the formal design and analysis of artificial intelligence-based systems,” in _International Conference on Computer Aided Verification_. Springer, 2019, pp. 432–442.
  * [16] T. Dreossi, A. Donzé, and S. A. Seshia, “Compositional falsification of cyber-physical systems with machine learning components,” _Journal of Automated Reasoning_ , vol. 63, no. 4, pp. 1031–1053, 2019.
  * [17] S. Ghosh, F. Berkenkamp, G. Ranade, S. Qadeer, and A. Kapoor, “Verifying controllers against adversarial examples with bayesian optimization,” in _2018 IEEE International Conference on Robotics and Automation (ICRA)_. IEEE, 2018, pp. 7306–7313.
  * [18] W. Ding, C. Xu, H. Lin, B. Li, and D. Zhao, “A survey on safety-critical scenario generation from methodological perspective,” _arXiv preprint arXiv:2202.02215_ , 2022.
  * [19] N. Webb, D. Smith, C. Ludwick, T. Victor, Q. Hommes, F. Favaro, G. Ivanov, and T. Daniel, “Waymo’s safety methodologies and safety readiness determinations,” _arXiv preprint arXiv:2011.00054_ , 2020.
  * [20] M. O’Kelly, A. Sinha, H. Namkoong, R. Tedrake, and J. C. Duchi, “Scalable end-to-end autonomous vehicle testing via rare-event simulation,” _Advances in neural information processing systems_ , vol. 31, 2018.
  * [21] D. Rempe, J. Philion, L. J. Guibas, S. Fidler, and O. Litany, “Generating useful accident-prone driving scenarios via a learned traffic prior,” in _Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition_ , 2022, pp. 17 305–17 315.
  * [22] O. Maler and D. Nickovic, “Monitoring temporal properties of continuous signals,” in _Formal Techniques, Modelling and Analysis of Timed and Fault-Tolerant Systems_. Springer, 2004, pp. 152–166.
  * [23] S. Greenhill, S. Rana, S. Gupta, P. Vellanki, and S. Venkatesh, “Bayesian optimization for adaptive experimental design: A review,” _IEEE Access_ , vol. 8, pp. 13 937–13 948, 2020.
  * [24] C. K. Williams and C. E. Rasmussen, _Gaussian processes for machine learning_. MIT press Cambridge, MA, 2006, vol. 2, no. 3.
  * [25] W. R. Thompson, “On the likelihood that one unknown probability exceeds another in view of the evidence of two samples,” _Biometrika_ , vol. 25, no. 3-4, pp. 285–294, 1933.
  * [26] C. M. Bishop and N. M. Nasrabadi, _Pattern recognition and machine learning_. Springer, 2006, vol. 4, no. 4.
  * [27] F. Pedregosa, G. Varoquaux, A. Gramfort, V. Michel, B. Thirion, O. Grisel, M. Blondel, P. Prettenhofer, R. Weiss, V. Dubourg, J. Vanderplas, A. Passos, D. Cournapeau, M. Brucher, M. Perrot, and E. Duchesnay, “Scikit-learn: Machine Learning in Python ,” _Journal of Machine Learning Research_ , vol. 12, pp. 2825–2830, 2011.
  * [28] U. Dahmen, T. Osterloh, and J. Roßmann, “Generation of virtual test scenarios for training and validation of ai-based systems,” in _2021 IEEE International Conference on Progress in Informatics and Computing (PIC)_ , 2021, pp. 64–71.
  * [29] “Asam opendrive,” https://www.asam.net/standards/detail/opendrive/, accessed: 2022-12-05.
  * [30] A. O. Griewank, “Generalized descent for global optimization,” _Journal of optimization theory and applications_ , vol. 34, no. 1, pp. 11–39, 1981.
  * [31] J. Rossmann, M. Schluse, C. Schlette, and R. Waspe, “Control by 3d simulation–a new erobotics approach to control design in automation,” in _International Conference on Intelligent Robotics and Applications_. Springer, 2012, pp. 186–197.
  * [32] L. Atorf, T. Cichon, and J. Roßmann, “Flexible data logging, management, and analysis of simulation results of complex systems for erobotics applications,” in _ESM 2015 (29th European Simulation and Modelling Conference)_ , 2015.
  * [33] L. Atorf and J. Roßmann, “Interactive analysis and visualization of digital twins in high-dimensional state spaces,” in _2018 15th International Conference on Control, Automation, Robotics and Vision (ICARCV)_. IEEE, 2018, pp. 241–246.