# Spike Stream Denoising via Spike Camera Simulation

Liwen Hu1,2, Lei Ma1,2, Zhaofei Yu1,2, Boxin Shi1,2,3, Tiejun Huang1,2,3
1National Engineering Research Center of Visual Technology (NERCVT), Peking
University
2Beijing Academy of Artificial Intelligence
3Institute for Artificial Intelligence, Peking University
Corresponding author.

###### Abstract

As a neuromorphic sensor with high temporal resolution, the spike camera shows
enormous potential in high-speed visual tasks. However, the high-speed
sampling of light propagation processes by existing cameras brings unavoidable
noise phenomena. Eliminating the unique noise in spike stream is always a key
point for spike-based methods. No previous work has addressed the detailed
noise mechanism of the spike camera. To this end, we propose a systematic
noise model for spike camera based on its unique circuit. In addition, we
carefully constructed the noise evaluation equation and experimental scenarios
to measure noise variables. Based on our noise model, the first benchmark for
spike stream denoising is proposed which includes clear (noisy) spike stream.
Further, we design a tailored spike stream denoising framework (DnSS) where
denoised spike stream is obtained by decoding inferred inter-spike intervals.
Experiments show that DnSS has promising performance on the proposed
benchmark. Eventually, DnSS can be generalized well on real spike stream.

Figure 1: Overview of real spike stream denoising which records a high-speed
disk with 2000 revolutions per minute. Based on spike camera simulation
(SCSim) where noise variables are estimated by using the spike-based noise
evaluation equation (SNEE), we systematically explore spike stream denoising
and a tailored framework (DnSS) is proposed.

##  1 Introduction

High-speed scenes have always been a major challenge for some vision
applications such as self-driving cars [21]. Images captured by conventional
cameras with low frame rate suffer from motion blur due to long-time
integration of light, which greatly limits the performance of frame-based
approaches. The emergence of spike camera [35, 7] provides a new perspective
for applications in high-speed scenes. The bio-inspired spike camera not only
has the advantage of high temporal resolution (40000 Hz), but also can report
per-pixel luminance intensity by firing spikes. Recently, it has shown
enormous potential for high-speed visual tasks such as reconstruction [32, 36,
34, 33, 37, 5], optical flow estimation [11], and depth estimation [30].
Different from event cameras [20, 6, 4, 24, 10] and traditional cameras, the
spike camera circuit can introduce unique noise into spike stream during the
sampling process. Eliminating noise in spike stream is a key issue for spike-
based methods [32, 33, 5].
At present, research on spike stream denoising is still in its infancy. One of
the main challenges is the frame-based and event-based methods cannot be
directly applied to the spike stream with unique data modality [35, 11, 33].
Although we can reconstruct clear image sequences first, and then synthesize
the spike stream with a noise-free spike camera model, such reconstruction-
based denoising pipeline depends on the performance of the reconstruction
methods. Traditional reconstruction methods can introduce extra noise and blur
while high-quality reconstruction methods would cost extra processing efforts.
Instead of adopting a reconstruction-based pipeline, we should fully utilize
features of spike stream to design a tailored denoising method. Another
challenge is the lack of datasets to evaluate the performance of denoising
methods. We can use spike camera simulator [11, 37, 33, 17] to generate noisy
spike stream and clean spike stream, but noise models in existing works are
unrealistic. In order to benchmark denoising methods for spike camera
accurately, we need to delve deeper into its distinct noise mechanism. In this
paper, as shown in Fig. 1, we propose a realistic simulation method for spike
cameras, SCSim, where its circuits-level noise is modeled. Further, we set up
the noise measurement experiment and evaluate statistics of the noise variable
by building the relationship between spike stream and luminance intensity.
Based on our noise model, the first synthetic benchmark is introduced to
evaluate denoising methods for spike camera. To remove the noise in spike
stream, we design a tailored spike stream denoising framework, DnSS, where
inter-spike interval (ISI [35]) in clean spike stream as the supervised signal
and the denoising spike stream can be obtained by decoding inferred ISI
sequences. We compare DnSS and the reconstruction-based denoising pipeline.
The experiment shows that lightweight DnSS achieves the best performance. Our
main contributions are summarized as follows:

  * $\bullet$

A spike camera simulation method: We analyze the basic principles of spike
camera circuit implementation in detail and model its distinct noise.

  * $\bullet$

A measurement method for spike camera noise: We propose the Spike-based Noise
Evaluation Equation (SNEE) to establish the relationship between noise and
spike stream, and we collect real spike stream and use SNEE to evaluate the
statistics of noise.

  * $\bullet$

Datasets for spike stream denoising: Based on SCSim and evaluated statistics
of noise variables, we propose two datasets for spike stream denoising.

  * $\bullet$

A spike stream denoising framework: In DnSS, the ISI in clean spike stream can
be inferred based on refine module (RF), and then we obtain the denoised spike
stream by using multi-stage update strategy (MUS) to decode the inferred ISI.

##  2 Related Work

###  2.1 Frame-based and Event-based Denoising

Image Denoising In the early years, hand-designed filters were used for image
denoising e.g., non-local means (NLM) algorithm [2], and the block-matching
based 3D filtering (BM3D) algorithm [19] which use a stack of non-local
similar patches. However , the difficulty of setting the hyper-parameters of
these methods will significantly affect the performance. With the development
of deep learning, dncnn[16] proposes first denoising neural network which
focus on learning a latent mapping from the noisy image to the clean image.
After that, numerous denoising networks are proposed to further improve the
performance [28, 15, 26, 27]. it is challenging and expensive to collect
plenty of training pairs for supervised denoising. To this end, some self-
supervised denoising models [13, 14, 1, 23, 29, 25, 12] are proposed. [13]
introduced Noise2Noise to train a deep denoiser with multiple noisy
observations of the same scenes. [1] and [14] use only one noisy observation
per scene to train deep denoisers. By generating sub-sampled paired images
with random neighbor sub-samplers as training image pairs, [12] can denoise
images more effienctly.

Event Denoising Since event stream have the unique data modality, event
denoising is treated as a classification problem i.e., distinguishing noise
event and real activity event. [22] proposed a filter to tackle the problem of
increased memory requirements for events denoising. [18] used another storage
technique for events and their timestamps to utilize less memory space. [9]
proposed a density matrix to denoising event stream in which each arriving
event is projected into its spatio-temporal region. [3] introduced a reliable
dataset, DVSNOISE20 by deriving an event probability mask using APS frames and
IMU motion data. Based on DVSNOISE20, [3] and [8] proposed a CNN and U-net
network to filter DVS noises, respectively.

###  2.2 Spike-based Applications and Simulation

The spike camera shows potential in high-speed visual tasks. Based on the
statistical characteristics of spike stream, [35] first reconstruct high-speed
scenes. [36, 34] and [33, 37] respectively use SNN and convolutional neural
network to reconstruct high-speed images from a spike stream, which greatly
improved the reconstruction quality. [31] present super-resolution framework
for the spike camera and recover external scenes with both high temporal and
high spatial resolution from spike stream. [11] present a deep learning
pipeline to estimate optical flow in high-speed scenes from spike stream. [30]
can clearly estimate the depth of high-speed scenes by introducing spike-based
transformer. High-speed scenes datasets for spike camera are mainly generated
through simulation. [33] first convert interpolated image sequences with high
frame rate into spike stream. [37, 17] add simple random noise to the ideal
spike camera model. [11] present the spike camera simulator (SPCS) combining
simulation function and rendering engine tightly.

##  3 Spike Camera Noise Simulation

Based on the circuit principle of spike camera, we model and measure the joint
interference of various noises for the generation of spike stream. We start by
introducing a noise-free spike camera model in Section 3.1. And then we
propose a practical model along with the noises in Section 3.2. To measure
statistics of noise, a reasonable experiment is designed. Finally, we evaluate
the spike camera simulation method (SCSim) in Section 3.3.

###  3.1 Ideal Spike Camera Model

The spike camera mimicking the retina fovea consists of an array of $H\times
W$ pixels and can report per-pixel luminance intensity by firing asynchronous
spikes. Specifically, as shown in Fig. 2, each pixel on the spike camera
sensor accumulates incoming light independently and persistently. At time $t$,
for pixel $(x,y)$, if the accumulated brightness arrives a fixed threshold
$\phi$ (as (1)), a spike is fired and then the accumulated brightness can be
reset to 0.

$\displaystyle{A}(x,y,t)=\int_{t_{x,y}^{\rm pre}}^{t}{I_{in}}(x,y,\tau)d\tau\
\geq\phi,$ (1)

where $x,y\in\mathbb{Z},x\leq H,y\leq W$, ${A}(x,y,t)$ is the accumulated
brightness at time $t$, ${I_{in}}(x,y,\tau)$ is the input and $t_{x,y}^{\rm
pre}$ expresses the last time when a spike is fired at pixel $(x,y)$ before
time $t$. If $t$ is the first time to send a spike, then $t_{x,y}^{\rm pre}$
is set as 0. In fact, due to the limitations of circuit technology, the spike
reading times are quantified. Hence, asynchronous spikes are read out
synchronously. Specifically, all pixels periodically judge the spike flag at
time $n\delta t,n\in\mathbb{Z}$, where $\delta t$ is a short interval of
microseconds. Therefore, the output of all pixels forms a $H\times W$ binary
spike frame. As time goes on, the camera would produce a sequence of spike
frames, i.e., a $H\times W\times N$ binary spike stream and can be
mathematically defined as,

$\displaystyle\begin{aligned} &{S}(x,y,n\delta t)=\\\ &\begin{cases}1&\mbox{
if $\exists t\in\left((n-1)\delta t,n\delta t\right]$, s.t.
${A}(x,y,t)\geq\phi$ },\\\ 0&\mbox{ if $\forall t\in\left((n-1)\delta
t,n\delta t\right]$, ${A}(x,y,t)<\phi$ }\\\ \end{cases}\end{aligned}$ (2)

Figure 2: The principle and circuit of spike camera. (a) is the work principle
of spike camera. (b) is the pixel circuit of spike camera. (c) is a specific
operation process. “Star” means that we have considered the corresponding
noise.

###  3.2 Spike Camera Model with Realistic Noise

We will systematically introduce the temporal noise and spatial noise in spike
camera according to its unique circuit.

Temporal Noise Temporal noise is a random variation in the signal that
fluctuates over time. The temporal noise in spike camera mainly includes shot
noise and thermal noise. The shot noise originates from randomness caused by
photon reception and, for pixel $(x,y)$, the probability that $n$ photons are
received between time $t$ and $t+\delta t$ is given by the Poisson probability
distribution, i.e.,

$P({ph}(x,y,t)=n)=\frac{\mu_{{ph}(x,y,t)}}{n!e^{\mu_{{ph}(x,y,t)}}}$ (3)

where $n\in\mathbb{N}$, ${ph}(x,y,t)$ is a random variable representing the
number of received photons from $t$ to $t+\delta t$, and $\mu_{{ph}(x,y,t)}$
is the expectation of ${ph}(x,y,t)$. The random number of photons can affect
the luminance intensity. Since $\delta t$ is enough small, we consider the
luminance intensity at time $t$ is proportional to the number of photons
between time $t$ and $t+\delta t$, i.e., ${L}(x,y,t)\propto{ph}(x,y,t)$.
Different from the previous simulator [11] which uses the ideal luminance
intensity $\mu_{{L}(x,y,t)}$ as the input, we sample ${L}(x,y,t)$ at time $t$
during simulation, i.e.,

${L}(x,y,t)=\mu_{{L}(x,y,t)}\dfrac{{ph}(x,y,t)}{\mu_{{ph}(x,y,t)}}.$ (4)

The threshold $\phi$ (as (1)) in the ideal spike camera model also fluctuates
over time. We start by writing the ideal threshold in the form of the circuit,
i.e,

$\phi=CV_{d}=C(V_{D}-V_{ref})$ (5)

The voltage fluctuation in the camera can be caused by the reset transistor,
which is affected by temperature where $V_{D}$ is reset voltage and $V_{ref}$
is reference voltage. As shown in Fig. 2(b), the reset transistor affected by
temperature can cause voltage fluctuation. We use the random variable
${V}^{T_{0}}(x,y,t)$ to describe fluctuating voltage in pixel $(x,y)$ at time
$t$ and ${V}^{T_{0}}(x,y,t)$ can be considered to obey Gaussian distribution,
i.e.,

${V}^{T_{0}}(x,y,t)\sim
N(0,(\sigma^{T_{0}})^{2}),\quad\sigma^{T_{0}}=\sqrt{\frac{kT_{0}}{C}},$ (6)

where $\sigma^{T_{0}}$ is the standard deviation of ${V}^{T_{0}}(x,y,t)$, $k$
is Boltzmann constant, $T_{0}$ is absolute temperature.

Spatial Noise Spatial noise i.e., fixed-pattern noise, is a random variation
in the signal that has nothing to do with time. The spatial noise in spike
camera mainly includes conversion rate, dark current, capacitor mismatch and
bias voltage. As shown in Fig. 2(a), we first develop the input $I_{in}$ based
on circuit principles, i.e., ${I_{in}}(x,y,\tau)=\alpha{L}(x,y,\tau)$ where
$\alpha$ is the photoelectric conversion rate, ${L}(x,y,\tau)$ refers to the
brightness of pixel $(x,y)$ at time $\tau$. The difference of photodiodes can
cause the conversion rate to fluctuate. We assume the conversion rate in pixel
$(x,y)$ obeys Gaussian distribution as,

