11institutetext:  Institut für Physik, Johannes Gutenberg-Universität Mainz,
D-55099 Mainz, Deutschland
Theoretische Physik, Universität des Saarlandes, D-66123 Saarbrücken,
Deutschland

First pacs description Second pacs description Third pacs description

# Injection and nucleation of topological defects in the quench dynamics of
the Frenkel-Kontorova model

Oksana Chelpanova 11 Shane P. Kelly 11 Giovanna Morigi 22 Ferdinand Schmidt-
Kaler 11 Jamir Marino 111122

###### Abstract

Topological defects have strong impact on both elastic and inelastic
properties of materials. In this article, we investigate the possibility to
controllably inject topological defects in quantum simulators of solid state
lattice structures. We investigate the quench dynamics of a Frenkel-Kontorova
chain, which is used to model discommensurations of particles in cold atoms
and trapped ionic crystals. The interplay between an external periodic
potential and the inter-particle interaction makes lattice discommensurations,
the topological defects of the model, energetically favorable and can tune a
commensurate-incommensurate structural transition. Our key finding is that a
quench from the commensurate to incommensurate phase causes a controllable
injection of topological defects at periodic time intervals. We employ this
mechanism to generate quantum states which are a superposition of lattice
structures with and without topological defects. We conclude by presenting
concrete perspectives for the observation and control of topological defects
in trapped ion experiments.

###### pacs:

nn.mm.xx

###### pacs:

nn.mm.xx

###### pacs:

nn.mm.xx

##  1 Introduction

The study of topological defects, such as solitons, is interdisciplinary,
permeating different areas of physics such as cold atoms [1, 2, 3, 4],
spintronics [5, 6, 7], and nano-friction in material science [8, 9, 10, 11,
12, 13]. They are also at the core of a variety of technological applications
encompassing quantum computing [14, 15, 16, 17] and atomic sensors [18, 19,
20]. In material science, they arise as lattice defects in solid state
crystals and are particularly important in determining the mechanical
properties such as rigidity [21].

In this article, we investigate the possibility to controllably inject and
eject topological defects in quantum simulators of solid state lattice
structures, with the perspective of preforming similar nano-scale control of
individual defects in real materials. To pursue this goal, we consider the
quench dynamics of the Frenkel-Kontorova model [22, 23], originally introduced
to understand the macroscopic and static properties of dislocations in
crystalline materials. This model can be simulated in a variety of quantum
platforms such as ultra cold atoms [24], and trapped ion crystals [25, 26, 10,
27, 28]; unlike in materials, microscopic control over these systems is
currently possible, allowing for experiments that can perform repeatable
quenches which are sensitive to the quantum nature of the constituents
particles.

The topological defects in this model are lattice discommensurations induced
by the interplay between an external periodic potential and inter-particle
interaction. By tuning the external potential, the lattice can undergo a
structural first-order transition from a (topologically trivial) commensurate
state, in which particles stay at the minima of an underlying periodic
potential, to an incommensurate state, where particles are dislocated from
their natural equilibria, and the local discommensuration of particles can be
mapped to a kink in a Frenkel-Kontorova model.

Our main result is that these topological defects can be injected from the
boundaries of the system by performing quenches which originate from the
commensurate phase. This mechanism occurs in a controlled fashion with atomic
discommensurations entering the chain at discrete intervals of time and
building a metastable state characterized by a temporal staircase of kinks,
with life-time scaling linearly with the chain size. We also explore the
impact of non-equilibrium fluctuations at the commensurate-incommensurate
transition, and show that they can mold novel states in the quantum Frenkel-
Kontorova model. By quenching the system close to the boundary of the
commensurate-incommensurate transition, quantum fluctuations lead to
statistical nucleation of kinks, which drive the chain in a quantum
superposition of configurations with and without topological defects.

We focus throughout the paper on parameters regimes proper of trapped ion
simulators, and we discuss in the concluding sections how to observe these
phenomena in low dimensional ultracold atomic wires, arguing for broad
applicability of our results to other quantum simulators.

##  2 Model

We consider $N$ atoms forming a one-dimensional chain with inter-particle
spacing $a_{0}$. The system is subject to a periodic external lattice
potential of period $a_{s}$. The physics of the discommensuration between the
atoms spacing and the lattice periodicity can be captured by the Frenkel-
Kontorova model [29, 28, 30, 23, 22]

$H=\sum_{n}\frac{p_{n}^{2}}{2}+{m_{K}^{2}}(1+\cos\phi_{n})+\sum_{r=\pm
1}\frac{(\phi_{n+r}-\phi_{n}-2\pi r\delta)^{2}}{2}$ (1)

where $p_{n}=\dot{\phi}_{n}.$ Here the first term describes the kinetic energy
of the atoms, the second one depicts the interaction with the lattice
potential, and the last term captures the inter-particle interaction. Here
$\phi_{n}/(2\pi)=x_{n}/a_{s}-n$ denotes a classical field, characterizing the
displacement of the $n$-th particle position $x_{n}$ from the coordinate of
the $n$-th minima of the lattice potential ($n=1,\ldots,N$). The dimensionless
parameter $m_{K}^{2}$ captures the amplitude of the periodic lattice potential
and it plays the role of the kink mass in the following. The ‘misfit
parameter’ $\delta=(a_{0}-a_{s})/a_{s}$ quantifies instead the degree of
discommensuration between $a_{0}$ and $a_{s}$.

Figure 1:  (a) Equilibrium phase diagram of the model (1) as a function of the
kink mass $m_{K}$ and the misfit parameter $\delta$. (b) Displacement,
$\phi_{n}$, of the $n$-th atom with respect to the $n$-th minimum of the
potential for the incommensurate phase (orange) and the commensurate phase
(yellow). (c) Equilibrium positions of the atoms in the incommensurate phase.
Some atoms sit at the minima of the lattice potential, while others in
discommensurate equilibrium positions which can be mapped into kinks of
$\phi_{n}$. (d) The atomic positions in the commensurate phase coincide with
the minima of the trapping potential. (e) Sketch of the energy spectrum of the
model as a function of momentum $q$ in the incommensurate phase. Blue crosses
correspond to collective sonic excitations, while red ones to parabolically
dispersing modes (optical branch).

By tuning the amplitude and misfit parameter, $m_{K}$ and $\delta$, of the
periodic lattice potential, the system at equilibrium undergoes a structural
phase transition between two different configurations for $\phi_{n}$. The
equilibrium phase diagram is sketched in Fig. 1(a). For
$\delta<\delta_{c}\approx 2m_{K}/\pi^{2}$ a strong lattice potential pins the
particles to the potential minima, forming a commensurate phase
($\phi_{n}=\pi$), see Fig. 1(b, d). Instead, for values of
$\delta\geq\delta_{c}$ atoms reside in a new incommensurate configuration
associated with kinks in $\phi_{n}$ placed at the location of the particles’
displacements, cf. Fig. 1(b, c). The kinks (anti-kinks) correspond to a local
distribution of holes (excess particles) in the chain [28, 30, 29, 31, 32,
30], and are represented by a ‘step’ in the displacement field, $\phi_{n}$,
equal to $2\pi$ ($-2\pi$). The number of kinks $Q=(\phi_{N}-\phi_{1})/(2\pi)$
is a parameter that allows us to distinguish between the commensurate and the
incommensurate phases [33]. This quantity estimates the number of ‘steps’ in
$\phi_{n}$, i.e., jumps in the value of the field $\phi$ from $\pi+2\pi m$ to
$\pi+2\pi(m+1),$ where $m$ is an integer number. For instance, the
incommensurate configuration in Fig. 1(b, c) has three ‘steps’ in $\phi_{n}$
which corresponds to $Q=3.$ For large values of the misfit parameter $\delta$,
the number of kinks, $Q$, is linearly proportional to $\delta$. Accordingly,
$Q$ is directly related to the spacing between kinks, $l_{\delta}\propto
1/\delta.$

In the following we investigate the dynamical injection of kicks and use the
width of the kink and the excitation spectrum as key parameters to interpret.
The width of the kink, $l_{m}$, is equal to the distance over which $\phi$
increases by $2\pi$, and it is roughly proportional to the inverse of the kink
mass, $l_{m}\propto 1/m_{K}$ [24, 33, 28]. For large $m_{K}$ the shape of the
kink is sharp, since any broadening of the incommensurate region requires
additional energy cost due to the lattice potential. For small $m_{K}$ the
lattice potential is almost negligible and even a single kink appears wide.
Following [33, 34] the excitation spectrum of the model (1) in the two phases
can be extracted by expanding the equations of motion for $\phi_{n}$ around
the classical solution $\phi_{n}^{(0)}$ in powers of small fluctuations
$\eta_{n}$. The corresponding equations of motion for the linear
displacements, $\eta_{q}$, in momentum space read

$\omega_{q}^{2}\eta_{q}=-{m_{K}^{2}}\cos(\phi^{0})\eta_{q}+c^{2}q^{2}\eta_{q}+\ldots,$
(2)

where $c$ is the speed of sound, which is equal to one for the model (1). The
spectrum in the incommensurate phase consists of two branches, cf. Fig. 1(e).
The first $Q$ modes correspond to gapless collective excitations, phonons,
describing the propagation of kinks through the lattice at the speed of sound
$c$ (blue line). The remaining $N-Q$ modes form an optical branch (red line)
and they describe small harmonic displacements of the atoms around the
potential minima: they are parabolically dispersing excitations with a gap set
by $m_{K}$. Within the commensurate phase, all the $N$ modes are optical.

##  3 Controlled Injection of Defects

We now show how to inject, and also eject, topological defects by keeping
$m_{K}$ fixed, and making an abrupt change of $\delta$, from $\delta_{i}$ to
$\delta_{f},$ which in terms of experimental parameters corresponds to the
sudden change of the inter-particle spacing $a_{0}$. We start by solving
numerically the equations of motion that govern the dynamics of the
topological defects in model (1), which read

$\displaystyle\ddot{\phi}_{1}=m_{K}^{2}\sin\phi_{1}+\left(\phi_{2}-\phi_{1}-{2\pi}\delta\right)$
(3) $\displaystyle\ddot{\phi}_{n=2:N-2}={m_{K}^{2}}\sin\phi_{n}+$
$\displaystyle+\left(\phi_{n+1}-\phi_{n}-{2\pi}\delta\right)+\left(\phi_{n-1}-\phi_{n}+{2\pi}\delta\right)$
$\displaystyle\ddot{\phi}_{N}={m_{K}^{2}}\sin\phi_{N}+\left(\phi_{N-1}-\phi_{N}+{2\pi}\delta\right).$

