# Supplementary Information:
Two-photon emission from superluminal and accelerating index perturbations

Jamison Sloan1†, Nicholas Rivera2, John D. Joannopoulos2, and Marin Soljačić2
1 Department of Electrical Engineering and Computer Science, Massachusetts
Institute of Technology, Cambridge, MA 02139, USA
2 Department of Physics, Massachusetts Institute of Technology, Cambridge, MA
02139, USA $\dagger$ Corresponding author e-mail<EMAIL_ADDRESS>
###### Contents

  1. I General theoretical formalism
    1. I.1 Mode expansion
    2. I.2 Macroscopic QED
  2. II Cherenkov emission into graphene plasmons
  3. III Supplementary concept figure
  4. IV Derivations for example systems
    1. IV.1 Derivation of quantum Cherenkov radiation in a homogeneous medium
    2. IV.2 Derivation of quantum Cherenkov radiation for plasmons
    3. IV.3 Derivation of two-photon Synchrotron radiation
  5. V Green’s functions in layered media

In this supplement, we provide additional details in support of the main text
of this work. In Section I, we outline the theoretical framework used to
obtain our results. In Section II, we provide additional results about the
production of graphene plasmon pairs from a fast index perturbation moving
above the surface. In Section IV, we outline the derivations used to obtain
results for specific example systems discussed in this work. Section V
provides a summary of properties of Maxwell Green’s functions in layered
media, which are used in the derivations for the two-plasmon Cherenkov effect
in graphene.

##  I General theoretical formalism

Here, we elaborate on the theoretical foundations of two-photon emission in
time-dependent media which are summarized in the Methods section of the main
text. More details about the theoretical foundations of two-photon emission in
general optical media can be found in sloan2020casimir.

###  I.1 Mode expansion

We begin by deriving an expression for the rate of two-photon emission using
an eigenmode expansion for the quantized electric field operator. We apply
perturbation theory to calculate elements of the $S$-matrix, which yields
emission probabilities, as well as the quantum state of the system after
emission. The Hamiltonian of the bare field is given by

$H=\sum_{n}\hbar\omega_{n}\left(a_{n}^{\dagger}a_{n}+\frac{1}{2}\right)$ (S1)

where $a_{n}$ is the annihilation operator for a mode with label $n$. In
space, the modes are assumed to be orthogonalized under the inner product
$\int
d^{3}r\,\mathbf{F}_{m}^{*}(\mathbf{r})\cdot\varepsilon_{\text{bg}}(\mathbf{r})\cdot\mathbf{F}_{n}(\mathbf{r})=\delta_{mn}$.
Here, $\varepsilon_{\text{bg}}(\mathbf{r})$ is the permittivity of the
background structure, which is assumed here to be lossless and nondispersive.
Then, the interaction picture electric field operator is given via a sum over
eigenmodes $\mathbf{F}_{n}(\mathbf{r})$

$\mathbf{E}(\mathbf{r},t)=\sum_{n}\sqrt{\frac{\hbar\omega}{2\epsilon_{0}}}\left(\mathbf{F}_{n}(\mathbf{r})a_{n}e^{-i\omega_{n}t}+\mathbf{F}_{n}^{*}(\mathbf{r})a_{n}^{\dagger}e^{i\omega_{n}t}\right).$
(S2)

The interaction Hamiltonian, written in the interaction picture, is given by

$V(t)=-\varepsilon_{0}\int
d^{3}r\,\Delta\varepsilon_{ij}(\mathbf{r},t)E_{i}(\mathbf{r},t)E_{j}(\mathbf{r},t),$
(S3)

where we have used index notation. We consider initial states which are the
vacuum, and final states which have two photons. Thus we write
$\ket{i}=\ket{0}$, and $\ket{f}=\ket{m,n}$, where $m,n$ are the mode labels.
If the interaction Hamiltonian were time-independent, such a process would not
be possible, as the initial and final states clearly do not conserve energy.
However, the time dependence of the medium lifts this constraint, and the
medium change $\Delta\varepsilon$ can act as a classical pump field. The
transition amplitude at leading order in perturbation theory is given by

$S_{fi}(t)=-\frac{i}{\hbar}\int_{-\infty}^{t}dt^{\prime}\braket{f}{V(t^{\prime})}{i}.$
(S4)

Plugging in the initial and final states defined above, we find that

$S_{fi}(t)=\frac{i}{2}\sqrt{\omega_{n}\omega_{m}}\int
d^{3}r\int_{-\infty}^{t}dt^{\prime}\,\Delta\varepsilon_{ij}(\mathbf{r},t^{\prime})e^{i(\omega_{m}+\omega_{n})t^{\prime}}F_{mj}^{*}(\mathbf{r})F_{ni}^{*}(\mathbf{r})$
(S5)

In order to study the emission of asymptotic states which would be detected
away from the source, we can take $t\to\infty$. This allows us to group the
time integral, the exponential, and $\Delta\varepsilon$ as the Fourier
transform of the change in susceptibility:

$\Delta\varepsilon(\mathbf{r},\omega)\equiv\int_{-\infty}^{\infty}dt\,e^{i\omega
t}\Delta\varepsilon(\mathbf{r},t).$ (S6)

For clarity, we note that this quantity is _not_ related to the dispersion of
the modulated material, which is assumed here not to be significant. With Eq.
S6, the matrix element as $t\to\infty$ simplifies to

$S_{fi}^{t\to\infty}=\frac{i}{2}\sqrt{\omega_{n}\omega_{m}}\int
d^{3}r\,\Delta\varepsilon_{ij}(\mathbf{r},\omega_{n}+\omega_{m})F_{mj}^{*}(\mathbf{r})F_{ni}^{*}(\mathbf{r}).$
(S7)

To find the total probability $P$ from this matrix element, we take the
modulus squared and sum over final photon states $m$ and $n$:

$P=\frac{1}{2}\sum_{m,n}|S_{fi}|^{2}=\frac{1}{8}\sum_{m,n}\omega_{m}\omega_{n}\left|\int
d^{3}r\,\Delta\varepsilon_{ij}(\mathbf{r},\omega_{n}+\omega_{m})F_{mj}^{*}(\mathbf{r})F_{ni}^{*}(\mathbf{r})\right|^{2},$
(S8)