${\alpha}(x,y)\sim N(\mu_{\alpha},(\sigma_{\alpha}^{S})^{2}),$ (7)

where $\mu_{\alpha}$ and $\sigma_{\alpha}$ contribute to the expectation and
variance of $\alpha(x,y)$, respectively. The dark current is also from
photodiodes. The dark current $I_{\rm{dark}}(x,y)$ in pixel $(x,y)$ obeys
Gaussian distribution as,

$I_{\rm{dark}}(x,y)\sim N(\mu_{\rm{dark}},(\sigma_{\rm{dark}}^{S})^{2}),$ (8)

where $\mu_{\rm{dark}}$ and $\sigma_{\rm{dark}}^{S}$ contribute to the
expectation and variance of $I_{\rm{dark}}(x,y)$. Besides, capacitance
difference in each pixel obeys Gaussian distribution and can be written as,

$C^{S}(x,y)\sim N(0,(\sigma_{C}^{S})^{2}),$ (9)

where $C^{S}(x,y)$ is the random variable describing the capacitance
nonuniformity at pixel $(x,y)$, $\sigma_{C}^{S}$ is the standard deviation of
$C^{S}(x,y)$. The bias voltage $V^{S}(x,y)$ mainly comes from reset voltage
nonuniformity in check module as shown in Fig. 2(c). We can assume the bias
voltage $V^{S}(x,y)$ at pixel $(x,y)$ obeys Gaussian distribution, i.e.,

$V^{S}(x,y)\sim N(0,(\sigma_{V}^{S})^{2}),$ (10)

where $\sigma_{V}^{S}$ is the standard deviation of $V^{S}(x,y)$.

The Sampling Process with Noise We can add the above noise model into the
generation of spike stream. At time $t$, for pixel $(x,y)$, the accumulated
brightness ${A}(x,y,t)$ in (1) can be rewritten as,

$\int_{t_{x,y}^{\rm
pre}}^{t}\\!\\!\\!\\!\\!\\!\alpha(x,y)L(x,y,\tau)+I_{\rm{dark}}(x,y)d\tau\geq\phi(x,y,t),$
(11)

where $\phi(x,y,t)$ is the threshold affected by noise. Further, $\phi(x,y,t)$
can be expressed as,

$\displaystyle(C+C^{S}(x,y))(V_{d}+V^{T_{0}}(x,y,t)+V^{S}(x,y)).$ (12)

The spike stream ${S}(x,y,n\delta t)$ in (2) can be rewritten as,

$\displaystyle\begin{cases}1\\!\\!\\!\\!&\mbox{ if $\exists
t\in\left((n-1)\delta t,n\delta t\right]$, s.t.
${A}(x,y,t)\geq\phi(x,y,t)$},\\\ 0\\!\\!\\!\\!&\mbox{ if $\forall
t\in\left((n-1)\delta t,n\delta t\right]$, ${A}(x,y,t)<\phi(x,y,t)$}.\\\
\end{cases}$ (13)

###  3.3 Noise Measurement

For simulating noise in spike camera accurately, we need to measure statistics
of noise, i.e., expectation and variance. Since statistics of temporal noise
is known, we only need to measure spatial noise, i.e., $\alpha(x,y)$,
$I_{\rm{dark}}(x,y)$, $V^{S}(x,y)$ and $C^{S}(x,y)$. We first establish the
equation spike-based noise evaluation equation (SNEE) according to the meaning
of spike. Further, we built the corresponding experimental scenes to capture
real spike stream. Finally, noise variables are estimated by solving the
equations of real data.

Spike-based Noise Evaluation Equation A spike fired by pixel $(x,y)$ in time
$t$ means that the accumulation at the pixel arrives threshold $\phi(x,y,t)$.
Hence, if accumulation would not be reset, the total accumulation at the pixel
$(x,y)$ in time $n\delta t$ can be estimated as,

$\sum_{S(x,y,i\delta t)=1}^{i}\\!\\!\\!\\!\\!\\!\\!\\!\phi(x,y,t_{i}),$ (14)

where $i\in\mathbb{Z}\cap[1,n]\;s.t.\;S(x,y,i\delta t)=1$, $t_{i}$ express the
time when the i-th spike is fired and $t_{i}\in((i-1)\delta t,i\delta t]$.
According to (11), the total accumulation at pixel $(x,y)$ in time $n\delta t$
also can be expressed as,

$\int_{0}^{n\delta t}\\!\\!\\!\alpha(x,y)L(x,y,\tau)+I_{\rm{dark}}(x,y)d\tau.$
(15)

Further, we have the spike-based noise evaluation equation (SNEE), i.e.,

$\displaystyle\sum_{S(x,y,i\delta
t)=1}^{i}\\!\\!\\!\\!\\!\\!\\!\\!\phi(x,y,t_{i})$ (16)
$\displaystyle=\\!\\!\\!\\!\\!\\!\\!\sum_{S(x,y,i\delta
t)=1}^{i}\\!\\!\\!\\!\\!\\!\\!\\!(C+C^{S}(x,y))(V_{d}+V^{T_{0}}(x,y,t_{i})+V^{S}(x,y))$
$\displaystyle=\int_{0}^{n\delta
t}\\!\\!\\!\alpha(x,y)L(x,y,\tau)+I_{\rm{dark}}(x,y)d\tau$

We can set experimental scenes to control $L(x,y,\tau)$ and eliminate temporal
noise by extending sampling time $n\delta t$. When experimental scene is
static and sampling time $n\delta t$ is long enough, (16) can be simplified
as,

$\displaystyle\sum_{S_{k}(x,y,i\delta
t)=1}^{i}\\!\\!\\!\\!\\!\\!\\!\\!(C+C^{S}(x,y))(V_{d}+V^{S}(x,y))$ (17)
$\displaystyle=\int_{0}^{n\delta
t}\\!\\!\\!\alpha(x,y)\mu_{k}+I_{\rm{dark}}(x,y)d\tau$

where $k\in\mathbb{Z}$, $\mu_{k}$ is the ideal luminance intensity in the k-th
static scene and $\mathbf{S}_{k}(x,y,t)$ is the spike stream captured from
k-th static scene.

Figure 3: Spike camera noise measurement setup. Figure 4: Comparison of spike
stream from [11], real spike camera and our model.(a) is the average number of
spike per sampling. (b) is histograms of Inter-spike interval (ISI) [35] in
spike stream. (c) is spike pattern where spike stream of $250\times 400\times
10$ is pressed onto a surface.

Spike camera noise measurement setup As shown in Fig. 3, we built the
experimental scences where the spike camera shot computer monitor under 25
different grayscale backgrounds respectively. Further, accroding to (17), we
can get a system of equations about $\alpha(x,y)$, $I_{\rm{dark}}(x,y)$,
$C^{S}(x,y)$ and $V^{S}(x,y)$ for each pixel $(x,y)$. We can solve the above
equations by optimization methods. Finally, the statics of spatial noise can
be evaluated from $\alpha(x,y)$, $I_{\rm{dark}}(x,y)$, $C^{S}(x,y)$ and
$V^{S}(x,y)$.

Evaluation of simulated spike stream We compare simulated spike stream and
captured real spike stream. Sampling scenes are solid backgrounds at 5 gray
levels which are beneficial to ensure consistency between real and synthetic
scenes. Sampling time is 0.1s and three kinds of spike stream are obtained
under each scene. Finally, we have $5\times 3$ spike stream. Fig. 4
illustrates reliability of SCSim.

##  4 Method

Different from image and event denoising, spike stream denoising is to infer
clean spike stream in a special spatio-temporal sampling from noisy spike
stream. Based on its unique data modality, we propose a tailored denoising
framework. In DnSS, denoised spike stream can be obtained by decoding inferred
ISI sequence from noisy spike stream. Next, we use
$\mathbf{S}^{t_{1},t_{2}}\in{\\{0,1\\}}^{H\cdot W\cdot(t_{2}-t_{1})}$ to
denote a spike stream from time $t_{1}$ to $t_{2}$.

###  4.1 Inter-spike Interval Plane

Spike stream is not suitable to be directly used as supervisory signal because
the encoded information is sparse (0 or 1). Instead of spike stream, we first
use its inter-spike interval plane, ISI [35], as supervision signal where the
information of spike stream in a temporal window across time $t$ is compressed
to a plane at time $t$ (written as $\mathbf{ISI}_{t}$).

###  4.2 Network Architecture

Figure 5: The architecture of proposed denoising method for spike camera
(DnSS). (a) Extract spatio-temporal denoising features of spike stream. (b)
Extract texture features in spike stream. (c) Fuse spatio-temporal denoising
features and texture features. (d) Infer ISI of clean spike stream based above
features. (e) Decode inferred ISI sequence to denoised spike stream. Figure 6:
The principle of multi-stage update strategy.

As shown in Fig. 5, in DnSS, We use the noisy spike stream in a temporal
window across time $t$, $\mathbf{S}_{n}^{t-\Delta t,t+\Delta t}$, to infer the
normalized spike interval of clean spike stream at time $t$,
$\tilde{\mathbf{ISI}}_{t}^{c}$ where elements are mapped to $0-1$ by taking
the inverse for each element in ${\mathbf{ISI}_{t}^{c}}$. Specifically, as
Fig. 5(a), the input spike stream $\mathbf{S}_{n}^{t-\Delta t,t+\Delta t}$
first is processed by 2-layer convolution where temporal noise is removed and
spike information is compressed. Then, we can get the denoising feature map
$\mathbf{F}_{dn}$ by using a U-Net and spatial attention module to remove
spatial noise.

In fact, all denoising methods are likely to introduce artifacts and loss of
details. Therefore, we design a refine module to restore high-frequency
details in spike stream where $\tilde{\mathbf{ISI}}_{t}^{n}$ is as input. As
Fig. 5(b), an spatial attention module is used to guide
$\tilde{\mathbf{ISI}}_{t}^{n}$ to pay more attention to regions with rich
textures. Then, we extract the texture feature map $\mathbf{F}_{tx}$ through
4-layer convolution. By introducing a fusion module, feature $\mathbf{F}_{tx}$
and feature $\mathbf{F}_{dn}$ can be fused as the fusion feature map
$\mathbf{F}_{fu}$ as Fig. 5(c). Finally, as Fig. 5(d), we use previous feature
maps to infer $\tilde{\mathbf{ISI}}_{t}^{c}$, i.e.,

$\displaystyle\tilde{\mathbf{ISI}}_{t}^{dn*}=f(\mathbf{F}_{dn},\mathbf{F}_{tx},\mathbf{F}_{fu}),$
(18)
$\displaystyle[{\Theta}_{t},\mathbf{1}-{\Theta}_{t}]=a([\tilde{\mathbf{ISI}}_{t}^{dn*},\tilde{\mathbf{ISI}}_{t}^{n}]),$
(19)
$\displaystyle\tilde{\mathbf{ISI}}_{t}^{dn}={\Theta}_{t}\otimes\tilde{\mathbf{ISI}}_{t}^{dn*}+(\mathbf{1}-\Theta_{t})\otimes\tilde{\mathbf{ISI}}_{t}^{n},$
(20)

where $f(\cdot)$ contributes the initial inference operation,
$\tilde{\mathbf{ISI}}_{t}^{dn*}$ is the normalized spike interval of denoising
spike stream at time $t$ which is first inferred, $a(\cdot)$ contributes is an
attention operation, $\otimes$ is product of elements and
$\tilde{\mathbf{ISI}}_{t}^{dn}$ is the normalized spike interval of denoising
spike stream at time $t$. The $\tilde{\mathbf{ISI}}_{t}^{dn}$ is supervised by
the normalized spike interval of clean spike stream at time $t$:

$\mathcal{L}=\|\tilde{\mathbf{ISI}}_{t}^{dn}-\tilde{\mathbf{ISI}}_{t}^{c}\|_{2}^{2},$
(21)

where $\mathcal{L}$ is our loss function and $\|\cdot\|_{2}$ denotes the
2-norm. After taking the inverse for each element in
$\tilde{\mathbf{ISI}}_{t}^{dn}$, we can obtain denoised spike interval
${\mathbf{ISI}}_{t}^{dn}$. Due to the limitation of temporal window, the
information in the input spike stream under dark light conditions may be
insufficient and the denoised interval ${\mathbf{ISI}}_{t}^{dn}$ would be
inaccurate. To this end, as shown in Fig. 6, we design a multi-stage update
strategy where, for pixel $(x,y)$,
$\\{{\mathbf{ISI}}_{t}^{dn}\\},t=t_{0},t_{0}+1,\dots,t_{0}+{\mathbf{ISI}}_{t_{0}}^{dn}(x,y)$
is used to update the final spike interval at time $t_{0}$,
$\overline{\mathbf{ISI}}_{t_{0}}^{dn}$, i.e.,

$\displaystyle\overline{\mathbf{ISI}}_{t_{0}}^{dn}(x,y)=t_{k}+{\mathbf{ISI}}_{t_{k}}^{dn}(x,y)-t_{0},$
(22)

where $k=\rm{argmax}_{m\in\mathbb{N^{*}}}\\{\forall r\in\mathbb{N},r\leq
m,s.t.,t_{r-1}<t_{r}<t_{r}+{\mathbf{ISI}}_{t_{r}}^{dn}(x,y)<t_{r-1}+{\mathbf{ISI}}_{t_{r-1}}^{dn}(x,y)\\}$.

Figure 7: Examples of noisy (clean) spike stream.

