# First-principles calculations of double resonance Raman spectra for
monolayer MoTe2

Jianqi Huang  Shenyang National Laboratory for Materials Science, Institute of
Metal Research, Chinese Academy of Sciences, School of Material Science and
Engineering, University of Science and Technology of China, Shenyang 110016,
P. R. China  Huaihong Guo<EMAIL_ADDRESS>College of Sciences, Liaoning
Petrochemical University, Fushun 113001, P. R. China  Lin Zhou  School of
Chemistry and Chemical Engineering, Frontiers Science Centre for
Transformative Molecules, Shanghai Jiao Tong University, Shanghai 200240, P.
R. China  Shishu Zhang  Center for Nanochemistry, Beijing Science and
Engineering Center for Nanocarbons, Beijing National Laboratory for Molecular
Sciences, College of Chemistry and Molecular Engineering, Peking University,
100871, Beijing, P. R. China  Lianming Tong  Center for Nanochemistry, Beijing
Science and Engineering Center for Nanocarbons, Beijing National Laboratory
for Molecular Sciences, College of Chemistry and Molecular Engineering, Peking
University, 100871, Beijing, P. R. China  Riichiro Saito  Department of
Physics, Tohoku University, Sendai, Miyagi 980-8578, Japan  Teng Yang
<EMAIL_ADDRESS>Shenyang National Laboratory for Materials Science,
Institute of Metal Research, Chinese Academy of Sciences, School of Material
Science and Engineering, University of Science and Technology of China,
Shenyang 110016, P. R. China  Zhidong Zhang  Shenyang National Laboratory for
Materials Science, Institute of Metal Research, Chinese Academy of Sciences,
School of Material Science and Engineering, University of Science and
Technology of China, Shenyang 110016, P. R. China

###### Abstract

Since double resonance Raman (DRR) spectra are laser-energy dependent, the
first-principles calculations of DRR for two-dimensional materials are
challenging. Here, the DRR spectrum of monolayer MoTe2 is calculated by home-
made program, in which we combine ab-initio density-functional-theory
calculations with the electron-phonon Wannier (EPW) method. Within the fourth-
order perturbation theory, we are able to quantify not only the electron-
photon matrix elements within the dipole approximation, but also the electron-
phonon matrix elements using the Wannier functions. The reasonable agreement
between the calculated and experimental Raman spectra is achieved, in which we
reproduce some distinctive features of transition metal dichalcogenides
(TMDCs) from graphene (for example, the dominant intervalley process involving
an electron or a hole). Furthermore, we perform an analysis of the possible
DRR modes over the Brillouin zone, highlighting the role of low-symmetry
points. Raman tensors for some DRR modes are given by first principles
calculations from which laser polarization dependence is obtained.

##  I introduction

Raman spectroscopy as a versatile probe tool has been widely used for
characterizing a broad range of physical properties including superconductive
[1], topological [2], ferroelectric [3] properties, magnetic ordering [4] and
phase transition [5], electronic interference effect [6, 7], phonon helicity
[8]. With the rise of two-dimensional (2D) van der Waals materials, Raman
spectroscopy plays an essential role in supplying information on the
heterostructure and intrinsically topological properties [9, 10, 11, 12]. In
particular, the second-order Raman spectra, which have been widely observed in
transition metal dichalcogenides (TMDCs) and other semiconducting 2D materials
[13, 14, 15, 16, 17, 18, 14, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28], host non
zone-centered two phonons, which are strongly dependent on the laser
excitation energy. Double resonance Raman (DRR) process is essential for
observing the two-phonon Raman spectra or defect-oriented Raman spectra, whose
intensity is comparable to or even larger than the resonant Raman process of a
zone-centered phonon [29, 30].

However, the first principles calculation of the DRR spectra has been a long-
term challenge before it is possible to accurately quantify electron-phonon
matrices with sufficiently dense grid under Wannier interpolation throughout
the Brillouin zone (BZ) [31]. With the developed electron-phonon Wannier
technique, Herziger et al. [32] and Torche et al. [33] have calculated the
double-resonant 2D mode of graphite by first-principles calculations. However,
compared to the 2D overtone band in graphite which is intensively studied and
has prior assignment due to empirical method, it remains challenging to
identify by first principles calculation the spectra of a system without prior
double-resonance information. It is still necessary to calculate the DRR
scattering amplitudes for all electron wavenumbers, $k$, for all phonon
wavevectors, $q$ with the given $k$, and for all combinations of two phonon
modes and electronic energy subbands as a function of laser excitation
energies.

So far, the DRR analysis in TMDC is largely semi-quantitative or even
qualitative, due to the lack of a quantitative treatment of electron-phonon
coupling. The earliest interpretations of DRR of TMDCs were merely by
comparing phonon frequencies with reference to some inelastic neutron
scattering results [13, 14], which therefore leads to some inconsistency in
mode assignments [15, 16, 17, 18, 14, 19, 20, 21, 22, 23]. Terrones et al.
[24] and Berkdemir et al. [25] assigned phonon modes in few-layer WSe2 by
referring to both the band structure and the high-symmetry point phonon. Guo
et al. [26] went further to specify the phonon wave vectors around the high-
symmetry M point for the DRR electron-phonon resonance in MoTe2. Carvalho et
al. [28] measured the second-order bands using more than twenty different
laser excitations in conjunction with DFT calculations. Livneh et al. [27]
made a comprehensive multiphonon spectral analysis in MoS2 based on group
theory. All the recent endeavors have been either focused on the high-symmetry
points in the reciprocal space, or treating the electron-phonon matrix element
as constant, which can hardly reveal the essential message on the interplay
between electron and phonon. Tatsumi et al. [8] made a computer program of
calculating first-order resonance Raman spectra by first principles
calculation in which electron-photon and electron-phonon matrix elements are
calculated by the wave function coefficient and Electron-Phonon-Wannier (EPW)
for zone-center phonon mode. Their concept can be extended for DDR spectra,
though a huge amount of calculations is required.

In this paper, we study theoretically the second-order Raman process of
monolayer MoTe2. We developed a computer program for calculating the DDR
spectra by first principles calculation in which we adopt Quantum-Espresso
(QE) and EPW for evaluating the electron-photon and electron-phonon matrix
elements, respectively, and compare with the experimental Raman spectra for
several excitation energies. The reasonable agreement between the calculated
and experimental double resonance Raman spectra is achieved, showing the
reliability of our method and revealing some distinctive features of TMDCs
from graphene. Furthermore, assignment of DRR modes on any possible
combinations of two-phonon modes at any random $k$ point is directly obtained.
Raman tensors for some DRR modes are obtained from which we can discuss
polarization dependence of two-phonon Raman spectra.

##  II methods

###  II.1 computational methods

Second-order double resonance Raman intensity with Raman shift $E_{RS}$ as a
function of the incident laser energy $E_{L}$ can be described as following
expression,

$\begin{split}I\left(E_{L}\right)\propto&\sum_{\textbf{q},\mu,\nu}\left|\bm{P}_{s}^{\dagger}\cdot\stackrel{{\scriptstyle\leftrightarrow}}{{R}}\left(\bm{q},\mu,\nu,E_{L}\right)\cdot\bm{P}_{i}\right|^{2}\\\
&\delta\left(E_{RS}\pm\hbar\omega_{\mu}\pm\hbar\omega_{\nu}\right),\end{split}$
(1)

where $\bm{P}_{i}$, $\bm{P}_{s}$ represent the polarization directions of
incident light and scattered light, respectively, $\omega_{\mu}$ and
$\omega_{\nu}$ denote the phonon frequencies corresponding to two phonon modes
$\mu$ and $\nu$. The Raman tensor takes a fourth-order perturbation form of

$\stackrel{{\scriptstyle\leftrightarrow}}{{R}}\left(\bm{q},\mu,\nu,E_{L}\right)=\sum_{\bm{k},i=f,n,n^{\prime},n^{\prime\prime}}\frac{\bm{D}_{fn^{\prime\prime}}\left(\bm{k}\right)\cdot
M_{n^{\prime\prime}n^{\prime}}^{ep,\nu_{-\bm{q}}}\left(\bm{k}+\bm{q}\right)\cdot
M_{n^{\prime}n}^{ep,\mu_{\bm{q}}}\left(\bm{k}\right)\cdot\bm{D}_{ni}^{\dagger}\left(\bm{k}\right)}{\left(E_{ni}-E_{L}-\mathrm{i}\gamma\right)\left(E_{n^{\prime}i}-E_{L}\pm\hbar\omega_{\mu}-\mathrm{i}\gamma\right)\left(E_{n^{\prime\prime}i}-E_{L}\pm\hbar\omega_{\mu}\pm\hbar\omega_{\nu}-\mathrm{i}\gamma\right)},$
(2)

where $\bm{D}_{fn^{\prime\prime}}$ is the electric dipole vector $\langle
f|\bm{D}|n^{\prime\prime}\rangle$ and
$M_{n^{\prime\prime}n^{\prime}}^{ep,\nu_{-\bm{q}}}$ is the electron-phonon
coupling matrix elements. $i$, $n$, $n^{\prime}$, $n^{\prime\prime}$ and $f$
denote, respectively, the initial state, the three intermediate states, and
the final state of an electron. Since the broadening due to electron-phonon
coupling is around 100 meV at room temperature, we adopt the value of damping
constant $\gamma$ in our simulation. Since the backscattering configuration is
set up in the experiment, the Raman spectra are calculated in the
backscattering configuration ($\bar{\textrm{Z}}$(XX)Z).

After obtaining the Raman spectrum, we take the following steps to assign the
combined modes that contribute to a specific Raman peak. With fixing the Raman
shift $E_{RS}$, we can obtain Raman intensity ($I_{\bm{q}}$) as a function of
$q$ in the Brillouin zone,

$\displaystyle I_{\bm{q}}$ $\displaystyle\propto$
$\displaystyle\sum_{\mu\nu}\left|\bm{P}_{s}^{\dagger}\cdot\stackrel{{\scriptstyle\leftrightarrow}}{{R}}\left(\bm{q},\mu,\nu\right)\cdot\bm{P}_{i}\right|^{2}\delta\left(E_{RS}\pm\hbar\omega_{\mu}\pm\hbar\omega_{\nu}\right)$
(3) $\displaystyle\equiv$ $\displaystyle\sum_{\mu,\nu}I^{\mu\nu}_{q}.$

Meanwhile, with fixing two vibration modes $\mu$ and $\nu$ at the same time,
we can obtain the Raman intensity ($I_{\bm{q}}^{\mu\nu}$) contributed by the
combined two phonon modes as a function of $q$. Considering the entire
Brillouin zone, the contribution of the fixed mode combination to total Raman
intensity is achieved by summation

$\zeta_{\mu\nu}=\sum_{\bm{q}}\frac{I_{\bm{q}}^{\mu\nu}}{I}.$ (4)

Now, we can extract the ones that contribute significantly to the total Raman
intensity by observing the value of $\zeta_{\mu\nu}$ of all possible mode
combinations.

We performed the electronic and phonon energy dispersion calculations on
monolayer MoTe2 by using first-principles density functional theory within the
local density approximation (LDA) as implemented in the QE code [34]. The
monolayer MoTe2 are separated by 25 Å from one another in a unit cell of the
calculation to eliminate the inter-few-layer interaction. We used norm-
conserving pseudopotentials (NCPP) within the local density approximation
(LDA) with a plane-wave cutoff energy of 120 Ry to describe the interaction
between electrons and ions. The spin-orbit split electronic band structures
were calculated by the relativistic pseudopotentials derived from an atomic
Dirac-like equation. The atomic coordinates were relaxed until the atomic
force was less than 10-5 Ry/Bohr. The Monkhorst-Pack scheme [35] was used to
sample the Brillouin zone over a 16$\times$16$\times$1 and 8$\times$8$\times$1
$k$-mesh for electronic and phonon energy calculation, respectively. The
phonon energy dispersion relations of MoTe2 were calculated by the density
functional perturbation theory [36]. On these basis, we calculated the
electrical dipole vector by using a modified version of the QE code. Further
by means of Wannier interpolation schemes as implemented in standard EPW [37,
38], we obtained the electron-phonon coupling matrix elements for each phonon
mode on a much fine grid of 45$\times$45$\times$1 $k$-mesh in the Brillouin
zone which is dense enough to achieve convincing results.

###  II.2 experimental method

Bulk crystals of MoTe2 were prepared through a chemical vapor transport method
[39]. Atomically thin crystals of MoTe2 were mechanically exfoliated from the
bulk crystals onto 90 nm SiO2/Si substrate. Raman spectroscopy for monolayer
MoTe2 was performed using 532, 633 and 785 nm excitation lasers for discussing
the observed phonon dispersion. The grating sizes were 1800 lines/mm for the
785, 633 and 532 nm laser excitation measurements. The magnification of the
objective lens was 100x. The accumulation times were 60-300 seconds. All
measurements were performed at room temperature in the backscattering
configuration. Typical Raman spectra of monolayer 2H-MoTe2 under different
laser excitation are analyzed and compared with theoretical calculations in
Fig. 1 below.