Here the effect of the discommensuration between $a_{0}$ and $a_{s}$ cancels
for particles in the bulk, as it produces opposite contributions from
interaction with right and left neighbours. However, the misfit parameter
affects directly dynamics of the boundary particles, because the contribution
from discommensuration is not compensated. In this way, after the quench of
$\delta$ only the edge is initially driven out of equilibrium such that the
boundary acts as a reservoir for topological defects (see also [35, 3, 36,
37]). Furthermore, due to their topological nature, once the defects enter the
bulk, they cannot be destroyed [38, 39, 40].

In order to see this concretely, we numerically initialize a system in the
equilibrium state with misfit parameter $\delta_{i}$, and then numerically
solve its evolution according to equation (3) with the misfit parameter
$\delta_{f}$. In Fig. 2(a) we plot a diagram of the possible dynamical
responses in the system using the long-time averaged number of kinks
$\bar{Q}$, where time averages are taken over several periods $\propto
2\pi/m_{K}$. In the white region of Fig. 2(a) no kinks or anti-kinks are
injected during dynamics. The red region in the phase diagram corresponds to
kinks’ injection after quenching the misfit parameter above a dynamical
critical threshold, $\delta_{c}^{d}$. Finally, the blue region in Fig. 2(a)
corresponds to quenches which decrease the total number of kinks. This occurs
by injecting anti-kinks into the system. Injection of a single kink (an anti-
kink) from the boundary is associated with the excitation of a collective
sonic mode.

Figure 2:  (a) Dynamical response of the system after a quench of the misfit
parameter from $\delta_{i}$ to $\delta_{f}$, with $m_{K}=0.5$ (here $N=80$).
We plot the time-averaged number of kinks, $\bar{Q}-Q_{(t=0)}$. The solid
lines mark the critical thresholds, $\delta_{c}^{d}$, corresponding to the
discrete injection of a certain number of (anti-)kinks. The dashed line
corresponds to $\delta_{i}=\delta_{f}$. The red region corresponds to the
injection of kinks into the system, while the blue one pertains to anti-kinks.
(b) Dynamics of the number of kinks in the system after the quench marked with
a $\color[rgb]{0,0,1}\star$ in panel (a). Here we fix $\delta_{i}=0$ and
$\delta_{f}=0.16>\delta_{c}^{d}$ and consider different system sizes $N$
(marked with different colors). The plot illustrates periodic injection of
kinks mirroring in a ‘step-like’ increase of $Q$ in time. The time scale
$\tau\propto N/c$ denotes the saturation time of $Q$ for $N=140$ ions after
which finite-size effects take over. The maximal number of kinks in the system
after the quench scales linearly with $N$ (see inset). (c) Energy of two kinks
as a function of the relative distance between their centers. The effective
kinks repulsion takes place at distances $r<l_{\delta}+l_{m}$. The energy
scale in panel (c) is normalized by the energy of two infinitely separated
kinks $E_{\infty}$.

Fig. 2(b) shows how the number of kinks in the system increases following the
quench from the commensurate to a target point of the incommensurate phase,
marked with $\color[rgb]{0,0,1}\star$ in Fig. 2(a). Dynamics start from $Q=0$
and exhibit a staircase structure in time due to the injection of kinks from
the boundaries up to time $\tau\propto N/c$ when finite-size effects become
relevant. The duration of the various plateaus of the staircase is
proportional to $l_{\delta_{f}}/c$. Arranging kinks at $r<l_{\delta}+l_{m}$
requires an additional energy cost. In order to show that, we calculate the
energy of two kinks as a function of the relative distance between their
centers $r$, cf. Fig. 2(c). If the kinks overlap, their energy is higher than
the energy of the configuration with two infinitely separated kinks, resulting
in short-distance repulsion [41, 29]. Thus, the average distance between kinks
after the quench is fixed to approximately $l_{m}+l_{\delta_{f}}$. As a
result, the number of kinks saturates at long times to $Q\approx
2N/(l_{m}+l_{\delta_{f}})$, which is basically given by the total number of
solitons fitting in the system, for a given size $N$ (cf. inset in Fig. 2(b)).
Note, that $\propto N/(l_{m}+l_{\delta_{f}})$ kinks enter the system from each
boundary. Indeed, up to times $\tau$, systems with a different number of
particles evolve with the same staircase temporal profile, since dynamics are
universally ruled by ${\delta_{f}}$ and $m_{K}$. Notice that the height of
each ‘step’ in Fig. 2(b) is equal to two, since we simultaneously inject one
kink from the left and another one from the right boundary, after the quench.
Finally, at late times $t>\tau$, the number of kinks, $Q(t)$, oscillates
around the steady state value, see Fig. 2. This behaviour results from both
sonic and optical modes. The first one alters the value of charge as a result
of scattering of the injected kinks against the boundaries. The second one
describes individual oscillations of particles around potential minima.

##  4 Dynamics close to the boundary of the commensurate-incommensurate
transition

An exciting feature of cold atom and trapped ion platforms is their
sensitivity to quantum fluctuations. Generally, the effect of quantum
fluctuations within a given phase only quantitatively modify the dynamics, but
the dynamics close to the boundary between two phases can have non-trivial
effects such as destroying order, generating new dynamical phases, or
provoking a quantum critical region in the system [42, 43, 44, 45, 46, 47, 48,
49, 50, 51]. As such, we investigate the effect of quantum fluctuations for
quenches close to the boundary of the commensurate-incommensurate transition,
and we introduce the quantum version of Hamiltonian (1):

$\hat{H}=\sum_{n}\frac{\hat{p}_{n}^{2}}{2}+m_{K}^{2}(1+\cos\hat{\phi}_{n})+\sum_{r=\pm
1}\frac{(\hat{\phi}_{n+r}-\hat{\phi}_{n}-2\pi r\delta)^{2}}{2}$ (4)

where we have promoted $p_{n}$ and $\phi_{n}$ to the dimensionless quantum
operators $\hat{p}_{n}$ and $\hat{\phi}_{n}$ with canonical commutation
relations $[\hat{\phi}_{n},\hat{p}_{m}]=ih_{\operatorname{eff}}\delta_{n,m}$.
Here we take $h_{\operatorname{eff}}$ as a dimensionless free parameter which
controls the strength of quantum fluctuations. Below we give an expression for
$h_{\operatorname{eff}}$ for the trapped ion implementation.

In the following, we take $h_{\operatorname{eff}}\ll 1$ and investigate the
impact of quantum fluctuations using the semi-classical approximation known as
the Truncated Wigner Approximation (TWA) [52, 53, 54]. We consider a quench
from the ground state of the Hamiltonian (4), that can be approximated [52] by
a Gaussian wave function with the Wigner function equal to

$W\left(\eta_{q}^{0},n_{q}^{0}\right)=\prod_{q}\exp\left[-{\left|\eta_{q}^{0}\right|^{2}}/{\sigma_{q}}-{\left|\psi_{q}^{0}\right|^{2}}/b_{q}\right],$
(5)

where $\eta^{0}_{q}=\phi^{0}_{q}-\langle\hat{\phi}_{q}\rangle_{0}$ and
$\psi^{0}_{q}=p^{0}_{q}-\left<\hat{p}_{q}\right>_{0}$ are the phase space
variables relative to the mean value of the operators in the initial state,
while
$\sigma_{q}=\langle\hat{\eta}_{q}^{2}\rangle_{0}=h_{\operatorname{eff}}/\omega_{q}$
and
$b_{q}=\langle\hat{\psi}_{q}^{2}\rangle_{0}=h_{\operatorname{eff}}\omega_{q}$
encode the variance of quantum fluctuations in the initial state. Finally, the
dispersion relation $\omega_{q}$ is determined by Eq. (2). At zeroth order in
$h_{\operatorname{eff}}$, the TWA approximation is simply the mean field
dynamics discussed above where the system evolves by the classical equations
of motion Eq. (3) with initial state taken as the mean field ground state. At
first order in $h_{\operatorname{eff}}$, the dynamics of phase space variables
are still described by Eq. 3, but the initial state is sampled from the
positive Wigner distribution in Eq. (5). Observables at later times are
calculated by averaging over the trajectories resulting from the initial
quantum noise. Taking into account higher orders quantum corrections can be
performed by including quantum jumps during dynamics in each trajectory, with
each jump introducing $O\left(h_{\operatorname{eff}}^{2}\right)$ corrections
[52, 55].

Figure 3:  (a) Separation of trajectories after a quench close to the boundary
between two phases $\delta_{i}=0\to\delta_{f}\lesssim\delta_{c}^{d}$ . Figure
shows different particle displacements, $\phi_{n}$, following the quench at a
time $t\propto N/c$. Different configurations correspond to different initial
conditions sampled from the distribution in Eq. (5), and color distinguishes
between configurations that remain in the commensurate phase (yellow) and
those that nucleated a kink (blue). The overall number of kinks in the
averaged value of $\phi_{n}$ (solid black line), at long times, is non-
vanishing. (b) Number of kinks as a function of time after the quench
calculated classically (yellow), semi-classically (red), and taking into
account next-order quantum corrections (blue). Dashed lines correspond to
time-averaged $Q$. The number of kinks at the classical level is zero and the
system remains in the commensurate phase, while for TWA dynamics and for its
first-order correction, the time-averaged $Q$ takes a non-integer value
between zero and one. (c) Energy of the system as a function of $\delta$
calculated numerically. For quenches $\delta_{f}<\delta_{c}^{d}$ entering the
incommensurate phase after the quench requires additional energy cost $\Delta
E=E(\delta_{c}^{d})-E(\delta_{f})$, which makes this transition impossible in
the classical case. Taking into account quantum fluctuations allows system to
tunnel this energy barrier and nucleate a kink. The inset of panel (a)
compares the probability of nucleating a kink through the energy barrier as a
function of the final misfit parameter, for TWA numerics (circles) and for the
WKB approximation (solid); fitting parameters are $P_{0}=0.81$ and $\alpha=5.$