where the factor of two has come from the indistinguishability of photons
labeled with $m$ and $n$. This is the result shown in Eq. 2 of the main text.

Critically, we see that the strength of the matrix element for two particular
modes $m$ and $n$ depends on $\Delta\varepsilon$ evaluated at the sum of
frequencies $\omega_{m}+\omega_{n}$. In this sense,
$\Delta\varepsilon(\mathbf{r},\omega)$ directly controls the emitted spectrum.
In our examples discussed in the main text, we see clearly how the spatial
profile of $\Delta\varepsilon(\mathbf{r},\omega)$ leads directly to frequency
and momentum phase-matching rules for the emitted photon pairs.

###  I.2 Macroscopic QED

While the mode expansion version detailed above works well in systems which
can be idealized as lossless and dispersionless, treating problems of emission
in vacuum in arbitrary environments requires the formalism of macroscopic
quantum electrodynamics (MQED) scheel2009macroscopic; rivera2020light. In this
framework, the Hamiltonian of the bare electromagnetic field is

$H_{\text{EM}}=\int_{0}^{\infty}d\omega\int
d^{3}r\,\hbar\omega\,\mathbf{f}^{\dagger}(\mathbf{r},\omega)\cdot\mathbf{f}(\mathbf{r},\omega),$
(S9)

where $\mathbf{f}^{(\dagger)}(\mathbf{r},\omega)$ is the annihilation
(creation) operator for a quantum harmonic oscillator at position $\mathbf{r}$
and frequency $\omega$. In such a medium, the electric field operator in the
interaction picture is given as

$\mathbf{E}(\mathbf{r},t)=i\sqrt{\frac{\hbar}{\pi\varepsilon_{0}}}\int_{0}^{\infty}d\omega\frac{\omega^{2}}{c^{2}}\int
d^{3}r^{\prime}\sqrt{\operatorname{Im}\varepsilon_{\text{bg}}(\mathbf{r}^{\prime},\omega)}\left(G(\mathbf{r},\mathbf{r}^{\prime},\omega)\mathbf{f}(\mathbf{r}^{\prime},\omega)e^{-i\omega
t}-\text{h.c.}\right).$ (S10)

Here, $G(\mathbf{r},\mathbf{r}^{\prime},\omega)$ is the electromagnetic
Green’s function of the background which satisfies
$\left(\nabla\times\nabla\times-\varepsilon_{\text{bg}}(\mathbf{r},\omega)\frac{\omega^{2}}{c^{2}}\right)G(\mathbf{r},\mathbf{r}^{\prime},\omega)=\delta(\mathbf{r}-\mathbf{r}^{\prime})I$,
where $I$ is the $3\times 3$ identity matrix. Additionally, we assume for
generality that the modulated material can also be dispersive, and is
described through a generalized linear response function
$\Delta\varepsilon(\mathbf{r},t,t^{\prime})$ (see sloan2020casimir for
details), so that the interaction Hamiltonian becomes

$V(t)=-\varepsilon_{0}\int
d^{3}r\,dt^{\prime}\,\Delta\varepsilon_{ij}(\mathbf{r},t,t^{\prime})E_{j}(\mathbf{r},t^{\prime})E_{i}(\mathbf{r},t).$
(S11)

We now take the initial and final states of the scattering problem to to be
$\ket{i}=\ket{0}$ and
$\ket{f}=\ket{m\mathbf{r}_{1}\omega_{1},n\mathbf{r}_{2}\omega_{2}}$. Then the
scattering matrix element is given as

$\displaystyle S_{fi}$
$\displaystyle=-\frac{i}{\hbar}\int_{-\infty}^{\infty}dt\braket{f}{V(t)}{i}$
(S12) $\displaystyle=-\frac{i(\omega_{1}\omega_{2})^{2}}{\pi c^{4}}\int
d^{3}x\int_{-\infty}^{\infty}dt\,dt^{\prime}\,e^{i\omega_{1}t^{\prime}}e^{i\omega_{2}t}\sqrt{\operatorname{Im}\varepsilon_{\text{bg}}(\mathbf{r}_{1},\omega_{1})\operatorname{Im}\varepsilon_{\text{bg}}(\mathbf{r}_{2},\omega_{1})}$
(S13) $\displaystyle\hskip
56.9055pt\times\Delta\varepsilon_{ij}(\mathbf{x},t,t^{\prime})G_{jm}^{*}(\mathbf{x},\mathbf{r}_{1},\omega_{1})G_{in}^{*}(\mathbf{x},\mathbf{r}_{2},\omega_{2}).$
(S14)

From the S-matrix, we can sum over final states to obtain the probability of
two-photon emission:

$\displaystyle P$ $\displaystyle=\frac{1}{2}\sum_{m,n}\int
d^{3}r_{1}\,d^{3}r_{2}\int_{0}^{\infty}d\omega_{1}\,d\omega_{2}\,|S_{fi}|^{2}$
(S15) $\displaystyle=\frac{1}{2}\sum_{m,n}\int
d^{3}r_{1}\,d^{3}r_{2}\,\int_{0}^{\infty}d\omega_{1}\,d\omega_{2}\Bigg{|}\frac{(\omega_{1}\omega_{2})^{2}}{\pi
c^{4}}\int
d^{3}x\,\sqrt{\operatorname{Im}\varepsilon_{\text{bg}}(\mathbf{r}_{1},\omega_{1})\operatorname{Im}\varepsilon_{\text{bg}}(\mathbf{r}_{2},\omega_{1})}$
(S16) $\displaystyle\hskip
85.35826pt\Delta\varepsilon_{ij}(\mathbf{x},\omega_{2},-\omega_{1})G_{jm}^{*}(\mathbf{x},\mathbf{r}_{1},\omega_{1})G_{in}^{*}(\mathbf{x},\mathbf{r}_{2},\omega_{2})\Bigg{|}^{2}.$
(S17)

Here, we have expressed $\Delta\varepsilon(\mathbf{r},t,t^{\prime})$ in terms
of its two-time Fourier transform, which we define as

$\Delta\varepsilon(\mathbf{r},\omega,\omega^{\prime})\equiv\int_{-\infty}^{\infty}dt\,dt^{\prime}\,e^{i\omega
t}\Delta\varepsilon(\mathbf{r},t,t^{\prime})e^{-i\omega^{\prime}t^{\prime}}.$
(S18)