##  III results and discussion

Figure 1: (Color online) Double Resonance Raman (DRR) process of monolayer
MoTe2. (a) Raman spectra of monolayer MoTe2 under three different laser
excitation energies, 2.33 eV (532 nm), 1.96 eV (633 nm) and 1.58 eV (785 nm),
showing the small DRR intensity peaks Pi (i = 1, 2, …, 7). The experiment and
simulated data are given in top and bottom panels by solid and dash lines,
respectively. The typical first-order Raman peaks of A${}^{\prime}_{1}$
(corresponding to A${}_{\textrm{1g}}$ in bulk MoTe2) and E′
(E${}_{\textrm{2g}}$) are only given in the experimental data. (b)
Decomposition of total simulated Raman spectra into different Raman scattering
process pathways, including the $ee$, $hh$, $eh$ and $he$ processes. (c)
Schema of the four typical DRR pathways. $eh$ and $he$ are not allowed for the
forbidden spin-flip in hole-phonon scattering.

Fig. 1(a) shows Raman spectra of monolayer MoTe2 under 2.33 eV (532 nm), 1.96
eV (633 nm) and 1.58 eV (785 nm) laser excitation energies from both
experiment (top panel in solid lines) and simulation (bottom panel in dashed
lines). The baseline correction of the raw experimental Raman spectra has been
performed to remove the fluorescence noise. The experimental Raman spectra
show two strong peaks, the in-plane E′ mode at $\sim$ 236 cm-1 (corresponding
to E${}_{\textrm{2g}}$ in bulk) and the out-of-plane A${}^{\prime}_{1}$ mode
at $\sim$ 171 cm-1 (A${}_{\textrm{1g}}$ in bulk) for monolayer MoTe2. These
spectra are assigned to the first-order Raman spectra of the $\Gamma$ point
phonons [40]. The calibration peak from Si at $\sim$ 300 cm-1 is assigned to
the 2TA mode [41]. Since the present simulation method is only applicable for
double resonance Raman peaks, the simulated Raman spectra can distinguish the
DRR spectrum from the first-order peaks. As shown in Fig. 1(a), seven peaks
are observed with relatively small intensities, which we denote as Pi (i = 1,
2, …, 7). The peak positions of each Pi are found be dispersive as a function
of excitation energy [42], either upshifting or downshifting by up to several
cm-1 by changing laser excitation energies, and were ascribed to the second-
order Raman process [26]. A reasonable agreement between experiment and our
calculation substantiates the double resonance origin of these peaks, as
evidenced in Fig. 1(a).

Double resonance Raman process usually consists of several pathways [43].
Depending on whether valence hole is involved and whether the optical
absorption/emission occur at the same $k$ point, four typical pathways are
$ee$ (only conduction electron involved), $hh$ (only valence hole involved),
$eh$ and $he$ (electron and hole both involved), as indicated in Fig. 1(c).
Our method can manage to decompose the contribution from each pathway to the
Raman intensity quantitatively. In Fig. 1(b) taking the laser excitation 1.96
eV, for example, we compare the Raman spectra from each pathway with the total
one. The $ee$ and $hh$ seem to contribute almost equally and dominantly to the
total intensity, in contrast to the negligible contribution from the $eh$ and
$he$ processes. Such a behavior is vastly different from what’s observed in
graphene [43], where $eh$ and $he$ dominate over $hh$ and $ee$. This is
sensible considering that there doesn’t exist electron-hole symmetry in MoTe2,
as indicated in the electronic band structure in Fig. 2(a), which strongly
hinders both electron and hole scattering simultaneously by intervalley
phonon. While in graphene with symmetric Dirac cone, intervalley scattering of
both conduction electron and valence hole by phonon can take place in parallel
[43]. Further, the inversion-symmetry breaking and strong spin-orbit coupling
(SOC) in monolayer MoTe2 give rise to a large spin splitting on the electronic
bands around the $K$ point, as shown in Fig. 2(a) and spin-valley locking
occurs around the valley. Such special band features in MoTe2 make the Raman
transition of $eh$ and $he$ pathways unallowable because of the forbidding
spin-flip in the hole-phonon scattering in such cases.

Figure 2: (Color online) (a) Electronic energy band structure and density of
states (DOS) of MoTe2. The vertical bar in black, red and blue represent the
position of electronic photo-excitation due to laser excitation at energy of
1.58 eV, 1.96 eV and 2.33 eV, respectively. (b) The phonon dispersion relation
of MoTe2. The point groups are listed above for different $q$ points with the
irreducible representations given on each colored band. Typical Raman-active
$\Gamma$-point phonons are visualized with their Raman tensor given.

In order to discuss the Raman active modes for two-phonon Raman scattering, we
first make analysis on group theoretical selection rules for two-phonon
scattering. Due to the dominant $ee$ and $hh$ processes for the DRR, the
$\mu_{th}$ and $\nu_{th}$ two phonons contributing to a single DRR process
should have the same wave vector with opposite sign to each other, qμ and -qν.
The origin of DRR which can arise from either combination, or subtraction, or
overtone of two phonon modes can therefore be analyzed at the same $q$ point
in the phonon dispersion relation in Fig. 2(b). The point groups along high-
symmetry line and $q$ points are put on top, such as point group
D${}_{\textrm{3h}}$ for the highest symmetric zone-center $\Gamma$ point,
C${}_{\textrm{3h}}$ for the zone corner $K$ point, C${}_{\textrm{2v}}$ for the
zone-edge middle $M$ point and between the $\Gamma$ and $M$ points, and
C${}_{\textrm{s}}$ for the rest of $q$ points (including between the $\Gamma$
(or $M$) and $K$ points) in the BZ. The irreducible representations are also
given to the corresponding phonon bands. Three typical zone-center Raman-
active modes (A${}^{\prime}_{1}$, E′ and E′′) are visualized on the right of
Fig. 2(b) with the Raman tensors. By the general methods of Birman [44, 45],
with the help of the Bilbao crystallographic server [46], we can determine
whether two-phonon modes for a $q$ point are Raman-active or not, by
correlating the irreducible representations of the combined modes of the group
at the $q$ point with the irreducible representations of the full space group
(D${}_{\textrm{3h}}$) at the $\Gamma$ point. To ascertain the Raman activity
of a two-phonon DRR mode, the reduced irreducible constituents from the
Kronecker products have to contain at least one of the three Raman-active
symmetries (A${}^{\prime}_{1}$, E′ and E′′) of zone-center point group
D${}_{\textrm{3h}}$. In the Appendix Table 2, we list all the possible
combinations of two phonons in the BZ and the reduction of symmetries, and
mark those Raman-active ones in blue color. First of all, the overtone modes,
even from two Raman-inactive symmetries, are all Raman active in the whole BZ,
since the decomposed symmetries include at least Raman-active
A${}^{\prime}_{1}$ symmetry. Second, for $q$ point beyond the zone-center
$\Gamma$ point (at $M$ or between $\Gamma$ and $M$ of C${}_{\textrm{2v}}$
group), especially with lower symmetry (C${}_{\textrm{s}}$ group), the
decomposed symmetries usually have two Raman-active symmetries
(A${}^{\prime}_{1}$ \+ E′) coexisting. Considering that A${}^{\prime}_{1}$ and
E′ have different circular polarization selectivity [8], we may anticipate
that DRR modes at these points should have non-zero Raman intensity in both
$\sigma_{+}\sigma_{+}$ and $\sigma_{+}\sigma_{-}$ configurations, which can be
used to test some possible two-phonon assignments.

Figure 3: (Color online) The assignment of the P3 DRR band. (a) Raman
intensity as a function of $q$ in the whole BZ. The main two-phonon
combinations are 2LA (overtone) and E′ \- TA (subtraction) whose Raman
intensities as a function of $q$ are given in (b) and (c), respectively. (d)
The percentage of the 2LA and E′ \- TA modes.

Next, let us assign each DRR mode. In Fig. 3, we take the P3 mode at laser
excitation energy 1.96 eV as an example, which was previously assigned to be
2LA(M) overtone mode [26]. Basically one specific Raman shift (205 cm-1 here)
with peak intensity is chosen, then we calculate the scattering cross section
as a function of $q$, Iq (see Eq. 3) from all possible two-phonon modes
including overtone, combination, subtraction ones. The Raman intensity Iq is
illustrated in Fig. 3(a), and the two dominating DRR modes are derived to be
the 2LA (overtone) and E′ \- TA (subtraction) as indicated on the right panel
of Fig. 3(a). The percentage $\zeta_{\mu\nu}$ (See Eq. 4) in Fig. 3(d) shows
that the 2LA dominates over the E′ \- TA. Such an assignment is consistent
with the previous group theory analysis made in the Appendix Table 2.

Table 1: A comparison of assignment for second-order Raman modes Pi (i = 1, 2,
…, 6) between the previous work [26] and the current work for EL = 1.96 eV.
The label (new) denotes new assignments by the present work. The percentage of
Raman intensity for each DRR mode is given.

mode |  Ref.[26, 47] |  this work |  $\zeta_{\mu\nu}$
---|---|---|---
P1 |  2TA(M) |  2TA(M) |  9.2%
|  E′(M) - LA(M) |  E′ \- LA |  84.9%
P2 |  E′(M) - TA(M) |  E′ \- TA |  18.4%
|  |  LA + TA (new) |  55.2%
|  |  2LA (new) |  13.9%
|  |  E′ \- LA (new) |  7.6%
P3 |  E′′(M) + TA(M) |  NA in $\bar{\textrm{Z}}$(XX)Z |
|  2LA(M) |  2LA(M) |  93.5%
|  |  E′ \- TA (new) |  3.1%
P4 |  A${}^{\prime}_{1}$(M) + LA(M) |  A${}^{\prime}_{1}$ \+ LA |  80.0%
|  |  E′ \+ LA (new) |  7.7%
P5 |  E′(M) + TA(M) |  E′ \+ TA |  27.4%
|  |  E′ \+ LA (new) |  58.2%
P6 |  E′(M) + LA(M) |  E′ \+ LA |  98.7%