Further, the final spike interval sequence can be decoded as denoised spike
stream $\mathbf{S}_{dn}^{0,T}$.

##  5 Experiments

###  5.1 Implementation Details

Datasets We propose noisy versions of SPIFT and PHM [11], $\text{SPIFT}_{N}$
and $\text{PHM}_{N}$. As shown in Fig. 7, our datasets include simulated spike
stream pairs by SCSim, i.e., clean (noisy) spike stream. We use
$\text{SPIFT}_{N}$ including 100 random high-speed scenes as our training set
while $\text{PHM}_{N}$ including 10 well-designed high-speed scenes as the
test set. Besides, to evaluate the generalization of DnSS in real spike
stream, we use the real spike stream [33, 36] captured from 5 high-speed
scenes i.e., “Car”, “Doll”, “Fan”, “Rotation1”, “Train”.

Figure 8: Imaging results of denoised (noisy) spike stream on proposed
dataset, $\text{PHM}_{N}$. Figure 9: Imaging results of denoised (noisy) spike
stream on dataset from [33, 36] .

Train setup We implemented our DnSS in PyTorch. In training, Adam optimizer is
adopted to optimize the networks and the initial learning rate is set to
0.0001. DnSS is trained for 50 epochs with a batch size of 32 on 1 NVIDIA
A100-SXM4-80GB GPU.

Metric Since there are no methods to measure the similarity of spike stream,
we cannot directly evaluate the performance of denoising methods. To this end,
we compare imaging quality of different spike stream, i.e., PSNR and SSIM,
where the state-of-the-art reconstruction method [33] is used.

###  5.2 Results

We compare our method with the spike stream denoising pipeline where denoised
image sequences are first obtained by reconstructing noisy spike stream, and
then the denoised spike stream is synthesized from denoised image sequences
with noise-free spike camera model. The performance of spike stream denoising
pipeline depends on reconstruction methods. Here, we use the traditional
method, TFP[35], the supervised method, Spk2ImgNet[33] and the self-supervised
method, SSML[5] respectively.

Table 1: Comparison of imaging results on $\text{PHM}_{N}$ dataset. | Method
---
| Param.
---
| PSNR
---
| SSIM
---
Noisy data | / | 27.70 | 0.698
TFP | / | 19.90 | 0.365
SSML | 2.38M | 27.08 | 0.723
Spk2ImgNet | 3.91M | 25.91 | 0.664
Ours | 0.27M | 29.61 | 0.847

Qualitative Evaluation on Synthetic Datasets. we compare the imaging
performance of noisy spike stream and denoised spike stream. As shown in
table. 1, the denoised spike stream from TFP and Spk2ImgNet cannot improve
image reconstruction quality and SSML can remove noise in spike stream to a
certain extent. Our denoising method gets the best average performance.
Moreover, the parameters of DnSS are far less than SSML and Spk2ImgNet. Fig. 8
shows the imaging results of denoised spike stream by different methods. TFP
can introduce a lot of motion blur during denoising. Although texture details
can be retained well by Spk2ImgNet, noise cannot be effectively eliminated.
Comparing the first two, SSML is able to retain partial textural details
during denoising. However, the denoising performance of SSML is limited. Our
method can effectively remove the noise in noisy spike stream without the loss
of texture details.

Table 2: Comparison of imaging results on dataset from [36, 33] Metric | Noise | TFP | SSML | Spk2ImgNet | Ours
---|---|---|---|---|---
NIQE($\downarrow$) | 4.07 | 4.29 | 4.10 | 4.32 | 3.76

Qualitative Evaluation on Real Datasets. We compare the imaging quality of
denoised spike stream by different methods are as input. Due to absence of GT
in spike camera datasets recording high-speed scenes, we use a no-reference
IQA metric, i.e., NIQE to evaluate different denoising methods in Table. 2.
Besides, the imaging results of different denoising methods are shown in Fig.
9. We can find our method can effectively remove noise introduced by hot
pixels which proves that our denoising method has pretty generalization on
real data.

###  5.3 Ablation

To investigate the effect of the proposed refine module (RF) and multi-stage
update strategy (MUS), we compare three baseline methods. Baseline1 is the
basic baseline without RF or MUS. Based on Baseline1, Baseline2 uses multi-
stage update strategy to reduce the inference error and Baseline3 adds refine
module to improve texture details.

Table 3: Evaluation of the proposed modules on $\text{PHM}_{N}$ dataset. Method | Param. | MUS | RF | PSNR | SSIM
---|---|---|---|---|---
Baseline1 | 0.20M | ✕ | ✕ | 28.67 | 0.822
Baseline2 | 0.20M | ✓ | ✕ | 29.13 | 0.828
Baseline3 | 0.27M | ✕ | ✓ | 29.11 | 0.841
Ours | 0.27M | ✓ | ✓ | 29.61 | 0.847
Figure 10: Ablation study of the proposed modules on $\text{PHM}_{N}$. Table 4: Evaluation of the proposed modules on dataset from [36, 33]. Metric | Baseline1 | Baseline2 | Baseline3 | Ours
---|---|---|---|---
NIQE($\downarrow$) | 4.15 | 4.14 | 3.80 | 3.76
Figure 11: Ablation study of the proposed modules on dataset from [36, 33].

Evaluation on Synthetic Datasets. As shown in Table. 3, we find that RF and
MUS can improve the performance of the basic baseline. By combining RF and
MUS, Our DnSS outperforms the basic baseline by around 1 dB. Fig. 10 shows the
imaging results of denoised spike stream.

Evaluation on Real Datasets. As shown in Table. 4, Our DnSS achieves the best
denoising performance on the real spike stream by introducing RF and MUS. The
imaging results in Fig. 11 also show advantages of RF and MUS intuitively.
Specifically, by using the MUS module into Baseline1, the wrong interval in
Baseline2 can be updated which makes the bad points in dark eliminated. By
using the RF module into Baseline1, the edge details become smoother.
Combining RF and MUS, Our DnSS can keep the details in the real scene better.

##  6 Conclusion

We propose a tailored and lightweight denoising framework for spike camera,
DnSS, where the ISI in clean spike stream can be inferred based on RF, and
then we obtain the denoised spike stream by using MUS to decode the inferred
ISI. To evaluate the performance of denoising methods, a spike camera
simulation method, SCSim, is proposed where we model and measure the noise in
spike camera refer to its unique circuit. Based on SCSim, we synthesize the
spike-based denoising datasets, $\text{SPIFT}_{N}$ and $\text{PHM}_{N}$. The
experiment demonstrates DnSS can denoise spike stream efficiently. Moreover,
DnSS can be generalized well on real spike stream.

## References

  * [1] Krull Alexander, Buchholz Tim-Oliver, and Jug Florian.  Noise2void - learning denoising from single noisy images.  In IEEE Conference on Computer Vision and Pattern Recognition (CVPR), page 2124–2132, 2019.
  * [2] Buades Antoni, Coll Bartomeu, and Morel J-M.  A non-local algorithm for image denoising.  In IEEE Conference on Computer Vision and Pattern Recognition (CVPR), page 60–65, 2005.
  * [3] R. Wes Baldwin, Mohammed Almatrafi, Vijayan Asari, and Keigo Hirakawa.  Event probability mask (epm) and event denoising convolutional neural network (edncnn) for neuromorphic cameras.  In IEEE Conference on Computer Vision and Pattern Recognition (CVPR), pages 1698–1707, 2020.
  * [4] Christian Brandli, Raphael Berner, Minhao Yang, Shih-Chii Liu, and Tobi Delbruck.  A 240 × 180 130 db 3 $\mu\rm s$ latency global shutter spatiotemporal vision sensor.  IEEE Journal of Solid-State Circuits (JSSC), 49(10):2333–2341, 2014\.
  * [5] Shiyan chen, Chaoteng Duan, Zhaofei Yu, Ruiqin Xiong, and Tiejun Huang.  Self-supervised mutual learning for dynamic scene reconstruction of spiking camera.  In Proceedings of the International Joint Conference on Artificial Intelligence, IJCAI, pages 2859–2866, 2022.
  * [6] Tobi Delbrück, Bernabe Linares-Barranco, Eugenio Culurciello, and Christoph Posch.  Activity-driven, event-based vision sensors.  IEEE International Symposium on Circuits and Systems (ISCAS), pages 2426–2429, 2010.
  * [7] Siwei Dong, Lin Zhu, Yonghong Tian, and Tiejun Huang.  An efficient coding method for spike camera using inter-spike intervals.  In Data Compression Conference (DCC), pages 437–437, 2017.
  * [8] Peiqi Duan, Zihao W. Wang, Xinyu Zhou, Yi Ma, and Boxin Shi.  Eventzoom: Learning to denoise and super resolve neuromorphic events.  In IEEE Conference on Computer Vision and Pattern Recognition (CVPR), pages 12824–12833, 2021.
  * [9] Yang Feng, Hengyi Lv, Hailong Liu, Yisa Zhang, Yuyao Xiao, and Chengshan Han.  Event density based denoising method for dynamic vision sensor.  Applied Sciences, 10(6):2024, 2020.
  * [10] Menghan Guo, Jing Huang, and Shoushun Chen.  Live demonstration: A 768× 640 pixels 200meps dynamic vision sensor.  IEEE International Symposium on Circuits and Systems (ISCAS), pages 1–1, 2017.
  * [11] Liwen Hu, Rui Zhao, Ziluo Ding, Lei Ma, Boxin Shi, Ruiqin Xiong, and Tiejun Huang.  Optical flow estimation for spiking camera.  In IEEE Conference on Computer Vision and Pattern Recognition (CVPR), pages 17844–17853, 2022.
  * [12] Tao Huang, Songjiang Li, Xu Jia, Huchuan Lu, and Jianzhuang Liu.  Neighbor2neighbor: Self-supervised denoising from single noisy images.  In IEEE Conference on Computer Vision and Pattern Recognition (CVPR), pages 14781–14790, 2021.
  * [13] Lehtinen Jaakko, Munkberg Jacob, Hasselgren Jon, Laine Samuli, Karras Tero, Aittala Miika, and Aila Timo.  Noisier2noise: Learning to denoise from unpaired noisy data.  In International Conference on Machine Learning (ICML), page 2965–2974, 2018.
  * [14] Batson Joshua and Royer Loic.  Noise2self: Blind denoising by self-supervision.  In International Conference on Machine Learning (ICML), page 524–533, 2019.
  * [15] Zhang Kai, Zuo Wangmeng, and Zhang Lei.  Ffdnet: Toward a fast and flexible solution for cnn-based image denoising.  27(9):4608–4622, 2017.
  * [16] Zhang Kai, Zuo Wangmeng, Chen Yunjin, Meng Deyu, and Zhang Lei.  Beyond a gaussian denoiser: Residual learning of deep cnn for image denoising.  page 3142–3155, 2017.
  * [17] Zhaodong Kang, Jianing Li, Lin Zhu, and Yonghong Tian.  Retinomorphic sensing: A novel paradigm for future multimedia computing.  In Proceedings of the ACM International Conference on Multimedia (ACMMM), page 144–152, 2021.
  * [18] Alireza Khodamoradi and Ryan Kastner.  $o(n)$-space spatiotemporal filter for reducing noise in neuromorphic vision sensors.  IEEE Transactions on Emerging Topics in Computing, 9(1):15–23, 2021\.
  * [19] Dabov Kostadin, Foi Alessandro, Katkovnik Vladimir, and Egiazarian Karen.  Image denoising by sparse 3-d transform-domain collaborative filtering.  page 2080–2095, 2007.
  * [20] Patrick Lichtsteiner, Christoph Posch, and Tobi Delbruck.  A 128 × 128 120 db 15 $\mu\rm s$ latency asynchronous temporal contrast vision sensor.  IEEE Journal of Solid-State Circuits (JSSC), 43(2):566–576, 2008\.
  * [21] Martin Litzenberger, Christoph Posch, D Bauer, Ahmed Nabil Belbachir, P Schon, B Kohn, and H Garn.  Embedded vision system for real-time object tracking using an asynchronous transient vision sensor.  In IEEE Digital Signal Processing Workshop, pages 173–178. IEEE, 2006.
  * [22] Hongjie Liu, Christian Brandli, Chenghan Li, Shih-Chii Liu, and Tobi Delbruck.  Design of a spatiotemporal correlation filter for event-based sensors.  In IEEE International Symposium on Circuits and Systems (ISCAS), pages 722–725, 2015.
  * [23] Moran Nick, Schmidt Dan, Zhong Yu, and Coady Patrick.  Noisier2noise: Learning to denoise from unpaired noisy data.  In IEEE Conference on Computer Vision and Pattern Recognition (CVPR), page 12061–12069, 2020.
  * [24] Christoph Posch, Daniel Matolin, and Rainer Wohlgenannt.  An asynchronous time-based image sensor.  IEEE International Symposium on Circuits and Systems (ISCAS), 49(10):2130–2133, 2008.
  * [25] Laine Samuli, Karras Tero, Lehtinen Jaakko, and Aila Timo.  High-quality self-supervised deep image denoising.  In Advances in Neural Information Processing Systems (NIPS), page 6970–6980, 2019.
  * [26] Guo Shi, Yan Zifei, Zhang Kai, Zuo Wangmeng, and Zhang Lei.  Toward convolutional blind denoising of real photographs.  In IEEE Conference on Computer Vision and Pattern Recognition (CVPR), page 1712–1722, 2019.
  * [27] Gu Shuhang, Li Yawei, Gool Luc, Van, and Timofte Radu.  Self-guided network for fast image denoising.  In IEEE International Conference on Computer Vision (ICCV), page 2511–2520, 2019.
  * [28] Lefkimmiatis Stamatios.  Universal denoising networks : A novel cnn architecture for image denoising.  In IEEE Conference on Computer Vision and Pattern Recognition (CVPR), page 3204–3213, 2018.
  * [29] Wu Xiaohe, Liu Ming, Cao Yue, Ren Dongwei, and Zuo Wangmeng.  Unpaired learning of deep image denoising.  In European Conference on Computer Vision (ECCV), 2020.
  * [30] Jiyuan Zhang, Lulu Tang, Zhaofei Yu, Jiwen Lu, and Tiejun Huang.  Spike transformer: Monocular depth estimation for spiking camera.  In European Conference on Computer Vision (ECCV), 2022.
  * [31] Jing Zhao, Jiyu Xie, Ruiqin Xiong, Jian Zhang, Zhaofei Yu, and Tiejun Huang.  Super resolve dynamic scene from continuous spike streams.  In IEEE International Conference on Computer Vision (ICCV), pages 2533–2542, 2021.
  * [32] Jing Zhao, Ruiqin Xiong, and Tiejun Huang.  High-speed motion scene reconstruction for spike camera via motion aligned filtering.  In International Symposium on Circuits and Systems (ISCAS), pages 1–5, 2020.
  * [33] Jing Zhao, Ruiqin Xiong, Hangfan Liu, Jian Zhang, and Tiejun Huang.  Spk2imgnet: Learning to reconstruct dynamic scene from continuous spike stream.  In IEEE Conference on Computer Vision and Pattern Recognition (CVPR), pages 11996–12005, 2021.
  * [34] Yajing Zheng, Lingxiao Zheng, Zhaofei Yu, Boxin Shi, Yonghong Tian, and Tiejun Huang.  High-speed image reconstruction through short-term plasticity for spiking cameras.  In IEEE Conference on Computer Vision and Pattern Recognition (CVPR), pages 6358–6367, 2021.
  * [35] Lin Zhu, Siwei Dong, Tiejun Huang, and Yonghong Tian.  A retina-inspired sampling method for visual texture reconstruction.  In IEEE International Conference on Multimedia and Expo (ICME), pages 1432–1437, 2019.
  * [36] Lin Zhu, Siwei Dong, Jianing Li, Tiejun Huang, and Yonghong Tian.  Retina-like visual image reconstruction via spiking neural model.  In IEEE Conference on Computer Vision and Pattern Recognition (CVPR), pages 1435–1443, 2020.
  * [37] Lin Zhu, Jianing Li, Xiao Wang, Tiejun Huang, and Yonghong Tian.  Neuspike-net: High speed video reconstruction via bio-inspired neuromorphic cameras.  In IEEE International Conference on Computer Vision (ICCV), pages 2400–2409, 2021.