Finally we make use of the following identity scheel2009macroscopic to
simplify the result:

$\frac{\omega^{2}}{c^{2}}\int
d^{3}x\operatorname{Im}\varepsilon_{\text{bg}}(\mathbf{x},\omega)G(\mathbf{r},\mathbf{x},\omega)G^{\dagger}(\mathbf{r}^{\prime},\mathbf{x},\omega)=\operatorname{Im}G(\mathbf{r},\mathbf{r}^{\prime},\omega)$
(S19)

R elabeling variables, we can write a final expression for the probability as

$\begin{split}P=\frac{1}{2\pi^{2}c^{4}}\int_{0}^{\infty}&d\omega\,d\omega^{\prime}\int
d^{3}r\,d^{3}r^{\prime}\,(\omega\omega^{\prime})^{2}\\\
&\times\operatorname{Tr}\left[\Delta\varepsilon(\mathbf{r},\omega,-\omega^{\prime})\operatorname{Im}G(\mathbf{r},\mathbf{r}^{\prime},\omega^{\prime})\Delta\varepsilon^{\dagger}(\mathbf{r}^{\prime},\omega,-\omega^{\prime})\operatorname{Im}G(\mathbf{r}^{\prime},\mathbf{r},\omega)\right].\end{split}$
(S20)

We note that in the limit of low loss in the Maxwell equations, the Imaginary
part of the Green’s function can be written in terms of the eigenmodes as

$\operatorname{Im}G(\mathbf{r},\mathbf{r}^{\prime},\omega)=\pi
c^{2}\sum_{n}\frac{1}{2\omega_{n}}\mathbf{F}_{n}^{*}(\mathbf{r})\mathbf{F}_{n}(\mathbf{r}^{\prime})\delta(\omega-\omega_{n}).$
(S21)

Using this correpsondence, Eq. S8 can be converted into Eq. S20, establishing
a firm connection between the eigenmode expansion and the MQED approach. We
briefly mention two important limits of Eq. S20. If the modulation to the
material takes place in a region which can be regarded as nondispersive, then
$\Delta\varepsilon(\mathbf{r},\omega,-\omega^{\prime})=\Delta\varepsilon(\mathbf{r},\omega+\omega^{\prime})$,
where the latter form is that described by Eq. S6. In this case, we have

$\begin{split}P=\frac{1}{2\pi^{2}c^{4}}\int_{0}^{\infty}&d\omega\,d\omega^{\prime}\int
d^{3}r\,d^{3}r^{\prime}\,(\omega\omega^{\prime})^{2}\\\
&\times\operatorname{Tr}\left[\Delta\varepsilon(\mathbf{r},\omega+\omega^{\prime})\operatorname{Im}G(\mathbf{r},\mathbf{r}^{\prime},\omega^{\prime})\Delta\varepsilon^{\dagger}(\mathbf{r}^{\prime},\omega+\omega^{\prime})\operatorname{Im}G(\mathbf{r}^{\prime},\mathbf{r},\omega)\right].\end{split}$
(S22)

This enables one to describe a situation where a nondispersive modulation
emits into a dispersive and lossy region nearby, which is the case for the
graphene plasmon system discussed in Section II. The other important limit is
when the modulation is assumed to be isotropic, in which case
$\Delta\varepsilon_{ij}(\mathbf{r},\omega,-\omega^{\prime})=\delta_{ij}\Delta\varepsilon(\mathbf{r},\omega,-\omega^{\prime})$.
In this case, we have

$\begin{split}P=\frac{1}{2\pi^{2}c^{4}}\int_{0}^{\infty}&d\omega\,d\omega^{\prime}\int
d^{3}r\,d^{3}r^{\prime}\,(\omega\omega^{\prime})^{2}\Delta\varepsilon(\mathbf{r},\omega,-\omega^{\prime})\Delta\varepsilon^{\dagger}(\mathbf{r}^{\prime},\omega,-\omega^{\prime})\\\
&\times\operatorname{Tr}\left[\operatorname{Im}G(\mathbf{r},\mathbf{r}^{\prime},\omega^{\prime})\operatorname{Im}G(\mathbf{r}^{\prime},\mathbf{r},\omega)\right].\end{split}$
(S23)

##  II Cherenkov emission into graphene plasmons

In a homogeneous medium of constant index $n$, interactions between the time-
variation $\Delta\varepsilon$ and the electric field $E$ are quite weak. The
coupling strength in the interaction Hamiltonian (Eq. S3) is, roughly
speaking, proportional to the magnitude of $\Delta\varepsilon$ as well as the
fluctuations of the vacuum electromagnetic field $\Braket{E^{2}}$ over the
region of index perturbation. A key paradigm in nanophotonics is the ability
to engineer electromagnetic excitations which have much higher field strengths
at the single photon level. These techniques have been applied to enhance
spontaneous emission purcell1946purcell, Casimir forces rodriguez2011casimir,
near-field radiative heat transfer volokitin2007near, and more. One
particularly interesting system for realizing these effects are thin materials
which support highly confined polaritons atwater2007promise;
basov2016polaritons; basov2017towards; low2017polaritons; iranzo2018probing;
ni2018fundamental; caldwell2013low; caldwell2014sub; caldwell2015low, thus
enabling strong light-matter interactions.

To demonstrate that this concept generalizes to our two-photon generation
mechanism, we show how a fast index perturbation near a sheet of graphene
jablan2009plasmonics; koppens2011graphene produces entangled pairs of graphene
plasmons with much higher efficiencies than in a homogeneous medium. While we
focus on the graphene example, this concept readily generalizes to other
planar systems (e.g. phonon-polaritons, exciton-polaritons, magnon-
polaritons). We assume that a substrate experiences a gaussian perturbation,
with a direction of propagation $x$ (Fig. S1a). Additionally, the soft wall is
given a width $W$ in the $y$ direction, which is assumed to be large compared
to the relevant wavelengths, and assumed to extend uniformly over $0<z<d$.
Mathematically, this is defined as

$\Delta\varepsilon(x,t)=\delta\varepsilon\,\exp\left[-\frac{(x-v_{0}t)^{2}}{2\sigma^{2}}\right],$
(S24)