Further, not only the mode assignment (2LA), but also the most probable phonon
vectors $q$ can be obtained from the calculated Raman spectra as shown in Fig.
3(a), which is distributed around the $M$ points with a ”8” shape, and agrees
reasonably with the results derived from equi-energy contour lines of band
structure by Guo et al. (see Fig. 2c in Ref. [26]. However, Guo et al. [26] in
theory and Caramazza et al. [47] in experiment also assigned the P3 band
partly to E′′(M) + TA(M), which is absent at M from our current analysis,
merely because we have used the backscattering configuration
($\bar{\textrm{Z}}$(XX)Z) in our calculation, in which E′′ is not Raman-
active.

The data of the remaining DRR modes Pi (i = 1, 2, 4-6) is given in Table 1 and
also in the Appendix. The mode assignments are compared directly with the
previous work by Guo et al. [26] and Caramazza et al. [47], as listed in Table
1. From Table 1, we confirm that most assignments in this work for all the six
DRR modes are consistent with that made by Guo et al. [26] and Caramazza et
al. [47]. Nevertheless, we also have non-zero $\zeta$ for other combination
modes which appear in the other regions of the $q$ space. One of the
advantages of the present calculation is to be able to supply the fraction of
each assignment when multiple scattering channels are possible, for example,
E′ \- LA ($\zeta_{\mu\nu}\approx$85%) is much more important than 2TA
($\zeta_{\mu\nu}\approx$9%) for P1, both P3 and P6 modes have one dominating
assignment (2LA and E′ \+ LA, respectively), all of which can not be extracted
merely from group theoretical analysis. Moreover, some additional assignments
can be unveiled from this work, for example, for P2, besides E′ \- TA, it has
LA + TA, 2LA, and E′ \- LA, in which LA + TA ($\zeta_{\mu\nu}\approx$55%) is
even more important than E′ \- TA ($\zeta_{\mu\nu}\approx$18%).

Figure 4: (Color online) The dispersion of DRR modes on laser energy EL. The
EL dependence of (a) Raman profile, (b) Raman shift, and (c) $q$-resolved
Raman intensity of the P3 and P6 bands.

What is more, as seen from the Raman intensity as a function of $q$ and laser
excitation energy EL in Fig. 4(c) and also in the Appendix Fig. 6-10, Raman
intensity can arise from general $q$ points with low symmetries, the $q$
points which contribute most to total Raman intensity appear either at or
between the high-symmetry points in BZ, which justifies the previous
theoretical treatment only on the high-symmetry points [15, 16, 17, 18, 14,
19, 20, 21, 22, 23, 24, 28, 27, 26, 47]. However, there is one DRR peak at
around 300 cm-1 which has not been explored before due to its overlap with the
Raman peak of Silicon. If looking more carefully into the experimental Raman
spectra in Fig. 1(a), one can see a satellite peak on the right of the Silicon
peak. This peak, which we designate as P7 band, is more visible from our
simulated Raman spectra, and can be assigned to be an overtone mode, namely
2E′′, at $q$ close to the $K$ point along $\Gamma$ and $K$ points, as analyzed
in Appendix Fig. 11.

Since the second-order Raman spectra are dispersive as a function of
excitation energy EL, we further substantiate our mode assignment by showing
the dispersion of two typical DRR bands P3 and P6 with EL in Fig. 4. From the
EL dependent band profile (Fig. 4(a)) and especially peak shift (with respect
to the frequency at 1.96 eV) (Fig. 4(b)), we can see the dispersive features
in the P3 and P6 DRR bands, both of which are blue-shifted with increasing EL,
with $\Delta\omega$/$\Delta$EL $\sim$ 25.6 and 3.33 cm-1/eV for P3 and P6,
respectively. The simulation results (in open squares and circles) agree
reasonably well with the experimental data (in solid squares and circles). The
dispersion difference between the P3 and P6 bands arises mainly from EL
dependent phonon wave vector $q$. As shown in Fig. 4(c), with increasing EL
from 1.96 eV to 2.33 eV, the $q$ changes from between $M$ and $K$, via $M$
point, to between $M$ and $\Gamma$ points for the P3 band. In contrast, the
$q$ for the P6 band keeps almost unchanged with increasing EL. The P3 band
dispersion can also be seen from the phonon dispersion relation in Fig. 2(b),
the band segment of the LA band around the $M$ point that is assigned to P3
(Table 1) is rather dispersive.

Figure 5: (Color online) The polarization selectivity of the DRR modes at EL =
1.96 eV. (a) The linearly polarized Raman spectra in the parallel
($\bar{\textrm{Z}}$(XX)Z) and perpendicular ($\bar{\textrm{Z}}$(XX)Z)
geometries in black and red, respectively. The intensity of the
($\bar{\textrm{Z}}$(XX)Z) configuration is multiplied by 5. The P4 band is
enlarged in inset. (b) The polar figure of Raman intensity for both co-
polarized and depolarized set-ups.

Finally, let us show the laser polarization dependence of the DRR modes. From
Table 1, we can see that most of the DRR modes have E′ involved, except for
the P4 band. The P4 band has two most probable assignments A${}^{\prime}_{1}$
\+ LA and E′ \+ LA, but the former of which has a much bigger fraction
($\zeta_{\mu\nu}\approx$80%) than the latter one
($\zeta_{\mu\nu}\approx$7.7%), suggesting a potentially strong linear
polarization dependence. This is indeed the case, as indicated in Fig. 5(a).
Figure 5(a) gives the calculated Raman spectra in both the parallel
($\bar{\textrm{Z}}$(XX)Z) and perpendicular ($\bar{\textrm{Z}}$(XX)Z)
geometries. The intensities in the two geometries are nonzero and almost equal
for all the DRR modes discussed here, except for the P4 band, which has zero
intensity in the perpendicular geometries. We also explore the angle
dependence of Raman intensity of the P4 band. The co-polarized and depolarized
geometries are both used and set-ups are schematized in Fig. 5(b). The polar
plot of both geometries is consistent with the above analysis. The isotropic
polar data are seen for the co-polarization geometry and comes from the in-
plane isotropy of monolayer MoTe2, while the depolarized geometry that the
incident and scattered laser polarizations have a relative angle $\theta$
gives a strong anisotropy ($\approx$ cos${}^{2}\theta$) arising from
A${}^{\prime}_{1}$, consistent with the group theoretical analysis in the
Appendix Table 2.

##  IV conclusion

In summary, we have calculated the second-order Raman spectra of MoTe2
monolayer, based on first-principles density functional calculation and time-
dependent perturbation theory. The non-empirical treatment of electron-phonon
interaction is performed, which allows us to quantify the contribution of all
possible two-phonon combinations to the double resonance Raman modes in every
single phonon wave vector. The Pi (i = 1, 2, …, 7) band assignments, which are
consistent with the DRR selective rule constructed based on group theory, show
some additional origin of two phonon modes which was not found in the previous
studies. The polarization dependence of the DRR modes are also investigated,
which should be observed experimentally. This study facilitates a deeper
understanding of the electron-phonon interaction and the second-order Raman
process in TMDC systems.

###### Acknowledgements.

We sincerely acknowledge Prof. Vitto Zheng Han for helping prepare the high-
quality figures. This project is supported by the National Natural Science
Foundation of China (52031014, 51702146) and the National Key R&D Program of
China (2017YFA0206301). H.G. acknowledges Department of Education of Liaoning
Province Grants No. LJKZ0391. L.Z. acknowledges the National Key Basic
Research Program of China (2021YFA1401400) and the National Natural Science
Foundation of China (52103344). L.T. acknowledges the National Natural Science
Foundation of China (21974004). R.S. acknowledges MEXT Grants No. JP18H01810.
The simulation work was carried out at National Supercomputer Center in
Tianjin, China, and the calculations were performed on TianHe-1(A).

##  Appendix A Raman selection rule analysis of two phonons

Table 2: Raman selection rule analysis of two phonons of MoTe2 with different
symmetries based on group theory. The light blue and brown color mark the
Raman active and inactive phonon combinations, respectively.

Product | Reduction | Product | Reduction | Product | Reduction
---|---|---|---|---|---
$\textrm{D}_{\textrm{3h}}$ | $\textrm{D}_{\textrm{3h}}$ | $\textrm{C}_{\textrm{3h}}$ | $\textrm{D}_{\textrm{3h}}$ | $\textrm{C}_{\textrm{3h}}$ | $\textrm{C}_{\textrm{2v}}$ | $\textrm{D}_{\textrm{3h}}$ | $\textrm{C}_{\textrm{2v}}$
$\textrm{A}^{\prime}_{1}\times\textrm{A}^{\prime}_{1}$ | $\textrm{A}^{\prime}_{1}$ | $\textrm{A}^{\prime}\times\textrm{A}^{\prime}$ | $\textrm{A}^{\prime}_{1}+\textrm{A}^{\prime}_{2}$ | $\textrm{A}^{\prime}$ | $\textrm{A}_{1}\times\textrm{A}_{1}$ | $\textrm{A}^{\prime}_{1}+\textrm{E}^{\prime}$ | $\textrm{A}_{1}+\textrm{B}_{1}$
$\textrm{A}^{\prime}_{1}\times\textrm{A}^{\prime}_{2}$ | $\textrm{A}^{\prime}_{2}$ | $\textrm{A}^{\prime}\times\textrm{A}^{\prime\prime}$ | $\textrm{A}^{\prime\prime}_{1}+\textrm{A}^{\prime\prime}_{2}$ | $\textrm{A}^{\prime\prime}$ | $\textrm{A}_{1}\times\textrm{A}_{2}$ | $\textrm{A}^{\prime\prime}_{1}+\textrm{E}^{\prime\prime}$ | $\textrm{A}_{2}+\textrm{B}_{2}$
$\textrm{A}^{\prime}_{1}\times\textrm{A}^{\prime\prime}_{1}$ | $\textrm{A}^{\prime\prime}_{1}$ | $\textrm{A}^{\prime}\times{}^{1}\textrm{E}^{\prime}$ | $\textrm{E}^{\prime}$ | ${}^{1}\textrm{E}^{\prime}$ | $\textrm{A}_{1}\times\textrm{B}_{1}$ | $\textrm{A}^{\prime}_{2}+\textrm{E}^{\prime}$ | $\textrm{A}_{1}+\textrm{B}_{1}$
$\textrm{A}^{\prime}_{1}\times\textrm{A}^{\prime\prime}_{2}$ | $\textrm{A}^{\prime\prime}_{2}$ | $\textrm{A}^{\prime}\times{}^{2}\textrm{E}^{\prime}$ | $\textrm{E}^{\prime}$ | ${}^{2}\textrm{E}^{\prime}$ | $\textrm{A}_{1}\times\textrm{B}_{2}$ | $\textrm{A}^{\prime\prime}_{2}+\textrm{E}^{\prime\prime}$ | $\textrm{A}_{2}+\textrm{B}_{2}$
$\textrm{A}^{\prime}_{1}\times\textrm{E}^{\prime}$ | $\textrm{E}^{\prime}$ | $\textrm{A}^{\prime}\times{}^{1}\textrm{E}^{\prime\prime}$ | $\textrm{E}^{\prime\prime}$ | ${}^{1}\textrm{E}^{\prime\prime}$ | $\textrm{A}_{2}\times\textrm{A}_{2}$ | $\textrm{A}^{\prime}_{1}+\textrm{E}^{\prime}$ | $\textrm{A}_{1}+\textrm{B}_{1}$
$\textrm{A}^{\prime}_{1}\times\textrm{E}^{\prime\prime}$ | $\textrm{E}^{\prime\prime}$ | $\textrm{A}^{\prime}\times{}^{2}\textrm{E}^{\prime\prime}$ | $\textrm{E}^{\prime\prime}$ | ${}^{2}\textrm{E}^{\prime\prime}$ | $\textrm{A}_{2}\times\textrm{B}_{1}$ | $\textrm{A}^{\prime\prime}_{2}+\textrm{E}^{\prime\prime}$ | $\textrm{A}_{2}+\textrm{B}_{2}$
$\textrm{A}^{\prime}_{2}\times\textrm{A}^{\prime}_{2}$ | $\textrm{A}^{\prime}_{1}$ | $\textrm{A}^{\prime\prime}\times\textrm{A}^{\prime\prime}$ | $\textrm{A}^{\prime}_{1}+\textrm{A}^{\prime}_{2}$ | $\textrm{A}^{\prime}$ | $\textrm{A}_{2}\times\textrm{B}_{2}$ | $\textrm{A}^{\prime}_{2}+\textrm{E}^{\prime}$ | $\textrm{A}_{1}+\textrm{B}_{1}$
$\textrm{A}^{\prime}_{2}\times\textrm{A}^{\prime\prime}_{1}$ | $\textrm{A}^{\prime\prime}_{2}$ | $\textrm{A}^{\prime\prime}\times{}^{1}\textrm{E}^{\prime}$ | $\textrm{E}^{\prime\prime}$ | ${}^{1}\textrm{E}^{\prime\prime}$ | $\textrm{B}_{1}\times\textrm{B}_{1}$ | $\textrm{A}^{\prime}_{1}+\textrm{E}^{\prime}$ | $\textrm{A}_{1}+\textrm{B}_{1}$
$\textrm{A}^{\prime}_{2}\times\textrm{A}^{\prime\prime}_{2}$ | $\textrm{A}^{\prime\prime}_{1}$ | $\textrm{A}^{\prime\prime}\times{}^{2}\textrm{E}^{\prime}$ | $\textrm{E}^{\prime\prime}$ | ${}^{2}\textrm{E}^{\prime\prime}$ | $\textrm{B}_{1}\times\textrm{B}_{2}$ | $\textrm{A}^{\prime\prime}_{1}+\textrm{E}^{\prime\prime}$ | $\textrm{A}_{2}+\textrm{B}_{2}$
$\textrm{A}^{\prime}_{2}\times\textrm{E}^{\prime}$ | $\textrm{E}^{\prime}$ | $\textrm{A}^{\prime\prime}\times{}^{1}\textrm{E}^{\prime\prime}$ | $\textrm{E}^{\prime}$ | ${}^{1}\textrm{E}^{\prime}$ | $\textrm{B}_{2}\times\textrm{B}_{2}$ | $\textrm{A}^{\prime}_{1}+\textrm{E}^{\prime}$ | $\textrm{A}_{1}+\textrm{B}_{1}$
$\textrm{A}^{\prime}_{2}\times\textrm{E}^{\prime\prime}$ | $\textrm{E}^{\prime\prime}$ | $\textrm{A}^{\prime\prime}\times{}^{2}\textrm{E}^{\prime\prime}$ | $\textrm{E}^{\prime}$ | ${}^{2}\textrm{E}^{\prime}$ | Product | Reduction
$\textrm{A}^{\prime\prime}_{1}\times\textrm{A}^{\prime\prime}_{1}$ | $\textrm{A}^{\prime}_{1}$ | ${}^{1}\textrm{E}^{\prime}\times{}^{1}\textrm{E}^{\prime}$ | $\textrm{A}^{\prime}_{1}+\textrm{A}^{\prime}_{2}$ | ${}^{2}\textrm{E}^{\prime}$ |  | $\textrm{D}_{\textrm{3h}}$ | $\textrm{C}_{\textrm{3h}}$
$\textrm{A}^{\prime\prime}_{1}\times\textrm{A}^{\prime\prime}_{2}$ | $\textrm{A}^{\prime}_{2}$ | ${}^{1}\textrm{E}^{\prime}\times{}^{2}\textrm{E}^{\prime}$ | $\textrm{E}^{\prime}$ | $\textrm{A}^{\prime}$ | $\textrm{C}_{\textrm{s}}$ | $\textrm{C}_{\textrm{2v}}$ | $\textrm{C}_{\textrm{s}}$
$\textrm{A}^{\prime\prime}_{1}\times\textrm{E}^{\prime}$ | $\textrm{E}^{\prime\prime}$ | ${}^{1}\textrm{E}^{\prime}\times{}^{1}\textrm{E}^{\prime\prime}$ | $\textrm{A}^{\prime\prime}_{1}+\textrm{A}^{\prime\prime}_{2}$ | ${}^{2}\textrm{E}^{\prime\prime}$ |  | $\textrm{A}^{\prime}_{1}+\textrm{A}^{\prime}_{2}+2\textrm{E}^{\prime}$ | $\textrm{A}^{\prime}+{}^{1}\textrm{E}^{\prime}+{}^{2}\textrm{E}^{\prime}$
$\textrm{A}^{\prime\prime}_{1}\times\textrm{E}^{\prime\prime}$ | $\textrm{E}^{\prime}$ | ${}^{1}\textrm{E}^{\prime}\times{}^{2}\textrm{E}^{\prime\prime}$ | $\textrm{E}^{\prime\prime}$ | $\textrm{A}^{\prime\prime}$ | $\textrm{A}^{\prime}\times\textrm{A}^{\prime}$ | $2\textrm{A}_{1}+2\textrm{B}_{1}$ | $2\textrm{A}^{\prime}$
$\textrm{A}^{\prime\prime}_{2}\times\textrm{A}^{\prime\prime}_{2}$ | $\textrm{A}^{\prime}_{1}$ | ${}^{2}\textrm{E}^{\prime}\times{}^{2}\textrm{E}^{\prime}$ | $\textrm{A}^{\prime}_{1}+\textrm{A}^{\prime}_{2}$ | ${}^{1}\textrm{E}^{\prime}$ |  | $\textrm{A}^{\prime\prime}_{1}+\textrm{A}^{\prime\prime}_{2}+2\textrm{E}^{\prime\prime}$ | $\textrm{A}^{\prime\prime}+{}^{1}\textrm{E}^{\prime\prime}+{}^{2}\textrm{E}^{\prime\prime}$
$\textrm{A}^{\prime\prime}_{2}\times\textrm{E}^{\prime}$ | $\textrm{E}^{\prime\prime}$ | ${}^{2}\textrm{E}^{\prime}\times{}^{1}\textrm{E}^{\prime\prime}$ | $\textrm{E}^{\prime\prime}$ | $\textrm{A}^{\prime\prime}$ | $\textrm{A}^{\prime}\times\textrm{A}^{\prime\prime}$ | $2\textrm{A}_{2}+2\textrm{B}_{2}$ | $2\textrm{A}^{\prime\prime}$
$\textrm{A}^{\prime\prime}_{2}\times\textrm{E}^{\prime\prime}$ | $\textrm{E}^{\prime}$ | ${}^{2}\textrm{E}^{\prime}\times{}^{2}\textrm{E}^{\prime\prime}$ | $\textrm{A}^{\prime\prime}_{1}+\textrm{A}^{\prime}_{2}$ | ${}^{1}\textrm{E}^{\prime\prime}$ |  | $\textrm{A}^{\prime}_{1}+\textrm{A}^{\prime}_{2}+2\textrm{E}^{\prime}$ | $\textrm{A}^{\prime}+{}^{1}\textrm{E}^{\prime}+{}^{2}\textrm{E}^{\prime}$
$\textrm{E}^{\prime}\times\textrm{E}^{\prime}$ | $\textrm{A}^{\prime}_{1}+\textrm{A}^{\prime}_{2}+\textrm{E}^{\prime}$ | ${}^{1}\textrm{E}^{\prime\prime}\times{}^{1}\textrm{E}^{\prime\prime}$ | $\textrm{A}^{\prime}_{1}+\textrm{A}^{\prime}_{2}$ | ${}^{2}\textrm{E}^{\prime}$ | $\textrm{A}^{\prime\prime}\times\textrm{A}^{\prime\prime}$ | $2\textrm{A}_{1}+2\textrm{B}_{1}$ | $2\textrm{A}^{\prime}$
$\textrm{E}^{\prime}\times\textrm{E}^{\prime\prime}$ | $\textrm{A}^{\prime\prime}_{1}+\textrm{A}^{\prime\prime}_{2}+\textrm{E}^{\prime\prime}$ | ${}^{1}\textrm{E}^{\prime\prime}\times{}^{2}\textrm{E}^{\prime\prime}$ | $\textrm{E}^{\prime}$ | $\textrm{A}^{\prime}$ | Raman active
$\textrm{E}^{\prime\prime}\times\textrm{E}^{\prime\prime}$ | $\textrm{A}^{\prime}_{1}+\textrm{A}^{\prime}_{2}+\textrm{E}^{\prime}$ | ${}^{2}\textrm{E}^{\prime\prime}\times{}^{2}\textrm{E}^{\prime\prime}$ | $\textrm{A}^{\prime}_{1}+\textrm{A}^{\prime}_{2}$ | ${}^{1}\textrm{E}^{\prime}$ | Raman inactive

##  Appendix B The double resonance Raman modes P1-P2, P4-P7

Figure 6: (Color online) The P1 double resonance Raman mode. (a) $q$-resolved
Raman intensity in the whole BZ. The main two-phonon combinations are E′ \- LA
(subtraction) and 2TA (overtone) whose $q$-dependent Raman intensities are
given in (b) and (c), respectively. (d) The percentage of the E′ \- LA and 2TA
modes. Figure 7: (Color online) The P2 double resonance Raman mode. (a)
$q$-resolved Raman intensity in the whole BZ. The main two-phonon combinations
are LA + TA, E′ \- TA, 2LA, E′ \- LA whose $q$-dependent Raman intensities are
given in (b-e), respectively. (f) The percentage of all the assigned modes.
Figure 8: (Color online) The P4 double resonance Raman mode. (a) $q$-resolved
Raman intensity in the whole BZ. The main two-phonon combinations are
A${}^{\prime}_{1}$ \+ LA and E′ \+ LA whose $q$-dependent Raman intensities
are given in (b) and (c), respectively. (d) The percentage of the
A${}^{\prime}_{1}$ \+ LA and E′ \+ LA modes. Figure 9: (Color online) The P5
double resonance Raman mode. (a) $q$-resolved Raman intensity in the whole BZ.
The main two-phonon combinations are E′ \+ LA and E′ \+ TA whose $q$-dependent
Raman intensities are given in (b-d), respectively. (e) The percentage of the
E′ \+ LA and E′ \+ TA modes. Figure 10: (Color online) The P6 double resonance
Raman mode. (a) $q$-resolved Raman intensity in the whole BZ. The dominating
two-phonon combination is from E′ \+ LA whose $q$-dependent Raman intensities
is given in (b). (c) The percentage of the E′ \+ LA. Figure 11: (Color online)
The P7 double resonance Raman mode. The dominating combination is from 2E′′
whose $q$ is distributed close the $K$ point between $\Gamma$ and $K$ points.

## References

  * Measson  _et al._ [2014] M. A. Measson, Y. Gallais, M. Cazayous, B. Clair, P. Rodiere, L. Cario, and A. Sacuto, Phys. Rev. B 89, 060503(R) (2014).
  * Kung  _et al._ [2017] H.-H. Kung, S. Maiti, X. Wang, S.-W. Cheong, D. L. Maslov, and G. Blumberg, Phys. Rev. Lett. 119, 136802 (2017).
  * Tenne  _et al._ [2006] D. A. Tenne, A. Bruchhausen, N. D. Lanzillotti-Kimura, A. Fainstein, R. S. Katiyar, A. Cantarero, A. Soukiassian, V. Vaithyanathan, J. H. Haeni, W. Tian, D. G. Schlom, K. J. Choi, D. M. Kim, C. B. Eom, H. P. Sun, X. Q. Pan, Y. L. Li, L. Q. Chen, Q. X. Jia, S. M. Nakhmanson, K. M. Rabe, and X. X. Xi, Science 313, 1614 (2006).
  * Huang  _et al._ [2020] B. Huang, J. Cenker, X. Zhang, E. L. Ray, T. Song, T. Taniguchi, K. Watanabe, M. A. McGuire, D. Xiao, and X. Xu, Nat. Nanotechnol. 15, 212 (2020).
  * Kim  _et al._ [2019] K. Kim, S. Y. Lim, J.-U. Lee, S. Lee, T. Y. Kim, K. Park, G. S. Jeon, C.-H. Park, J.-G. Park, and H. Cheong, Nat. Commun. 10, 345 (2019).
  * Miranda  _et al._ [2017] H. P. C. Miranda, S. Reichardt, G. Froehlicher, A. Molina-Sanchez, S. Berciaud, and L. Wirtz, Nano Lett. 17, 2381 (2017).
  * Zhang  _et al._ [2022] S. Zhang, J. Huang, Y. Yu, S. Wang, T. Yang, Z. Zhang, L. Tong, and J. Zhang, Nature Communications 13, 1254 (2022).
  * Tatsumi and Saito [2018] Y. Tatsumi and R. Saito, Phys. Rev. B 97, 115407 (2018).
  * Chen  _et al._ [2016] S.-Y. Chen, T. Goldstein, D. Venkataraman, A. Ramasubramaniam, and J. Yan, Nano Letters 16, 5852 (2016).
  * Zhang  _et al._ [2016] K. Zhang, C. Bao, Q. Gu, X. Ren, H. Zhang, K. Deng, Y. Wu, Y. Li, J. Feng, and S. Zhou, Nat. Commun. 7, 13552 (2016).
  * Liu  _et al._ [2018] Y. Liu, Q. Gu, Y. Peng, S. Qi, N. Zhang, Y. Zhang, X. Ma, R. Zhu, L. Tong, J. Feng, Z. Liu, and J.-H. Chen, Adv. Mater. 30, 1706402 (2018).
  * Zhang  _et al._ [2020] K. Zhang, X. Pang, T. Wang, F. Han, S.-L. Shang, N. T. Hung, A. R. T. Nugraha, Z.-K. Liu, M. Li, R. Saito, and S. Huang, Phys. Rev. B 101, 014308 (2020).
  * Wakabayashi  _et al._ [1975] N. Wakabayashi, H. G. Smith, and R. M. Nicklow, Phys. Rev. B 12, 659 (1975).
  * Sourisseau  _et al._ [1991] C. Sourisseau, F. Cruege, M. Fouassier, and M. Alba, Chem. Phys. 150, 281 (1991).
  * Sekine  _et al._ [1984] T. Sekine, K. Uchinokura, T. Nakashizu, E. Matsuura, and R. Yoshizaki, J. Phys. Soc. Jpn. 53, 811 (1984).
  * Stacy and Hodul [1985] A. M. Stacy and D. T. Hodul, J. Phys. Chem. Sol. 46, 405 (1985).
  * Chen and Wang [1974] J. M. Chen and C. S. Wang, Solid State Commun. 14, 857 (1974).
  * Sourisseau  _et al._ [1989] C. Sourisseau, M. Fouassier, M. Alba, A. Ghorayeb, and O. Gorochov, Mater. Sci. Eng.: B 3, 119 (1989).
  * Feldman  _et al._ [1996] Y. Feldman, G. L. Frey, M. Homyonfer, V. Lyakhovitskaya, L. Margulis, H. Cohen, G. Hodes, J. L. Hutchison, and R. Tenne, J. Am. Chem. Soc. 118, 5362 (1996).
  * Frey  _et al._ [1999] G. L. Frey, R. Tenne, M. J. Matthews, M. S. Dresselhaus, and G. Dresselhaus, Phys. Rev. B 60, 2883 (1999).
  * Windom  _et al._ [2011] B. C. Windom, W. Sawyer, and D. W. Hahn, Tribol. Lett. 42, 301 (2011).
  * Li  _et al._ [2012] H. Li, Q. Zhang, C. C. R. Yap, B. K. Tay, T. H. T. Edwin, A. Olivier, and D. Baillargeat, Adv. Funct. Mater. 22, 1385 (2012).
  * Chakraborty  _et al._ [2013] B. Chakraborty, H. S. S. R. Matte, A. K. Sood, and C. N. R. Rao, J. Raman Spectrosc. 44, 92 (2013).
  * Terrones  _et al._ [2014] H. Terrones, E. Del Corro, S. Feng, J. M. Poumirol, D. Rhodes, D. Smirnov, N. R. Pradhan, Z. Lin, M. A. Nguyen, A. L. Elías, T. E. Mallouk, L. Balicas, M. A. Pimenta, and M. Terrones, Sci. Rep. 4, 4215 (2014).
  * Berkdemir  _et al._ [2013] A. Berkdemir, H. R. Guti$\acute{e}$rrez, A. R. Botello-M$\acute{e}$ndez, N. Perea-L$\acute{o}$pez, A. L. Elías, C.-I. Chia, B. Wang, V. H. Crespi, F. L$\acute{o}$pez-Urías, J.-C. Charlier, H. Terrones, and M. Terrones, Sci. Rep. 3, 1755 (2013).
  * Guo  _et al._ [2015] H. Guo, T. Yang, M. Yamamoto, L. Zhou, R. Ishikawa, K. Ueno, K. Tsukagoshi, Z. Zhang, M. S. Dresselhaus, and R. Saito, Phys. Rev. B 91, 205415 (2015).
  * Livneh and Spanier [2015] T. Livneh and J. E. Spanier, 2D Materials 2, 035003 (2015).
  * Carvalho  _et al._ [2017] B. Carvalho, Y. Wang, S. Mignuzzi, D. Roy, M. Terrones, C. Fantini, V. H. Crespi, L. M. Malard, and M. A. Pimenta, Nat. Commun. 8, 14670 (2017).
  * Ferrari and Basko [2013] A. Ferrari and D. Basko, Nat. Nanotechnol. 8, 235 (2013).
  * Liu  _et al._ [2015] H.-L. Liu, H. Guo, T. Yang, Z. Zhang, Y. Kumamoto, C.-C. Shen, Y.-T. Hsu, L.-J. Li, R. Saito, and S. Kawata, Phys. Chem. Chem. Phys. 17, 14561 (2015).
  * Giannozzi  _et al._ [2017] P. Giannozzi, O. Andreussi, T. Brumme, O. Bunau, M. Buongiorno Nardelli, M. Calandra, R. Car, C. Cavazzoni, D. Ceresoli, M. Cococcioni, N. Colonna, I. Carnimeo, A. Dal Corso, S. de Gironcoli, P. Delugas, R. A. DiStasio, A. Ferretti, A. Floris, G. Fratesi, G. Fugallo, R. Gebauer, U. Gerstmann, F. Giustino, T. Gorni, J. Jia, M. Kawamura, H.-Y. Ko, A. Kokalj, E. Küçükbenli, M. Lazzeri, M. Marsili, N. Marzari, F. Mauri, N. L. Nguyen, H.-V. Nguyen, A. Otero-de-la Roza, L. Paulatto, S. Poncé, D. Rocca, R. Sabatini, B. Santra, M. Schlipf, A. P. Seitsonen, A. Smogunov, I. Timrov, T. Thonhauser, P. Umari, N. Vast, X. Wu, and S. Baroni, Journal of Physics: Condensed Matter 29, 465901 (2017).
  * Herziger  _et al._ [2014] F. Herziger, M. Calandra, P. Gava, P. May, M. Lazzeri, F. Mauri, and J. Maultzsch, Phys. Rev. Lett. 113, 187401 (2014).
  * Torche  _et al._ [2017] A. Torche, F. Mauri, J.-C. Charlier, and M. Calandra, Phys. Rev. Materials 1, 041001(R) (2017).
  * Giannozzi  _et al._ [2009] P. Giannozzi, S. Baroni, N. Bonini, M. Calandra, R. Car, C. Cavazzoni, D. Ceresoli, G. L. Chiarotti, M. Cococcioni, I. Dabo, A. D. Corso, S. D. Gironcoli, S. Fabris, G. Fratesi, R. Gebauer, U. Gerstmann, C. Gougoussis, A. Kokalj, M. Lazzeri, L. Martin-Samos, N. Marzari, F. Mauri, R. Mazzarello, S. Paolini, A. Pasquarello, L. Paulatto, C. Sbraccia, S. Scandolo, G. Sclauzero, A. P. Seitsonen, A. Smogunov, P. Umari, and R. M. Wentzcovitch, J. Phys.: Condens. Matter 21, 395502 (2009).
  * Monkhorst and Pack [1976] H. J. Monkhorst and J. D. Pack, Phys. Rev. B 13, 5188 (1976).
  * Baroni  _et al._ [2001] S. Baroni, S. de Gironcoli, A. Dal Corso, and P. Giannozzi, Rev. Mod. Phys. 73, 515 (2001).
  * Noffsinger  _et al._ [2010] J. Noffsinger, F. Giustino, B. D. Malone, C.-H. Park, S. G. Louie, and M. L. Cohen, Computer Physics Communications 181, 2140 (2010).
  * Ponc$\acute{e}$  _et al._ [2016] S. Ponc$\acute{e}$, E. Margine, C. Verdi, and F. Giustino, Computer Physics Communications 209, 116 (2016).
  * Lieth [1977] R. M. A. Lieth,  _Preparation and Crystal Growth of Materials with Layered Structures_ (Springer: Berlin, 1977).
  * Yamamoto  _et al._ [2014] M. Yamamoto, S.-T. Wang, M.-Y. Ni, Y.-F. Lin, S.-L. Li, S. Aikawa, W.-B. Jian, K. Ueno, K. Wakabayashi, and K. Tsukagoshi, ACS Nano 8, 3895 (2014).
  * Spizzirri [2010] P. G. Spizzirri, in  _Microscopy: Science, Technology, Applications and Education_ , Vol. 2, edited by A. Méndez-Vilas and J. Díaz (Formatex Research Center, 2010) pp. 1389–1396.
  * Saito  _et al._ [2001] R. Saito, A. Jorio, A. G. Souza Filho, G. Dresselhaus, M. S. Dresselhaus, and M. A. Pimenta, Phys. Rev. Lett. 88, 027401 (2001).
  * Venezuela  _et al._ [2011] P. Venezuela, M. Lazzeri, and F. Mauri, Phys. Rev. B 84, 035433 (2011).
  * Birman [1962] J. L. Birman, Phys. Rev. 127, 1093 (1962).
  * Birman [1963] J. L. Birman, Phys. Rev. 131, 1489 (1963).
  * Aroyo  _et al._ [2006] M. I. Aroyo, A. Kirov, C. Capillas, J. M. Perez-Mato, and H. Wondratschek, Acta Cryst. A 62, 115 (2006).
  * Caramazza  _et al._ [2018] S. Caramazza, A. Collina, E. Stellino, F. Ripanti, P. Dore, and P. Postorino, Eur. Phys. J. B 91, 35 (2018).


# Comment on “Weak values and the past of a quantum particle”

Lev Vaidman  Raymond and Beverly Sackler School of Physics and Astronomy
Tel-Aviv University, Tel-Aviv 6997801, Israel

(May 2023)

###### Abstract

In a recent paper, Hance, Rarity and Ladyman [Phys. Rev. Res. 5, 023048
(2023)] criticized recent proposals connecting weak values and the past of a
quantum particle. I argue that their conclusion follows from a conceptual
error in understanding the approach to the past of the particle they discuss.

Hance, Rarity and Ladyman (HRL) HRL  discuss the connection between the
presence of a quantum particle in the past and weak values, the topic I
introduced in 2013 past . They claim to analyze it according to my definition
and also according to an alternative approach Cheshire . In this Comment I
argue that there is a conceptual error in the HRL paper in presentation of
these approaches and consequently, their conclusion

> “that these approaches specifically are not useful for helping identify the
> past path of quantum particles”

misses the target.

The conceptual error of the HRL paper is that according to their presentation,
the discussed approaches argue for the existence of an independent ontological
concept of the presence of a pre- and postselected particle. Perhaps I should
refrain from discussing the alternative approach Cheshire , but I can say that
this is definitely not true for my approach. The definition of the “presence
of the particle” in past  is operational: the particle was where it left a
(weak) trace. Therefore, to “identify the past path of quantum particles” is
to find the locations where they left a trace. The weak values of the local
operators are a useful tool for calculating these local traces. There is a
controversy about the faithfulness of this method, but HRL mainly criticise
the connection between weak values and hypothetical “particle presence” and
not between weak values and weak traces, which can be calculated using
standard quantum mechanics. HRL write

> “These approaches simply assume the particle was present wherever the weak
> value of an operator containing the spatial projection operator is nonzero.”

The approach past  defines that the particle was present where it left a
trace. The purpose of this Comment is to clarify the approach by pointing out
several misconceptions in the HRL paper.

HRL write (in the Introduction) about “weak values only being defined over
ensembles”. I, as co-author of the original paper which introduced weak values
AAV , disagree with this statement. It is correct that we usually need an
ensemble to observe a weak value, but nothing prevents us from defining it for
a single system beyond .

HRL also discuss disturbance in weak measurements. Apparently, they attach a
weak value to a weak measurement. We do need to know the results of the
preselection measurement and the postselection measurement, but the discussion
of the presence of the particle between these measurements does not require
weak measurement: the environment “measures” the weak value by being
disturbed. The particle is also disturbed by the environment; the weak values
are then modified and it is a subtle issue when this can or cannot be
neglected PNAS .

The definition of presence based on a weak trace requires the existence of all
possible types of local interactions with the environment. These interactions
must be non-vanishing but can be arbitrarily small. Their purpose is to serve
as a reference to the trace left on the environment in the discussed
experiments relative to a hypothetical experiment with a well-localized
particle in the same location. Then I do agree with

> “the existence of at least one operator formed from the product of the
> spatial projection operator for a location and some other operator, with a
> nonzero weak value, is both a necessary and a sufficient condition for
> particle presence at that location”.

It follows from the fact that nonvanishing local interactions ensure the
first-order trace in the environment.

Note that in optical interferometric experiments we always have a finite
interaction of the photon with mirrors. (The HRL energy exchange estimate
$10^{-33}$, has to be replaced by considering a much larger momentum exchange
of every photon bouncing off a mirror. The amplitude of the orthogonal
component of the quantum state of the mirror due to the bouncing photon in the
same experiment is of the order $10^{-17}$.)

I agree with the HRL view that

> “any attempt to form a definition of presence for quantum particles should
> correspond to our intuitions about classical presence, unless we have a good
> reason for it to deviate from this.
>
> The classical conception of a particle presence—being present at a certain
> place at a certain time—can be characterized as follows:
>
> (i) Every particle is located in space at all times.
>
> (ii) Particles cannot be on more than one path simultaneously.
>
> (iii) Particle trajectories are continuous (or at least as continuous as
> space is) so particles cannot get from one place to another without passing
> through the space in between.
>
> (iv) Particles interact with other objects and/or fields local to their
> location.
>
> (v) If a particle is on a path at a given time, and that path is within some
> region, then the particle is also located in that region at that time.
>
> (vi) If a particle’s property is at a location, the particle must be at that
> location too.”

My attempt for a new definition came exactly when I found “a good reason”. In
the nested Mach-Zehnder interferometer V07  there is a contradiction between
(iii) and (iv) . The traces left on the environment that provide evidence of
particle interactions have disconnected parts. We cannot retain all the
classical characterizations of presence in the quantum world, so in past  I
abandoned (iii) and adopted (iv) as the definition.

My definition allows for keeping all other properties, although there is a
very subtle and paradoxical situation about property (ii). In the nested Mach-
Zehnder interferometer the particle is present in every one of the three arms
at the same time, but it is not present in any two (or three) arms
simultaneously. To be present in a particular location is to leave a trace
there, and the quantum states of the environment in all three arms are
changed, so the particle was in three places. However, the traces in the
environment are entangled such that orthogonal components of the local states
of the environment, which provide evidence of the presence in every location,
are entangled with undisturbed states in all other locations. Thus, there is
no trace corresponding to simultaneous presence in different locations and one
can claim that (ii) is fulfilled. A similar paradoxical situation arises in
the Hardy paradox Hardy ; HardyV ; Product  which describes a pre- and
postselected system of two particles: an electron was present in one arm and a
positron was present in another, but the particles were not present in these
arms simultaneously.

There is no similar difficulty with (v), since the definition is that if
anywhere in the region there is a non-vanishing trace, the particle was in
this region. Note that due to the unavoidable momentum exchange of the photon
with mirrors, there is no such thing as “undisturbed inner interferometer”
discussed by HRL. Contrary to the HRL claim, the fact that traces might have
various properties (e.g. sign) is not neglected in the two-state vector
formalism and asking specific questions like: Was the photon in two places
together? leads to paradoxical answers, as in the discussion of (ii), see
ChinaCom .

I am puzzled by the HRL claim about the inconsistency of Peled . Why “would we
expect a particle to necessarily have a non-zero weak value for the spatial
projection operator for any path along which it travels” when interactions
with other degrees of freedom lead to the trace?

I want to repeat that I disagree with Section VI of HRL, my weak value is
defined for a single system. The fact that the weak value of the velocity of a
particle can be larger than the speed of light (see Sec. VIII of AV90 ) does
not contradict the special theory of relativity. The experiments involve
postselection and their low probability of success prevents a superluminal
change in the probability of finding a quantum particle.

I also disagree with the claim of Sec. VII of the HRL paper, according to
which the weak value approach is intended to show that some quantum protocols
“are not as ‘spooky’ as they appear”. The weak value approach helps to find
quantum protocols which are “spooky” if analyzed in classical terms. My papers
based on the weak value approach cited by HRL 57 ; 58 ; 59 ; 60 ; 62  do not
try to remove paradoxical features. Instead, these papers try to correct
erroneous claims about alleged counterfactual communication. In particular,
HRL are correct in their weak values analysis of the protocol described in
Laws  and shown on their Fig. 2. The photons reaching detector $D_{0}$ were
not present at Bob’s site according to the weak trace criterion (all weak
values of local operators on Bob’s site vanish). However, there is no
contradiction with the approach because, in this case, Bob’s communication
with Alice fails. The click at $D_{0}$ means that the photon did not perform
any test of the presence or absence of Bob’s shutter because the probability
of this click does not depend on Bob’s actions, see Counterport .

Finally, let me comment on the concluding sentence of HRL

> “we have shown that weak value approaches to the path of a particle do not
> contribute any new physics—the assumption of a connection between particle
> presence and weak values does not give us anything testable.”

First, weak values, as all other concepts and results of the two-state vector
formalism are fully consistent with the standard formalism of quantum
mechanics, so neither new physics, i.e., a deviation from the Schrödinger
equation, nor introducing some new ontology, is proposed. My approach
introduces new concepts (which I believe are useful), in particular, the local
presence of a pre- and post-selected particle defined by the local trace it
leaves on the environment. The formalism predicts that these traces can be
found based on finite weak values of local operators, and this statement is
definitely testable.

This work has been supported in part by the U.S.-Israel Binational Science
Foundation (Grant No. 735/18) and by the Israel Science Foundation Grant No.
2064/19.

## References

  * [1] J.R. Hance, J. Rarity and J. Ladyman, Weak values and the past of a quantum particle, Phys. Rev. Res. 5, 023048 (2023).
  * [2] L. Vaidman, Past of a quantum particle, Phys. Rev. A 87, 052104 (2013).
  * [3] Y. Aharonov, S. Popescu, D. Rohrlich, and P. Skrzypczyk, Quantum Cheshire cats, New J. Phys. 15, 113015 (2013).
  * [4] Y. Aharonov D. Alabert, and L. Vaidman, How the result of a measurement of a component of the spin of a spin-1/2 particle can turn out to be 100, Phys. Rev. Lett. 60, 1351 (1988).
  * [5] L. Vaidman, A. Ben-Israel, J. Dziewior, L.Knips, M. Weißl, J. Meinecke, C. Schwemmer, R. Ber, and H. Weinfurter, Weak value beyond conditional expectation value of the pointer readings, Phys. Rev. A 96, 032114 (2017).
  * [6] J. Dziewior, L. Knips, D. Farfurnik, K. Senkalla, N. Benshalom, J. Efroni, J. Meinecke, S. Bar-Ad, H. Weinfurter, and L. Vaidman, Universality property of local weak interactions and its application for interferometric alignment, Proc. Natl. Acad. Sci. USA, 116, 2881 (2019).
  * [7] L. Vaidman, Impossibility of the counterfactual computation for all possible outcomes, Phys. Rev. Lett. 98, 160403 (2007).
  * [8] L. Hardy, Quantum mechanics, local realistic theories, and Lorentz-invariant realistic theories, Phys. Rev. Lett. 68, 2981-2984 (1992).
  * [9] L. Vaidman, Lorentz-invariant “elements of reality” and the joint measurability of commuting observables Phys. Rev. Lett. 70, 3369-3372 (1993).
  * [10] X-Y. Xu, W-W Pan, Q-Q. Wang, J. Dziewior, L. Knips, Y. Kedem, Kai Sun, J-S. Xu, Y-J. Han, C-F. Li, G-C Guo, and L. Vaidman, Measurements of nonlocal variables and demonstration of the failure of the product rule for a pre- and postselected pair of photons, Phys. Rev. Lett. 122, 100405 (2019).
  * [11] A. Ben-Israel, L. Knips, J. Dziewior, J. Meinecke, A. Danan, H.Weinfurter and L. Vaidman, An improved experiment to determine the ‘past of a particle’ in the nested Mach-Zehnder interferometer, Chin. Phys. Lett. 34, 020301 (2017).
  * [12] U. Peleg and L. Vaidman, Comment on “Past of a quantum particle revisited”, Phys. Rev. A 99, 026103 (2019).
  * [13] Y. Aharonov and L. Vaidman, Properties of a quantum system during the time interval between two measurements, Phys. Rev. A 41, 11 (1990).
  * [14] L. Vaidman, Comment on “Protocol for direct counterfactual quantum communication”, Phys. Rev. Lett. 112, 208901 (2014).
  * [15] L. Vaidman, Counterfactuality of ‘counterfactual’ communication, J. Phys. A: Math. Theor. 48, 465303 (2015).
  * [16] L. Vaidman, “Counterfactual” quantum protocols, Int. J. Quantum Inf. 14, 1640012 (2016).
  * [17] L. Vaidman, Analysis of counterfactual communication protocols, Phys. Rev. A 99, 052127 (2019).
  * [18] Y. Aharonov and L. Vaidman, Modification of Counterfactual communication protocols that eliminates weak particle traces, Phys. Rev. A 99, 010103(R) (2019).
  * [19] H. Salih, W. McCutcheon, J. R. Hance, and J. Rarity, The laws of physics do not prohibit counterfactual communication, NPJ Quantum Inf. 8, 60 (2022).
  * [20] J. Dressel, G. Reznik, L. Vaidman, Counterportation and the two-state vector formalism, arXiv:2303.08962 (2023).


# Basic Engagement Zones

Alexander Von Moll111Aerospace Engineer and Isaac Weintraub222Electronics
Engineer, AIAA Senior Member Air Force Research Laboratory, 2210 8th St.
WPAFB, OH 45433 USA

##  1 Introduction

Navigation in threat-laden environments is a fundamental problem for many
applications and missions [1, pp. III-1]. †† This paper is based on work
performed at the Air Force Research Laboratory (AFRL) Control Science Center.
Distribution Statement A: Approved for Public Release; Distribution is
Unlimited. PA# AFRL-2023-5574. For example for a supply airplane to reach its
destination in a reasonable time, it may have to fly over or near unfriendly
areas containing one or more ground-based enemy defensive assets, which could
fire upon the airplane should it come too close. Ideally, the airplane and its
crew strives to reach a designated goal as quickly as possible, while avoiding
aggressive maneuvers to avoid interception, should one of the defensive assets
decide to engage. While there is a wealth of literature concerning navigation
in obstacle-laden environments (c.f., e.g., [2]), there is a gap concerning
the case when the “obstacles” are not only dynamic but also adversarial.

When it comes to modeling adversarial engagements between mobile agents, a
common paradigm has been differential games (c.f., e.g., [3]). If a scenario
can be described wherein two players have a diametrically opposed objective
that can be distilled to a single scalar value, one can formulate a zero-sum
differential game. The solution of a differential game may be comprised of the
Value function that describes the equilibrium value of the cost/reward as a
function of the initial condition as well as the equilibrium control policies
for the two players (often, the Nash equilibrium is the desired solution
concept). Two particular classes of differential game problems are
particularly relevant to the adversarial navigation problem described above:
pursuit-evasion and target guarding (which also appears in the literature as
reach-avoid). In the former, a Pursuer generally seeks to capture an Evader in
minimum time or minimum control effort, while in the latter, a Defender
generally seeks to capture an Attacker before the latter is able to reach some
target or goal region; many of these problems and the literature surrounding
them are described in [4] . However most of these solutions, in themselves,
are not well-suited to address the navigation problem. With pursuit-evasion
solutions, for example, it’s difficult to incorporate objectives for the
Evader outside of merely evading without completely changing the problem
formulation. Target Guarding solutions, on the other hand, typically delineate
where breach or capture can be guaranteed but generally do not account for
possible constraints on the Defender.

Some notable and relevant exceptions include the following works. In [5] the
authors explicitly formulated a target guarding game with a fixed time limit.
This could be seen as a constraint on the Defender’s onboard fuel, for
example. However, the solution only answers the question of what is the
outcome (breach or capture) and equilibrium headings of the agents as a
function of the initial positions assuming that the engagement (i.e., the
playout of the differential game) begins immediately. In other words, it does
not directly address the strategy of either agent before or after the
engagement. Another example is the work in [6]. There, a zero-sum stochastic
differential game is formulated in which one player desires to keep the state
of the system in the safe set for some specified duration (akin to a range or
fuel constraint) while the other desires to cause the system to enter the
unsafe set. Like the previous work, this solution describes what the agents
should do once the clock starts but does not directly address how the Evader
should navigate to avoid (a priori) getting into a situation in which it could
lose this game.

In order to begin addressing what may happen before or after an engagement
takes place we introduce the following.

######  Definition 1 (Engagement Zone).

Given a Mobile Agent, $A$, Threat, $T$, their respective dynamic models,
$\dot{\mathbf{x}}_{A}$, $\dot{\mathbf{x}}_{T}$, and an Agent strategy,
$u_{A}(t)$, an Engagement Zone (EZ) is a region of the state space in which it
is possible for the Threat to neutralize the Mobile Agent if the latter does
not deviate from its current strategy.

Note this is just one particular definition for EZ. The interpretation is
this: if the Agent remains outside of the EZ associated with a particular
Threat that the latter may not have any incentive to actively engage the Agent
as 1) neutralization (e.g., capture) cannot be guaranteed and 2) the Agent
does not even need to actively maneuver/evade to avoid the Threat. The latter
point distinguishes the EZ from the less conservative approach of win/lose
regions obtained from solutions of differential games. Navigating near the
win/lose region boundary could necessitate an aggressive maneuver by the Agent
should the Threat decide to engage.