# Robustness of entanglement in Hawking radiation for optical systems immersed
in thermal baths

Ivan Agullo<EMAIL_ADDRESS>Department of Physics and Astronomy, Louisiana
State University, Baton Rouge, LA 70803, U.S.A.  Anthony J. Brady
<EMAIL_ADDRESS>Department of Electrical and Computer Engineering,
University of Arizona, Tucson, Arizona 85721, USA  Dimitrios Kranas
<EMAIL_ADDRESS>Department of Physics and Astronomy, Louisiana State
University, Baton Rouge, LA 70803, U.S.A.

###### Abstract

Entanglement is the quantum signature of Hawking’s particle pair-creation from
causal horizons, for gravitational and analog systems alike. Ambient thermal
fluctuations, ubiquitous in realistic situations, strongly affects the
entanglement generated in the Hawking process, completely extinguishing it
when the ambient temperature is of the same order as the Hawking temperature.
In this work, we show that optical analog systems have a built-in robustness
to thermal fluctuations which are at rest in the laboratory. In such systems,
horizons move relative to the laboratory frame at velocities close to the
speed of light. We find that a subtle interplay between this relative velocity
and dispersion protects the Hawking-generated entanglement—allowing ambient
temperatures several orders of magnitude larger than the Hawking temperature
without significantly affecting entanglement.

##  I Introduction

That causal horizons spontaneously emit pairs of entangled particles is a
deeply fascinating prediction [1, 2] that subtly combines various aspects of
physics: the kinematics of causal barriers, thermodynamics, and quantum
mechanics. As a whole, the Hawking process (as we generically call this
phenomenon, regardless of the physical system under consideration) is agreed
to be a genuinely quantum process. However, which particular aspects of this
process are quantum, and which are not? In a classical universe, there are two
aspects of the Hawking effect which would not exist: (i) particle-creation
starting from an initial vacuum (i.e., spontaneous emission) and (ii) quantum
mechanically entangled Hawking pairs.

The first aspect is difficult to recreate, since it is challenging to
completely isolate event horizons; thermal fluctuations and other sources of
stochastic noise are ubiquitous. This is certainly true for astrophysical
black holes, which are immersed in the cosmic microwave background radiation
(amongst other source, such as the stochastic background of gravitational
waves, the cosmic background of neutrinos, etc.). When embedded in a populated
environment, a horizon emits induced or stimulated Hawking radiation, and it
can be difficult to distinguish between spontaneous emission (originating from
vacuum fluctuations) and induced emission.

What about the second aspect—entanglement? It turns our that background noise
also deteriorates the Hawking-generated entanglement, even causing it to
completely disappear for sufficiently intense noise. Mathematically, the core
of Hawking’s pair production is a process of two-mode squeezing (see, e.g.,
Refs. [3, 4]). In other words, the time evolution induces the following
transformation, which relates creation and annihilation operators of the
output modes to those of the input

$\displaystyle\hat{a}^{\rm(in)}_{1}$ $\displaystyle\to$
$\displaystyle\hat{a}^{\rm(out)}_{1}=\cosh
r\,\hat{a}^{\rm(in)}_{1}+e^{i\phi}\,\sinh r\,\hat{a}^{\rm(in)\,\dagger}_{2},$
(I.1) $\displaystyle\hat{a}^{\rm(in)}_{2}$ $\displaystyle\to$
$\displaystyle\hat{a}^{\rm(out)}_{2}=\cosh
r\,\hat{a}^{\rm(in)}_{2}+e^{i\phi}\,\sinh r\,\hat{a}^{\rm(in)\,\dagger}_{1}.$
(I.2)

where $r$ and $\phi$ are the squeezing intensity and squeezing angle,
respectively. In this expression, $\hat{a}^{\rm(in)}_{i}$ ($i=1,2$) are the
annihilation operators for the progenitors of Hawking pairs, while
$a^{\rm(out)}_{1}$ labels the outgoing Hawking quanta and $a^{\rm(out)}_{2}$
their entangled partners.

Under the above transformation, the vacuum evolves to a two-mode squeezed
vacuum, in which case each subsystem is in a mixed thermal state, but there
exist quantum correlations between the subsystems that purify the global
state. If one replaces the initial vacuum by a thermal state of equal quanta
$n_{\rm env}$ in each of the two input modes, it is easy to check that the
entanglement in the final state degrades and disappears for $n_{\rm env}\geq
e^{r}\sinh r$ (see, e.g., Appendix A of Ref. [5] for details). Hence, when
thermal fluctuations dominate, the squeezer is incapable of entangling the
modes, and the final state is separable. This is a generic fact about two-mode
squeezing processes—including the Hawking process of gravitational and analog
systems.

The goal of this article it to provide a quantitative analysis for the impact
that thermal noise has on Hawking-generated entanglement. For concreteness, we
focus on an optical analog system containing a white-black hole pair that is
moving at a finite velocity (a fraction of the speed of light) with respect to
the lab frame. One of the main messages of our work is that this relative
velocity introduces a large blue-shift, which in turn significantly increases
the _threshold temperature_ (as measured in the lab) at which entanglement
vanishes. The robustness of entanglement to high lab temperatures originates
from (i) dispersive effects of the medium and (ii) the Lorentz boost between
the lab frame—in which the equilibrium temperature of the thermal bath is
naturally defined—and the comoving frame—in which the Hawking temperature is
naturally defined. We argue that this is an interesting advantage for optical
systems, absent in other analog systems for which horizons are at rest in the
lab. We also argue that the robustness to thermal fluctuations has an
interesting analogy to gravitational black holes.

Regarding previous works, the authors of Ref. [6] identified the fragility of
quantum entanglement to ambient noise and computed a threshold temperature for
the analog Hawking effect in Bose-Einstein condensates (BECs). A similar
analysis was performed in Refs. [7, 5] for optical analog systems assuming the
rather unphysical situation of a thermal bath at rest with the propagating
horizons. In either of these configurations, one does not find the enhanced
robustness that we elucidate in this paper because there is no relative boost
between the thermal bath and the horizons. That a relative boost can add
robustness to some aspects of the Hawking process in optical analogs was
pointed out qualitatively in Ref. [8], albeit with no reference to
entanglement.

The rest of this article is organized as follows. Section II contains a brief
summary of the Hawking effect in optical systems and includes references to
more extended treatments where the interested reader can find additional
details. Section III contains the main analysis of the paper—specifically, a
quantitative study of the complex entanglement spectrum generated in the
Hawking process and how such deteriorates in noisy environments. Section IV
contains a qualitative comparison with the physics of gravitational black
holes. Conclusions and take-home points are gathered in Section V. Throughout,
we adopt natural units in which $c=\hbar=1$.

##  II Optical analog white-black hole horizons

This section provides a summary of the Hawking process in optical analog
horizons [8]. The reader is referred to [8, 9] for further details.

Figure 1: Illustration of the structure of $in$, $int$, and $out$ modes for an
optical analog white-black hole in the comoving frame (adapted from [7]). The
analog white-black hole is generated by a strong electromagnetic pulse via the
Kerr effect. There are four ingoing modes (three arriving at the black hole
horizon and one at the white hole horizon), and four outgoing-modes. There are
two real, propagating interior modes (int-modes) between horizons; the other
two modes within this region are evanescent (i.e., they have complex
wavenumber, and their amplitudes are exponentially suppressed).

###  II.1 The system

Optical analogs [8, 10, 11, 12] rely on the Kerr effect, whereby a strong
electromagnetic pulse propagating in an optical medium modifies the local
refractive index. In turn, weak probe waves propagating thereon experience the
perturbed refractive index near the pulse. Probes initially faster than the
strong pulse will slow down when trying to overtake it and, for a strong
enough pulse, its rear end becomes an impenetrable barrier [8, 13, 14]. This
is the analog of a white hole horizon. Similarly, an analog black hole horizon
forms in the front end of the pulse. We thus have a pair of white- and black-
hole event horizons that share an interior. We consider here that the material
is a thin fiber and assume symmetry in the transverse $y$-$z$ plane.
Restricting further to one fixed polarization, e.g. along the $y$ axis, the
problem becomes effectively one dimensional [9].

From the perspective of an observer at rest in the laboratory, the pair white-
black hole propagates with the strong pulse, which moves at group velocity
$u$. The effective refractive index, therefore, depends on space and time,
$n_{\rm eff}=n+\delta n(x,t)$, where $n$ is the refractive index of the
material in absence of any pulse, and $\delta n(x,t)$ is the perturbation
generated by the strong pulse. It is more convenient to move to the frame
comoving with the pulse because, in this frame, the pulse is static and the
refractive index depends only on the spatial coordinate $\chi=\gamma(x-ut)$,
where $\gamma=1/\sqrt{1-u^{2}}$ is the standard Lorentz factor. This, in turn,
implies that frequencies of weak probes are conserved. We will denote
frequencies in the comoving frame by $\omega$ and the comoving wavenumber in
the direction of propagation as $k(\omega)$. Due to the spatial dependence of
the refrative index, the wavenumer $k$ is not conserved, and modes of same
frequency $\omega$ but different wavenumbers $k_{i}(\omega)$ can mix or
scatter; see Fig. 1 for an illustration.

Figure 2: Graphical solution to the dispersion relation in the comoving frame
for fixed comoving frequency $\omega$; see, e.g., Ref. [9] for explicit
details. Allowed comoving wavenumbers $k_{i}$ are labelled. The central
vertical and horizontal black lines indicate $k=0$ and $\omega^{\rm lab}=0$,
respectively. The symplectic norm of a mode is determined by the sign of its
lab frame frequency $\omega_{\rm lab}$ (shown in the vertical axis); thus, the
$k_{1}$ mode is the only “negative norm” mode in this model.

For concreteness, we consider the microscopic model proposed in Ref. [9]
(building upon previous work [15]), which describes a dielectric medium with
one Sellmeier pole (i.e., with a single absorption resonance; see, e.g., Ref.
[16] for an introduction to optical properties of dielectric media). Figure 2
illustrates a graphical resolution of the dispersion relation for a fixed
comoving frequency $\omega$ far away from the strong pulse. There are four
real solutions, which we denote by $k_{i}$, $i=1,2,3,4$. The modes $k_{1}$ and
$k_{4}$ are short-wavelength modes, while $k_{2}$ and $k_{3}$ are of the same
order as $\omega$. The modes $k_{1}$, $k_{2}$ and $k_{4}$ are left-movers in
the comoving frame (negative group velocity, using the convention that the
white hole is on the left of the black hole; see Fig. 1). The mode $k_{3}$
moves to the right. One important property of the modes is the sign of their
lab-frame frequency ${\rm sgn}({\omega}^{(\rm lab)})$, where ${\omega}^{(\rm
lab)}=\gamma(\omega+uk)$, because it determines the sign of the symplectic
norm of the modes.II.1II.1II.1In general, for all complex solutions
$f(\vec{x},t)$ to the equations of motion, the symplectic norm is defined from
the symplectic product via $(f_{1},f_{2})\equiv
i\,\int_{\Sigma}d^{3}\Sigma\,n^{\mu}\,(\bar{f}_{1}\nabla_{\mu}f_{2}-\nabla_{\mu}\bar{f}_{1}\,f_{2})$,
where $d^{3}\Sigma$ is the volume element of the equal-time hypersurface
$\Sigma$, and $n^{\mu}$ its future-directed unit normal (see, for instance
[17, 18]). The mode $k_{1}$ is the only mode with negative ${\omega}^{(\rm
lab)}$, as it can be seen directly from Fig. 2. This is crucial to understand
the origin of spontaneous particle creation.