for $0<z<d$, and zero otherwise. The width of the wall in the propagation
direction is denoted by $\sigma$, similarly to the free-space Cherenkov
emission discussed in the main text. While we depict the perturbation to be
free standing, this would in practice occur in a nonlinear substrate or
superstrate, leading to effects similar in magnitude. This index perturbation
then travels parallel to the graphene surface at some fraction the speed of
light $\beta$. In the case of a uniform medium discussed in the main text, the
“moving wall” of permittivity variation emits photon pairs into the
surrounding space. Here, the relevant electromagnetic excitations are surface
plasmon polaritons. Consequently, an index perturbation with sufficiently high
velocity will emit pairs of entangled plasmons into the graphene sheet. Since
the velocity of graphene plasmons can easily go below $c/100$, the
superluminal condition is easier to reach than in a homogeneous medium. Using
our formalism, we obtain a general expression (see Supplementary Section IV.2
for details) for the two-plasmon emission rate per unit frequency $\omega$ and
in-plane angle $\theta$:

$\begin{split}\frac{1}{W}\frac{d\Gamma}{d\omega d\theta}&=\frac{\beta
c}{(2\pi)^{3}}\frac{\delta\varepsilon^{2}}{\omega_{0}^{2}}\int
d\omega^{\prime}\int
dq\,q^{2}q^{\prime}\,e^{-2(\omega+\omega^{\prime})^{2}/\omega_{0}^{2}}\\\
&\times[g(q+q^{\prime})]^{2}(\hat{q}\cdot\hat{q}^{\prime}-1)^{2}\operatorname{Im}r_{p}(\omega,q)\operatorname{Im}r_{p}(\omega^{\prime},q^{\prime}).\end{split}$
(S25)

Here, the geometry (Fig. S1a) is parameterized as
$q^{\prime}=\sqrt{q_{\perp}^{\prime 2}+q_{\parallel}^{\prime 2}}$, where
$q_{\perp}^{\prime}=-q_{\perp}$,
$q_{\parallel}^{\prime}=(\omega+\omega^{\prime})/\beta c-q_{\parallel}$, and
$(q_{\parallel},q_{\perp})=(q\cos\theta,q\sin\theta)$. Additionally,
$g(q)=(1-e^{-qd})/q$ is a geometric factor which comes from the finite
thickness of the film in the $z$-direction. The frequency scale is set by
$\omega_{0}\equiv\sqrt{2}v/\sigma$, using the same definition as for the free-
space Cherenkov emission. Finally, $\operatorname{Im}r_{p}(\omega,q)$ is the
imaginary part of the p-polarized reflection coefficient for the interface
above the material novotny2012principles. This expression can be easily
modified to account for a different distribution of index perturbation in the
$z$-direction. The result of Eq. S25 can be separated into interpretable
pieces. The exponential factor sets $\omega_{0}$ as the scale of frequencies
which can be excited by the pulse, similarly to the case of two-photon
Cherenkov radiation in a homogeneous medium. The factor
$(\hat{q}\cdot\hat{q}^{\prime}-1)^{2}$ is responsible for controlling the
angular relation between the two plasmons. The imaginary part of the
reflectivity peaks where modes are present. There are two factors of
$\operatorname{Im}r_{p}$ — one for each plasmon — so that the total expression
is peaked for values of $(\omega,q)$ and $(\omega^{\prime},q^{\prime})$ which
both lie along the dispersion. In the lossless limit,
$\operatorname{Im}r_{p}(\omega,q)\propto q\delta\left(q-q(\omega)\right)$, so
both points lie exactly on the dispersion $q(\omega)$. With loss present, this
condition is relaxed so that the two frequencies and wavevectors instead lie
near the dispersion, set by its linewidth. Models for the optical properties
of graphene are taken from koppens2011graphene; jablan2009plasmonics;
christensen2017classical.

Figure S1: Emission of entangled plasmon pairs by moving index perturbations
near graphene. (a) Schematic of an index perturbation $\Delta\varepsilon$
moving at constant velocity $v_{0}$ near a sheet of graphene. (b) Angular
correlation of photon pairs emitted in lossless Drude model. Contours show
values of constant $\hbar\omega\beta/(4\alpha E_{\text{F}})$, where $\alpha$
is the fine structure constant. The ratio between frequencies
$\omega/\omega^{\prime}$ is shown on the colorbar below, indicating at any
point on a contour how one frequency relates to the other. (c) Frequency
distribution $(1/W)d\Gamma/d\omega$ of emitted plasmon pairs for an index
perturbation $\delta\varepsilon/\varepsilon=10^{-3}$ of velocity $\beta c$
which has dimensions $\sigma=400$ nm, $d=1\,\mu$m. Graphene with three
different doping levels $E_{\text{F}}=\\{0.2,0.5,0.8\\}$ eV are shown, and
each with three different conductivity models. A relaxation constant
$\tau=10^{-14}$ s is used, and dispersion of $\tau$ may be taken into account
using the methods provided in jablan2009plasmonics. The electron-hole pair
excitation region at $\hbar\omega=2E_{\text{F}}$ is marked with a dotted line.
Nonlocal corrections are small when the frequency spectrum does not
substantially enter this region. At high velocity and low doping levels (right
panel for 0.2 eV), the nonlocal corrections are substantial, and account for
the production of electron-hole pairs.

Similarly to a homogeneous medium, phase matching and momentum conservation
place constraints on the emitted angles and frequencies of the two plasmons.
Fig. S1b shows the kinematic constraints between the two emitted plasmons
assuming a simplified lossless Drude model of the surface conductivity (i.e.
assuming that the plasmon has a dispersion $\omega\propto\sqrt{q}$). The two
axes show the angle of emission of each of the two plasmons measured with
respect to the direction of propagation. If the first plasmon is emitted at an
angle which satisfies, $0\leq\theta\leq\pi$, then the second must satisfy
$\pi\leq\theta^{\prime}\leq 2\pi$ in order to obey total momentum conservation
in the direction transverse to propagation. The contours show the kinematic
constraints for different values of the dimensionless constant
$\hbar\omega\beta/(4\alpha E_{\text{F}})$. With higher frequency and velocity,
or lower Fermi level, the contours are pushed toward the line
$\theta^{\prime}=\theta+\pi$. Additionally, the colorbar shows the ratio of
plasmon frequencies $\omega/\omega^{\prime}$, where $\omega$ and
$\omega^{\prime}$ are respectively the lower and higher of the frequencies,
indicating the frequency-angle entanglement between the two plasmons.