Some earlier papers by the authors have begun to address navigation around
EZs. In [7] the single-Agent, single-EZ navigation problem was introduced.
There, path plans were devised for which the Agent avoided entering into the
EZ and others which allowed some penetration into the EZ to reduce overall
travel time or reach the goal at a specified arrival time. Then the results
were extended in [8] to address navigation of a single vehicle around (or
through) two EZs. In both of these works, a notional EZ model based on a
cardioid shape was assumed. This EZ was not tied to a particular model of the
Threat itself. In practice, for a particular Agent and Threat model, one may
compute EZs via simulation, in which case the EZ may be data-based and not
analytic which could make planning a path more computationally expensive. This
paper makes the following two contributions: 1) a set of EZs that are based on
first-principle models of the Agent and Threat 2) a direct comparison of EZ-
based navigation and path-planning to the nominal approach of
circumnavigation.

The remainder of the paper is organized as follows. Sections 2 and 3 model the
Threat as a Pursuer and a Turret, respectively, and derive the associated EZs.
Section 4 highlights path planning as a potential application of EZs. Section
5 concludes the paper.

##  2 Pursuit-Evasion

###  2.1 Basic Pursuit-Evasion Model

In this section, the Threat from Definition 1 is specialized to be a Pursuer.
Perhaps the most basic dynamic model to consider for pursuit-evasion is that
of both agents moving with simple motion in the two-dimensional plane (i.e.,
$A,P\in\mathbb{R}^{2}$). That is, both agents move with constant speeds and
have control over their instantaneous heading [9]:

$\dot{\mathbf{x}}=\begin{bmatrix}\dot{\mathbf{x}}_{P}\\\
\dot{\mathbf{x}}_{A}\end{bmatrix}=\begin{bmatrix}\dot{x}_{P}\\\ \dot{y}_{P}\\\
\dot{x}_{A}\\\ \dot{y}_{A}\end{bmatrix}=\begin{bmatrix}v_{P}\cos ψ_{P}\\\
v_{P}\sin ψ_{P}\\\ v_{A}\cos ψ_{A}\\\ v_{A}\sin
ψ_{A}\end{bmatrix},\qquad\mathbf{x}(0)=\begin{bmatrix}x_{P_{0}}\\\
y_{P_{0}}\\\ x_{A_{0}}\\\ y_{A_{0}}\end{bmatrix}$ (1)

where the Mobile Agent’s control is $u_{A}\equiv ψ_{A}$. For this scenario,
neutralization is said to occur if $A$ comes within $r\geq 0$ distance of $P$
(i.e., $A$ is captured by $P$):

$\mathcal{N}=\left\\{\mathbf{x}\mid\overline{AP}\leq r\right\\},$ (2)

where $\mathcal{N}$ is the neutralization set, and the notation
$\overline{AP}$ is used to indicate the distance between two points, e.g., $A$
and $P$. Define the speed ratio of Mobile Agent and Pursuer speeds as
$μ=\tfrac{v_{A}}{v_{P}}>0$. Note that since the EZ is based on the Mobile
Agent implementing zero control capture may be possible, depending on the
initial conditions, even when $μ>1$. It is common to consider such a pursuit-
evasion scenario as-is without any additional constraints. Alternatively, some
work has been done on similar scenarios taking place in the presence of
obstacles [10], in bounded domains [11], and under a variety of other
geometric and integral constraints [12]. This paper considers the Pursuer to
be range-limited with maximum range $R$ as in [9]. Thus the Pursuer is capable
of reaching any point within the disk of radius $R$ centered at its initial
position, and thus its reachability region is