If the system is illuminated with a (ingoing) wave packet sharply centered at
frequency $\omega$ and which has contribution from only one of the four
possible wave-numbers $k_{i}$, the scattered (outgoing) wave will be a wave
packet centered at the same frequency $\omega$ but will have contributions
from all possible wavenumbers. When the dynamics mixes modes with the same
sign of their symplectic norms, we have a simple scattering process; i.e.,
intensities of the incoming modes get distributed among the outgoing modes,
but the total intensity is conserved. However, when modes with different signs
of their symplectic norm mix (e.g., a positive norm mode and a negative norm
mode), the process classically corresponds to an amplifier, in the sense that
the scattered waves generically are more intense than the incident ones (with
energy supplied by the strong pulse, similar to how the mass sources Hawking
radiation in semi-classical black holes). Quantum mechanically, a scattering
process between two modes is described by a beam splitter, while an
amplification process corresponds to a two-mode squeezer. Two-mode squeezing
creates entangled quanta. Thus, for the optical analogs, any scattering
process involving the $k_{1}$ mode (the only mode with negative norm) will
generally lead to creation of entangled pairs.

###  II.2 Formalism

In a previous work [7], we have solved the dynamics of the system
numerically—again, using the model of Ref. [9]—and derived the evolution
matrix describing the dynamics (see Refs. [14, 13, 19, 20, 21, 22, 23, 24, 25,
26, 27] for numerical efforts in similar systems). The evolution takes a
simple form in the Heisenberg picture, as we now describe. Since comoving
frequency $\omega$ is conserved, it suffices to focus on one frequency at a
time. For each frequency $\omega$, there are eight total asymptotic modes—four
ingoing modes moving towards the horizons which scatter to the four outgoing
modes moving away from the horizons; see Fig. 1.

Let us denote by $\hat{a}_{k_{i}}^{\rm in}$ and $\hat{a}_{k_{i}}^{{\rm
in}\,\dagger}$ the annihilation and creation operators for wave-packet modes
peaked at each of the four ingoing modes. From them, we can construct
canonical pairs, $\hat{Q}_{k_{i}}\equiv\frac{1}{\sqrt{2}}\,(a_{k_{i}}^{\rm
in}+\hat{a}_{k_{i}}^{{\rm in}\,\dagger})$ and
$\hat{P}_{k_{i}}\equiv-\frac{i}{\sqrt{2}}\,(a_{k_{i}}^{\rm
in}-\hat{a}_{k_{i}}^{{\rm in}\,\dagger})$. We collect these operators in a
(column) vector

$\hat{\bm{R}}^{\rm(in)}\equiv\\\
\left(\hat{Q}_{k_{1}}^{\rm(in)},\hat{P}_{k_{1}}^{\rm(in)},\hat{Q}_{k_{2}}^{\rm(in)},\hat{P}_{k_{2}}^{\rm(in)},\hat{Q}_{k_{3}}^{\rm(in)},\hat{P}_{k_{3}}^{\rm(in)},\hat{Q}_{k_{4}}^{\rm(in)},\hat{P}_{k_{4}}^{\rm(in)}\right)^{\top},$
(II.1)

If we denote by $\hat{\bm{R}}^{\rm(out)}$ the vector similarly constructed
from the outgoing modes, Heisenberg evolution can be recasted as an expression
for $\hat{\bm{R}}^{\rm(out)}$ in terms of $\hat{\bm{R}}^{\rm(in)}$. Although
the existence of optical horizons originates in non-linear optics, the
evolution of weak probes is well approximated by linear equations, with the
non-linearities induced by the strong pulse all encoded in the optical
properties of the medium, in analogy with the quantum field theory in curved
spacetimes used in Hawking’s original derivation. This implies that the
relation we are looking for can be written in matrix form as

$\hat{\bm{R}}^{\rm(out)}=\bm{S}\hat{\bm{R}}^{\rm(in)}$ (II.2)

where the $8\times 8$ matrix $\bm{S}$ contains all the information about the
dynamics of the system. Furthermore, due to the linearity of the system, the
matrix $\bm{S}$ can be obtained by solving the classical equations of motion.
In other words, $\bm{S}$ contains the Bogoliubov coefficients between in and
out complex solutions to the equations of motion.

In what follows, we shall focus on a special class of states, Gaussian quantum
states (e.g., ground states of quadratic Hamiltonians, coherent states etc.),
as these are easy to generate and manipulate in the lab. Moreover, it can be
shown that linear dynamics [i.e., Eq. (II.2)] maps Gaussian states to Gaussian
states. We use many efficient techniques for dealing with Gaussian states and
refer to the reader to our previous work [7, 5] and Refs. [28, 29] for further
details about these techniques.

Gaussian states are completely characterized by their first and second moments
(similar to multivariate Gaussian probability distributions). The first
moments are codified in the vector
$\bm{\mu}\equiv\expectationvalue*{\hat{\bm{R}}}$, while the second moments are
conveniently encoded in the covariance matrix

$\displaystyle\bm{\sigma}^{ij}$
$\displaystyle\equiv\expectationvalue{\left\\{\hat{\bm{R}}^{i}-\bm{\mu}^{i},\hat{\bm{R}}^{j}-\bm{\mu}^{j}\right\\}},$
(II.3)

where the expectation value is computed with respect to the quantum state
under consideration (either pure or mixed), and $\\{\cdot,\cdot\\}$ denotes
the symmetric anti-commutator. For a Gaussian state, rather than working with
a density operator $\hat{\rho}$, we can equivalently describe the state by a
pair $(\bm{\mu},\bm{\sigma})$. For linear systems, evolution preserves the
Gaussian character of the state, and the relation between the in and out state
is given by the matrix $\bm{S}$ discussed above

$\displaystyle\bm{\mu}^{\rm(out)}$
$\displaystyle=\bm{S}\cdot\bm{\mu}^{\rm(in)},$ (II.4)
$\displaystyle\bm{\sigma}^{\rm(out)}$
$\displaystyle=\bm{S}\cdot\bm{\sigma}^{\rm(in)}\cdot\bm{S}^{\top}.$ (II.5)

Our primary focus in this work is on entanglement in optical analogs. To
quantify such, we use the logarithmic negativity (LN) [30, 31]. This is an
easily-computable measure of entanglement for pure states and mixed states.
For Gaussian quantum states (and other types of quantum states), the value of
the LN has an operational meaning as the exact costII.2II.2II.2The cost refers
to exhange of a ‘currency’—which is entangled bits, or ebits—where 1 ebit is
equal to the entanglement contained in 1 Bell pair. that is required to
prepare or simulate the quantum state under consideration [32, 33].
Importantly, when restricting to Gaussian states and when one of the two
subsystems is made of a single mode, the LN is greater than zero if and only
if the state is entangled—regardless of the size of the other subsystem—and is
larger for states with more entanglement. In other words, the LN is a faithful
quantifier of entanglement for 1 versus $N$ mode bipartitions [28, 29]. (We
later apply this technical fact to optical event horizons to find threshold
lab temperatures at which entanglement between subsystems strictly vanishes.)

Given a system $A$ of $N$ modes and system $B$ of $M$ modes occupying a
Gaussian quantum state $\rho_{AB}$ with covariance matrix $\bm{\sigma}_{AB}$,
the LN can be computed from the symplectic eigenvalues $\tilde{\nu}_{I}$ of
the partially transposed covariance matrix
${\widetilde{\bm{\sigma}}_{AB}\equiv(\bm{I}_{A}\oplus\bm{T}_{B})\bm{\sigma}_{AB}(\bm{I}_{A}\oplus\bm{T}_{B})}$,
where $\bm{T}_{B}=\bigoplus_{k=1}^{M}\big{(}\begin{smallmatrix}1&0\\\
0&-1\end{smallmatrix}\big{)}$ (see [28, 29] for further details),

${\rm
LN}(\rho_{AB})=\sum_{J=1}^{M+N}\max\left[0,-\log_{2}(\tilde{\nu}_{J})\right].$
(II.6)

###  II.3 Analytical description of dynamics

To finish this summary, we consider a convenient way of describing the
dynamics of the optical system under consideration. In previous works [7, 5],
we showed that the scattering process of the four ingoing modes to the four
outgoing modes can be accurately described by a specific combination of two-
mode squeezers and beam splitters; see Fig. 3. Each squeezer is associated
with Hawking particle production from either the black hole horizon or the
white hole horizon, which originates from scattering of the modes $k_{1}$ and
$k_{4}$ to the modes $k_{1}$ and $k_{3}$. Since $k_{1}$ has negative
symplectic norm, this scattering process is described by a two-mode squeezer.
Each beam splitter describes the scattering between the modes $k_{2}$ and
$k_{3}$ (and results in greybody factors). In realistic situations, such
mixing is extremely weak and can be neglected, but we maintain this scattering
in our analysis for completeness. The mixing between modes $k_{2}$ and $k_{1}$
(and likewise for modes $k_{2}$ and $k_{4}$) is negligibly small, due to the
large differences in their wavelengths (see, e.g., Table 1).

Figure 3: Symplectic-circuit of the Hawking process for an optical white-black
hole pair (adapted from [7]). The two-mode squeezer (blue box) and beam
splitter (gray rectangle) on the left correspond to the black hole, while the
elements to the right correspond to the white hole horizon.

Taking all this into account, the evolution from the ingoing to the outgoing
modes can be described by the circuit of Fig. 3—we call this a symplectic
circuit because it is constructed from simple symplectic transformations [7,
5]. From the symplectic circuit, it is straightforward to write an analytical
expression for the evolution matrix of the system, ${\bm{S}}$, in terms of the
parameters describing the squeezers and beam splitters; see Appendix A for
details. Such parameters can be extracted from numerical simulations. The
circuit in Fig. 3 and the analytical expression for $\bm{S}$ is extremely
useful to understand many aspects of the physics of the system. The full
scattering matrix describing the transformation from the four ingoing modes to
the four outgoing modes is

$\bm{S}_{\rm WB}=\\\
\begin{pmatrix}(1+\cos^{2}\theta\sinh^{2}r_{H})\bm{I}_{2}&\cos\theta\sin\theta\sinh
r_{H}\bm{\sigma}_{z}&-\cos\theta\sinh
r_{H}\bm{\sigma}_{z}&\cos^{2}\theta\cosh{r_{H}}\sinh{r_{H}}\bm{\sigma}_{z}\\\
-\cos\theta\sin\theta\sinh{r_{H}}\bm{\sigma}_{z}&\cos^{2}\theta\bm{I}_{2}&\sin\theta\bm{I}_{2}&-\cos\theta\sin\theta\cosh{r_{H}}\bm{I}_{2}\\\
\cos\theta\sinh{r_{H}}\bm{\sigma}_{z}&\sin\theta\bm{I}_{2}&0&\cos\theta\cosh{r_{H}}\bm{I}_{2}\\\
-\cos^{2}\theta\cosh{r_{H}}\sinh{r_{H}}\bm{\sigma}_{z}&-\cos\theta\sin\theta\cosh{r_{H}}\bm{I}_{2}&\cos\theta\cosh{r_{H}}\bm{I}_{2}&(\sin^{2}\theta-\cos^{2}\theta\sinh^{2}r_{H})\bm{I}_{2}\end{pmatrix},$
(II.7)

where the subscript WB refers to “white-black” hole, $\bm{I}_{2}$ is the
$2\times 2$ identity matrix, and $\bm{\sigma}_{z}$ is the Pauli-Z matrix. The
squeezing intensity $r_{H}$ is a function of frequency and, for physical
setups where the analogy with the Hawking effect is on firm ground,
$r_{H}(\omega)={\rm arctanh}\,e^{-\omega/(2T_{H})}$, where $T_{H}$ is the
Hawking temperature associated to the event horizons. The angle $\theta$ is
also frequency-dependent and is related to the greybody factor $\Gamma$
(typically introduced to describe classical scattering in the Hawking effect
of black holes) via $\Gamma=\cos^{2}\theta$. The scattering matrix
$\bm{S}_{\rm WB}$ assumes that the strong pulse generating the horizons is
symmetric, which amounts to say that their Hawking temperatures and greybody
factors are the same. (It is straightforward to drop this assumption though.)

###  II.4 Physical setup