Let us now revisit the quench of the misfit parameter within the commensurate
phase, $\delta_{i}=0\to\delta_{f}\lesssim\delta_{c}^{d}.$ At the classical
level this excites the displacement field $\phi_{n}^{(0)}$ in an optical mode.
At small times the dynamics of fluctuations are governed by
$\ddot{\eta}_{q}\approx
m_{K}^{2}\cos\left(\phi^{0}(t)\right)\eta_{q}-c^{2}q^{2}\eta_{q}$ (cf. with
Eq. (2)), with a periodically driven mass term set by the motion of
$\phi_{n}^{(0)}(t)$. In cosmology [56] or cold atoms [57], a similar structure
for dynamics of the low-energy modes results in an instability that can
trigger non-linear effects. In our system the parametric amplification of low-
energy modes can indeed cause the separation of quantum trajectories on times
of order $\log{h_{\operatorname{eff}}}$. Fig. 3(a) shows final configurations
of the field for different samplings of the Wigner distribution associated to
the pre-quench state. The yellow trajectories remain within the commensurate
phase and contain no kinks, while the blue trajectories are the ones which
have tunnelled into the incommensurate region, and in fact they display kinks
in the profile of the displacement field, $\phi_{n}$. Such a semi-classical
ensemble represents a quantum state that is a superposition of an atomic
lattice in the commensurate configuration and a configuration with a kink. We
check that this result remains qualitatively unaltered after adding higher
order quantum corrections in Fig. 3(b).

Classically, a quench below the dynamical critical threshold,
$\delta_{i}=0\to\delta_{f}\lesssim\delta_{c}^{d}$, keeps the system within the
commensurate phase since the injection of a single kink would require an
additional energy cost $\Delta E=E(\delta_{c}^{d})-E(\delta_{f})$. By adding
quantum fluctuations, however, the system can tunnel through this barrier and
nucleate a kink at the boundary. The fraction of trajectories that nucleated a
kink after the quench of $\delta$ can be estimated via WKB approximation [58].
In this approximation, nucleation of a kink can be understood as tunneling of
a quantum particle through the energy barrier separating a commensurate from
an incommensurate configuration, see Fig. 3(c). The transition probability
according to WKB reads

$P=P_{0}\exp\left(-{\alpha}{h_{\operatorname{eff}}^{-1}}\int_{\delta_{f}}^{\delta_{c}^{d}}\sqrt{E(\delta_{c}^{d})-E(\delta)}d\delta\right),$
(6)

which we fit in good agreement with numerical results in the inset of Fig.
3(a). Here the effective Planck constant $h_{\operatorname{eff}}$ determines
the width of the region in $\delta_{c}^{d}-\delta_{f}$ where nucleation of
kinks can be observed, namely, $\sqrt{\Delta
E}(\delta_{c}^{d}-\delta_{f})\lesssim h_{\operatorname{eff}}\alpha^{-1}$.
Quenches away from this region do not produce states with the above quantum
superposition. Instead, quantum noise provokes dephasing among trajectories
and can result in the broadening of the width of a single kink such that atoms
will be harder to resolve.

Similar arguments hold for quenches deep into incommensurate region. As shown
in the portrait of dynamical responses (Fig. 2(a)) there is a set of critical
thresholds $\delta_{c_{n}}^{d}$ that correspond to the injection of multiple
($n=\pm 2,\pm 4,\ldots$) solitons into the system. Therefore, by quenching
$\delta$ in between $\delta_{c_{n}}^{d}$ and $\delta_{c_{n+2}}^{d}$ we can
generate states which are in a superposition of $n$ and $n+1$ kinks.

Notice that if the initial state was instead a Gibbs thermal state,
transitions into the incommensurate state could also occur by thermal
activation [59], and it would give rise qualitatively to the same results of
Fig. 3. The major difference is that the system would be in a mixed state of
crystal configurations with different values of $Q$ and with quantum coherence
decreasing with increasing temperature.

Figure 4: Number of kinks in a chain after the quench of the misfit parameter
from the commensurate to the incommensurate phase in a long-range interacting
model (cf. Eq. (7)) for different system sizes $N.$ Here we fix $m_{K}=0.5$
and $\delta_{f}=0.11.$ Inset shows the scaling of the maximal number of kinks,
injected after the quench as a function of the system size.

##  5 Trapped Ion Realization

We now discuss how the quench dynamics presented in the core of the paper can
find application in trapped ion chain [25, 26, 10, 27]. An ionic crystal is
formed by trapping ions in a 1D Paul trap whose strength sets the inter-
particle distance $a_{0}$. The external lattice potential is provided by an
optical lattice with period $a_{s}$, and amplitude, $\epsilon$, that sets the
mass $m_{K}^{2}=(2\pi)^{2}{\left(4\pi\epsilon_{0}\right)\hbar\epsilon
a_{0}^{3}}/{\left(4q^{2}a_{s}^{2}\right)}$. Here $q$ is an electron charge and
$\epsilon_{0}$ is a vacuum permittivity. The effective Plank constant takes
the form
$\hbar_{\operatorname{eff}}={\hbar(2\pi)^{2}}\sqrt{{(4\pi\epsilon_{0})a_{0}^{3}}/({2q^{2}M})}/{a_{s}^{2}},$
where $M$ is the mass of a single 40Ca ion. In a trapped ion implementation
the parameters of model (1) can be changed in the range $0<m_{K}<1$,
$-0.5\leq\delta\leq 0.5$ and $h_{\operatorname{eff}}\leq 10^{-2}.$ The
discussion above holds except for an unscreened Coulomb interaction that
creates a long range modification to the harmonic interaction between
neighboring ions. In this section we explore how the long range interaction
modifies the ability to controllably inject topological defects into the ionic
lattice by considering the dynamics under the modified Hamiltonian

${H}=\sum_{n}\frac{p_{n}^{2}}{2}+m_{K}^{2}(1+\cos\phi_{n})+\sum_{r>0}\frac{(\phi_{n+r}-\phi_{n}-2\pi
r\delta)^{2}}{2|r|^{3}}$ (7)

where $r$ is an integer and runs over all neighbours.

As the inter-particle interaction increases, the form of a single kink
modifies and it becomes less sharp on the tails [28, 60]. As a result, the
effective interaction between kinks is also modified. Furthermore, in case of
the long-range interaction all ions are sensitive to the presence of
dislocations, and thus, injection of each kink into the system contributes
immediately into dynamics of the ions in the middle of the chain. Nonetheless,
as we show below, our previous results remain qualitatively similar.

Similarly to the model with short range interactions, quenching the misfit
parameter can drive commensurate-incommensurate transition in the long-range
interacting model (7). Fig. 4 shows the number of kinks in the system $Q$ as a
function of time in this case, and again shows kinks entering chain one by one
from the boundaries in a periodic time intervals. Quantitatively, the duration
of the plateau is shortened as the form of a single kink is broadened by the
Coulomb repulsion. Furthermore, in the short-range interacting model injection
of kinks takes place independently from left and right boundaries and it
continues up to the time when first injected kink reaches the opposite
boundary. However, in the long-range interacting model, dynamics of all atoms
in the chain is immediately affected by the presence of the topological
defects. Thus, the saturation time for a long-range interacting system, is
approximately two times shorter than for its short-range interacting
counterpart, $\tau\approx N/2c.$ It results from the fact that at this time
$N/2(l_{m}+l_{\delta})$ kinks entered the chain from right boundary and
$N/2(l_{m}+l_{\delta})$ kinks entered the chain from the left and the given
density of kinks already minimizes the total energy of the system.

The role of Coulomb repulsion can be effectively suppressed by increasing
$m_{K}$, whose inverse sets its relative strength with respect to the lattice
potential. In this case, a staircase structure in $Q(t)$ can be observed more
clearly. Even in a long-range interacting model the maximal number of kinks in
the system scales linearly with the system size, sf. inset in Fig. 4, which is
also the case for the nearest-neighbour interactions.

In order to observe these dynamics in experiments, it is necessary to resolve
the commensurate and incommensurate configurations. One possibility would be
to image the ions positions by applying non-destructive far-field measurements
[61] using the diffractive pattern obtained from the scattered light.
Alternatively, one could resort to AC-Stark shift methods [62] in order to
distinguish the positions of the ions more accurately. The above results
pertain to relevant parameters regime of a typical trapped ion experiment [25,
10], and experimental details such as the effect of the external Paul trap,
and an estimate of the temperature required to observe quantum suppositions of
kinks in trapped ions experiments are elaborated in [65].

##  6 Outlook

The effective long-wavelength description of Eq. (1), is a sine-Gordon model
with $m_{K}$ controlling the strength of the non-linearity ($\sim
m_{K}^{2}\cos\phi$), and the misfit parameter $\delta$ coupling to the
gradient of the field ($\sim\delta\partial_{x}\phi$). This model is known as
Pokrovsky-Talapov field theory [33, 66, 67], and it finds numerous
applications in ultra-cold gases [68, 24, 57, 69, 70, 71, 72], trapped ions
applied to nano-friction [10, 8, 9], magnetism [30, 73, 29, 74] and
tangentially also strongly correlated systems [75, 76]. This naturally
suggests broad applicability of our results to different platforms.

In addition to the trapped ion simulators discussed above, another example
occurs in the cold atoms implementation of Ref. [24], where the Pokrovsky-
Talapov model describes the dynamics of the phase difference between two Raman
tunnel-coupled one-dimensional Bose quasi-condensates. The wave-vector
difference of the Raman beams plays the role of the solitons chemical
potential ($\delta$ in the model of Eq. (1)). Quantum quench dynamics would
therefore become accessible through interferometric measurements of the
correlation functions of the phase difference [77, 78].

These two examples represent just a few of the several interdisciplinary
outreaches of our results which establish one more step towards merging the
fields of topological defects and dynamical phase transitions in low-
dimensional quantum simulators. In the far future one might imagine nano-scale
control of topological defects in real materials by quenching high field
optical potentials onto material crystals to controllablly affect material
properties. Experiments in trapped ion quantum simulators may offer the first
steps in such a direction.

###### Acknowledgements.

We thank R.J. Valencia-Tortora and V. Vuletić for useful discussions. J. M.
acknowledges E. Demler and V. Kasper for previous collaborations on related
topics. This project has been supported by the Deutsche Forschungsgemeinschaft
(DFG, German Research Foundation) – Project-ID 429529648 – TRR 306 QuCoLiMa
(“Quantum Cooperativity of Light and Matter”), and by the Dynamics and
Topology Centre funded by the State of Rhineland Palatinate and Topology
Centre funded by the State of Rhineland Palatinate. The authors gratefully
acknowledge the computing time granted on the supercomputer MOGON 2 at
Johannes Gutenberg-University Mainz (hpc.uni-mainz.de).