$\mathcal{R}_{P}=\left\\{(x,y)\mid\left(x-x_{P_{0}}\right)^{2}+\left(y-y_{P_{0}}\right)^{2}\leq
R^{2}\right\\}.$ (3)

Note that one possible interpretation of the range-limit is that the Pursuer
may have arbitrary endurance as long as it stays within $R$ of its initial
position. Under this interpretation, the range-limit may represent a maximum
communication range or sensing range for some base station located at the
Pursuer’s initial position. Another interpretation is that the Pursuer may
only traverse a total distance $R$ along its trajectory, which does not
necessarily preclude the Pursuer from stopping. Based on the constant speed
model of the Pursuer, however, stopping would require, e.g., modulating its
heading, $ψ_{P}$, infinitely fast. In this work, such a maneuver, although
mathematically possible, is considered to consume the Pursuer’s remaining
range at the same rate as if it were moving. Therefore, the range-limit is
analogous to a time-limit – the Pursuer may be equipped with, for example, a
solid rocket motor which, once started, has a fixed time at which fuel will
run out (barring changes in environmental conditions or changes in various
forces experienced along the trajectory).

######  Remark 1.

If the Pursuer has additional, unmodeled dynamics (e.g., a bounded turn rate)
then this analysis is conservative.

###  2.2 Pursuit-Evasion EZ