In our numerical simulations, we model the perturbation of the refractive
index by a sech profile of the form ${\delta n(x,t)=\delta
n_{0}\,\sech^{2}\left(\frac{t-x/u}{D}\right)}$. This is a common and
convenient choice [8, 14, 10]. In this expression, $u$ is the group velocity
of the perturbation caused by the strong pulse, and $x$ and $t$ are space-time
coordinates in the lab frame. This profile depends on two real positive
numbers, $\delta n_{0}$ and $D$, which determine the amplitude and width of
the perturbation, respectively. References [7, 34] have explored values of
$\delta n_{0}\in[0.01,0.1]$ and $D\in[2,10]$ fs. Away from these ranges, one
does not have a horizon for all frequencies or dispersive effects are large
and jeopardize the analogy with the Hawking effect. For the calculations in
this paper, we will choose the concrete values $\delta n_{0}=0.05$ and
$D=6\,\rm{fs}$. Other choices within the ranges specified above do not change
our conclusions. Our tools permit to extend the analysis beyond these ranges,
but analyzing such lacks motivation since we would be away from the analogy
with the Hawking effect.

For $\delta n_{0}=0.05$ and $D=6\,\rm{fs}$, the velocity of the horizons is
$u=0.41$, and we obtain a value of $T_{H}=0.00046$ PHz ($T_{H}=3.51$ K if we
restore the values of $\hbar$ and Boltzman’s constant $k_{B}$) for the Hawking
temperature. We extract this temperature form numerical simulations [7,
34].II.3II.3II.3Reference [9] contains an analytical expression for $T_{H}$ in
terms of the parameters of the microphysical model, which provides an
excellent approximation. Table 1 contains the values of $k_{i}$ and lab
frequencies for three representative values of the comoving frequency
$\omega$, namely $\omega/T_{H}=10^{-1},\,1\,,10$. Observe the large blue-shift
that the comoving frequency $\omega$ suffers when translated to the laboratory
frame. For the modes $k_{1}$, $k_{3}$ and $k_{4}$ (the three modes primarily
involved in the Hawking process), we have $\omega^{\rm lab}_{i}\gg\omega$. As
a consequence, one would need temperatures in the lab frame much larger than
the Hawking temperature $T_{H}$ in order to significantly populate these modes
with thermal quanta.

Although it is tempting to associate the large blue-shift exclusively with the
Lorentz boost, dispersive effects play a key role. Notice that, for $u=0.41$,
we have $\gamma=1.1$, which is close to unity. Consequently, if the system had
the vacuum dispersion relation $k=\pm\omega$, the expression $\omega^{\rm
lab}_{k_{i}}=\gamma\,(\omega+u\,k_{i})$ would imply that $\omega^{\rm
lab}_{k_{i}}$ is of the same order as the comoving frequency $\omega$. For the
optical system that we are considering, we find instead that $k_{1}$, $k_{3}$
and $k_{3}$ are all much larger than $\omega$ (see Table 1), due to the
bending of the dispersion relation away from the straight line $k=\pm\omega$;
see Fig. 2. In turn, this causes $\omega^{\rm lab}_{i}$ to be much larger than
$\omega$ for $i=1,3,4$. On the other hand, the mode $k_{2}$ lives in a region
in which the dispersion is close to vacuum dispersion and, consequently, its
lab frequency is similar to its comoving frequency. In other words, there is
no blue-shift for the $k_{2}$ mode; but the contribution of this mode to the
Hawking process is faint. The results that we obtain in the next section
originate from the blue-shift of $\omega^{\rm lab}_{i}$ (compared to $\omega$)
for $i=1,3,4$.

$\omega$ | $(k_{1},\omega^{\rm lab}_{k_{1}})$ | $(k_{2},\omega^{\rm lab}_{k_{2}})$ | $(k_{3},\omega^{\rm lab}_{k_{3}})$ | $(k_{4},\omega^{\rm lab}_{k_{4}})$
---|---|---|---|---
$10^{-1}$ |  $(-1.62,-.737)\times 10^{4}$ | $(-1.41,.459)\times 10^{-1}\,\,$ | $(1.15,.530)\times 10\,\,$ | $(1.62,.737)\times 10^{4}$
$1$ |  $(-1.63,-.739)\times 10^{4}$ |  $(-1.41,.459)$ |  $(1.15,.530)\times 10^{2}$ |  $(1.62,.735)\times 10^{4}$
$10$ |  $(-1.68,-.761)\times 10^{4}$ |  $(-1.41,.459)\times 10$ |  $(1.15,.532)\times 10^{3}$ |  $(1.56,.709)\times 10^{4}$
Table 1: Comoving wavenumbers $k_{i}$ and lab frequencies
$\omega_{k_{i}}^{(\rm lab)}$ for three representative values of the comoving
frequency $\omega$. All quantities in units of the Hawking temperature
$T_{H}=4.6\times 10^{-4}$ PHz (in Kelvin, $T_{H}=3.51$ K). The numbers here
are for a refractive index perturbation with $\delta n_{0}=0.05$ and
$D=6\,\rm{fs}$, for which $u=0.41$.

##  III Thermal noise in the laboratory frame

We evaluate the impact that initial thermal fluctuations have on entanglement.
We focus on thermal fluctuations since they are a common noise source in
experimental settings. As explained above, we focus on situations in which the
thermal bath is at rest in the lab frame and, thus, at rest with the optical
medium. Due to the form of the dispersion relation, each mode $k_{i}$
($i=1,2,3,4$) has a different population of thermal quanta; for a simplified
analysis with isotropic noise, corresponding to a thermal bath at rest in the
comoving frame, see Ref. [5].

###  III.1 Ingoing and outgoing state

In general, a thermal state for a quantum system of $N$ modes is a mixed state
with density operator given by $\hat{\rho}=\,e^{-\beta\,\hat{H}}/Z$, where
$Z={\rm Tr}e^{-\beta\hat{H}}$ is a normalization factor (the partition
function), $\hat{H}$ is the Hamiltonian of the system (which, here, is
quadratic in the canonical variables), and $\beta=1/T$ the inverse
temperature. This is, in fact, a Gaussian state and can be simply described by
its first moments and its covariance matrix. The first moments of a thermal
state are zero, and its covariance matrix is $\bm{\sigma}_{\rm
th}=\bigoplus_{i=1}^{N}\,(1+2\,n_{\rm env;i})\,\bm{I}_{2}$, where $n_{\rm
env;i}$ is the mean number of quanta for the mode $i$. For the optical system
under consideration, the ingoing state of the modes is described by

$\displaystyle\bm{\mu}_{\rm th}^{\rm(in,\,lab)}$ $\displaystyle=\bm{0},$
(III.1) $\displaystyle\bm{\sigma}_{\rm th}^{(\rm in,\,lab)}$
$\displaystyle=\bigoplus_{i=1}^{4}\,(1+2\,n^{(\rm lab)}_{\rm
env;i})\bm{I}_{2}$ (III.2)

in the laboratory frame, where $n^{(\rm lab)}_{\rm env;i}$ are given by the
Bose-Einstein formulaIII.1III.1III.1We write $|\omega^{\rm lab}_{k_{i}}|$
rather than just $\omega^{\rm lab}_{k_{i}}$ because the mode $k_{1}$ has
negative lab frequency for positive comoving frequency $\omega$. To compute
the number of physical quanta, $|\omega^{\rm lab}_{k_{i}}|$ must be used in
this expression.

$n^{\rm(lab)}_{\rm env;i}=\frac{1}{e^{|\omega^{\rm lab}_{k_{i}}|/T_{\rm
env}^{\rm(lab)}}-1}\,.$ (III.3)

The evolution problem of the last section was formulated in the comoving
frame, for which frequencies $\omega$ are conserved. Therefore, we must
translate the thermal state from the lab to the comoving frame. It is well
known that the absence of correlations makes this an easy task; all we must do
is transform $n^{\rm(lab)}_{\rm env;i}$ to the comoving frame. In doing so,
note that $n^{\rm(lab)}_{\rm env;i}$ is the number of quanta per volume in
phase space. In other words, $n^{\rm(lab)}_{\rm env;i}\,{\rm d}^{3}x\,{\rm
d}^{3}k$ is the number of quanta within ${\rm d}^{3}x$ and with wavenumber
inside a volume ${\rm d}^{3}k$ centered at $\vec{k}$. Volumes in phase space
are Lorentz invariant, which automatically implies that $n^{\rm(lab)}_{{\rm
env};i}$ is a Lorentz scalar (see, e.g., [35]). Hence, the number of quanta in
the comoving frame in the mode $k_{i}$ is given by $n^{\rm(co)}_{\rm
env;i}=n^{\rm(lab)}_{\rm env;i}$. We can thus describe the quantum state in
the comoving frame (via $\bm{\mu}^{\rm(in,\,co)}_{\rm th}$ and
$\bm{\sigma}^{\rm(in,\,co)}_{\rm th}$) by simply relabeling $n_{\rm
env;i}^{\rm(lab)}$ as $n_{\rm env;i}^{\rm(co)}$ in Eqs. (III.1) and (III.2).

Observe that the number of noisy quanta in the mode $k_{i}$, for a thermal
bath at rest in the lab frame, is determined solely by the ratio
$|\omega^{\rm(lab)}_{i}|/T^{\rm(lab)}_{\rm env}$. This is interesting because,
although the modes $k_{i}$ have the same comoving frequency $\omega$, they
have very different lab frequencies (see Table 1). This in turn implies that
the modes are initially populated with noisy quanta in an exponentially
different manner. In particular, the values of lab frequencies written in
Table 1 result in the hierarchy $n^{\rm(co)}_{\rm env;1}\approx
n^{\rm(co)}_{\rm env;4}\ll n^{\rm(co)}_{\rm env;3}\ll n^{\rm(co)}_{\rm
env;2}$. This hierarchy affects the entanglement between different subsystems,
as we discuss in the next section.

Figure 4: Entanglement versus lab temperature $T_{\rm env}^{(\rm lab)}$ for
various bipartitions evaluated at $\omega/T_{H}\approx 1$. Here
$T_{H}=4.\times 10^{-4}$ PHz (in Kelvin, $T_{H}=3.51$ K). Solid curves
correspond to numerically obtained values. Dotted curves correspond to the
analytic circuit approximation. The agreement is excellent. Vertical dashed
lines label the threshold lab temperature $T_{\rm env}^{\star\,(\rm lab)}$ at
which entanglement vanishes for that bipartition; see second row of Table 2.

Finally, given an initial thermal state in the lab frame with temperature
$T_{\rm env}^{(\rm lab)}$ and using the general transformations (II.4) and
(II.5), the outgoing state, in the comoving frame, is characterized by

$\displaystyle\bm{\mu}^{\rm(out,\,co)}$ $\displaystyle=\bm{S}_{\rm
WB}\cdot\bm{\mu}^{\rm(in\,co)}_{\rm th}=\bm{0},$ (III.4)
$\displaystyle\bm{\sigma}^{\rm(out,\,co)}$ $\displaystyle=\bm{S}_{\rm
WB}\cdot\bm{\sigma}^{\rm(in,\,co)}_{\rm th}\cdot\bm{S}_{\rm WB}^{\top}.$
(III.5)

The covariance matrix $\bm{\sigma}^{(\rm out,\,co)}$ encodes information about
all four outgoing modes; though, we often focus on particular subsystems by
looking at the reduced moments—i.e., sub-matrices of $\bm{\sigma}^{(\rm
out,co)}$. From the output covariance matrix, we can extract information about
the number of particles created (see Appendix B for analytical expressions for
these quantities), the amount of entanglement generated in the scattering
process (see below), etc. We remark that, although entanglement calculations
in next section are done in the comoving frame, results do not change by going
back to the lab frame because a Lorentz boost does not mix the modes.

###  III.2 Entanglement degradation

$\omega$ |  $(k_{2},k_{3},k_{4}|\,k_{1})$ |  $(k_{4},k_{1},k_{2}|\,k_{3})$ |  $(k_{3},k_{4}|\,k_{1})$ |  $(k_{4},k_{1}|\,k_{3})$ |  $(k_{1}|\,k_{2})$ |  $(k_{1}|\,k_{3})$ |  $(k_{1}|\,k_{4})$
---|---|---|---|---|---|---|---
$10^{-1}$ |  $1.41\times 10^{6}$ |  $1.40\times 10^{5}$ |  $9.51\times 10^{3}$ |  $620$ |  $3.86\times 10^{-3}$ |  $2.30$ |  $1.74\times 10^{3}$
1 |  $2.22\times 10^{4}$ |  $1.41\times 10^{4}$ |  $8.06\times 10^{3}$ |  $4.42\times 10^{3}$ |  $4.42\times 10^{-2}$ |  $1.20\times 10^{2}$ |  $1.64\times 10^{3}$
10 |  $1.50\times 10^{3}$ |  $1.50\times 10^{3}$ |  $1.41\times 10^{3}$ |  $1.41\times 10^{3}$ |  $1.27$ |  $1.40\times 10^{3}$ |  $340$
Table 2: Threshold lab temperatures $T_{\rm env}^{\star\,(\rm lab)}$ at which
entanglement for various bipartitions vanishes for three representative values
of the comoving frequency $\omega$. All quantities in units of the Hawking
temperature $T_{H}=4.6\times 10^{-4}$ PHz (in Kelvin $T_{H}=3.51$ K). The
label $(A|B)$ represents a bipartite cut across the modes $A$ and modes $B$
for which entanglement is evaluated.