Figure S1c shows the width-normalized angular spectrum $(1/W)d\Gamma/d\omega$
for 2-plasmon emission for two different velocities $\beta=\\{0.4,0.8\\}$. We
have computed these results for three different optical conductivity models
for graphene (Drude, local relaxation time approximation, and nonlocal) in
order to gain insights as to how dissipation, dispersion, intraband
transitions, and nonlocality impact the production of plasmon pairs. Different
doping levels $E_{\text{F}}$ for the graphene sample are also considered. For
the more highly doped samples considered (i.e. $E_{\text{F}}=0.5,0.8$ eV), the
Drude conductivity model is seen to provide a good approximation to the more
physically accurate conductivity models. However, intraband transitions and
nonlocal effects become important in samples with less doping. Notably, for
sufficiently high velocities ($\beta=0.8$), and relatively low doping levels
($E_{\text{F}}=0.2$ eV), nonlocal corrections to the conductivity become
important. In this regime, particle-hole pairs in the region
$\hbar\omega>2E_{\text{F}}$ are produced. We now comment on the total rate of
this process, making a comparison to the free-space case discussed in the main
text. To compare the 3D and 2D processes on equal footing, we consider the
pair production rate _per unit volume_ of dielectric which has its
permittivity changed. Accounting for the substrate height and pulse width
$\sigma$ in both cases, we have $\Gamma/V\approx 10^{4}\,$ s${}^{-1}\,\mu$m-3
for the homogeneous medium, and $\Gamma/V\approx 2\times 10^{8}\,$
s${}^{-1}\,\mu$m-3 for plasmons. This corresponds to a four order of magnitude
increase in efficiency, showing that 2D polaritonic platforms such as graphene
may hold promising opportunities for the observation and utilization of this
two-photon process. The ability to electrically tune the dispersion to change
the angular and frequency correlations of emission is also an attractive
feature of graphene-based systems.

##  III Supplementary concept figure

Figure S2: Concept figure for for colliding ultrashort light pulses with
nonlinear structures to create accelerating index perturbations. In this
particular example, we show a spring-shaped structure which is made of a
nonlinear material. When a short pulse collides with the structure, the
intersection creates an index perturbation which follows the helical
trajectory defined by the structure. Such an index perturbation can travel
superluminally, and realize acceleration with a periodicity as described in
the “synchrotron” geometry (main text Fig. 3), thus emitted entangled photon
pairs.

##  IV Derivations for example systems

###  IV.1 Derivation of quantum Cherenkov radiation in a homogeneous medium

Here, we provide a derivation of main text Eq. 4 using the mode expansion
method outlined in the first section. We use the index perturbation

$\Delta\varepsilon(\mathbf{r},t)=\delta\varepsilon\,\exp\left[-\frac{(z-v_{0}t)^{2}}{2\sigma^{2}}\right].$
(S26)

To compute the emission rate, we compute the Fourier transform associated with
this index change, defined as
$\Delta\varepsilon(z,\omega)=\int_{-\infty}^{\infty}dt\,e^{i\omega
t}\Delta\varepsilon(z,t)$. This is a Gaussian integral, given as

$\Delta\varepsilon(z,\omega)=\delta\varepsilon\,e^{-z^{2}/2\sigma^{2}}\sqrt{\frac{\pi}{(v_{0}^{2}/2\sigma^{2})}}\exp\left[\frac{(i\omega+v_{0}z/\sigma^{2})^{2}}{4(v_{0}^{2}/2\sigma^{2})}\right]=\delta\varepsilon\sqrt{\frac{4\pi}{\omega_{0}^{2}}}e^{-\omega^{2}/\omega_{0}^{2}}e^{i(\omega/v_{0})z},$
(S27)

where in the second step, we have simplified, and defined
$\omega_{0}^{2}\equiv 2v_{0}^{2}/\sigma^{2}$. We see that in frequency space,
this perturbation is a plane wave, much like the current distribution of a
free electron which produces Cherenkov radiation. In a homogeneous medium of
constant nondispersive background index $n=\sqrt{\varepsilon}$, the eigenmodes
are expressed as

$\mathbf{F}_{k}(\mathbf{r})=\frac{\hat{\epsilon}_{k}}{\sqrt{Vn^{2}}}e^{i\mathbf{k}\cdot\mathbf{r}},$
(S28)

where $\hat{\epsilon}_{k}$ is a polarization, and $V$ is the mode
normalization volume . These modes obey the dispersion relation
$\omega_{k}=ck/n$. Then, following Eq. S7, the $S$-matrix element which
connects the initial ground state $\ket{i}=\ket{0}$ to a 2-photon final state
$\ket{f}=\ket{k,k^{\prime}}$ is

$\begin{split}S_{fi}=\frac{2}{V}\left(\frac{\delta\varepsilon}{\varepsilon}\right)(2\pi)^{3}\frac{\sqrt{\pi\omega_{k}\omega_{k^{\prime}}}}{\omega_{0}}&\exp\left[-\frac{(\omega_{k}+\omega_{k^{\prime}})^{2}}{\omega_{0}^{2}}\right]\left(\hat{\epsilon}_{k}^{*}\cdot\hat{\epsilon}_{k^{\prime}}^{*}\right)\\\
&\times\delta^{(2)}(\mathbf{k}_{\perp}+\mathbf{k}_{\perp}^{\prime})\delta\left(\frac{\omega_{k}+\omega_{k^{\prime}}}{v_{0}}-(k_{z}+k_{z}^{\prime})\right).\end{split}$
(S29)

Here, the delta functions make momentum conservation and phase matching
manifest. The probability of the process is given by
$P=\frac{1}{2}\sum_{k,k^{\prime}}|S_{fi}|^{2}$, where the factor of $1/2$
accounts for the indistinguishability of the two photons. Defining the
transition rate as the probability per unit time, we can write the transition
rate $\Gamma$ per unit area $A$ of the wall as