######  Lemma 1.

Under the model Equation 1 if the Mobile Agent, starting from $A_{0}$ and
moving with $ψ_{A}$ cannot be neutralized by $P$ under collision course
guidance then the point $A_{0}$ is outside of the EZ.

###### Proof.

Collision course, wherein the Pursuer takes a straight-line path along the
line connecting its initial position to the Mobile Agent’s position at the
time of neutralization, yields the minimum time trajectory for $P$ to
intercept $A$ [9]. Therefore, if the collision course trajectory results in
$P$ travelling a distance greater than $R$, then the Mobile Agent cannot be
neutralized under any guidance law. Thus, from Definition 1, the initial
position of $A$ must be outside the EZ. ∎

######  Corollary 1.

The boundary of the EZ is the locus of all Mobile Agent initial positions
corresponding to neutralization by the Pursuer via a collision course
trajectory of length equal to its maximum range, $R$.

This result is useful in obtaining the boundary of the EZ as follows. Define
the quantity $ρ\equiv\overline{P_{0}A_{0}}$ to be the initial distance between
the agents resulting in the Pursuer taking a collision course trajectory of
length $R$. Then the Engagement Zone is defined mathematically by

$\mathcal{Z}=\left\\{A_{0}\mid\overline{P_{0}A_{0}}\leq ρ(ξ;μ,R,r)\right\\},$
(4)

where $ξ$, the aspect angle, is the angle the Agent’s heading makes with the
initial line-of-sight angle. Define a set of coordinate axes
$(\hat{x},\hat{y})$ whose origin is $P_{0}$ in which the $\hat{x}$ vector is
aligned with the heading of the Mobile Agent.

######  Remark 2.

This model and the resulting EZ may be the simplest way of directly accounting
for the effect of aspect angle and closing speed on the engagement which was
the main physical phenomena under investigation in [7].

####  2.2.1 Fast Pursuer

For the case of a fast Pursuer ($μ\leq 1$), the Law of Cosines applied to the
triangle $\triangle P_{0}E_{0}E_{f}$ may be used to obtain an expression for
$ρ$ (see Figure 1):

$ρ(ξ;μ,R,r)=μR\left[\cos
ξ+\sqrt{\cos^{2}ξ-1+\frac{(R+r)^{2}}{μ^{2}R^{2}}}\right],\qquad μ\leq 1,\
ξ\in\left\\{-π,π\right\\}.$ (5)

Figure 1 shows an example EZ for the fast Pursuer case along with relevant
geometry. The symbol $\partial$ is used to denote the boundary of a set or
region.

Figure 1:  Range-limited, simple motion pursuit-evasion Engagement Zone for
fixed Agent heading and fast Pursuer (speed ratio $μ=0.7$) with finite capture
radius ($r=0.25$) and a maximum range of $R=1$.

The paper by the authors [7] contains an expression for an EZ that is similar
to Equation 5 but is not based on a particular model of the engagement between
the Mobile Agent and the Threat. That expression (with changed notation) is

$\tilde{ρ}=\frac{\cos
ξ+1}{2}\left(\tilde{ρ}_{\max}-\tilde{ρ}_{\min}\right)+ρ_{\min}$ (6)

where $\tilde{ρ}_{\max}\equiv\tilde{ρ}(ξ=0)$ and
$\tilde{ρ}_{\min}\equiv\tilde{ρ}(ξ=π)$. The two models for the distance of the
EZ boundary from the Threat’s initial position may be compared by setting
$\tilde{ρ}_{\max}=\max_{ξ}ρ=\left(1+μ\right)R+r$ and
$\tilde{ρ}_{\min}=\min_{ξ}ρ=\left(1-μ\right)R+r$. Figure 2 shows a comparison
of Equations 5 and 6 for a variety of $μ<1$. Note that the biggest deviation
between the two models occurs for higher $μ$ and for aspect angles near the
center of the domain $\lvert ξ\rvert\in\left[0,π\right]$.

Figure 2:  Comparison of the Pursuit-Evasion EZ model presented in this paper
with a previously developed generic EZ model from [7].

This section is focused on the particular interpretation of the EZ that is
based on fixing $ψ_{A}$ and varying $ξ$ to construct the region $\mathcal{Z}$.
Another valid interpretation of the EZ is based on fixing $ξ$ and varying
$ψ_{A}$. The same geometry depicted in Figure 1 applies for the example shown.
Based on Equation 5 it is clear that if $ξ$ is held constant $ρ$ is constant
as it does not depend directly on $ψ_{A}$. Thus, the EZ for a fixed aspect
angle is a circle centered at $P_{0}$, as shown by the dashed green circle in
Figure 1. These two interpretations of the EZ boundary yield the same $A_{0}$
for the points corresponding to the $ξ$ and $ψ_{A}$ values used.

####  2.2.2 Slow Pursuer

When the Pursuer is relatively slow ($μ>1$) then there are aspect angles ($ξ$
values) for which neutralization is possible only if the Mobile Agent started
on or within the capture disk of the Pursuer. This is because of the fact that
in order for the Pursuer to move and neutralize the Mobile Agent their
separation distance must be non-increasing at the time of neutralization. Thus
the optimal Pursuer heading must be $\cos ψ_{P}^{*}\leq\cos^{-1}\tfrac{1}{μ}$.
This is akin to the Usable Part of the terminal surface in the verbiage of
Isaacs [3]. An analogous result appears in the related scenario covered in
[13] in which this condition is used to delineate the positions from which a
slow Observer may establish contact with a fast Target. The locus of initial
Mobile Agent positions corresponding to collision course with
$ψ_{P}^{*}=\cos^{-1}\tfrac{1}{μ}$ (with the Pursuer travelling $R$ distance or
less) forms a portion of the EZ boundary. From the Law of Sines, the Pursuer’s
travel distance may be expressed as

$\overline{P_{0}P_{f}}=\frac{r\sin(ξ+ψ_{P})}{μ\sin ξ-\sin(ξ+ψ_{P})}.$ (7)

Invoking the Law of Sines once more and substituting in the above with
$ψ_{P}=\cos^{-1}\tfrac{1}{μ}$ gives

$ρ(ξ;μ,R,r)=\frac{r\sqrt{μ^{2}-1}}{μ\sin\lvert ξ\rvert-\sin\left(\lvert
ξ\rvert+\cos^{-1}\tfrac{1}{μ}\right)}.$ (8)

It can be verified that the aspect angle which satisfies both Equations 5 and
8 is

$ξ_{c}=\pm\sin\left(\frac{\left(R+r\right)\sqrt{μ^{2}-1}}{μR\sqrt{μ^{2}-1+\frac{r^{2}}{R^{2}}}}\right).$
(9)

Of course, the last possible $ξ$ value in this family of trajectories
corresponds to the one where $\overline{P_{0}P_{f}}=0$ which gives
$ξ_{\max}=π-\cos^{-1}\tfrac{1}{μ}$. Finally, the complete expression for $ρ$
in the case of a slow Pursuer is given by

$ρ(ξ;μ,R,r)=\begin{cases}μR\left[\cos
ξ+\sqrt{\cos^{2}ξ-1+\frac{\left(R+r\right)^{2}}{μ^{2}R^{2}}}\right]&\text{ if
}\lvert ξ\rvert\in\left[0,ξ_{c}\right]\\\ \frac{r\sqrt{μ^{2}-1}}{μ\sin\lvert
ξ\rvert-\sin\left(\lvert ξ\rvert+\cos^{-1}\frac{1}{μ}\right)}&\text{ if
}\lvert ξ\rvert\in\left(ξ_{c},π-\cos^{-1}\frac{1}{μ}\right]\\\ r&\text{
otherwise}\end{cases},\qquad μ>1,\ ξ\in\left[-π,π\right].$ (10)

Figure 3 shows an example for the slow Pursuer case along with relevant
geometry.

Figure 3:  Range-limited, simple motion pursuit-evasion Engagement Zone for
fixed Agent heading and slow Pursuer (speed ratio $μ=1.5$) with finite capture
radius ($r=0.25$) and a maximum range of $R=1$. Two example trajectories are
shown: (blue/red) the Agent starts on the boundary of the EZ wherein the
capture configuration is general, and (light blue/pink) the Agent starts on
the EZ boundary wherein the Pursuer-Agent distance rate is zero at capture.

##  3  Turret-Evasion

###  3.1 Basic Turret Model

In this section, the Threat from Definition 1 is specialized to be a Turret,
denoted by $T$. The Turret is a stationary agent with a look angle, $θ$,
bounded turn rate, $ω\in\left[-\bar{\omega},\bar{\omega}\right]$, and finite
range, $R$. As in the previous section, the Mobile Agent, moves with constant
speed and controls its instantaneous heading. A similar model has been used,
e.g., in [14, 15, 16]. The kinematics of the system in polar form are

$\dot{\mathbf{x}}=\begin{bmatrix}\dot{\mathbf{x}}_{T}\\\
\dot{\mathbf{x}}_{A}\end{bmatrix}=\begin{bmatrix}\dot{θ}\\\ \dot{x}_{A}\\\
\dot{y}_{A}\end{bmatrix}=\begin{bmatrix}ω\\\ v_{A}\cos ψ_{A}\\\ v_{A}\sin
ψ_{A}\end{bmatrix},\qquad\dot{\mathbf{x}}(0)=\begin{bmatrix}θ_{0}\\\
x_{A_{0}}\\\ y_{A_{0}}\end{bmatrix}$ (11)

where, again, the Mobile Agent’s control is $u_{A}\equiv ψ_{A}$.
Neutralization is said to occur if $\cos θ=\tfrac{x_{A}}{\overline{AT}}$,
$\sin θ=\tfrac{y_{A}}{\overline{AT}}$, and $\overline{AT}\leq R$. Define the
ratio of the Mobile Agent’s speed and Turret’s maximum turn rate as
$μ=\tfrac{v_{A}}{\bar{ω}}$. Note that, without the range constraint on the
Turret that neutralization is possible from any initial condition
(eventually). Similar to the range-limited Pursuer considered in the previous
section, all the points within a distance of $R$ from the $T$ are considered
to be reachable by the Turret, i.e.,

$\mathcal{R}_{T}=\left\\{(x,y)\mid\left(x-x_{T}\right)^{2}+\left(y-y_{T}\right)^{2}\leq
R\right\\}.$ (12)

###  3.2 Turret EZ

The EZ constructed for the turret-evasion model is fundamentally different
from the one constructed for the Pursuit-Evasion model in the previous
section. That is due to the fact that in the former it was assumed that $P$
was removed from the scenario once its range/fuel had been exhausted, whereas
here the Turret is always “active” even if it is not moving. For example, if
the Mobile Agent begins very far away and is aimed at the Turret’s position it
will eventually collide with the Turret’s beam thereby becoming neutralized.
Based on Definition 1, that point, no matter how far away, is considered to be
in the EZ. Therefore, the construction of the turret-evasion EZ is only based
upon terminal Agent positions which are exiting the Turret’s reachable set,
$\mathcal{R}_{T}$ (i.e., the right half-circle in the $(\hat{x},\hat{y})$
coordinate system). For a particular terminal Agent position, $A_{f}$, it is
assumed that $T$ turned in the shortest direction to align with the Agent.
Define the angle traversed by the Turret during an engagement as $γ$, thus
from this assumption $γ\in\left[-π,π\right]$. Moreover, because only those
$A_{f}$ positions which are exiting $\mathcal{R}_{T}$ matter the range for $γ$
is even further restricted:

$γ\in\begin{cases}\left[-\frac{π}{2}-θ_{0},\ 0\right]\cup\left[0,\
\frac{π}{2}-θ_{0}\right]&\text{ if }\cos θ_{0}>0\\\ \left[-π,\
\frac{π}{2}-θ_{0}\right]\cup\left[\frac{3π}{2}-θ_{0},\ π\right]&\text{ if
}\cos θ_{0}<0,\ \sin θ_{0}>0\\\ \left[-π,\
-\frac{3π}{2}-θ_{0}\right]\cup\left[-\frac{π}{2}-θ_{0},\ π\right]&\text{ if
}\cos θ_{0}<0,\ \sin θ_{0}<0\end{cases}$ (13)

This range for $γ$ gives results in the final Turret look angle to be in the
range $θ_{f}\in\left[-π/2,\ π/2\right]$.

######  Proposition 1.

Under the model Equation 11 the most limiting scenario (in terms of initial
Agent positions, $A_{0}$) is when $T$ turns at its maximum angular speed and
neutralizes $A$ exactly when it is exiting $\mathcal{R}_{T}$, i.e., $\cos
θ_{f}=\tfrac{x_{A_{f}}}{R}$ and $\sin θ_{f}=\tfrac{y_{A_{f}}}{R}$ where
$θ_{f}\equiv θ_{0}+γ$.