## References

  * [1] Pitaevskii L. Stringari S. Bose-Einstein Condensation and Superfluidity International series of monographs on physics (Oxford University Press) 2016.
  * [2] Ramanathan A. other Phys. Rev. Lett.1062011130401.
  * [3] Wright K. C. et al. Phys. Rev. Lett.1102013025302.
  * [4] Damski B. Zurek W. H. Phys. Rev. Lett.1042010160404.
  * [5] Göbel B., Mertig I. Tretiakov O. A. Physics Reports89520211.
  * [6] Wieder B. J. et al. Nat Rev Mater72022196.
  * [7] Vedmedenko E. Y. et al. J. Phys. D: Appl. Phys.532020453001.
  * [8] Liu L., Zhou M., Jin L., Li L., Mo Y., Su G., Li X., Zhu H. Tian Y. Friction72019199.
  * [9] Vanossi A., Bechinger C. Urbakh M. Nat Commun1120201.
  * [10] Gangloff D. A., Bylinskii A. Vuletić V. Phys. Rev. Research22020013380.
  * [11] Fogarty T. et al. Phys. Rev. Lett.1152015233602.
  * [12] Bohlein T. et al. Nat. Mater.112012126.
  * [13] Timm L. et al. Phys. Rev. Research22020033198.
  * [14] Nayak C., Simon S. H., Stern A., Freedman M. Das Sarma S. Rev. Mod. Phys.8020081083.
  * [15] Tokura Y. et al. Nature Reviews Physics12019126.
  * [16] Sarma S. D. et al. Physics today59200632.
  * [17] Price H., Chong Y., Khanikaev A., Schomerus H. et al. Journal of Physics: Photonics42022032501.
  * [18] Amico L., Boshier M., Birkl G. et al. AVS Quantum Science32021039201.
  * [19] Bland T. et al. Phys. Rev. Res.42022043171.
  * [20] Eckel S. et al. Nature5062014200.
  * [21] Chaikin P. M. Lubensky T. C. Topological defects (Cambridge University Press) 1995 p. 495–589.
  * [22] Frenkel Y. I. Kontorova T. Zh. Eksp. Teor. Fiz8193889.
  * [23] Frank F. C., van der Merwe J. H. Mott N. F. Proceedings of the Royal Society of London. Series A. Mathematical and Physical Sciences1981949205.
  * [24] Kasper V., Marino J. et al. Phys. Rev. B1012020224102.
  * [25] von Lindenfels D., Gräb O. other Phys. Rev. Lett.1232019080602.
  * [26] Linnet R. B. et al. Phys. Rev. Lett.1092012233005.
  * [27] Cetina M. et al. New J. Phys152013053001.
  * [28] Landa H., Cormick C. Morigi G. Condensed Matter52020.
  * [29] Braun O. M. Kivshar Y. S. The Frenkel-Kontorova model: concepts, methods, and applications (Springer) 2004.
  * [30] Bak P. Reports on Progress in Physics451982587.
  * [31] Pruttivarasin T. et al. New J. Phys132011075012.
  * [32] García-Mata I., Zhirov O. Shepelyansky D. Eur. Phys. J. D412007325.
  * [33] Aristov D. N. Luther A. Phys. Rev. B652002165412.
  * [34] Marcovitch S. Reznik B. Phys. Rev. A782008052303.
  * [35] Yakimenko A. I. et al. Phys. Rev. A912015033607.
  * [36] Kochetov B. A. other Phys. Rev. E1002019062202.
  * [37] Ma X. et al. Nat Commun1120201.
  * [38] Rajaraman R. Solitons and instantons. an introduction to solitons and instantons in quantum field theory. 1. reprint as a paperback (Jan 1987).
  * [39] Dauxois T. Peyrard M. Physics of Solitons (Cambridge University Press) 2006.
  * [40] Brox J. et al. Phys. Rev. Lett.1192017153602.
  * [41] Carretero-González R. et al. Communications in Nonlinear Science and Numerical Simulation1092022106123.
  * [42] Sachdev S. Physics world12199933.
  * [43] Fernandes R. M. et al. Annual Review of Condensed Matter Physics102019133.
  * [44] Shimshoni E., Morigi G. Fishman S. Phys. Rev. Lett.1062011010401.
  * [45] Gornostyrev Y. N. et al. Phys. Rev. B6019991013.
  * [46] Shimshoni E., Morigi G. Fishman S. Phys. Rev. A832011032308.
  * [47] Lazarides A., Tieleman O. Morais Smith C. Phys. Rev. B802009245418.
  * [48] Lerose A., Marino J., Žunkovič B., Gambassi A. Silva A. Phys. Rev. Lett.1202018130603.
  * [49] Lerose A., Žunkovič B., Marino J., Gambassi A. Silva A. Phys. Rev. B992019045128.
  * [50] Berdanier W., Marino J. Altman E. Phys. Rev. Lett.1232019230604.
  * [51] Marino J., Eckstein M., Foster M. S. Rey A. M. Reports on Progress in Physics852022116001.
  * [52] Polkovnikov A. Annals of Physics32520101790.
  * [53] Bistritzer R. Altman E. Proceedings of the National Academy of Sciences10420079955.
  * [54] Horváth D. X. et al. Phys. Rev. A1002019013613.
  * [55] Lancaster J., Gull E. Mitra A. Phys. Rev. B822010235124.
  * [56] Greene P. B., Kofman L. Starobinsky A. A. Nuclear Physics B5431999423.
  * [57] Neuenhahn C., Polkovnikov A. Marquardt F. Phys. Rev. Lett.1092012085304.
  * [58] Berry M. V. Mount K. Reports on Progress in Physics351972315.
  * [59] Reichl L. E. A modern course in statistical physics (1999).
  * [60] Braun O. M., Kivshar Y. S. Zelenskaya I. I. Phys. Rev. B4119907118.
  * [61] Wolf S. et al. Phys. Rev. Lett.1162016183002.
  * [62] Schmiegelow C. T. et al. Phys. Rev. Lett.1162016033002.
  * [63] Knap M. et al. Phys. Rev. Lett.1112013147205.
  * [64] Baltrusch J. D. et al. Phys. Rev. A842011063821.
  * [65] Chelpanova O., Kelly S., Schmidt-Kaler F., Morigi G. Marino J. In preparation (2023).
  * [66] Pokrovsky V. Talapov A. Publ., New York1984.
  * [67] Pokrovsky V. L. Talapov A. L. Phys. Rev. Lett.42197965.
  * [68] Lovas I. et al. Phys. Rev. B1062022075426.
  * [69] Gritsev V., Polkovnikov A. Demler E. Phys. Rev. B752007174511.
  * [70] Whitlock N. K. Bouchoule I. Phys. Rev. A682003053609.
  * [71] Büchler H. P., Blatter G. Zwerger W. Phys. Rev. Lett.902003130401.
  * [72] Haller E., Hart R., Mark M. J., Danzl J. G. et al. Nature4662010597.
  * [73] Caux J.-S., Essler F. H. L. Löw U. Phys. Rev. B682003134431.
  * [74] Chepiga N. arXiv preprint arXiv:2209.103902022.
  * [75] Aristov D. N., Cheianov V. V. Luther A. Phys. Rev. B662002073105.
  * [76] Lebwohl P. Stephen M. J. Phys. Rev.1631967376.
  * [77] Gring M., Kuhnert M., Langen T. et al. Science33720121318.
  * [78] Schweigler T. et al. Nature5452017.


#  Adaptive Testing for Connected and Automated Vehicles with
Sparse Control Variates in Overtaking Scenarios

Jingxuan Yang*  —, Honglin He*  —, Yi Zhang*  —, , Shuo Feng*  —,
and Henry X. Liu*  — This work is supported by National Key Research and
Development Program under Grant 2021YFB2501200 and National Natural Science
Foundation of China under Grant 62133002. (Corresponding author: Shuo
Feng.)Jingxuan Yang and Honglin He are with the Department of Automation,
Tsinghua University, Beijing 100084, China (email: {yangjx20,
hehl21}@mails.tsinghua.edu.cn).Yi Zhang is with the Department of Automation,
Beijing National Research Center for Information Science and Technology
(BNRist), Tsinghua University, Beijing 100084, China (e-mail:
zhyi@tsinghua.edu.cn).Shuo Feng is with the University of Michigan
Transportation Research Institute, Ann Arbor, MI 48109, USA (e-mail:
fshuo@umich.edu).Henry X. Liu is with the Department of Civil and
Environmental Engineering, University of Michigan, Ann Arbor, MI 48109, USA
(e-mail: henryliu@umich.edu).

###### Abstract

Testing and evaluation is a critical step in the development and deployment of
connected and automated vehicles (CAVs). Due to the black-box property and
various types of CAVs, how to test and evaluate CAVs adaptively remains a
major challenge. Many approaches have been proposed to adaptively generate
testing scenarios during the testing process. However, most existing
approaches cannot be applied to complex scenarios, where the variables needed
to define such scenarios are high dimensional. Towards filling this gap, the
adaptive testing with sparse control variates method is proposed in this
paper. Instead of adaptively generating testing scenarios, our approach
evaluates CAVs’ performances by adaptively utilizing the testing results.
Specifically, each testing result is adjusted using multiple linear regression
techniques based on control variates. As the regression coefficients can be
adaptively optimized for the CAV under test, using the adjusted results can
reduce the estimation variance, compared with using the testing results
directly. To overcome the high dimensionality challenge, sparse control
variates are utilized only for the critical variables of testing scenarios. To
validate the proposed method, the high-dimensional overtaking scenarios are
investigated, and the results demonstrate that our approach can further
accelerate the evaluation process by about 30 times.

###### Index Terms:

Adaptive testing, connected and automated vehicles, sparse control variates,
overtaking scenarios

##  I Introduction