$\frac{\Gamma}{A}=\frac{\pi}{2}\left(\frac{\delta\varepsilon}{\varepsilon}\right)^{2}\frac{(2\pi)^{3}v_{0}}{\omega_{0}^{2}V^{2}}\sum_{k,k^{\prime}}|\hat{\epsilon}\cdot\hat{\epsilon}^{\prime}|^{2}\omega_{k}\omega_{k^{\prime}}e^{-\frac{2(\omega_{k}+\omega_{k^{\prime}})^{2}}{\omega_{0}^{2}}}\delta^{(2)}(\mathbf{k}_{\perp}+\mathbf{k}_{\perp}^{\prime})\delta\left(\frac{\omega_{k}+\omega_{k^{\prime}}}{v_{0}}-(k_{z}+k_{z}^{\prime})\right),$
(S30)

where we have used $\delta(\omega=0)=T/2\pi$ to define the time of
interaction, and $\delta(\mathbf{k}_{\perp}=0)=A/(2\pi)^{2}$ to define the
area of the semi-infinite moving wall. We also use the shorthand
$\hat{\epsilon}$ to denote $\hat{\epsilon}_{\k}$, and likewise for the primed
version. The polarization factor in the rate indicates that the photons prefer
to be oriented in the same direction, and cannot be emitted orthogonally. The
mode sums are converted to a continuum as

$\sum_{k,k^{\prime}}\to\frac{1}{V^{2}}\sum_{\epsilon,\epsilon^{\prime}}\int\frac{d^{3}k}{(2\pi)^{3}}\frac{d^{3}k^{\prime}}{(2\pi)^{3}}.$
(S31)

Then, the polarization sum is handled as

$\sum_{\epsilon,\epsilon^{\prime}}|\hat{\epsilon}\cdot\hat{\epsilon}^{\prime}|^{2}=1+(\hat{k}\cdot\hat{k}^{\prime})^{2}.$
(S32)

For the transverse delta function, we can write

$\delta^{(2)}(\mathbf{k}_{\perp}+\mathbf{k}_{\perp}^{\prime})=\frac{1}{k_{\perp}}\delta(k_{\perp}-k_{\perp}^{\prime})\delta(\phi-\phi^{\prime}-\pi).$
(S33)

We also have $k_{z}=k\cos\theta$ and likewise for the primed version. This
allows us to write

$\begin{split}\frac{\Gamma}{A}=v_{0}\frac{\pi}{2}\left(\frac{\delta\varepsilon}{\varepsilon}\right)^{2}&\frac{1}{\omega_{0}^{2}(2\pi)^{2}}\int_{0}^{\infty}dk\,k^{2}\,dk^{\prime}\,k^{\prime
2}\int_{-1}^{1}d(\cos\theta)\,d(\cos\theta^{\prime})e^{-\frac{2(\omega_{k}+\omega_{k^{\prime}})^{2}}{\omega_{0}^{2}}}\omega_{k}\omega_{k^{\prime}}\\\
&\times\left(1+(\hat{k}\cdot\hat{k}^{\prime})^{2}\right)\frac{\delta(k\sin\theta-k^{\prime}\sin\theta^{\prime})}{k\sin\theta}\delta\left(\frac{\omega_{k}+\omega_{k^{\prime}}}{v_{0}}-(k\cos\theta+k^{\prime}\cos\theta^{\prime})\right)\end{split}$
(S34)

The two delta-function conditions must be satisfied simultaneously. These
equations (after using the dispersion relation $\omega=ck/n$ to convert $k$’s
to $\omega$’s) are

$\displaystyle\frac{\omega+\omega^{\prime}}{n\beta}$
$\displaystyle=\omega\cos\theta+\omega^{\prime}\cos\theta$ (S35)
$\displaystyle\omega\sin\theta$
$\displaystyle=\omega^{\prime}\sin\theta^{\prime}$ (S36)

There is a trivial solution which corresponds to the onset of the radiation.
This occurs when $\cos\theta=\cos\theta^{\prime}=1/n\beta$. Consequently,
$\omega=\omega^{\prime}$ here. Thus this solution corresponds to 2 photons of
the same frequency produced at the classical Cherenkov angle. The pair of
equations, however, is singular in $\omega,\omega^{\prime}$. Taking the
determinant of the system and setting to zero reveals an underlying constraint
on the angles, which is

$\sin\theta+\sin\theta^{\prime}=n\beta\sin(\theta+\theta^{\prime})$ (S37)

Performing the delta function integrals, we find a final differential decay
rate per unit frequency and angle (and normalized to the area $A$) as

$\frac{1}{A}\frac{d\Gamma}{d\omega
d\theta}=\frac{v_{0}}{8\pi\omega_{0}^{2}}\left(\frac{\delta\varepsilon}{\varepsilon}\right)^{2}\left(\frac{n}{c}\right)^{3}(\omega\omega^{\prime})^{2}e^{-2(\omega+\omega^{\prime})^{2}/\omega_{0}^{2}}\frac{\left(1+(\cos\theta\cos\theta^{\prime}-1)^{2}\right)}{\left(\frac{\cot\theta}{n\beta}-\frac{1}{\sin\theta}\right)},$
(S38)

which is the result quoted in main text Eq. 4.

###  IV.2 Derivation of quantum Cherenkov radiation for plasmons

Here, we provide details of the derivation of Eq. S25. We consider an index
perturbation

$\Delta\varepsilon(x,t)=\delta\varepsilon\,\exp\left[-\frac{(x-v_{0}t)^{2}}{2\sigma^{2}}\right],$
(S39)

for $0<z<d$, and zero otherwise. The Fourier transform follows a form
identical to that of Eq. S27. Then, implementing Eq. S20 with the form of the
Green’s function shown in Eq. S59, we obtain

$\begin{split}\frac{dP}{d\omega
d\omega^{\prime}}=\frac{\delta\varepsilon^{2}}{\omega_{0}^{2}}e^{-2(\omega+\omega^{\prime})^{2}/\omega_{0}^{2}}(2\pi)WL&\int\frac{d^{2}q}{(2\pi)^{2}}\frac{d^{2}q^{\prime}}{(2\pi)^{2}}\left(\frac{1-e^{-(q+q^{\prime})d}}{q+q^{\prime}}\right)^{2}qq^{\prime}\operatorname{Im}r_{p}(\omega,q)\operatorname{Im}r_{p}(\omega^{\prime},q^{\prime})\\\
&\times\operatorname{Tr}\left[\hat{\epsilon}_{q^{\prime}}^{*}\otimes\hat{\epsilon}_{q^{\prime}}\cdot\hat{\epsilon}_{q}\otimes\hat{\epsilon}_{q}^{*}\right]\delta(q_{\perp}+q_{\perp}^{\prime})\delta\left(\frac{\omega_{q}+\omega_{q^{\prime}}}{v_{0}}-(q_{\parallel}+q_{\parallel}^{\prime})\right).\end{split}$
(S40)