(a) $θ_{0}=\tfrac{π}{6}$

(b) $θ_{0}=\tfrac{5π}{6}$

Figure 4:  Range-limited turret-evasion Engagement Zone for fixed Agent
heading and Turret (speed ratio $μ=0.5$) with maximum range of $R=1$.

##  4 Application to Path Planning

One possible application of these results is to plan paths which stay outside
of the EZ. Perhaps the most basic instantiation of such a scenario is to
specify initial and goal positions for a vehicle, place a range-limited
Pursuer in between, and specify a desire for the vehicle to reach the goal
position in minimum time. This usage of the EZ was the subject of [7, 8],
although, in the former, considerations for acceptance of some entry into the
EZ were included. Mathematically speaking, we wish to solve the following
problem:

$\displaystyle\min_{ψ(t)}\qquad$ $\displaystyle t_{f}\qquad$ (14) subject to
$\displaystyle\mathbf{x}_{A}(0)=\begin{bmatrix}x_{A_{0}}&y_{A_{0}}\end{bmatrix}^{\top},\qquad$
$\displaystyle\mathbf{x}_{A}(t_{f})=\begin{bmatrix}x_{A_{f}}&y_{A_{f}}\end{bmatrix}^{\top},\qquad$
$\displaystyle\left(\mathbf{x}_{A}(t),\ ψ(t)\right)\notin\mathcal{Z},\qquad$
$\displaystyle\forall t\in\left[0,\ t_{f}\right],$

where $\mathcal{Z}$ is given by Equation 4. Note that the instantaneous EZ
shape depends on the Agent’s position as well as its heading. Equation 14 is
solved by first discretizing the trajectory and then solving via a nonlinear
program. Specifically, even collocation is used (wherein the path constraint
of staying outside the EZ is imposed at the discrete nodes) and the Ipopt
solver [17] is used within the JuMP package [18] for the Julia programming
language.

In order to illustrate some of the potential benefits of the EZ-based path
planning approach several circumnavigation-based nominal trajectories are
described in the following. Each circumnavigation-based trajectory has an
associated circle, centered on $P$ with a particular radius, which the Agent
navigates around in minimum time. In particular, the Agent heads from its
initial position towards the tangent on the associated circle, remains on the
circle for some portion of its trajectory, and then departs the circle
tangentially to reach the specified goal location. Therefore, the associated
travel time is given by

$t_{f}^{\bigcirc}=\frac{1}{μ}\left(\sqrt{\overline{PA_{0}}-\hat{R}^{2}}+\sqrt{\overline{PA_{f}}-\hat{R}^{2}}+\hat{R}\lvert
θ_{2}-θ_{1}\rvert\right),$ (15)

where $\hat{R}$ is the radius of the circle and $θ_{1}$, $θ_{2}$ are given by

$θ_{1}=\operatorname*{atan2}\left(y_{A_{0}},\
x_{A_{0}}\right)-\cos^{-1}\left(\frac{\hat{R}}{\overline{PA_{0}}}\right),\qquad
θ_{2}=\operatorname*{atan2}\left(y_{A_{f}},\
x_{A_{f}}\right)+\cos^{-1}\left(\frac{\hat{R}}{\overline{PA_{f}}}\right),$
(16)

where $\operatorname*{atan2}$ is the two-argument inverse tangent function.
Three different circumnavigation radii are considered and summarized in Table
1. The last column in the table gives the improvement (negative) or loss
(positive) of the EZ-based path plan w.r.t. the associated circumnavigation
path plan in terms of time, i.e.,

$\text{EZ \%
difference}=\left(\frac{t_{f}^{\textrm{EZ}}-t_{f}^{\bigcirc}}{t_{f}^{\bigcirc}}\right)\cdot
100.$ (17)

The trajectories are shown in Figure 5.

Table 1: Circumnavigation Trajectories Label | $\hat{R}$ | Time | EZ % difference
---|---|---|---
Reach | $R$ | 7.46 | -1.61%
Worst | $(1+μ)R+r$ | 8.44 | -13.0%
Apol | $(1-μ)R+r$ | 7.04 | 4.32%
Figure 5:  Example path planning problem with $μ=0.9$, $r=0.2$, and
$R=\tfrac{2-r}{μ+1}\approx 0.95$, roughly corresponding to the example in [7,
Fig. 4]. The highlighted portion of the green EZ trajectory indicates
activation of the constraint, i.e., $\left(\mathbf{x}_{A}(t),\
ψ(t)\right)\in\partial\mathcal{Z}$.

The circumnavigation trajectory labeled Reach is based on navigating around
the circle corresponding to the reachable region of $P$ accounting for its
capture radius, i.e. the capturability region, hence $\hat{R}=R+r$. If the
Agent follows this trajectory it need not react to Pursuer if the latter were
to begin moving since capture is not possible for any of the points along the
trajectory (depending on whether the boundary of the capturability region is
considered capture or not). The EZ-based path plan results in a time savings
of 1.61% w.r.t. this path plan due to the fact that the former passes through
the capturability region. For those points inside said region, if $P$ were to
begin moving $A$ would be guaranteed to exit the region before capture could
occur simply by maintaining its heading (by way of these points being outside
the EZ). Thus $P$ has little to no incentive to begin moving in the first
place.

Next, the circumnavigation trajectory labeled Worst is based on the distance
$A$ would have to be from $P$ such that if $A$ were heading directly to $P$
when the latter began moving that $A$ would just barely avoid capture, hence
$\hat{R}=(1+μ)R+r$. This trajectory is overly conservative because at no point
along the trajectory is $A$ ever pointed towards $P$. Nonetheless, this
trajectory may be employed in practice due to the relative ease in estimating
this worst-case distance as opposed to computing a more realistic EZ shape
(e.g., in higher fidelity scenarios where the EZ can only be computed via a
computationally expensive numerical procedure). Additionally, safety is
guaranteed irrespective of the Agent’s heading. This is in contrast with the
Reach trajectory since the Agent’s heading must point tangentially or outside
the capturability region when it is on the boundary. The EZ-based path plan
results in a time savings of 13% w.r.t. this path plan.

Finally, the circumnavigation trajectory labeled Apol is based on the distance
$A$ would have to be from $P$ such that $A$ could just barely escape from $P$,
hence $\hat{R}=(1-μ)R+r$. This condition is described in detail in [9].
Navigating along this trajectory is inherently riskier as it may require $A$
to perform an extreme maneuver in order to avoid being captured by $P$ (unlike
any of the other trajectories considered). In fact, when $A$’s position is
such that $\overline{AP}=(1-μ)R+r$ there is only one heading which can
guarantee escape: heading directly away from $P$. Depending on where $A$ is
located when $P$ begins moving, performing this evasive maneuver could
lengthen $A$’s trajectory considerably. This is a high price to pay
considering that the EZ-based path plan is only 4.32% longer. The degree of
lengthening that the Pursuer can achieve on the Evader’s path is out of the
scope of this study but will be considered in future work. One instance of
that problem has been considered, for example, in [19].

##  5 Conclusion

This paper has established a more formal definition for an Engagement Zone and
derived some basic EZs associated with fundamental engagement models
associated with pursuit-evasion and turret-evasion. The basic EZs presented in
this paper capture the most salient aspects of the Pursuer-Agent and Turret-
Agent engagements: namely the geometry of the aspect angle and the relative
differences in capability (i.e., maximum speeds, range, etc.). One of the main
advantages of utilizing EZs for path planning is that they encode an overall
desire for Agent to go somewhere without requiring an aggressive maneuver or
active evasion should the Pursuer or Turret begin its pursuit. It was shown
that there is some advantage in terms of time savings in EZ-based navigation
around a single range-limited Pursuer as compared with circumnavigating the
capturability region. Future work will focus on investigating worst-case path
lengthening in the event that the Threat decides to engage with the Agent
while it is en route.

## References

  * Scott [2017] Scott, K. D., “Countering Air and Missile Threats,” Tech. rep., Joint Chiefs of Staff, United States of America, 2017.
  * Patle et al. [2019] Patle, B., Babu L, G., Pandey, A., Parhi, D., and Jagadeesh, A., “A review: On path planning strategies for navigation of mobile robot,” _Defence Technology_ , Vol. 15, 2019, pp. 582–606.  10.1016/j.dt.2019.04.011.
  * Isaacs [1965] Isaacs, R., _Differential Games: A Mathematical Theory with Applications to Optimization, Control and Warfare_ , Wiley, New York, 1965.
  * Weintraub et al. [2020] Weintraub, I. E., Pachter, M., and Garcia, E., “An Introduction to Pursuit-evasion Differential Games,” _2020 American Control Conference (ACC)_ , IEEE, Denver, CO, USA, 2020, pp. 1049–1066.  10.23919/acc45564.2020.9147205.
  * Chen et al. [2022] Chen, X., Yu, J., Yang, D., and Niu, K., “A geometric approach to reach-avoid games with time limits,” _IET Control Theory & Applications_, 2022.  10.1049/cth2.12374.
  * Patil et al. [2023] Patil, A., Zhou, Y., Fridovich-Keil, D., and Tanaka, T., “Risk-Minimizing Two-Player Zero-Sum Stochastic Differential Game via Path Integral Control,” 2023.
  * Weintraub et al. [2022] Weintraub, I. E., Von Moll, A., Carrizales, C., Hanlon, N., and Fuchs, Z., “An Optimal Engagement Zone Avoidance Scenario in 2-D,” _Scitech_ , AIAA, San Diego, 2022.  10.2514/6.2022-1587.
  * Dillon et al. [2023] Dillon, P. M., Zollars, M. D., Weintraub, I. E., and Von Moll, A., “Optimal Trajectories for Aircraft Avoidance of Multiple Weapon Engagement Zones,” _Journal of Aerospace Information Systems_ , Vol. 20, 2023, pp. 520–525.  10.2514/1.I011224.
  * Weintraub et al. [2023a] Weintraub, I., Von Moll, A., and Pachter, M., “Range-Limited Pursuit-Evasion,” _2023 National Aerospace and Electronics Conference (NAECON)_ , Dayton, OH, 2023a.  Accepted.
  * Oyler et al. [2016] Oyler, D. W., Kabamba, P. T., and Girard, A. R., “Pursuit-evasion games in the presence of obstacles,” _Automatica_ , Vol. 65, 2016, pp. 1–11.  10.1016/j.automatica.2015.11.018.
  * Zhou et al. [2016] Zhou, Z., Zhang, W., Ding, J., Huang, H., Stipanović, D. M., and Tomlin, C. J., “Cooperative pursuit with Voronoi partitions,” _Automatica_ , Vol. 72, 2016, pp. 64–72.  10.1016/j.automatica.2016.05.007.
  * Ibragimov [1998] Ibragimov, G. I., “A game of optimal pursuit of one object by several,” _Journal of Applied Mathematics and Mechanics_ , Vol. 62, 1998, pp. 187–192.  10.1016/S0021-8928(98)00024-0.
  * Weintraub et al. [2023b] Weintraub, I., Von Moll, A., Garcia, E., Casbeer, D., and Pachter, M., “Surveillance of a Faster Fixed-Course Target,” _Transactions on Aerospace & Electronic Systems_, 2023b.  10.1109/TAES.2023.3237129, none.
  * Von Moll et al. [2023] Von Moll, A., Fuchs, Z., Shishika, D., Maity, D., Dorothy, M., and Pachter, M., “Turret Escape Differential Game,” _Journal of Dynamics and Games_ , 2023.  10.3934/jdg.2023012, presented at the 19th ISDG.
  * Galyaev and Maslov [2013] Galyaev, A. A., and Maslov, E. P., “Evading a rotating detection zone on a plane,” _Journal of Computer and Systems Sciences International_ , Vol. 52, 2013, pp. 377–385.  10.1134/s1064230713030076.
  * Ivanov and Maslov [1993] Ivanov, M., and Maslov, E., “A problem of avoidance of a rotating segment,” _Computers & Mathematics with Applications_, Vol. 26, No. 6, 1993, pp. 67–75.  10.1016/0898-1221(93)90118-F.
  * Wächter and Biegler [2006] Wächter, A., and Biegler, L. T., “On the implementation of an interior-point filter line-search algorithm for large-scale nonlinear programming,” _Mathematical Programming_ , Vol. 106, 2006, pp. 25–57.  10.1007/s10107-004-0559-y.
  * Lubin et al. [2023] Lubin, M., Dowson, O., Garcia, J. D., Huchette, J., Legat, B., and Vielma, J. P., “JuMP 1.0: Recent improvements to a modeling language for mathematical optimization,” _Mathematical Programming Computation_ , 2023\.  10.1007/s12532-023-00239-3.
  * Zhang et al. [2022] Zhang, Y., Zhang, P., Wang, X., Song, F., Li, C., and Hao, J., “An open loop Stackelberg solution to optimal strategy for UAV pursuit-evasion game,” _Aerospace Science and Technology_ , Vol. 129, 2022, p. 107840.  10.1016/j.ast.2022.107840.