Testing and evaluation are major challenges for the development and deployment
of connected and automated vehicles (CAVs). The past few years have witnessed
increasingly rapid advances in the field of testing scenario library
generation (TSLG) [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12]. The goal of TSLG is
usually to purposely generate safety-critical testing scenarios that can
improve the evaluation efficiency of CAVs while ensuring the evaluation
unbiasedness. As CAVs are usually black boxes, to evaluate the criticality
values of different testing scenarios, surrogate models (SMs) are usually
constructed by leveraging prior knowledge of CAVs. Due to the various types of
CAVs, however, the performance dissimilarities between SMs and the CAVs under
test usually exist, which may compromise the effectiveness of the testing
scenarios and decrease the evaluation efficiency. Therefore, how to adaptively
test different types of CAVs becomes a critical problem.

Towards addressing this problem, several adaptive testing methods have been
proposed [13, 14, 15, 16]. The basic idea of these existing methods is to
adaptively adjust the testing scenarios by leveraging the testing results of
CAVs during the testing process. With more testing results of CAVs, more
posteriori knowledge of CAVs can be obtained, and therefore the testing
scenarios can be more customized and optimized for the CAVs under test, which
can adaptively improve the testing efficiency [16]. However, most existing
adaptive testing methods have limitations in dealing with complex scenarios,
because of the “Curse of Dimensionality” (CoD) problem. For example, the CAV
overtaking scenarios with two background vehicles could have 12-dimensional
states (lateral and longitudinal positions and velocities of all three
vehicles) and 3-dimensional actions (i.e., accelerations). Although the
overtaking scenarios seem simple, if they last for 10 seconds at a frequency
of 10 Hz, the dimension of the scenarios could exceed 103, which cannot be
handled by most existing adaptive testing methods.

The goal of this paper is to develop adaptive testing methods for high-
dimensional scenarios using control variates (CVs). Instead of adaptively
generating testing scenarios, we evaluate CAVs’ performances by adaptively
utilizing the testing results. The CVs are some random variables with means
known [17], which usually correlate the performance index of the event of
interest. Through the regression of control parameters associated with the
CVs, the testing results can be adaptively adjusted into a much narrower
interval, which can greatly reduce the estimation variance. The worse the
prior knowledge of CAVs is, the worse the testing results are, and the more
they could be improved by the CVs. However, there also exists the CoD problem
if we directly apply traditional CVs, since the number of control parameters
would increase exponentially with the dimension of the scenarios.

To address the CoD problem, we propose the adaptive testing with sparse
control variates (ATSCV) method in this paper to apply sparse control variates
(SCVs) to the testing results, which only controls the variates associated
with the critical variables of testing scenarios. To identify the critical
variables, we applied the naturalistic and adversarial driving environment
(NADE) [18], which can generate high-dimensional testing scenarios by
adjusting the critical variables of the naturalistic driving environment
(NDE). Specifically, the critical variables are the scenes where NADE makes
adversarial adjustments to the maneuvers of background vehicles, the dimension
of which is much smaller than that of the scenarios. By generating control
variates for the sparse critical variables, the SCVs can address the CoD
problem for the adaptive testing method. Then, the SCVs are determined by
multiple linear regression (MLR) techniques [19, 20], which can essentially
adjust the testing results. As the MLR techniques can leverage the posteriori
knowledge of CAVs, the adjusted testing results can better evaluate the
performance of the CAVs under test, resulting in an adaptive testing method.
We note that this method is complementary to the adaptive testing methods that
can adaptively generate testing scenarios during the testing process. To
validate the proposed ATSCV method, the overtaking scenarios are investigated.
Comparing with the estimation efficiency in NADE, the new adaptive testing
method can further accelerate the evaluation process by about 30 times.

The rest of this paper is organized as follows. Section II introduces the
overtaking scenarios and presents the testing scenario library generation
methods in NDE and NADE. In Section III, the adaptive testing with sparse
control variates method is proposed. Section IV evaluates the accuracy and
efficiency of the ATSCV method with a case study on overtaking scenarios.
Finally, Section V concludes the paper.

##  II Problem Formulation

###  II-A Overtaking Scenarios

Figure 1: Illustration of the overtaking scenarios.

The overtaking scenarios are shown in Fig. 1. The key to handle the TSLG
problem is to model the overtaking scenarios as Markov decision processes
(MDPs). In overtaking scenarios, the leading vehicle (LV) in the left lane
runs at a constant speed, whose state (including longitudinal position,
lateral position and longitudinal velocity) and action (i.e., acceleration)
are denoted as
$s_{\mathrm{LV}}=(x_{\mathrm{LV}},y_{\mathrm{LV}},v_{\mathrm{LV}})$ and
$a_{\mathrm{LV}}=0$, respectively. The background vehicle (BV) follows LV,
whose state is denoted as
$s_{\mathrm{BV}}=\left(x_{\mathrm{BV}},y_{\mathrm{BV}},v_{\mathrm{BV}}\right)$.
Let $\mathcal{A}$ be the action space (i.e., acceleration space), and
$\mathcal{A}^{+}\triangleq\mathcal{A}\cup\\{\Rsh\\}$ be the total action
space, where $\Rsh$ represents right lane change. The action of BV
$a_{\mathrm{BV}}\in\mathcal{A}^{+}$. Specifically,
$a_{\mathrm{BV}}\in\mathcal{A}$ before BV cuts in, $a_{\mathrm{BV}}={\Rsh}$
when BV cuts in and after that $a_{\mathrm{BV}}=0$. The automated vehicle (AV)
in the right lane runs at a constant speed before BV cuts in and follows BV
after that, whose state is denoted as
$s_{\mathrm{AV}}=\left(x_{\mathrm{AV}},y_{\mathrm{AV}},v_{\mathrm{AV}}\right)$.
The action of AV $a_{\mathrm{AV}}=0$ before BV cuts in, and
$a_{\mathrm{AV}}\in\mathcal{A}$ after that.

The state of an overtaking scenario is the collection of three vehicles’
states $(s_{\mathrm{LV}},s_{\mathrm{BV}},s_{\mathrm{AV}})$. Since these state
variables have interdependencies, we can formulate the state as

$s\triangleq\left(v_{\mathrm{BV}},R_{1},\dot{R}_{1},R_{2},\dot{R}_{2}\right)\in\mathcal{S},$
(1)

where $\mathcal{S}$ is the set of all feasible states,

$\displaystyle R_{1}$
$\displaystyle=x_{\mathrm{LV}}-x_{\mathrm{BV}},\quad\dot{R}_{1}=v_{\mathrm{LV}}-v_{\mathrm{BV}},$
(2) $\displaystyle R_{2}$
$\displaystyle=x_{\mathrm{BV}}-x_{\mathrm{AV}},\quad\dot{R}_{2}=v_{\mathrm{BV}}-v_{\mathrm{AV}}.$

The action of the overtaking scenario is defined as the action of BV, that is,
$a\triangleq a_{\mathrm{BV}}\in\mathcal{A}^{+}$. Then an overtaking scenario
is given by

$x=\left(s_{0},a_{0},s_{1},a_{1},\dots,s_{m},a_{m}\right)\in\mathcal{X},$ (3)

where $\mathcal{X}$ is the set of all feasible scenarios,

$s_{k}=\left(v_{\mathrm{BV},k},R_{1,k},\dot{R}_{1,k},R_{2,k},\dot{R}_{2,k}\right),~{}k=0,1,\dots,m$
(4)

is the state and $a_{k}=a_{\mathrm{BV},k}$ is the action of BV at
$k^{\mathrm{th}}$ time step, $0\leqslant m\leqslant M$, where $M$ is the
maximum time step. When AV rear-ends BV after BV cuts in or BV does not cut in
until AV passes BV, the overtaking scenarios terminate.

Suppose that the overtaking scenarios are simulated for 10 seconds at a
frequency of 10 Hz, then the dimension of the scenario will exceed 103,
resulting in the CoD problem. We note that the overtaking scenarios are
significantly different from the simple scenarios that are usually studied
such as car-following scenarios [10, 12] and cut-in scenarios [8, 10, 11],
where the accident events rely heavily on the initial states. On the contrary,
the accident event of overtaking scenarios is more stochastic and complicated,
because the right lane-changing action of BV is probabilistic and BV may have
many chances of time steps to cut in, which essentially results in different
cut-in and car-following scenarios. That is also the reason why the cut-in
scenarios are usually low dimensional, while the overtaking scenarios are much
more high dimensional. Moreover, the riskiest scenarios for the upcoming Level
3 automated lane keeping system (ALKS) will arguably be the overtaking
scenarios [21]. Therefore, the overtaking scenarios are selected for the case
study.

###  II-B Testing Scenario Library Generation

Let $\Omega=\mathcal{X}$ be the sample space, including all feasible
overtaking scenarios. Consider the probability space
$(\Omega,\mathcal{F},\mathbb{P})$, where $\mathcal{F}$ is a $\sigma$-algerbra
of subsets of $\Omega$ and $\mathbb{P}$ is a probability measure on
$\mathcal{F}$. Since the scenarios are discretized, we choose
$\mathcal{F}=\mathcal{P}(\Omega)$ to be the power set of $\Omega$. Then an
event $A\subset\Omega$ is a subset of scenarios in the sample space $\Omega$.
Let $X:x\mapsto x$, $x\in\mathcal{X}$ be the random variable of scenarios. To
verify and validate AV, the accident event is usually the event of interest,
which is defined as $A=\left\\{x\in\mathcal{X}\,\,|R_{1,m}\leqslant
d_{\mathrm{accid}}\right\\}$, where $d_{\mathrm{accid}}$ is the distance
threshold for accident. The accident rate $\mu=\mathbb{P}(A)$ is utilized as
the performance index of AV, which can also be represented as an expectation,
i.e.,

$\mu=\mathbb{E}_{p}[\mathbb{I}_{A}(X)]=\sum_{x\in\mathcal{X}}\mathbb{P}(A|x)p(x),$
(5)

where $\mathbb{I}_{A}(X)$ is the indicator function of the event $A$,

$\mathbb{I}_{A}(X)=\begin{cases}1,&\mathrm{if}~{}X\in A,\\\
0,&\mathrm{if}~{}X\notin A,\\\ \end{cases}$ (6)

and $p$ is the naturalistic joint distribution of $x$. Assuming the Markovian
property, the joint distribution can be factorized as

$p(x)=p(s_{0})\prod_{k=0}^{m}p(a_{k}|s_{k}).$ (7)

The essence of testing AV in naturalistic driving environment (NDE) is to
estimate the performance index $\mu$ by Monte Carlo simulation, i.e.,