here, $L$ is a length factor which comes from the delta function in the
parallel direction. It is related to the time factor $T$ by $L=\beta cT$.
Performing the integral $\int d^{2}q^{\prime}=\int
dq_{\parallel}^{\prime}\,dq_{\perp}^{\prime}$, and using

$\operatorname{Tr}\left[\hat{\epsilon}_{q^{\prime}}^{*}\otimes\hat{\epsilon}_{q^{\prime}}\cdot\hat{\epsilon}_{q}\otimes\hat{\epsilon}_{q}^{*}\right]=\frac{1}{4}\left(\hat{q}\cdot\hat{q}^{\prime}-1\right)^{2}$
(S41)

we obtain the general result

$\begin{split}\frac{1}{W}\frac{d\Gamma}{d\omega d\theta}=\frac{\beta
c}{(2\pi)^{3}}\frac{\delta\varepsilon^{2}}{\omega_{0}^{2}}\int
d\omega^{\prime}\int
dq\,q^{2}q^{\prime}\,e^{-2(\omega+\omega^{\prime})^{2}/\omega_{0}^{2}}&\left(\frac{1-e^{-(q+q^{\prime})d}}{q+q^{\prime}}\right)^{2}(\hat{q}\cdot\hat{q}^{\prime}-1)^{2}\\\
&\times\operatorname{Im}r_{p}(\omega,q)\operatorname{Im}r_{p}(\omega^{\prime},q^{\prime}),\end{split}$
(S42)

where we have $q^{\prime}=\sqrt{q_{\perp}^{\prime 2}+q_{\parallel}^{\prime
2}}$, where $q_{\perp}^{\prime}=-q_{\perp}$,
$q_{\parallel}^{\prime}=(\omega+\omega^{\prime})/\beta c-q_{\parallel}$, and
$(q_{\parallel},q_{\perp})=(q\cos\theta,q\sin\theta)$. The relation between
$q^{\prime}$ and $q$ in terms of emission frequencies and angles is thus given
by

$q^{\prime 2}=q^{2}\sin^{2}\theta+\left(\frac{\omega+\omega^{\prime}}{\beta
c}-q\cos\theta\right)^{2}.$ (S43)

For numerical implementation, we also note that the overlap between wave
vectors can be expressed in terms of components as

$\hat{q}\cdot\hat{q}^{\prime}=\frac{q_{\parallel}q_{\parallel}^{\prime}+q_{\perp}q_{\perp}^{\prime}}{qq^{\prime}}.$
(S44)

###  IV.3 Derivation of two-photon Synchrotron radiation

Here, we provide details of the derivation of Eq. 6 of the main text, for two-
photon synchrotron radiation. For an index perturbation which follows a
circular trajectory, we use

$\Delta\varepsilon(\mathbf{r},t)=\delta\varepsilon\,e^{-\frac{1}{2\sigma^{2}}(\rho^{2}+z^{2}+R^{2})}e^{\frac{R\rho}{\sigma^{2}}\cos(\theta-\Omega
t)},$ (S45)

where $(\rho,\theta,z)$ are cylindrical coordinates, $\Omega$ is the angular
frequency of precession, $\sigma$ is the Gaussian width, and $R$ is the radius
of the circle. Throughout this derivation, we will make use of the Jacobi-
Anger expansion, shown in two forms here:

$\displaystyle e^{iz\cos\theta}$
$\displaystyle=\sum_{n=-\infty}^{\infty}e^{in\theta}i^{n}J_{n}(z)$ (S46)
$\displaystyle e^{z\cos\theta}$
$\displaystyle=\sum_{n=-\infty}^{\infty}e^{in\theta}I_{n}(z)$ (S47)

Here, $I_{n}(x)=i^{-n}J_{n}(ix)$ is a modified Bessel function of the first
kind, in terms of the Bessel function of the first kind $J_{n}(x)$. Using the
second form, we can take the Fourier transform of the perturbation as

$\Delta\varepsilon(\mathbf{r},\omega)=\delta\varepsilon\,e^{-\frac{1}{2\sigma^{2}}(\rho^{2}+z^{2}+R^{2})}\sum_{n=-\infty}^{\infty}e^{-in\theta}I_{n}\left(\frac{R\rho}{\sigma^{2}}\right)2\pi\delta(\omega-n\Omega),$
(S48)

Then, to compute the scattering matrix element, we have that

$\begin{split}\int
d^{3}r\,\Delta\varepsilon(\mathbf{r},\omega_{k}+\omega_{k^{\prime}})\mathbf{F}_{k}^{*}(\mathbf{r})\cdot\mathbf{F}_{k^{\prime}}^{*}(\mathbf{r})=\frac{\delta\varepsilon}{\varepsilon
V}(\hat{\epsilon}_{k}^{*}\cdot\hat{\epsilon}_{k^{\prime}}^{*})&(2\pi)^{5/2}\sigma^{3}e^{-(\mathbf{k}+\mathbf{k}^{\prime})^{2}\sigma^{2}/2}\\\
&\times\sum_{m=-\infty}^{\infty}J_{m}(KR)\delta(\omega_{k}+\omega_{k^{\prime}}-m\Omega).\end{split}$
(S49)

Here, $K=\sqrt{(k_{x}+k_{x}^{\prime})^{2}+(k_{y}+k_{y}^{\prime 2})^{2}}$ is
the net in-plane wavevector of emission. As an intermediate step, we have made
use of the integral

$\int_{0}^{\infty}d\rho\,\rho
J_{m}(K\rho)I_{m}\left(\frac{R\rho}{\sigma^{2}}\right)e^{-\frac{\rho^{2}}{2\sigma^{2}}}=\sigma^{2}e^{-\frac{\sigma^{2}K^{2}}{2}}e^{\frac{R^{2}}{2\sigma^{2}}}J_{m}(KR)$
(S50)