We study how ambient thermal noise in the lab frame degrades the Hawking-
generated entanglement and, in particular, determine the threshold temperature
above which entanglement is extinguished. This is a rich problem because the
threshold temperature depends on the comoving frequency $\omega$. This
dependence is easy to understand because the intensity of the Hawking
squeezers, $r_{H}={\rm arctanh}\,e^{-\omega/(2T_{H})}$, falls exponentially
with $\omega$; hence, entanglement is weaker for modes with large $\omega$ and
is consequently more fragile. Moreover, entanglement and its sensitivity to
noise also depends on which bipartitions and subsystems of the outgoing state
we consider. We report results for three representative values of the comoving
frequency: $\omega/T_{H}=10^{-1},\,1,\,10$. For completeness, we show results
obtained from both numerical resolution of the equations of the model for the
dielectric medium [9], as well as from the analytical expressions derived from
Eq. (II.7). As shown in Fig. 4, both approaches agree well when computing the
entanglement.

In Fig. 4, we plot the LN versus the lab temperature $T_{\rm env}^{\rm(lab)}$
at $\omega/T_{H}\approx 1$ for various subsystems of modes, where, e.g., we
use the notation $(k_{1}|k_{4})$ to represent a bipartite cut across the
reduced subsystem containing only the Hawking pairs emitted by the white hole
horizon; similarly, $(k_{1},k_{2},k_{4}|k_{3})$ represents the cut
partitioning the white hole exterior and black hole exterior, etc. The
bipartition with the largest amount of entanglement is
$(k_{2},k_{3},k_{4}|k_{1})$, which partitions the positive norm modes
$(k_{2},k_{3},k_{4})$ from the negative norm mode $k_{1}$. The subystem with
the least amount of entanglement is $(k_{1}|k_{2})$, which we do not show on
the plot as the entanglement in this subsystem is five orders of magnitude
smaller than the rest.

For each subsystem, there exists a threshold lab temperature $T_{\rm
env}^{\star\,\rm(lab)}$ at which point entanglement vanishes within that
subsystem. The values of $T_{\rm env}^{\star\,\rm(lab)}$ are indicated on Fig.
4 by dashed vertical lines. We also give explicit values for $T_{\rm
env}^{\star\,\rm(lab)}$ (in units of $T_{H}=3.51$ K) in Table 2 for three
representative comoving frequencies $\omega/T_{H}=10^{-1},1,10$. The threshold
temperatures drastically vary across the different subsystems and for
different frequencies, which highlights the richness of the entanglement
spectrum. Generally though, entanglement is more robust for low frequencies.
The lowest threshold temperature is found for the subsystem $(k_{1}|k_{2})$;
the reason being that the entanglement in this subsystem is the weakest, as
such entanglement does not originate from the Hawking effect. For various
subsystems whose entanglement originates from the Hawking process, and at
various frequencies (low and high), we observe threshold temperatures $T_{\rm
env}^{\star\,\rm(lab)}/T_{H}\sim\order{10^{3}}$ or larger.

##  IV Comparison with gravitational black holes

Figure 5: Illustration depicting the correspondence between reference frames
for an astrophysical black hole (ff frame and $\mathcal{J}^{+}$) and an
optical analog (lab frame and comoving frame). Blue-shift between ingoing and
outgoing frequencies is also shown. For the ff frame, the black hole event
horizon is in motion; similarly, for the optical analog, the white-black hole
(i.e., the pulse) is in motion with respect to the lab frame. For observers at
$\mathcal{J}^{+}$, the black hole is stationary; similarly, for the optical
analog, the white-black hole is stationary in the frame comoving with the
pulse.

The robustness to thermal fluctuations is reminiscent of what happens for
astrophysical black holes. In the astrophysical case, the progenitors of the
Hawking radiation are modes supported on ultra-high frequencies as measured by
observers in the distant past (at past null infinity, $\mathcal{J}_{-}$).
Therefore, any thermal fluctuations in the past would be incapable of
significantly affecting the outgoing Hawking spectrum (measured by observers
at future null infinity, $\mathcal{J}_{+}$) or its entanglement with the black
hole, unless those thermal fluctuations had an ultra-high temperature (larger
than the Planck temperature) in the natural frame at $\mathcal{J}_{-}$. To
better elucidate these points, in both the astrophysical and analog settings,
we examine the scenario from the viewpoint of preferred reference frames.

For astrophysical black holes, there are two preferred frames—the ingoing and
outgoing frames. The ingoing frame is the natural frame to set up the initial
state; it is equivalent to a freely falling (ff) frame (see, e.g., Refs. [36,
37, 38]) defined by an observer which crosses the horizon at the same instance
the surface of the collapsing object does. The outgoing frame at future null
infinity, $\mathcal{J}_{+}$ is the natural one to evaluate the Hawking
radiation. There is an exponential blue-shift between the ingoing and outgoing
frequencies (defined in the ingoing and outgoing frames respectively),
$\omega^{\rm(in)}\sim e^{\kappa u}\,\omega^{\rm(out)}$, where $\kappa$ is the
surface gravity of the black hole (proportional to $T_{H}$), and $u$ is the
retarded time along $\mathcal{J}_{+}$; see Fig. 5 for an
illustration.IV.1IV.1IV.1Recall that Hawking quanta reach $\mathcal{J}_{+}$
for asymptotically large $u$. The blue-shift between these frames has a
gravitational origin; wave packets reaching $\mathcal{J}_{+}$ at late $u$,
when propagated back in time towards $\mathcal{J}_{-}$, come arbitrarily close
to the horizon, blue-shifting as they approach the horizon. The trip from the
horizon to $\mathcal{J}_{-}$ does not undo the blue-shift. Thus, it is the
time asymmetry induced by gravitational collapse that causes an exponentially
large blue-shift.

The analogs to the ingoing and outgoing frames in the optical systems
discussed in this article are the lab and comoving frames, respectively; see
Fig. 5 for a visual. For the optical system, there is a large blue-shift
between the lab frame and comoving, similar to the astrophysical black hole.
The difference from the astrophysical case is, of course, the physical origin
of the blue-shift. For optical systems, the blue-shift originates from the
combined effects of a Lorentz boost and dispersion. Moreover, contrary to the
astrophysical case, the blue-shift between lab and comoving frames is not
exceedingly (exponentially) large. However, for all practical purposes, it is
still large enough to protect the Hawking effect from spurious thermal
fluctuations (see Tables 1 and 2).

##  V Conclusions

Optical systems present an intriguing setting to recreate the physics of the
Hawking effect. They have attractive advantages, such as current experimental
capabilities to manipulate individual photons and to measure their
entanglement. These systems also have the peculiar feature that the analog
horizons move with respect to the lab frame at a fraction of the speed of
light. We have shown in this paper that moving horizons actually provides an
advantage: it naturally protects entanglement generated during the Hawking
process against background thermal noise.

The robustness of entanglement to harmful thermal fluctuations derives from an
interplay between dispersive effects and the Lorentz boost between the lab
frame and the frame comoving with the horizons. Modes with comoving frequency
$\omega$ that are of the order of the Hawking temperature $T_{H}$ appear to be
high frequency modes as seen in the lab frame due to the boost. Consequently,
the lab temperatures needed to appreciably affect Hawking-generated
entanglement are several orders of magnitude larger than $T_{H}$.

Such robustness is absent for systems with analog horizons that are at rest in
the laboratory frame, like BECs. For BEC analogs, both temperatures $T_{H}$
and $T_{\rm env}$ are defined in the same frame and, consequently, $T_{\rm
env}\approx T_{H}$ suffices to decohere the outgoing radiation and eliminate
all quantum traces of the Hawking effect [6]; in which case, the final state
can then be accounted for by a classical process of amplification of thermal
radiation. One is required to work at low temperatures $T_{\rm env}\lesssim
T_{H}$ for quantum features to survive.

In this sense, optical systems share some analogy with astrophysical black
holes. For the latter, the entanglement produced during Hawking’s evaporation
is shielded against background thermal fluctuations (such as the cosmic
microwave background radiation) for which $T_{\rm env}\gg T_{H}$ for black
holes of astrophysical origin. One needs an exponentially high background
temperature to noticeably populate the Hawking progenitors and thus decohere
astrophysical Hawking particles from their partners. The difference is that,
for astrophycial black holes, the protection originates entirely from
gravitational blue-shifting; while for optical analogs, it originates from an
interplay of blue-shifting (due to a Lorentz boost) and dispersive effects.

We conclude that entanglement generated by optical analog event horizons is
extremely robust to background thermal fluctuations in the lab frame. This
observation—together with the fact that precise manipulation of quantum states
of light is a staple of modern quantum technologies [39, 40]—indicates that
optical analogs [8, 10, 11, 12] are promising candidates for observing the
genuinely quantum nature of the Hawking effect in a practical setting,
hopefully in the near future.

###### Acknowledgements.

We have benefited from discussions with A. Bhardwaj, D. Bermudez, M. Jacquet,
D. Seehy and J. Wilson. I.A. and D.K. are supported by the NSF grant
PHY-2110273, and by the Hearne Institute for Theoretical Physics. A.J.B.
acknowledges discussion with Quntao Zhuang and support from the Defense
Advanced Research Projects Agency (DARPA) under Young Faculty Award (YFA)
Grant No. N660012014029.

##  Appendix A Symplectic matrix for white-black hole

To construct the symplectic matrix for the white-black hole, $\bm{S}_{\rm WB}$
of Eq. (II.7), we can regard the entire scattering process as a succession of
four distinct events. [For justification of this, see Refs. [7, 5].]

High-frequency progenitor modes $k_{1}^{\rm(in)}$ and $k_{4}^{\rm(in)}$
approach the black hole side of the strong pulse (affected negligibly by a
dispersion-induced potential barrier, due to their high frequencies) and
participate in a two mode squeezing process converting into Hawking pairs
$k_{3}^{(\rm out)}$ and $k_{1}^{(\rm int)}$, where $k_{3}^{(\rm out)}$ denotes
the outgoing Hawking radiation and $k_{1}^{(\rm int)}$ denotes the entangled
Hawking partner that falls into the white-black hole. This pair creation event
is described by the symplectic matrix

$\bm{S}_{1}=\begin{pmatrix}\cosh{r_{H}}\,\bm{I}_{2}&0&0&\sinh{r_{H}}\,\bm{\sigma}_{z}\\\
0&\bm{I}_{2}&0&0\\\ 0&0&\bm{I}_{2}&0\\\
\sinh{r_{H}}\,\bm{\sigma}_{z}&0&0&\cosh{r_{H}}\,\bm{I}_{2}\end{pmatrix},$
(A.1)

where the ordering of elements follows the convention in Eq. (II.1).

As the outgoing Hawking radiation, $k_{3}^{(\rm out)}$ (a lower frequency
mode), propagates outward, it meets a potential barrier where it mixes with
the low-frequency mode $k_{2}^{\rm(in)}$, partly getting scattered back into
the white-black hole as $k_{2}^{\rm(int)}$. This scattering process is
described by a symplectic matrix

$\displaystyle\bm{S}_{2}$ $\displaystyle=\begin{pmatrix}\bm{I}_{2}&0&0&0\\\
0&\cos{\theta_{\omega}}\bm{I}_{2}&0&-\sin\theta_{\omega}\bm{I}_{2}\\\
0&0&\bm{I}_{2}&0\\\
0&\sin\theta_{\omega}\bm{I}_{2}&0&\cos{\theta_{\omega}}\bm{I}_{2}\end{pmatrix},$
(A.2)

where $\cos^{2}\theta$ is the greybody factor of the potential barrier (i.e.,
its transmission probability). Taken together, $\bm{S}_{2}\bm{S}_{1}$ denotes
the dynamics induced by the black hole event horizon.

Following the scattering processes induced by black hole, we have the white
hole dynamics. Intuitively, the white hole dynamics follows the black hole
dynamics because the black hole pair creation mechanism actually sources the
white hole scattering processes. Indeed, the back-scattered mode
$k_{2}^{\rm(int)}$ traverses through the interior of the white-black hole
unscathed. However, as it exits the white hole event horizon, it meets the
$k_{3}^{\rm(in)}$ mode (the time-reverse of the outgoing Hawking radiation) at
a potential barrier and scatters. This scattering process is described by the
matrix

$\bm{S}_{3}=\begin{pmatrix}\bm{I}_{2}&0&0&0\\\
0&\cos{\theta_{\omega}}\bm{I}_{2}&\sin\theta_{\omega}\bm{I}_{2}&0\\\
0&-\sin\theta_{\omega}\bm{I}_{2}&\cos{\theta_{\omega}}\bm{I}_{2}&0\\\
0&0&0&\bm{I}_{2}\end{pmatrix},$ (A.3)

Observe the sign reversal of the angle $\theta$, which is due to the fact
that, for symmetric pulses, the white hole can be understood as the time-
reverse of the black hole.

The portion of ingoing mode $k_{3}^{\rm(in)}$ that makes it to the white-hole
event horizon then interacts with the $k_{1}^{\rm(int)}$ from the black hole
pair creation process by a two-mode squeezing interaction,A.1A.1A.1Locally,
the $k_{1}^{\rm(int)}$ mode occupies a thermal state at temperature $T_{H}$;
thus, the white hole pair creation process is stimulated by the Hawking
partners generated by the black hole.

$\bm{S}_{4}=\begin{pmatrix}\cosh{r_{H}}\bm{I}_{2}&0&-\sinh{r_{H}}\bm{\sigma}_{z}&0\\\
0&\bm{I}_{2}&0&0\\\ 0&0&0&\bm{I}_{2}\\\
-\sinh{r_{H}}\bm{\sigma}_{z}&0&\cosh{r_{H}}\bm{I}_{2}&0\end{pmatrix}.$ (A.4)