$\hat{\mu}_{n}=\frac{1}{n}\sum_{i=1}^{n}\mathbb{P}(A|X_{i}),\quad X_{i}\sim
p.$ (8)

The NDE suffers from the inefficiency problem, as the accident event is
usually a rare event. To improve the efficiency of Monte Carlo simulation,
importance sampling (IS) is adopted to test AV with scenarios generated from a
different probability distribution $q$, which is called the importance
function and needs to satisfy $q(x)>0$ whenever $\mathbb{P}(A|x)p(x)>0$ [11,
9, 10]. Since

$\displaystyle\mu$
$\displaystyle=\sum_{x\in\mathcal{X}}\mathbb{P}(A|x)p(x)=\sum_{x\in\mathcal{X}}\frac{\mathbb{P}(A|x)p(x)}{q(x)}q(x)$
(9)
$\displaystyle=\mathbb{E}_{q}\left[\frac{\mathbb{P}(A|X)p(X)}{q(X)}\right],$

the performance index can be estimated as

$\displaystyle\hat{\mu}_{q}$
$\displaystyle=\frac{1}{n}\sum_{i=1}^{n}\frac{\mathbb{P}(A|X_{i})p(X_{i})}{q(X_{i})},\quad
X_{i}\sim q$ (10)
$\displaystyle=\frac{1}{n}\sum_{i=1}^{n}\mathbb{P}(A|X_{i})\frac{p(s_{0,i})}{q(s_{0,i})}\prod_{k=0}^{m_{i}}\frac{p(a_{k,i}|s_{k,i})}{q(a_{k,i}|s_{k,i})},$

where $X_{i}=(s_{0,i},a_{0,i},\dots,s_{m_{i},i},a_{m_{i},i})$, $i=1,\dots,n$.

###  II-C Naturalistic and Adversarial Driving Environment Generation

The IS method is able to overcome the inefficiency problem with properly
selected importance function and in fact, there exists optimal importance
function whose estimation variance is zero [9]. However, the IS method faces
the CoD problem if the testing scenarios are high-dimensional. To address both
the inefficiency problem and the CoD problem, the naturalistic and adversarial
driving environment (NADE) has been proposed to only sample critical variables
with IS, while other variables remain its naturalistic distribution [18].

Denote $x=(x_{c},x_{-c})$, where $x_{c}=\\{x_{c_{1}},\dots,x_{c_{l}}\\}$,
$l=0,1,\dots,L$ is the set of critical variables, $c_{1},\dots,c_{l}$ are
called the critical moments, $L$ is the maximum number of control steps among
all scenarios $x\in\mathcal{X}$ and $x_{-c}$ is the set of other variables.
Let $X_{c}:x\mapsto x_{c}$ be the random variable of critical variables and
$X_{-c}:x\mapsto x_{-c}$ be the random variable of other variables. Then we
have $X=(X_{c},X_{-c})$. The random variable $X_{c}$ controls how to identify
critical variables, and thus it’s also called the control policy. The
importance function can then be formulated as $q(x)=q(x_{c})p(x_{-c})$, and
therefore the performance index can be estimated in NADE as

$\displaystyle\tilde{\mu}_{q}$
$\displaystyle=\frac{1}{n}\sum_{i=1}^{n}\frac{\mathbb{P}(A|X_{i})p(X_{c,i})}{q(X_{c,i})},\quad
X_{i}\sim q$ (11)
$\displaystyle=\frac{1}{n}\sum_{i=1}^{n}\mathbb{P}(A|X_{i})\frac{p(s_{0,i})}{q(s_{0,i})}\prod_{k\in\mathcal{T}_{c,i}}\frac{p(a_{k,i}|s_{k,i})}{q(a_{k,i}|s_{k,i})},$

where $X_{c,i}$ is the random variable of critical variables of $X_{i}$, and
$\mathcal{T}_{c,i}$ is the set of critical moments of the $i^{\mathrm{th}}$
test.

It can be shown that with appropriate control policies and importance
functions, NADE is able to greatly address the CoD problem and increase the
evaluation efficiency [18].

##  III Adaptive Testing with Sparse Control Variates

###  III-A Control Variates

Control variates (CVs) can be usefully combined with the mixture importance
sampling. In mixture IS, we sample $X_{1},\dots,X_{n}$ from the mixture
importance function $q_{\alpha}=\sum_{j=1}^{J}\alpha_{j}q_{j}$, where
$\alpha_{j}\geqslant 0$, $\sum_{j=1}^{J}\alpha_{j}=1$ and the $q_{j}$ are
importance functions. A control variate is a vector

$h(x)=(h_{1}(x),\dots,h_{J}(x))^{\top}$ (12)

for which $\sum_{x\in\mathcal{X}}h(x)=\theta$, where $\theta$ is a known
value. Using the importance functions as control variates provides at least as
good a variance reduction as we get from ordinary importance sampling [17].
When we combine mixture IS from $q_{\alpha}$ with control variates based on
the component densities, the estimated performance index is

$\hat{\mu}_{q_{\alpha},\beta}=\frac{1}{n}\sum_{i=1}^{n}\frac{\mathbb{P}(A|X_{i})p(X_{i})-\sum_{j=1}^{J}\beta_{j}q_{j}(X_{i})}{\sum_{j=1}^{J}\alpha_{j}q_{j}(X_{i})}+\sum_{j=1}^{J}\beta_{j}$
(13)

for $X_{i}\sim q_{\alpha}$, where $\beta_{j}$ are control parameters.

We can compare the variance of mixture importance sampling to that of
importance sampling with the individual mixture components $q_{j}$. Let
$\beta^{*}$ be any minimizer over $\beta$ of
$\mathrm{Var}_{q_{\alpha}}(\hat{\mu}_{q_{\alpha},\beta})$. It can be proved
[22] that

$\mathrm{Var}_{q_{\alpha}}(\hat{\mu}_{q_{\alpha},\beta^{*}})\leqslant\min_{1\leqslant
j\leqslant J}\frac{\sigma_{q_{j}}^{2}}{n\alpha_{j}},$ (14)

where

$\sigma_{q_{j}}^{2}=\mathrm{Var}_{q_{j}}\left(\frac{\mathbb{P}(A|X)p(X)}{q_{j}(X)}\right),~{}j=1,\dots,J.$
(15)

If any one of the $q_{j}$ is optimal then we will get a zero variance
estimator of $\mu$. In practice, the optimal value $\beta^{*}$ is not known,
and we can estimate it by multiple linear regression (MLR). Letting
$Y_{i}=\mathbb{P}(A|X_{i})p(X_{i})/q_{\alpha}(X_{i})$, $i=1,\dots,n$ and
$Z_{ij}=q_{j}(X_{i})/q_{\alpha}(X_{i})-1$, $i=1,\dots,n$, $j=1,\dots,J-1$,
then the $\beta^{*}$ can be estimated as the coefficients obtained from MLR of
$Y_{i}$ on $Z_{ij}$.

###  III-B CoD of Control Variates

Considering the Markov chain structure of overtaking scenarios, the mixture
importance function of states is

$q_{\alpha}(s)=\sum_{j=1}^{J}{\alpha_{j}q_{j}(s)},~{}\forall s\in\mathcal{S},$
(16)

and the mixture importance function of actions given certain state
$s\in\mathcal{S}$ is

$q_{\alpha}(a|s)=\sum_{j=1}^{J}{\alpha_{j}q_{j}(a|s)},~{}\forall
a\in\mathcal{A}^{+}.$ (17)

Therefore, the mixture importance function of the overtaking scenarios is

$q_{\alpha}(x)=q_{\alpha}(s_{0})\prod_{k=0}^{m}q_{\alpha}(a_{k}|s_{k}),~{}\forall
x\in\mathcal{X},$ (18)

which is the product of $m+2$ individual mixture importance functions.
According to Eq. (13), the control variates are

$q_{j_{0},\dots,j_{m+1}}(x)=q_{j_{0}}(s_{0})q_{j_{1}}(a_{0}|s_{0})\cdots
q_{j_{m+1}}(a_{m}|s_{m}),$ (19)

where $j_{0},\dots,j_{m+1}=1,\dots,J$. Denote the corresponding control
parameters as $\beta_{j_{0},\dots,j_{m+1}}$.

The number of control parameters is $J^{m+2}$, which will increase
exponentially with the dimension of scenarios, leading to the CoD problem when
solving optimal control parameters via least squares method.

###  III-C Sparse Control Variates

To address the CoD problem of CVs, we propose the adaptive testing with sparse
control variates (ATSCV) method to apply sparse control variates (SCVs) for
testing results in NADE. The SCVs are constructed by the importance functions
of critical variables of testing scenarios. We use “sparse” to describe these
CVs because the critical variables are sparse in NADE. The number of critical
variables can be limited to a much smaller number than the whole dimension of
scenarios, and so can the number of SCVs, which will greatly address the CoD
problem of CVs.

Using scenarios sampled from $q_{\alpha}$, the performance index can be
estimated in NADE as

$\tilde{\mu}_{q_{\alpha}}=\frac{1}{n}\sum_{i=1}^{n}\frac{\mathbb{P}(A|X_{i})p(X_{c,i})}{q_{\alpha}(X_{c,i})},\quad
X_{i}\sim q_{\alpha}.$ (20)

Let $\mathcal{X}_{l}=\\{x\in\mathcal{X}:|x_{c}|=l\\}$, $l=0,1,\dots,L$ be the
set of scenarios that are controlled $l$ steps, satisfying
$\bigcup_{l=0}^{L}\mathcal{X}_{l}=\mathcal{X}$. Denote the performance index
of scenarios over $\mathcal{X}_{l}$ as

$\mu_{l}\triangleq\mathbb{E}_{p}[\mathbb{I}_{A}(X)\mathbb{I}_{\mathcal{X}_{l}}(X)],~{}l=0,1,\dots,L.$
(21)

Similar to Eq. (20), the estimation of $\mu_{l}$ can be written as

$\tilde{\mu}_{l,q_{\alpha}}=\frac{1}{n}\sum_{i=1}^{n}\frac{\mathbb{P}(A|X_{i})\mathbb{I}_{\mathcal{X}_{l}}(X_{i})p(X_{c,i})}{q_{\alpha}(X_{c,i})},$
(22)

then we have

$\mu=\sum_{l=0}^{L}\mu_{l},$ (23)

and