By squaring the matrix element, and summing over final two-photon states, the
probability of two-photon emission from this system per unit time $T$ can be
written as

$\begin{split}\frac{P}{T}=\left(\frac{\delta\epsilon}{\epsilon}\right)\frac{\sigma^{6}}{32\pi^{2}}\sum_{m=-\infty}^{\infty}\int\frac{d^{3}k}{(2\pi)^{3}}\frac{d^{3}k^{\prime}}{(2\pi)^{3}}&\omega_{k}\omega_{k^{\prime}}\left(1+(\hat{\mathbf{k}}\cdot\hat{\mathbf{k}}^{\prime})^{2}\right)\\\
&\times
e^{-\sigma^{2}\left(\mathbf{k}+\mathbf{k}^{\prime}\right)^{2}}J_{m}^{2}(KR)\delta(\omega_{k}+\omega_{k^{\prime}}-m\Omega)\end{split}$
(S51)

After some algebra, we find that the differential rate of down-conversion per
unit angle $\theta$ and frequency $\omega$ into each harmonic $m$ is given as:

$\frac{d\Gamma_{m}}{d\omega
d\theta}=\left(\frac{\delta\varepsilon}{\varepsilon}\right)^{2}\frac{\sigma^{6}n^{6}\omega^{3}(m\Omega-\omega)^{3}\sin\theta}{16\pi
c^{6}}\int
d\Omega^{\prime}\,\left(1+(\hat{\mathbf{k}}\cdot\hat{\mathbf{k}}^{\prime})^{2}\right)e^{-\sigma^{2}\left(\mathbf{k}+\mathbf{k}^{\prime}\right)^{2}}J_{m}^{2}(KR).$
(S52)

The total differential rates can be found by summing over many harmonics $m$.

##  V Green’s functions in layered media

In the derivation of Cherenkov-like radiation into graphene plasmons, we used
expressions for Green’s functions for polaritons on layered surfaces. We
provide a summary of this method here. In a dielectric with permittivity
$\varepsilon(\mathbf{r},\omega)$, the Maxwell dyadic Green’s function
$G(\mathbf{r},\mathbf{r}^{\prime},\omega)$ is defined so that

$\left(\nabla\times\nabla\times-\frac{\omega^{2}}{c^{2}}\varepsilon(\mathbf{r},\omega)\right)G(\mathbf{r},\mathbf{r}^{\prime},\omega)=\delta^{(3)}(\mathbf{r}-\mathbf{r}^{\prime})I,$
(S53)

where $I$ is the $3\times 3$ identity matrix. In a 2D translation-invariant
medium (such as in a geometry of planar slabs), the Green’s function can be
written as a Fourier transform in the invariant directions (denoted by $\rho$)
as

$G(\mathbf{r},\mathbf{r}^{\prime},\omega)=\int\frac{d^{2}q}{(2\pi)^{2}}e^{i\mathbf{q}\cdot(\rho-\rho^{\prime})}g_{\mathbf{q}}(z,z^{\prime},\omega),$
(S54)

where $g_{\mathbf{q}}(z,z^{\prime},\omega)$ is now a reduced Green’s function
for the $z$-direction only. Substituting this expansion into Eq. S53, we find
that the reduced Green’s function satisfies the equation

$\left((\partial_{z}\hat{z}+i\mathbf{q})\times(\partial_{z}\hat{z}+i\mathbf{q})\times-\frac{\omega^{2}}{c^{2}}\varepsilon(z,\omega)\right)g_{\mathbf{q}}(z,z^{\prime},\omega)=\delta(z-z^{\prime})I$
(S55)

If the final interface is at $z=0$, and we have evanescently decaying waves in
the electrostatic limit, for points $\mathbf{r},\mathbf{r}^{\prime}$ above the
interface, we have the solution

$g_{\mathbf{q}}(z,z^{\prime},\omega)=e^{-q(z+z^{\prime})}\frac{c^{2}q}{\omega^{2}}r_{p}(\omega,q)\hat{\epsilon}_{q}\otimes\hat{\epsilon}_{q}^{*},$
(S56)

where $\hat{\epsilon}_{q}=(\hat{q}+i\hat{z})/\sqrt{2}$, and $r_{p}(\omega,q)$
is the $p$-polarized reflectivity of the surface novotny2012principles.

We can now use this expression, along with properties of $g_{\mathbf{q}}$, to
obtain a useful form for
$\operatorname{Im}G(\mathbf{r},\mathbf{r}^{\prime},\omega^{\prime})$. From the
reciprocity of $G$, it follows that
$g_{-\mathbf{q}}(z,z^{\prime},\omega)=g_{\mathbf{q}}^{T}(z^{\prime},z,\omega)$.
This allows us to write

$\displaystyle\operatorname{Im}G(\mathbf{r},\mathbf{r}^{\prime},\omega)$
$\displaystyle=\operatorname{Im}\int\frac{d^{2}q}{(2\pi)^{2}}e^{i\mathbf{q}\cdot(\rho-\rho^{\prime})}g_{\mathbf{q}}(z,z^{\prime},\omega)$
(S57)
$\displaystyle=\int\frac{d^{2}q}{(2\pi)^{2}}e^{i\mathbf{q}\cdot(\rho-\rho^{\prime})}\left(\frac{g_{\mathbf{q}}(z,z^{\prime},\omega)-g_{\mathbf{q}}^{\dagger}(z^{\prime},z,\omega)}{2i}\right).$
(S58)

Taking the dagger of $g_{\mathbf{q}}$ as defined in Eq. S56 leaves the
polarization factor unchanged, and results in conjugation only on the
reflection coefficient. This allows us to write the form

$\operatorname{Im}G(\mathbf{r},\mathbf{r}^{\prime},\omega)=\int\frac{d^{2}q}{(2\pi)^{2}}e^{i\mathbf{q}\cdot(\rho-\rho^{\prime})}e^{-q(z+z^{\prime})}\frac{c^{2}q}{\omega^{2}}\operatorname{Im}r_{p}(\omega,q)\hat{\epsilon}_{q}\otimes\hat{\epsilon}_{q}^{*},$
(S59)

which is used to derive Eq. S25