Again, the relative sign between $S_{1}$ and $S_{4}$ is due to the inverse
character of the white hole relative to the black hole. The product
$\bm{S}_{3}\bm{S}_{4}$ denotes the white hole scattering process.

Altogether, $\bm{S}_{\rm WB}=\bm{S}_{4}\bm{S}_{3}\bm{S}_{2}\bm{S}_{1}$, which
reduces to Eq. (II.7); see Fig. 1 for a visual aid and Fig. 3 for a circuit
description of the processes.

##  Appendix B Mean number of quanta

When analyzing particle pair production in analog gravity or quantum fields in
curved spacetime, one often computes the mean number of quanta produced. For
completeness, we do so here for the optical white-black hole pair and an
thermal initial state.

Quite generally, starting from the output covariance matrix
$\bm{\sigma}^{\rm(out)}$, we can compute the outgoing number of quanta in the
mode $k_{i}^{\rm(out)}$ via

$\expectationvalue{\hat{n}^{\rm(out)}_{i}}=\frac{1}{4}\text{Tr}\\{\bm{\sigma}_{i}^{\rm(red)}\\}+\frac{1}{2}\bm{\mu}^{\rm(red)\,\top}\bm{\mu}^{\rm(red)}-\frac{1}{2}\,.$
(B.1)

where $\bm{\mu}_{i}^{\rm(red)}$ and $\bm{\sigma}_{i}^{\rm(red)}$ are the sub-
vector and sub-matrix of $\bm{\mu}^{\rm(out)}$ and $\bm{\sigma}^{\rm(out)}$
for the mode $k_{i}$.B.1B.1B.1The relation between the mean number of quanta
and the first and second moments [Eq. (B.1)] holds for any quantum state, not
just Gaussian states. Though the state under consideration is indeed a
Gaussian state. Assuming that the modes are initially populated with thermal
quanta and utilizing Eqs. (III.4) and (III.5), it is straightforward to obtain
the following analytical expressions for the outgoing quanta in each mode.
First, for brevity, define $N_{\rm env;k}\equiv(1+2n_{\rm env;k})$. Then,

$\expectationvalue{\hat{n}^{\rm(out)}_{1}}=\frac{1}{2}\bigg{[}-1+N_{\rm
env;1}\cosh^{4}r_{H}+\sinh^{2}r_{H}\Big{[}\cosh^{2}r_{H}\Big{(}N_{\rm
env;4}\cos^{4}\theta-2N_{\rm env;1}\sin^{2}\theta\Big{)}\,\\\
+\cos^{2}\theta\,\Big{(}N_{\rm env;3}+N_{\rm
env;2}\sin^{2}\theta\Big{)}\Big{]}+N_{\rm
env;1}\sin^{4}\theta\sinh^{4}r_{H}\bigg{]},$ (B.2)
$\expectationvalue{\hat{n}^{\rm(out)}_{2}}=\frac{1}{2}\bigg{[}-1+N_{\rm
env;2}\cos^{4}\theta+N_{\rm
env;3}\sin^{2}\theta+\cos^{2}\theta\sin^{2}\theta\Big{(}N_{\rm
env;4}\cosh^{2}r_{H}+N_{\rm env;1}\sinh^{2}r_{H}\Big{)}\bigg{]},$ (B.3)
$\expectationvalue{\hat{n}^{\rm(out)}_{3}}=\frac{1}{2}\bigg{[}-1+N_{\rm
env;2}\sin^{2}\theta+\cos^{2}\theta\Big{(}N_{\rm env;4}\cosh^{2}r_{H}+N_{\rm
env;1}\sinh^{2}r_{H}\Big{)}\bigg{]},\quad\quad\quad\quad\,\,$ (B.4)
$\expectationvalue{\hat{n}^{\rm(out)}_{4}}=\frac{1}{2}\bigg{[}-1+N_{\rm
env;4}\cosh^{4}r_{H}\sin^{4}\theta+N_{\rm
env;4}\sinh^{4}r_{H}+\cosh^{2}r_{H}\bigg{(}\cos^{2}\theta\,\Big{[}N_{\rm
env;3}+N_{\rm env;2}\sin^{2}\theta\Big{]}\\\ +\sinh^{2}r_{H}\Big{[}N_{\rm
env;1}\cos^{4}\theta-2N_{\rm env;4}\sin^{2}\theta\Big{]}\bigg{)}\bigg{]}.$
(B.5)

If we substitute $n_{\rm env;i}=0$ ($N_{\rm env;i}=1$) into this expression,
such that the initial state is vacuum, we obtain the radiation spontaneously
emitted by the system. Terms proportional to $n_{\rm env;i}$ thus correspond
to induced or stimulated radiation.

## References

  * Hawking [1975] S. W. Hawking, Particle Creation by Black Holes, Commun. Math. Phys. 43, 199 (1975), [Erratum: Commun.Math.Phys. 46, 206 (1976)].
  * Unruh, William George [1981] Unruh, William George, Experimental black-hole evaporation?, Physical Review Letters 46, 1351 (1981).
  * Frolov and Novikov [2012] V. Frolov and I. Novikov,  _Black hole physics: Basic concepts and new developments_ , Vol. 96 (Springer Science & Business Media, 2012).
  * Jacobson [2005] T. Jacobson, Introduction to quantum fields in curved spacetime and the hawking effect, in  _Lectures on Quantum Gravity_ (Springer, 2005) pp. 39–89.
  * Brady  _et al._ [2022] A. J. Brady, I. Agullo, and D. Kranas, Symplectic circuits, entanglement, and stimulated Hawking radiation in analog gravity (2022), arXiv:2209.11317 [gr-qc] .
  * Bruschi, David Edward and Friis, Nicolai and Fuentes, Ivette and Weinfurtner, Silke [2013] Bruschi, David Edward and Friis, Nicolai and Fuentes, Ivette and Weinfurtner, Silke, On the robustness of entanglement in analogue gravity systems, New Journal of Physics 15, 113016 (2013).
  * Agullo  _et al._ [2022] I. Agullo, A. J. Brady, and D. Kranas, Quantum Aspects of Stimulated Hawking Radiation in an Optical Analog White-Black Hole Pair, Physical Review Letters 128, 091301 (2022).
  * Philbin, Thomas G and Kuklewicz, Chris and Robertson, Scott and Hill, Stephen and König, Friedrich and Leonhardt, Ulf [2008] Philbin, Thomas G and Kuklewicz, Chris and Robertson, Scott and Hill, Stephen and König, Friedrich and Leonhardt, Ulf, Fiber-optical analog of the event horizon, Science 319, 1367 (2008).
  * Linder, Malte F and Schützhold, Ralf and Unruh, William G [2016] Linder, Malte F and Schützhold, Ralf and Unruh, William G, Derivation of Hawking radiation in dispersive dielectric media, Physical Review D 93, 104010 (2016).
  * Drori, Jonathan and Rosenberg, Yuval and Bermudez, David and Silberberg, Yaron and Leonhardt, Ulf [2019] Drori, Jonathan and Rosenberg, Yuval and Bermudez, David and Silberberg, Yaron and Leonhardt, Ulf, Observation of stimulated Hawking radiation in an optical analogue, Physical Review Letters 122, 010404 (2019).
  * Rosenberg, Yuval [2020] Rosenberg, Yuval, Optical analogues of black-hole horizons, Philosophical Transactions of the Royal Society A 378, 20190232 (2020).
  * Aguero-Santacruz and Bermudez [2020] R. Aguero-Santacruz and D. Bermudez, Hawking radiation in optics and beyond, Phil. Trans. Roy. Soc. Lond. A 378, 20190223 (2020), arXiv:2002.07907 [gr-qc] .
  * Petev, Mike and Westerberg, Niclas and Moss, Daniel and Rubino, Elenora and Rimoldi, C and Cacciatori, SL and Belgiorno, F and Faccio, D [2013] Petev, Mike and Westerberg, Niclas and Moss, Daniel and Rubino, Elenora and Rimoldi, C and Cacciatori, SL and Belgiorno, F and Faccio, D, Blackbody emission from light interacting with an effective moving dispersive medium, Physical Review Letters 111, 043902 (2013).
  * Rubino, Elenora and Lotti, A and Belgiorno, F and Cacciatori, SL and Couairon, Arnaud and Leonhardt, Ulf and Faccio, D [2012] Rubino, Elenora and Lotti, A and Belgiorno, F and Cacciatori, SL and Couairon, Arnaud and Leonhardt, Ulf and Faccio, D, Soliton-induced relativistic-scattering and amplification, Scientific Reports 2, 1 (2012).
  * Belgiorno, F and Cacciatori, SL and Dalla Piazza, F [2015] Belgiorno, F and Cacciatori, SL and Dalla Piazza, F, Hawking effect in dielectric media and the Hopfield model, Physical Review D 91, 124063 (2015).
  * Fox [2002] M. Fox, Optical properties of solids (2002).
  * Wald [1995] R. M. Wald,  _Quantum Field Theory in Curved Space-Time and Black Hole Thermodynamics_ , Chicago Lectures in Physics (University of Chicago Press, Chicago, IL, 1995).
  * Fabbri, Alessandro and Navarro-Salas, José [2005] Fabbri, Alessandro and Navarro-Salas, José,  _Modeling black hole evaporation_ (World Scientific, 2005).
  * Gaona-Reyes and Bermudez [2017] J. L. Gaona-Reyes and D. Bermudez, The theory of optical black hole lasers, Annals Phys. 380, 41 (2017), arXiv:1701.05655 [gr-qc] .
  * Moreno-Ruiz and Bermudez [2020] A. Moreno-Ruiz and D. Bermudez, Hawking temperature in dispersive media: Analytics and numerics, Annals Phys. 420, 168268 (2020), arXiv:1908.02368 [gr-qc] .
  * Macher and Parentani [2009] J. Macher and R. Parentani, Black/white hole radiation from dispersive theories, Physical Review D 79, 124008 (2009).
  * Busch  _et al._ [2014] X. Busch, I. Carusotto, and R. Parentani, Spectrum and entanglement of phonons in quantum fluids of light, Phys. Rev. A 89, 043819 (2014), arXiv:1311.3507 [cond-mat.quant-gas] .
  * Finazzi and Parentani [2011a] S. Finazzi and R. Parentani, Spectral properties of acoustic black hole radiation: broadening the horizon, Phys. Rev. D 83, 084010 (2011a), arXiv:1012.1556 [gr-qc] .
  * Finazzi and Parentani [2011b] S. Finazzi and R. Parentani, On the robustness of acoustic black hole spectra, J. Phys. Conf. Ser. 314, 012030 (2011b), arXiv:1102.1452 [gr-qc] .
  * Finazzi and Parentani [2012] S. Finazzi and R. Parentani, Hawking radiation in dispersive theories, the two regimes, Phys. Rev. D 85, 124027 (2012), arXiv:1202.6015 [gr-qc] .
  * Michel and Parentani [2014] F. Michel and R. Parentani, Probing the thermal character of analogue Hawking radiation for shallow water waves?, Phys. Rev. D 90, 044033 (2014), arXiv:1404.7482 [gr-qc] .
  * Michel  _et al._ [2016] F. Michel, J.-F. Coupechoux, and R. Parentani, Phonon spectrum and correlations in a transonic flow of an atomic Bose gas, Phys. Rev. D 94, 084027 (2016), arXiv:1605.09752 [cond-mat.quant-gas] .
  * Weedbrook  _et al._ [2012] C. Weedbrook, S. Pirandola, R. García-Patrón, N. J. Cerf, T. C. Ralph, J. H. Shapiro, and S. Lloyd, Gaussian quantum information, Reviews of Modern Physics 84, 621 (2012).
  * Serafini [2017] A. Serafini,  _Quantum continuous variables: a primer of theoretical methods_ (CRC press, 2017).
  * Peres [1996] A. Peres, Separability criterion for density matrices, Physical Review Letters 77, 1413 (1996).
  * Plenio [2005] M. B. Plenio, Logarithmic negativity: A full entanglement monotone that is not convex, Physical Review Letters 95, 090503 (2005).
  * Wang and Wilde [2020a] X. Wang and M. M. Wilde, Cost of quantum entanglement simplified, Phys. Rev. Lett. 125, 040502 (2020a).
  * Wang and Wilde [2020b] X. Wang and M. M. Wilde, $\alpha$-logarithmic negativity, Phys. Rev. A 102, 032416 (2020b).
  * [34] I. Agullo, A. Brady, and D. Kranas, In preparation.
  * Weinberg [2008] S. Weinberg,  _Cosmology_ (2008).
  * Unruh, William G [1976] Unruh, William G, Notes on black-hole evaporation, Physical Review D 14, 870 (1976).
  * Fredenhagen and Haag [1990] K. Fredenhagen and R. Haag, On the Derivation of Hawking Radiation Associated with the Formation of a Black Hole, Commun. Math. Phys. 127, 273 (1990).
  * Agullo  _et al._ [2010] I. Agullo, J. Navarro-Salas, G. J. Olmo, and L. Parker, Hawking Radiation by Kerr Black Holes and Conformal Symmetry, Phys. Rev. Lett. 105, 211305 (2010).
  * Lvovsky and Raymer [2009] A. I. Lvovsky and M. G. Raymer, Continuous-variable optical quantum-state tomography, Reviews of Modern Physics 81, 299 (2009).
  * Bachor and Ralph [2019] H.-A. Bachor and T. C. Ralph,  _A Guide to Experiments in Quantum Optics_ (John Wiley & Sons, 2019).