$\displaystyle\tilde{\mu}_{q_{\alpha}}$
$\displaystyle=\sum_{l=0}^{L}\frac{1}{n}\sum_{i=1}^{n}\frac{\mathbb{P}(A|X_{i})\mathbb{I}_{\mathcal{X}_{l}}(X_{i})p(X_{c,i})}{q_{\alpha}(X_{c,i})}=\sum_{l=0}^{L}\tilde{\mu}_{l,q_{\alpha}}.$
(24)

Let $q_{j_{1},\dots,j_{l}}(x_{c})=q_{j_{1}}(x_{c_{1}})\cdots
q_{j_{l}}(x_{c_{l}})$ be the importance functions of critical variables that
sample $x_{c_{1}},\dots,x_{c_{l}}$ from $q_{j_{1}},\dots,q_{j_{l}}$
respectively, where $j_{1},\dots,j_{l}=1,\dots,J$, $l=1,\dots,L$. We deal with
the CoD problem by applying $q_{j_{1},\dots,j_{l}}$ instead of
$q_{j_{0},\dots,j_{m+1}}$ as SCVs to the estimation
$\tilde{\mu}_{l,q_{\alpha}}$. If the control steps $l\ll m$, then the number
of control variates can be reduced considerably, which will greatly address
the CoD problem. For a good enough control policy, this assumption can be
guaranteed. Since the numbers of SCVs with different control steps are
different, we group the testing results by their control steps and apply the
MLR within each group.

Denote the summation of the product of control variates and corresponding
control parameters as

$h_{l}(x_{c})\triangleq\sum_{j_{1},\dots,j_{l}}\beta_{j_{1},\dots,j_{l}}q_{j_{1},\dots,j_{l}}(x_{c}),~{}l=1,\dots,L,$
(25)

and the summation of $h_{l}(x_{c})$ over $x\in\mathcal{X}_{l}$ as

$\theta_{l}\triangleq\sum_{x\in\mathcal{X}_{l}}h_{l}(x_{c}),~{}l=1,\dots,L.$
(26)

With the SCVs $q_{j_{1},\dots,j_{l}}$, the estimation
$\tilde{\mu}_{l,q_{\alpha}}$ can be adaptively evaluated as

$\tilde{\mu}_{l,q_{\alpha},\beta_{l}}=\frac{1}{n}\sum_{i=1}^{n}\frac{\mathbb{P}(A|X_{i})p(X_{c,i})-h_{l}(X_{c,i})}{q_{\alpha}(X_{c,i})}\mathbb{I}_{\mathcal{X}_{l}}(X_{i})+\theta_{l}$
(27)

for $l=1,\dots,L$. Note that for $l=0$, the scenarios are sampled from
naturalistic distribution $p$, and therefore we do not apply SCVs to its
estimation, i.e.,

$\tilde{\mu}_{0,q_{\alpha},\beta_{0}}=\frac{1}{n}\sum_{i=1}^{n}\mathbb{P}(A|X_{i})\mathbb{I}_{\mathcal{X}_{0}}(X_{i}).$
(28)

Therefore, the performance index estimated by the ATSCV method is

$\tilde{\mu}_{q_{\alpha},\beta}=\sum_{l=0}^{L}\tilde{\mu}_{l,q_{\alpha},\beta_{l}}.$
(29)

##  IV Simulation Analysis

###  IV-A Generation of NDE

The goal of NDE is to generate human-like driving behaviors with probability
distributions consistent with the naturalistic driving data (NDD). Since the
overtaking scenarios are modeled with MDPs, after generating the initial
state, the actions of vehicles are determined by current states, which can be
sampled from empirical distributions obtained from NDD. The decision process
can be represented by a decision tree, where each node denotes vehicle states
and each path denotes a specific realization of vehicle behaviors. For
simplicity, as a replacement of empirical distributions, we use intelligent
driver model (IDM) [23] to model car-following scenarios and stochastic
minimizing overall braking induced by lane changes (MOBIL) model [18] to
govern lane-changing behaviors. All vehicles select actions independently and
simultaneously for each time step (0.1 s).

The initial state of overtaking scenarios is set as

$s_{0}=[8,R_{1},-5,5,-5],~{}R_{1}\sim\mathcal{U}(30,32),$ (30)

where $\mathcal{U}$ is the uniform distribution. The overtaking scenario has
two stages separated by the time when BV cuts in. Before cutting in, BV is
controlled by IDM and stochastic MOBIL model, while LV and AV run at a
constant speed. After cutting in, AV is controlled by IDM, while LV and BV run
at a constant speed. The cut-in maneuver is set completed within one time
step. The simulation continues until AV rear-ends BV or AV passes BV.

###  IV-B Generation of NADE

The key of NADE generation is to construct new distributions as the
replacement of the naturalistic distributions in NDE, for sampling actions of
BV [18]. In NADE, we only twist the action distribution of BV at critical
moments, while keeping naturalistic distribution as in NDE at other time
steps. To construct the importance function, at each time step, the action of
BV will be evaluated by maneuver criticality, which is computed as a
multiplication of exposure frequency $p(a|s)$ and maneuver challenge
$\mathbb{P}(A|s,a)$, i.e.,

$V(a|s)=\mathbb{P}(A|s,a)p(a|s).$ (31)

The exposure frequency represents the probability of the action $a$ given the
state $s$ in NDE, while the maneuver challenge measures the accident
probability of AV given the state-action pair $(s,a)$.

Let $p_{R}$ be the probability of right lane change induced from stochastic
MOBIL model, then the exposure frequency is

$p(a|s)=\begin{cases}p_{R},&\mathrm{if}~{}a={\Rsh},\\\
1-p_{R},&\mathrm{if}~{}a=\mathrm{IDM}(s)\in\mathcal{A},\\\
0,&\mathrm{otherwise},\\\ \end{cases}$ (32)

where $\mathrm{IDM}(s)$ is the acceleration calculated by IDM given the state
$s$. Since the AV model is usually a black-box, we use surrogate models (SMs)
to approximate the maneuver challenge. In this paper, IDM and full velocity
difference model (FVDM) [23] are selected as SMs with adjusted parameters as
follows: 1) IDM. 2) FVDM-I. 3) FVDM-II with $\kappa=6$.

Let $S_{j}$, $j=1,\dots,J$ denotes the accident event between BV and the
$j^{\mathrm{th}}$ SM, where $J=3$, then the maneuver challenge can be
approximated by $\mathbb{P}(S_{j}|s,a)$, $j=1,\dots,J$. If BV follows LV with
$a\in\mathcal{A}$, then we have

$\mathbb{P}(S_{j}|s,a)=\sum_{a\in\mathcal{A}^{+}}{\mathbb{P}(S_{j}|s^{\prime},a)p(a|s^{\prime})},$
(33)

where $s^{\prime}$ is the next state of the overtaking scenario given the
state-action pair $(s,a)$. The maneuver challenge $\mathbb{P}(S_{j}|s,a)$ can
be solved by recursion with Eq. (33) and the following two terminal
conditions.

  1. 1.

If BV gets past by AV, then there will be no accident. Therefore,
$\mathbb{P}(S_{j}|s,a)=0$ for $R_{2}<0$.

  2. 2.

After BV cuts in, whether the accident will occur is predicted by the
$j^{\mathrm{th}}$ SM and the state after cutting in. Thus, for $a={\Rsh}$,
$\mathbb{P}(S_{j}|s,a)=1$ if AV rear-ends BV, and otherwise
$\mathbb{P}(S_{j}|s,a)=0$.

Input: $p$, $q_{\alpha}$, $X_{c,i}$, and $\mathbb{P}(A|X_{i})$, $i=1,\dots,n$

Output: $\tilde{\mu}_{q_{\alpha},\hat{\beta}}$

1 initialize $Y_{l}$ and $Z_{l}$ as empty arrays, $l=0,\dots,L$;

2 for  _$i\leftarrow 1$ to $n$_ do

3 $l\leftarrow$ number of control steps of $X_{c,i}$;

4 if  _$l=0$_ then

5 append $Y_{l}$ with $\mathbb{P}(A|X_{i})$;

6 append $Z_{l}$ with 0;

7

8 else

9 append $Y_{l}$ with $\mathbb{P}(A|X_{i})p(X_{c,i})/q_{\alpha}(X_{c,i})$;

10 append $Z_{l}$ with
$\mathrm{vec}(q_{j_{1},\dots,j_{l}}(X_{c,i})/q_{\alpha}(X_{c,i}))$,
$j_{1},\dots,j_{l}=1,\dots,J-1$;

11 /* $\mathrm{vec}(\cdot)$ is the vectorization operator */

12 end if

13

14 end for

15for  _$l\leftarrow 0$ to $L$_ do

16 $Z_{l}\leftarrow Z_{l}-\mathrm{average}(Z_{l})$;

17 $\mathrm{MLR}\leftarrow$ multiple linear regression of $Y_{l}$ on $Z_{l}$;

18 $\hat{\beta}_{l}\leftarrow$ estimated coefficients from $\mathrm{MLR}$;

19 $\hat{\eta}_{l}\leftarrow$ estimated intercept from $\mathrm{MLR}$;

20
$\tilde{\mu}_{l,q_{\alpha},\hat{\beta}_{l}}\leftarrow\mathrm{length}(Y_{l})\hat{\eta}_{l}/n$;

21

22 end for

23return
$\tilde{\mu}_{q_{\alpha},\hat{\beta}}\leftarrow\sum_{l=0}^{L}\tilde{\mu}_{l,q_{\alpha},\hat{\beta}_{l}}$;

Algorithm 1 Adaptive testing with sparse control variates by multiple linear
regression

The criticality of BV can then be calculated as the summation of maneuver
criticality over all actions, i.e.,

$C_{j}(s)=\sum_{a\in\mathcal{A}^{+}}V_{j}(a|s),~{}j=1,\dots,J,$ (34)

where $V_{j}(a|s)=\mathbb{P}(S_{j}|s,a)p(a|s)$. If any $C_{j}(s)>0$, then the
moment of state $s$ is defined as the critical moment and the action of BV
will be adjusted by sampling from the mixture importance function
$q_{\alpha}(a|s)=\sum_{j=1}^{J}\alpha_{j}q_{j}(a|s)$, where $\alpha_{j}=1/3$,
$j=1,\dots,J$ is the weight of the $j^{\mathrm{th}}$ SM,

$q_{j}(a|s)=\begin{cases}\epsilon
p(a|s)+(1-\epsilon)\frac{V_{j}(a|s)}{C_{j}(s)},&\mathrm{if}~{}C_{j}(s)>0,\\\
p(a|s),&\mathrm{otherwise},\end{cases}$ (35)

where $\epsilon=0.1$ is the weight of the naturalistic distribution.

The ATSCV method can be applied while evaluating the testing results in NADE.
Algorithm 1 describes this process. The key is to identify the critical
variables of testing scenarios to formulate the sparse control variates for
each testing result. Then the performance index of different control steps can
be estimated as the intercept obtained from multiple linear regression of
weighted testing results on corresponding sparse control variates.

###  IV-C Evaluation Results

The accuracy and efficiency of AV evaluation in NDE and NADE are validated in
our simulation. Fig. 2 shows the evaluation results of the accident rate per
test in NDE and NADE. The gray line represents the testing results in NDE, and
the bottom $x$-axis indicates the number of tests. The blue line represents
the testing results in NADE, and the top $x$-axis for the number of tests. It
can be seen that NADE obtains the same accident rate estimation with NDE by a
much smaller number of tests.

Figure 2: Accident rate of AV in NDE and NADE. Figure 3: RHW of AV evaluation
in NDE and NADE.

The relative half-width (RHW) [11] is selected to measure the evaluation
precision. With the confidence level $100(1-\gamma)\%$, the RHW of
$\tilde{\mu}_{q_{\alpha},\beta}$ is defined as

$l_{r}=z_{\gamma}\frac{\sqrt{\mathrm{Var}(\tilde{\mu}_{q_{\alpha},\beta})}}{\tilde{\mu}_{q_{\alpha},\beta}},$
(36)

where $z_{\gamma}=\Phi^{-1}(1-\gamma/2)$, $\gamma=0.1$, and $\Phi^{-1}$ is the
inverse cumulative distribution function of the standard normal distribution
$\mathcal{N}(0,1)$. Fig. 3 presents the minimum number of tests for reaching a
predetermined precision threshold (RHW is 0.1) in NDE and NADE. It can be
found that NADE requires only 385 number of tests, while NDE requires 5.5
$\times$ 104 number of tests, resulting in an acceleration of evaluation for
about 143 times.

Figure 4: Accident rate of AV in NADE with and without SCVs. Figure 5: RHW of
AV evaluation in NADE with and without SCVs.

To investigate the influences of the ATSCV method, we compare the efficiency
of AV evaluation in NADE with and without SCVs. As shown in Fig. 4, NADE and
ATSCV can converge to the same accident rate after a sufficient number of
tests. To reach the 0.1 relative half-width, the total required number of
tests is 385 and 12, respectively, as shown in Fig. 5. Therefore, the ATSCV
method can further accelerate the evaluation by about 30 times comparing with
the evaluation results in NADE. Fig. 6 presents the adjusted testing results
in NADE with and without SCVs, where each blue point is the value

$\frac{\mathbb{P}(A|X_{i})p(X_{c,i})}{q_{\alpha}(X_{c,i})},~{}i=1,\dots,n$
(37)

in Eq. (20), and each red point represents the value

$\frac{\mathbb{P}(A|X_{i})p(X_{c,i})-h_{l}(X_{c,i})}{q_{\alpha}(X_{c,i})}+\theta_{l}$
(38)

in Eq. (27) for $X_{i}\in\mathcal{X}_{l}$, $l=1,\dots,L$ or
$\mathbb{P}(A|X_{i})$ for $X_{i}\in\mathcal{X}_{0}$. It can be seen that the
ATSCV method is able to adjust the testing results into a much narrower
interval, resulting in a considerable reduction of the estimation variance.

To further verify the reliability of the ATSCV method, the simulation is
repeated 100 times. Results show that the average accelerated rate is 28.34
and the standard deviation is 6.27. The required numbers of tests in all 100
experiments are less than the evaluation results in NADE.

Figure 6: Adjusted testing results in NADE with and without SCVs.

##  V Conclusion

In this paper, the adaptive testing with sparse control variates (ATSCV)
method is proposed to adaptively evaluate the testing results of connected and
automated vehicles in complex scenarios. The major idea is to apply SCVs of
critical variables for the evaluation of scenarios in NADE. With scenarios
modeled as MDPs, the mixture importance functions are the summation of
weighted defensive importance functions obtained using SMs. The SCVs are
constructed as the importance functions of critical variables, each of which
is the product of a combination of importance functions obtained by different
SMs at critical moments. The overtaking scenarios are investigated for safety
testing. Comparing with the evaluation results in NDE and NADE, the proposed
ATSCV method is always more efficient. In the near future, the theoretical
analysis with rigorous proofs of the ATSCV method will be developed and more
realistic cases with large-scale naturalistic driving data will be studied.

## References

  * [1] A. Li, S. Chen, L. Sun, N. Zheng, M. Tomizuka, and W. Zhan, “Scegene: Bio-inspired traffic scenario generation for autonomous driving testing,” _IEEE Transactions on Intelligent Transportation Systems_ , 2021.
  * [2] J. Wang, A. Pun, J. Tu, S. Manivasagam, A. Sadat, S. Casas, M. Ren, and R. Urtasun, “Advsim: Generating safety-critical scenarios for self-driving vehicles,” in _Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition_ , 2021, pp. 9909–9918.
  * [3] 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.
  * [4] Y. Tian, K. Pei, S. Jana, and B. Ray, “Deeptest: Automated testing of deep-neural-network-driven autonomous cars,” in _Proceedings of the 40th International Conference on Software Engineering_ , 2018, pp. 303–314.
  * [5] L. Li, W.-L. Huang, Y. Liu, N.-N. Zheng, and F.-Y. Wang, “Intelligence testing for autonomous vehicles: A new approach,” _IEEE Transactions on Intelligent Vehicles_ , vol. 1, no. 2, pp. 158–166, 2016.
  * [6] L. Li, Y.-L. Lin, N.-N. Zheng, F.-Y. Wang, Y. Liu, D. Cao, K. Wang, and W.-L. Huang, “Artificial intelligence test: A case study of intelligent vehicles,” _Artificial Intelligence Review_ , vol. 50, no. 3, pp. 441–465, 2018.
  * [7] L. Li, X. Wang, K. Wang, Y. Lin, J. Xin, L. Chen, L. Xu, B. Tian, Y. Ai, J. Wang _et al._ , “Parallel testing of vehicle intelligence via virtual-real interaction,” _Science Robotics_ , 2019.
  * [8] S. Feng, Y. Feng, X. Yan, S. Shen, S. Xu, and H. X. Liu, “Safety assessment of highly automated driving systems in test tracks: A new framework,” _Accident Analysis & Prevention_, vol. 144, p. 105664, 2020.
  * [9] S. Feng, Y. Feng, C. Yu, Y. Zhang, and H. X. Liu, “Testing scenario library generation for connected and automated vehicles, part i: Methodology,” _IEEE Transactions on Intelligent Transportation Systems_ , vol. 22, no. 3, pp. 1573–1582, 2021.
  * [10] S. Feng, Y. Feng, H. Sun, S. Bao, Y. Zhang, and H. X. Liu, “Testing scenario library generation for connected and automated vehicles, part ii: Case studies,” _IEEE Transactions on Intelligent Transportation Systems_ , vol. 22, no. 9, pp. 5635–5647, 2021.
  * [11] D. Zhao, H. Lam, H. Peng, S. Bao, D. J. LeBlanc, K. Nobukawa, and C. S. Pan, “Accelerated evaluation of automated vehicles safety in lane-change scenarios based on importance sampling techniques,” _IEEE Transactions on Intelligent Transportation Systems_ , vol. 18, no. 3, pp. 595–607, 2016.
  * [12] D. Zhao, X. Huang, H. Peng, H. Lam, and D. J. LeBlanc, “Accelerated evaluation of automated vehicles in car-following maneuvers,” _IEEE Transactions on Intelligent Transportation Systems_ , vol. 19, no. 3, pp. 733–744, 2017.
  * [13] G. E. Mullins, P. G. Stankiewicz, R. C. Hawthorne, and S. K. Gupta, “Adaptive generation of challenging scenarios for testing and evaluation of autonomous vehicles,” _Journal of Systems and Software_ , vol. 137, pp. 197–215, 2018\.
  * [14] M. Koren, S. Alsaif, R. Lee, and M. J. Kochenderfer, “Adaptive stress testing for autonomous vehicles,” in _2018 IEEE Intelligent Vehicles Symposium (IV)_. IEEE, 2018, pp. 1–7.
  * [15] J. Sun, H. Zhou, H. Xi, H. Zhang, and Y. Tian, “Adaptive design of experiments for safety evaluation of automated vehicles,” _IEEE Transactions on Intelligent Transportation Systems_ , 2021.
  * [16] S. Feng, Y. Feng, H. Sun, Y. Zhang, and H. X. Liu, “Testing scenario library generation for connected and automated vehicles: an adaptive framework,” _IEEE Transactions on Intelligent Transportation Systems_ , vol. 23, no. 2, pp. 1213–1222, 2022.
  * [17] A. B. Owen, “Monte carlo theory, methods and examples,” 2013.
  * [18] S. Feng, X. Yan, H. Sun, Y. Feng, and H. X. Liu, “Intelligent driving intelligence test for autonomous vehicles with naturalistic and adversarial environment,” _Nature Communications_ , vol. 12, no. 1, pp. 1–14, 2021.
  * [19] K. A. Marill, “Advanced statistics: linear regression, part ii: multiple linear regression,” _Academic Emergency Medicine_ , vol. 11, no. 1, pp. 94–102, 2004.
  * [20] D. J. Olive, “Multiple linear regression,” in _Linear Regression_. Springer, 2017, pp. 17–83.
  * [21] UNECE, “Uniform provisions concerning the approval of vehicles with regard to automated lane keeping systems,” 2021. [Online]. Available: https://unece.org/sites/default/files/2021-03/R157e.pdf
  * [22] A. Owen and Y. Zhou, “Safe and effective importance sampling,” _Journal of the American Statistical Association_ , vol. 95, no. 449, pp. 135–143, 2000\.
  * [23] J. W. Ro, P. S. Roop, A. Malik, and P. Ranjitkar, “A formal approach for modeling and simulation of human car-following behavior,” _IEEE Transactions on Intelligent Transportation Systems_ , vol. 19, no. 2, pp. 639–648, 2017.