\makeatletter
\@ifundefined{HCode}
{
\documentclass[CRPHYS,Unicode,screen,Thematic]{cedram}
\newenvironment{noXML}{}{}
\newenvironment{Table}{\begin{table}}{\end{table}}
\def\thead{\noalign{\relax}\hline}
\def\endthead{\noalign{\relax}\hline}
\def\tbody{\noalign{\relax}\hline}
\def\tabnote#1{\vskip4pt\parbox{.98\linewidth}{#1}}
\def\tsup#1{$^{{#1}}$}
\def\tsub#1{$_{{#1}}$}
\def\sfrac#1#2{{#1}/{#2}}
\def\stfrac#1#2{({#1})/{#2}}
\def\stofrac#1#2{(({#1})/{#2})}
\def\ensuremath#1{$#1$}
\usepackage{amssymb}
\RequirePackage{etoolbox}
\usepackage{amsmath}
\usepackage{marvosym}
\usepackage{upgreek,bm}
\csdef{Seqnsplit}{\\}
\def\rmpi{\uppi}
\def\tbody{\noalign{\relax}\hline}
\def\xmorerows#1#2{#2}
\def\jobid{crphys20240461}
%\graphicspath{{/tmp/\jobid_figs/web/}}
\graphicspath{{./figures/}}
\def\xsection{}
\newcounter{runlevel}
\let\MakeYrStrItalic\relax
\def\refinput#1{}
\def\bnabla{\bm{\nabla}}
\def\rmmu{\upmu}
\def\ndash{\text{--}}
\def\0{\phantom{0}}
\def\botline{\\\hline}
\def\sfrac#1#2{{#1}/{#2}}
\def\stfrac#1#2{({#1})/{#2}}
\def\stofrac#1#2{(({#1})/{#2})}
\def\back#1{}
\DOI{10.5802/crphys.212}
\datereceived{2024-06-03}
\daterevised{2024-08-02}
\dateaccepted{2024-10-01}
\ItHasTeXPublished
\def\rmpi{\uppi}
\usepackage{braket}
\usepackage{leftindex}
\usepackage{amsmath}
\usepackage{amssymb}
}
{\documentclass[crphys]{article}
\usepackage{upgreek}
\def\CDRdoi{10.5802/crphys.212}
\let\newline\break
\def\xmorerows#1#2{\morerows{#1}{#2}}
%%%%%%%%%%%%
\def\sfrac#1#2{{#1}/{#2}}
\def\stfrac#1#2{({#1})/{#2}}
\def\stofrac#1#2{(({#1})/{#2})}
%%%%%%%%%%%%
\usepackage[T1]{fontenc} 
\usepackage{marvosym}
\makeatletter
\usepackage{braket}
\usepackage{leftindex}
\usepackage{amsmath}
\usepackage{amssymb}
}
\makeatother

\begin{document}

\def\braket#1{\langle#1\rangle}

\begin{noXML}

%\makeatletter
%\def\TITREspecial{\relax}
%\def\cdr@specialtitle@english{Simulating gravitational problems with condensed matter analog models: a special issue in memory of Renaud Parentani (1962--2020)}
%\def\cdr@specialtitle@french{Simuler des probl\`emes gravitationnels avec des mod\`eles analogues en mati\`ere condens\'ee : un num\'ero sp\'ecial en m\'emoire de Renaud Parentani (1962--2020)}
%\makeatother

\CDRsetmeta{articletype}{review}

\title{Resonant analogue configurations in atomic condensates}

\alttitle{Configurations analogues r\'{e}sonantes dans les condensats
atomiques}

\shortrunauthors

\author{\firstname{Juan Ram\'on} \lastname{Mu\~noz de Nova}\CDRorcid{0000-0001-6229-9640}\IsCorresp}
\address{Departamento de F\'isica de Materiales, Universidad Complutense de Madrid, 28040 Madrid, Spain}
\email[J. R. Mu\~noz de Nova]{jrmnova@fis.ucm.es}

\author{\firstname{Pablo} \lastname{Fern\'andez Palacios}\CDRorcid{0000-0002-0622-6794}}
\address{Instituto de Energ\'ia Solar, Universidad Polit\'ecnica de Madrid, 28040 Madrid, Spain}

\author{\firstname{Pedro~} \lastname{Alc\'azar~Guerrero}}
\address{Catalan Institute of Nanoscience and Nanotechnology (ICN2), CSIC and BIST, 
Campus UAB, Bellaterra, 08193 Barcelona, Spain}
\address{Department of Physics, Campus UAB, Bellaterra, 08193 Barcelona, Spain}

\author{\firstname{Ivar} \lastname{Zapata}}
\address{mrHouston Tech Solutions, 28002 Madrid, Spain}

\author{\firstname{Fernando} \lastname{Sols}\CDRorcid{0000-0002-0947-286X}}
\addressSameAs{1}{Departamento de F\'isica de Materiales, Universidad Complutense de
Madrid, 28040 Madrid, Spain}

\thanks{European Union's Horizon 2020 research and innovation programme
(Grant Agreement No. 847635), Spain's Agencia Estatal de
Investigaci\'on (Grant No. PID2022-139288NB-I00), Universidad
Complutense de Madrid (Grant No. FEI-EU-19-12)}

\begin{abstract}
As a contribution to a memorial volume, we provide a comprehensive
discussion of resonant configurations in analogue gravity, focusing on
its implementation in atomic condensates and combining review features
with original insights and calculations. In particular, we jointly
analyze the analogues of the Andreev and Hawking effects using a
microscopic description based on the Bogoliubov approximation. We
perform a detailed study of the thermality of the Andreev and Hawking
spectra for canonical black-hole solutions, finding that both can be
described by a gray-body distribution to a very good approximation. We
contemplate several resonant scenarios whose efficiency to enhance
anomalous scattering processes is compared to that of non-resonant
setups. The presence of quantum signatures in analogue configurations,
such as the violation of Cauchy--Schwarz inequalities or entanglement,
is analyzed, observing that resonant configurations highly increase the
entanglement signal, especially for the Andreev effect. We also discuss
how these results have served as inspiration for the rapidly expanding
field of quantum information in high-energy colliders. Finally, we
study the physics of black-hole lasers as further examples of resonant
analogue structures, distinguishing three stages in its time evolution.
For short times, we compute the linear and non-linear spectrum for
different models. For intermediate times, we generalize the current
analysis of the BHL--BCL crossover. For long times, we discuss the
emerging concept of spontaneous Floquet state and its potential
implications.
\end{abstract}

\begin{altabstract}
En guise de contribution \`{a} ce volume comm\'{e}moratif, nous
proposons une discussion approfondie des configurations r\'{e}sonnantes
en gravit\'{e} analogue, en nous concentrant sur leur mise en \oe{}uvre
dans les condensats atomiques et en combinant une revue de la
litt\'{e}rature avec des analyses et des calculs originaux. En
particulier, nous analysons conjointement les analogues des effets
Andreev et Hawking en utilisant une description microscopique bas\'{e}e
sur l'approximation de Bogoliubov. Nous r\'{e}alisons une \'{e}tude
d\'{e}taill\'{e}e de la thermalit\'{e} des spectres d'Andreev et de
Hawking pour les solutions canoniques de trous noirs, en constatant que
les deux peuvent \^{e}tre d\'{e}crits par une distribution de corps
gris avec une tr\`{e}s bonne approximation. Nous envisageons plusieurs
sc\'{e}narios r\'{e}sonants dont l'efficacit\'{e} pour am\'{e}liorer
les processus de diffusion anormaux est compar\'{e}e \`{a} celle de
configurations non r\'{e}sonantes. La pr\'{e}sence de signatures
quantiques dans les configurations analogues, telles que la violation
des in\'{e}galit\'{e}s de Cauchy--Schwarz ou l'intrication, est
analys\'{e}e. Nous observons que les configurations r\'{e}sonantes
augmentent fortement le signal d'intrication, en particulier pour
l'effet Andreev. Nous discutons \'{e}galement de la fa\c{c}on dont ces
r\'{e}sultats ont servi d'inspiration pour le domaine en pleine
expansion de l'information quantique dans les collisionneurs de haute
\'{e}nergie. Enfin, nous \'{e}tudions la physique des lasers \`{a}
trous noirs comme autres exemples de structures analogues
r\'{e}sonantes, en distinguant trois \'{e}tapes dans leur \'{e}volution
temporelle. Pour les temps courts, nous calculons le spectre
lin\'{e}aire et non lin\'{e}aire pour diff\'{e}rents mod\`{e}les. Pour
les temps interm\'{e}diaires, nous g\'{e}n\'{e}ralisons l'analyse
actuelle du croisement BHL--BCL. Pour les temps longs, nous discutons du
concept \'{e}mergent d'\'{e}tat de Floquet spontan\'{e} et de ses
implications potentielles.
\end{altabstract}

\keywords{\kwd{Analog gravity}\kwd{Quantum gases}\kwd{Andreev
processes}\kwd{Quantum information}\kwd{High-energy colliders}
\kwd{Time-crystals}}

\altkeywords{\kwd{Gravit\'{e} analogique}\kwd{Gaz quantiques}
\kwd{Processus d'Andreev}\kwd{Information quantique}
\kwd{Collisionneurs \`{a} haute \'{e}nergie}\kwd{Cristaux temporels}}

\maketitle

\end{noXML}

\section{Introduction}\label{sec:Intro}

In the 70s, Stephen Hawking~\cite{Hawking1974} predicted, on the basis
of a semiclassical calculation (in which fields are quantized but the
background spacetime is treated as classical), that black holes
spontaneously emit radiation with a Planckian spectrum. This became
known as Hawking radiation, one of the most celebrated predictions of
modern theoretical physics. As it involves thermodynamics, quantum
mechanics, and general relativity, understanding Hawking radiation is
regarded as the first step towards a quantum theory of gravity.
However, its detection in an astrophysical scenario is quite unlikely
because its effective temperature of emission, the Hawking temperature,
is of the order of $T_{\mathrm{H}}\sim 10^{-8}~\mathrm{K}$ for a black hole of
several solar masses, much lower for instance than the temperature of
the cosmic microwave background, $T_{\mathrm{CMB}}\simeq 2.7~\mathrm{K}$.

It was later noted by Unruh~\cite{Unruh1981} that the equations of
motion describing the fluctuations of an irrotational flow are formally
analogue to those of a massless scalar field in a curved spacetime
described by the so-called acoustic metric. This connection established
the rich field of analogue gravity, where tabletop experiments are used
to replicate gravitational effects in a controlled setup. A number of
vastly different systems have been proposed for implementing analogue
experiments, including atomic Bose--Einstein
condensates~\cite{Garay2000,Lahav2010}, water
waves~\cite{Schutzhold2002,Weinfurtner2011}, non-linear optical
fibers~\cite{Philbin2008,Drori2019}, ion
rings~\cite{Horstmann2010,Wittemer2019}, quantum fluids of
light~\cite{Carusotto2013,Nguyen2015}, and even superconducting
transmon qubits~\cite{Shi2023}. 

In the specific case of Hawking radiation, the subsonic/supersonic
regions of a flowing fluid are akin to the exterior/interior of a black
hole, since sound cannot travel upstream in a supersonic flow, in the
same way as light is trapped inside a black hole. As a result, the
analogue of the event horizon is provided by the subsonic/supersonic
interface, and Hawking radiation is associated with the spontaneous
emission of correlated phonons into the subsonic and supersonic
regions. A classical stimulated version of this emission has been
observed in hydraulic~\cite{Weinfurtner2011,Euve2016} and
optical~\cite{Drori2019} analogues. Due to their intrinsic quantumness,
well-controlled behavior, and low temperature, Bose--Einstein
condensates are the most promising candidates to exhibit the genuine
Hawking effect, of a quantum nature. Indeed, the study of quantum field
theory in curved spacetimes using atomic condensates can be regarded as
another paradigm of quantum simulation~\cite{Viermann2022}. 

From the theoretical point of view, the initial proposal of
Ref.~\cite{Garay2000} was subsequently expanded in a number of
works~\cite{Leonhardt2003a,Balbinot2008,Carusotto2008,Macher2009a,Recati2009,Zapata2011,Larre2012,deNova2014,Finazzi2014,Busch2014,deNova2015,Michel2016}. 
From the experimental point of view, the first analogue black hole in
an atomic condensate was achieved by the Technion group in
2010~\cite{Lahav2010}. After progressively improving the
setup~\cite{Shammass2012,Schley2013,Steinhauer2014}, the Technion
experiment provided the first claimed detection of the Hawking effect
in 2016~\cite{Steinhauer2016}, which was later confirmed by a precise
agreement with the theoretical
predictions~\cite{deNova2019,Kolobov2021}, including the measurement of
the Hawking temperature and the observation of the stationarity of the
spontaneous emission of Hawking radiation. The detection of the Hawking
effect is far from being the end of the road for the field, and
analogues of the dynamical Casimir effect~\cite{Jaskula2012}, Sakharov
oscillations~\cite{Hung2013}, superradiance~\cite{Torres2017},
inflation~\cite{Eckel2018}, Unruh effect~\cite{Hu2019}, quasi-normal
ringdown~\cite{Torres2020}, backreaction~\cite{Patrick2021} or
cosmological particle creation~\cite{Steinhauer2022a}, have also been
observed in the laboratory.  

Interestingly, subsonic/supersonic interfaces in condensates provide
yet another insightful conceptual connection, since it was
shown~\cite{Zapata2009a} that there also takes place the analogue of
the Andreev reflection in superconductors~\cite{Andreev1964}, where an
electron (hole) incident on a normal/superconductor interface from the
normal side is reflected as a hole (electron). The Andreev reflection
has already been observed in
superconductors~\cite{Bozhko1982,Benistant1983}, as well as in
superfluid $^3$He~\cite{Enrico1993,Okuda1998}, but there is not yet
experimental evidence of this behavior in condensates. 

During the quest for the detection of the Hawking effect, it was
proposed~\cite{Zapata2011} that setups involving multiple internal
reflections, similar to those occurring in a Fabri--Perot
interferometer, might be advantageous due to the non-thermal character
of the resulting spectrum. The precedent existed of the black-hole
laser (BHL)~\cite{Corley1999}, where multiple internal reflections on a
pair of horizons result in an increasing unstable output of Hawking
radiation. Moreover, it was later shown that resonant configurations
may enhance the quantum contribution of the Hawking
effect~\cite{deNova2014,deNova2015}.


In this work, we provide a comprehensive discussion of resonant
analogue configurations in atomic condensates and their most important
features, including an original unified analysis of the Andreev and
Hawking effects. We will devote special attention to key contributions
from Renaud Parentani, highlighted throughout the article.

We begin by providing a general introduction to gravitational analogues
in atomic condensates in Section~\ref{sec:HawkingCondensates}. In this
respect, the work by Macher and
Parentani~\cite{Macher2009a,Macher2009}, along with that of Recati, 
Pavloff and Carusotto~\cite{Recati2009}, represented a milestone,
since it established the microscopic analogy of the Hawking effect
without resorting to any effective metric. This microscopic approach is
the framework upon which we will develop our results. In particular, we
draw an interesting analogy between the Bogoliubov--de Gennes (BdG)
equations and the relativistic Klein--Gordon equation, which provides an
insightful way to address the quantization of the problem. We discuss
how the bosonic version of the Andreev reflection emerges in the
interior of an analogue black hole~\cite{Zapata2009a} and gives rise to
a spontaneous emission of quasiparticles, to which we refer as the
Andreev effect, in analogy with the Hawking effect. We assess the
thermality of canonical analogue configurations, finding that their
Andreev spectrum is also Planckian to a very good approximation due to
its universal scaling at low frequencies~\cite{deNova2014}, an
unnoticed result in the\break literature.

Section~\ref{sec:ResonantHawking} is devoted to the study of
Andreev--Hawking processes in resonant analogue structures. Resonant
Andreev scenarios have been also considered in normal/superconductor
junctions~\cite{Prada2005}. Here, we examine a double barrier structure
with alternating supersonic and subsonic segments in the region between
the two barriers, in whose original proposal Parentani directly
contributed~\cite{Zapata2011}. Another setup analyzed in this section
is a resonant flat-profile structure, involving a stationary
homogeneous flowing condensate with a piecewise modulated speed of
sound, first discussed in Ref.~\cite{deNova2014}.  Our results show
that the Andreev signal is also highly enhanced by resonant peaks in
the spectrum, and it may even become larger than the Hawking signal, in
contrast with typical analogue configurations.

It was noticed that an optical lattice, with its characteristic
multiple barrier structure, might enhance some of the properties of a
resonant structure~\cite{deNova2014a,deNova2017b}. The second part of
Section~\ref{sec:ResonantHawking} addresses the long-lasting
quasi-stationary outcoupling of a condensate through an optical
lattice, which may have an ideal uniform or a Gaussian shape, with a
characteristic time scale much longer than that of conventional
black-hole configurations.

Section~\ref{sec:QuantumHR} analyzes the  quantumness of the Andreev
and Hawking effects by borrowing concepts and techniques from quantum
optics~\cite{Walls2008}. This approach has also been recently advocated
in the original gravitational scenario~\cite{Scully2022}. With the help
of these tools, we jointly characterize the Andreev and Hawking effects
as the spontaneous production of hybrid Andreev--Hawking modes from the
non-degenerate parametric amplification of the vacuum. The quantumness
of the Andreev--Hawking effect is evaluated through different quantum
correlations such as the violation of Cauchy--Schwarz (CS) inequalities
or entanglement. Originally, the violation of CS inequalities by
quantum Hawking radiation was first discussed in
Ref.~\cite{deNova2014}. In an almost parallel effort, Busch and
Parentani~\cite{Busch2014}, as well as Finazzi and
Carusotto~\cite{Finazzi2014}, analyzed the entanglement of Hawking
radiation in condensates using the generalized Peres--Horodecki
criterion~\cite{Simon2000}. The inspiring visit of Renaud Parentani to
our group in Madrid in November 2013 motivated us to unify both
approaches within a common framework in Ref.~\cite{deNova2015}. All
these techniques were later employed in the detection of the
entanglement of Hawking radiation in the 2016\break
experiment~\cite{Steinhauer2016}. 

In the last part of Section~\ref{sec:QuantumHR}, we briefly discuss how
the study of quantum correlations in the Andreev--Hawking effect has
motivated the research on quantum information in high-energy
colliders~\cite{Afik2021}. This has rapidly become an emergent field of
research by itself, which has already achieved its first milestone with
the pioneering observation of entanglement in quarks by the ATLAS and
CMS collaborations at the Large Hadron Collider
(LHC)~\cite{ATLAS2024,CMS2024}, representing
also the highest-energy detection of entanglement ever. A pedagogical
introduction to this fascinating topic for a readership outside the
high-energy field is presented in Ref.~\cite{Afik2022}. 

Section~\ref{sec:BHL} addresses the emergence of a black-hole laser in
resonant configurations. In an atomic condensate, the BHL effect arises
because of its superluminal dispersion relation, which allows the
radiation reflected at the inner horizon to travel back to the outer
one, further stimulating the production of Andreev--Hawking
radiation~\cite{Leonhardt2003,Barcelo2006,Jain2007,Coutant2010,Finazzi2010,Bermudez2018,Burkle2018}. 
Other analogue setups have been proposed to observe the BHL
effect~\cite{Faccio_2012,Peloquin2016,RinconEstrada2021,Katayama2021}.
The work by Parentani and collaborators~\cite{Coutant2010,Finazzi2010}
was instrumental in determining the properties of a BHL, including the
first full microscopic BdG computation of the spectrum of dynamical
instabilities, in analogy with the microscopic derivation of the
Hawking effect~\cite{Recati2009,Macher2009a}. Parentani and Michel also
pioneered the study of the non-linear regime of a
BHL~\cite{Michel2013}, achieved once the initial instability has grown
up to saturation, and the numerical study of its
dynamics~\cite{Michel2015}, in parallel to the work by Mu\~noz\break 
de Nova, Finazzi and Carusotto~\cite{deNova2016}.

Our discussion of the black-hole laser further extends the original
work of Parentani in all the stages of its time evolution. At short
times, using the protocol to construct BHL solutions of
Ref.~\cite{deNova2017a}, we compute the linear and non-linear spectrum
for different BHL models, including that of a double barrier structure,
original of the present work. Our results confirm all the trends
anticipated in Ref.~\cite{Michel2013}.

At intermediate times in the evolution of a BHL, one has to take into
account the Bogoliubov--Cherenkov--Landau (BCL) mode present in the
supersonic region~\cite{Carusotto2006}, which is analogous to the
undulation in hydraulic setups~\cite{Coutant2012}. Because of its
zero-frequency, the BCL mode is resonantly excited by any obstacle in
the flow and overshadows the BHL effect in real
experiments~\cite{Steinhauer2014,Kolobov2021,Steinhauer2022}. The
BHL--BCL problem has attracted a number of studies in the theoretical
literature~\cite{Kolobov2021,Steinhauer2022,Tettamanti2016,Steinhauer2017,Wang2016,Wang2017,Llorente2019,Tettamanti2021,deNova2023}, 
and the observation of the BHL effect still remains a major challenge
in the analogue field. Here, we generalize the discussion in
Ref.~\cite{deNova2023} of the BHL--BCL crossover, originally based on a
flat-profile model, underlining the crucial role played by the
$\mathbb{Z}_2$ symmetry of a quantum BHL, first predicted by Michel and
Parentani~\cite{Michel2013}.

For sufficiently long times, the BHL displays a dynamical phase diagram
where it can only reach two states~\cite{deNova2016,deNova2021}: the
true non-linear ground state or the so-called Continuous Emission of
Solitons (CES) state, which represents a realization of a spontaneous
Floquet state~\cite{deNova2022}. The original conception of spontaneous
Floquet state was heavily influenced by richful discussions with Renaud
Parentani during the visit of one of us (JRMdN) to Orsay in 2015. Here
we analyze in detail the CES state arising from a flat-profile BHL
solution, and discuss interesting implications of spontaneous Floquet
states, including a specific and tangible realization of time operator
in quantum mechanics.

The inspiration of Renaud Parentani, and in some cases his direct
involvement, is a common thread of the work discussed in this article.

\section{Andreev and Hawking effects in atomic condensates}\label{sec:HawkingCondensates}

\subsection{Gross--Pitaevskii and Bogoliubov--de Gennes equations}\label{subsec:GPBdG}

We begin by reviewing the basic concepts and techniques for the study
of atomic condensates. We consider the following general
second-quantization Hamiltonian for interacting
bosons~\cite{Fetter2003}:
{\begin{equation}\label{eq:HamiltonianManyBody}
\hat{H}=\int\mathrm{d}\mathbf{x}~\hat{\Psi}^{\dagger}(\mathbf{x})\left[-\frac{\hbar^2}{2m}\nabla^2+V(\mathbf{x},t)+\frac{g}{2}\hat{\Psi}^{\dagger}(\mathbf{x})\hat{\Psi}(\mathbf{x})\right]\hat{\Psi}(\mathbf{x}),
\end{equation}}\unskip
where $m$ is the mass of the atoms, $V(\mathbf{x},t)$ is some external
potential, and the bosons interact through the contact pseudopotential
$W(\mathbf{x}-\mathbf{x}')=g\delta(\mathbf{x}-\mathbf{x}')$~\cite{Pitaevskii2003}. 
The field operator $\hat{\Psi}(\mathbf{x})$ satisfies the canonical
commutation relation $[\Psi(\mathbf{x}),\Psi^\dagger
(\mathbf{x}')]=\delta(\mathbf{x}-\mathbf{x}')$, which leads to the
Heisenberg equation of motion
{\begin{equation}\label{eq:HeisenbergEquationOfMotion}
    \mathrm{i}\hbar\partial_t\hat{\Psi}(\mathbf{x},t)=\left[-\frac{\hbar^2}{2m}\nabla^2+V(\mathbf{x},t)+g\hat{\Psi}^{\dagger}(\mathbf{x},t)\hat{\Psi}(\mathbf{x},t)\right]\hat{\Psi}(\mathbf{x},t).
\end{equation}}\unskip
Close to $T=0$, the condensate can be described by a coherent state,
characterized by a macroscopic wavefunction $\Psi(\mathbf{x},t)$ that
is normalized to the total particle number,
{\begin{equation}\label{eq:Normalization}
\int\mathrm{d}\mathbf{x}~{|\Psi(\mathbf{x},t)|}^2=N.
\end{equation}}\unskip
Quantum fluctuations around the condensate are accounted by expanding
the field operator around its coherent expectation value (see
Section~\ref{subsec:QuantumOptics} for a thorough discussion on
coherent states) as 
{\begin{equation}
   \hat{\Psi}(\mathbf{x},t)=\Psi(\mathbf{x},t)+\hat{\varphi}(\mathbf{x},t).
\end{equation}}\unskip
Plugging this expansion into
Equation~(\ref{eq:HeisenbergEquationOfMotion}) yields,  at leading
order, the \textit{time-dependent} Gross--Pitaevskii (GP) equation,
{\begin{equation}\label{eq:TDGP}
    \mathrm{i}\hbar\partial_t \Psi(\mathbf{x},t)=\left[-\frac{\hbar^2}{2m}\nabla^2+V(\mathbf{x},t)+g{|\Psi(\mathbf{x},t)|}^2\right]\Psi(\mathbf{x},t),
\end{equation}}\unskip
and, at linear order in the quantum fluctuations, the
\textit{time-dependent} Bogoliubov--de Gennes (BdG) equations
{\begin{equation}\label{eq:TDBdG} \mathrm{i}\hbar\partial_t\hat{\Phi}=M(t)\hat{\Phi},\quad \hat{\Phi}=
\left[\begin{array}{@{}c@{}}\hat{\varphi}\\ \hat{\varphi}^\dagger\end{array}\right],\quad M(t)=\left[\begin{array}{@{}cc@{}} N(t) & A(t)\\
-A^*(t) &-N(t)\end{array}\right],
\end{equation}}\unskip
where
{\begin{equation}\label{eq:TDBdGMatrix}
    N(t)=-\dfrac{\hbar^2}{2m}\nabla^2+V(\mathbf{x},t)+2g{|\Psi(\mathbf{x},t)|}^2,\quad A(t)=\Psi^2(\mathbf{x},t).
\end{equation}}\unskip
The GP equation is thus a non-linear Schr\"odinger equation that
describes the condensate dynamics, where the non-linearity stems from
the interactions between the condensate atoms, while the linear
dynamics of the quantum fluctuations is governed by the BdG equations.
Notice that the BdG equations also describe the linear dynamics of  the
fluctuations of the GP wavefunction $\Psi'(\mathbf{x},t)$ around a
reference solution $\Psi(\mathbf{x},t)$,
$\Psi'(\mathbf{x},t)=\Psi(\mathbf{x},t)+\varphi(\mathbf{x},t)$,
resulting in the substitution $\hat{\varphi}\to \varphi$ in
Equation~(\ref{eq:TDBdG}).

We focus on time-independent configurations,
$V(\mathbf{x},t)=V(\mathbf{x})$, and consider stationary condensates,
which are accounted by 
{\begin{equation}\label{eq:StationaryFieldExpansion}
   \hat{\Psi}(\mathbf{x},t)=\left[\Psi_0(\mathbf{x})+\hat{\varphi}(\mathbf{x},t)\right]\mathrm{e}^{-\mathrm{i}\mu t/\hbar},
\end{equation}}\unskip
$\mu$ being the chemical potential. This results in the
\textit{time-independent} GP equation 
{\begin{equation}\label{eq:TIGP}
    \mu \Psi_0(\mathbf{x})=\left[-\frac{\hbar^2}{2m}\nabla^2+V(\mathbf{x})+g{|\Psi_0(\mathbf{x})|}^2\right]\Psi_0(\mathbf{x}),
\end{equation}}\unskip
and the \textit{stationary} BdG equations
{\begin{equation}\label{eq:TIBdG} \mathrm{i}\hbar\partial_t\hat{\Phi}=M_0\hat{\Phi},\quad\hat{\Phi}=\left[\begin{array}{@{}c@{}}\hat{\varphi}\\ \hat{\varphi}^\dagger\end{array}\right],\quad M_0=\left[\begin{array}{@{}cc@{}} N_0 & A_0\\
-A_0^* &-N_0\end{array}\right],
\end{equation}}\unskip
where now
{\begin{equation}\label{eq:TIBdGMatrix}
    N_0=-\dfrac{\hbar^2}{2m}\nabla^2+V(\mathbf{x})+2g{|\Psi_0(\mathbf{x})|}^2-\mu,\quad A_0=\Psi^2_0(\mathbf{x}).
\end{equation}}\unskip
Since $M_0$ is time-independent, any solution to the stationary BdG
equations can be expanded in terms of a complete set of eigenmodes
{\begin{equation}\label{eq:BdGEigenmode}
    M_0z_n=\epsilon_n z_n,\quad z_n\equiv \left[\begin{array}{@{}c@{}}u_n\\ v_n\end{array}\right].
\end{equation}}\unskip
The matrix operator $M_0$ is non-Hermitian and it can possess complex
eigenvalues, representing dynamical instabilities which grow
exponentially in time. For the present moment, we assume that the
system is dynamically stable and ignore the presence of Nambu--Goldstone
modes; we will come back later to this issue in Section~\ref{sec:BHL}. 

Even though $M_0$ is non-Hermitian, the BdG eigenmodes do form an
orthonormal basis under the inner product
{\begin{equation}\label{eq:BdGProduct}
    (z_n|z_m)\equiv\braket{z_n|\sigma_z| z_m} =\int\mathrm{d}\mathbf{x}~[u^*_nu_m-v^*_nv_m], 
\end{equation}}\unskip
with $\braket{z_n|z_m}$ the standard scalar product for two spinors and
$\sigma_i$ the usual Pauli matrices. This is because $\Lambda\equiv
\sigma_z M_0$ is indeed Hermitian, and thus $M_0$ is pseudo-Hermitian,
i.e.,
{\begin{equation}\label{eq:PseudoHermitian}
    (z_n|M_0z_m)=\braket{z_n|\Lambda z_m}=\braket{\Lambda z_n|z_m}=(M_0z_n|z_m),
\end{equation}}\unskip
which implies the conservation of the inner product between solutions
of the BdG equations and the orthogonality between eigenmodes,
{\begin{equation}\label{eq:EigenOrto}
    (\epsilon_m-\epsilon^*_n)(z_n|z_m)=0.
\end{equation}}\unskip
Actually, the conservation of the norm for any solution $z$ of the BdG
equations, $\mathrm{i}\hbar\partial_t z=M_0z$, can be derived from a continuity
equation, in analogy with the Schr\"odinger equation (see
Equation~(\ref{eq:SchrodingerCurrent})),
{\begin{eqnarray}\label{eq:QuasiparticleCurrent}
\begin{array}{c}
\displaystyle    \partial_t(z^\dagger\sigma_z z)+\boldsymbol{\nabla}\cdot \mathbf{j}=0,\vspace*{2.5pt}\\
\displaystyle  \mathbf{j}=-\frac{\mathrm{i}\hbar}{2m}\left[u^*\boldsymbol{\nabla} u-u\boldsymbol{\nabla}u^*+v^*\boldsymbol{\nabla} v-v\boldsymbol{\nabla}v^*\right],
\end{array}
\end{eqnarray}}\unskip
where $u,v$ are the components of the spinor $z$ and $\mathbf{j}$ is
the quasiparticle current. 

However, in contrast to the Schr\"odinger case, both $M_0$ and the
inner product are not positive definite. Indeed, by noticing that
$\sigma_x M^*_0 \sigma_x=-M_0$ and $\sigma_x \sigma_z
\sigma_x=-\sigma_z$, we can define a conjugate mode as $\bar{z}_n\equiv
\sigma_x z^*_n$, which has opposite eigenvalue and norm  
{\begin{equation}\label{eq:ConjugateRelations}
    M_0\bar{z}_n=-\epsilon^*_n \bar{z}_n,\quad (z_n|z_m)=-(\bar{z}_n|\bar{z}_m)^*.
\end{equation}}\unskip
This symmetry stems from that of the field spinor $\hat{\Phi}$, which
is self-conjugate, $\hat{\bar{\Phi}}=\hat{\Phi}$. Unless otherwise
stated, the modes $z_n$ are chosen with positive norm, $(z_n|z_n)=1$.

The above properties of the inner product suggest that the correct
analogy for the BdG equations should rather be established with the
Klein--Gordon (KG) equation for an Hermitian scalar field,
{\begin{equation}\label{eq:KleinGordon}
    \left[\square-\frac{m^2c^2}{\hbar^2}\right]\hat{\phi}=0,\quad \square\equiv \partial_\mu \partial^\mu=\nabla^2-\frac{1}{c^2}\partial^2_t, 
\end{equation}}\unskip
where we take the Minkowski metric as
$\eta_{\mu\nu}=\mathrm{diag}[-1,1,1,1]$. By invoking its canonical
momentum $\hat{\Pi}(\mathbf{x})$, which satisfies the equation of
motion $\hat{\Pi}=\hbar \partial_t\hat{\phi}$ and the commutation
relation
$[\hat{\phi}(\mathbf{x}),\hat{\Pi}(\mathbf{x}')]=\mathrm{i}\delta(\mathbf{x}-\mathbf{x}')$, 
the KG equation can be recasted as the BdG equations (\ref{eq:TIBdG}),
{\begin{equation}\label{eq:TIKG} \mathrm{i}\hbar\partial_t\hat{\Phi}=M_{0}\hat{\Phi},\quad\hat{\Phi}=\left[\begin{array}{@{}c@{}}\hat{\phi}\\ \mathrm{i}\hat{\Pi}\end{array}\right],\quad M_{0}=\left[\begin{array}{@{}cc@{}} 0 & 1\\
H_0 &0 \end{array}\right], 
\end{equation}}\unskip
with $H_0\equiv -(\hbar c \nabla)^2+m^2c^4$. The KG modes are derived
from the eigenvalue problem 
{\begin{equation}
    M_0z_n=\epsilon_n z_n,\quad z_n\equiv \left[\begin{array}{@{}c@{}}\phi_n\\ \chi_n\end{array}\right],
\end{equation}}\unskip
equivalent to the more usual equation $\epsilon^2_n\phi_n=H_0\phi_n$.
Notice that $M_0$ is again non-Hermitian, while $\Lambda=\sigma_x M_0$
is, which implies the conservation of the KG inner product 
{\begin{equation}\label{eq:KGProduct}
    (z_n|z_m)\equiv\braket{z_n|\sigma_x| z_m} =\int\mathrm{d}\mathbf{x}~[\phi^*_n\chi_m+\chi^*_n\phi_m].
\end{equation}}\unskip
Conjugate solutions are defined now as $\bar{z}_n\equiv \sigma_z
z^*_n$, since $\sigma_z M^*_0 \sigma_z=-M_0$ and $\sigma_z \sigma_x
\sigma_z=-\sigma_x$, so Equation~(\ref{eq:ConjugateRelations}) is
satisfied. Moreover, the field spinor is also self-conjugate; this
property can be directly traced back here to the Hermitian character of
the field $\hat{\phi}$. 

The field spinor can be expanded in terms of the complete set of
eigenmodes $\{z_n,\bar{z}_n\}$. In both BdG and Hermitian KG cases, the
self-conjugate character of $\hat{\Phi}$ implies that this expansion is
of the form
{\begin{equation}\label{eq:QuantumFieldFluctuations}
   \hat{\Phi}(\mathbf{x},t)=\sum_{n}\hat{a}_{n}(t)z_{n}(\mathbf{x})+\hat{a}^{\dagger}_{n}(t)\bar{z}_{n}(\mathbf{x}),
\end{equation}}\unskip
where $\hat{a}_n(t)\equiv (z_n|\hat{\Phi}(t))$ is the quantum amplitude
of the mode $z_n$. The canonical commutation rules for the field spinor
can be expressed in matrix form as 
{\begin{equation}
    [\hat{\Phi}(\mathbf{x}),\hat{\Phi}^\dagger(\mathbf{x}')]=\sigma_i\delta(\mathbf{x}-\mathbf{x}'),
\end{equation}}\unskip
with $\sigma_i$ the Pauli matrix characterizing the corresponding inner
product, Equations (\ref{eq:BdGProduct}), (\ref{eq:KGProduct}). Using
this relation, it is straightforward to prove that the quantum
amplitudes $\hat{a}_{n}$ behave as bosonic annhihilation operators,
{\begin{equation}\label{eq:Aniquilacion}
    [\hat{a}_n,\hat{a}^\dagger_m]=[(z_n|\hat{\Phi}),(\hat{\Phi}|z_m)]=(z_n|z_m)=\delta_{nm},
\end{equation}}\unskip
whose equation of motion is simply
{\begin{equation}\label{eq:BogoliubovQuantum}
    \mathrm{i}\hbar\partial_t \hat{a}_n=(z_n|M_0\hat{\Phi})=\hbar\omega_n\hat{a}_n\Longrightarrow \hat{a}_n(t)=\hat{a}_n\mathrm{e}^{-\mathrm{i}\omega_nt},
\end{equation}}\unskip
$\hbar \omega_n=\epsilon_n$ being the frequency of the mode. Hence, we
arrive at the usual result
{\begin{equation}\label{eq:QuantumFieldFluctuationsExplicit}
   \hat{\Phi}(\mathbf{x},t)=\sum_{n}\hat{a}_{n}z_{n}(\mathbf{x})\mathrm{e}^{-\mathrm{i}\omega_n t}+\hat{a}^{\dagger}_{n} \bar{z}_{n}(\mathbf{x})\mathrm{e}^{\mathrm{i}\omega_n t}.
\end{equation}}\unskip
Remarkably, this expansion allows to diagonalize the KG Hamiltonian in
an elegant and compact way:
{\begin{eqnarray}
    \hat{H}_{\mathrm{KG}}=\frac{1}{2}\int\mathrm{d}\mathbf{x}~\hat{\Pi}^2+(\hbar c)^2{|\boldsymbol{\nabla}\hat{\phi}|}^2+m^2c^4\hat{\phi}^2    
 =\frac{1}{2}(\hat{\Phi}|M_0\hat{\Phi})=\sum_{n} \epsilon_n\left(\hat{a}^{\dagger}_{n}\hat{a}_{n}+\frac{1}{2}\right).
\end{eqnarray}}\unskip

In the BdG case, the field expansion diagonalizes the grand-canonical
Hamiltonian\break $\hat{K}=\hat{H}-\mu \hat{N}$, with
{\begin{equation}\label{eq:ParticleNumberOperator}
\hat{N}=\int\mathrm{d}\mathbf{x}~\hat{\Psi}^{\dagger}(\mathbf{x})\hat{\Psi}(\mathbf{x})
\end{equation}}\unskip
the particle-number operator. After expanding up to quadratic order, in
the spirit of the Bogoliubov approximation, one obtains:
{\begin{equation}\label{eq:GrandCanonicalEnergy}
    \hat{K}\simeq K_0+K'_V+\tfrac{1}{2}(\hat{\Phi}|M_0\hat{\Phi})=K_0+K_V+\sum_{n} \epsilon_n \hat{a}^{\dagger}_{n}\hat{a}_{n},
\end{equation}}\unskip
where $K_0\equiv K[\Psi_0]$ is the mean-field energy of the condensate,
obtained by replacing $\hat{\Psi}$ by $\Psi_0$, and
{\begin{eqnarray}
 K'_V&=&\frac{1}{2}\int\mathrm{d}\mathbf{x}~[\hat{\varphi}^\dagger,N_0\hat{\varphi}
 +A_0\hat{\varphi}^\dagger]=-\frac{1}{2}\sum_{n}\epsilon_n\braket{z_n|z_n}, \nonumber\\
     K_V&=&K'_V+\frac{1}{2}\sum_{n}\epsilon_n=-\sum_{n}\int\mathrm{d}\mathbf{x}~\epsilon_n{|v_n|}^2 
\end{eqnarray}}\unskip
are $c$-number contributions arising from the zero-point motion of the
quasiparticles. Notice that the grand-canonical Hamiltonian $\hat{K}$
is the one governing the dynamics instead of $\hat{H}$, since we have
extracted the global phase $\mathrm{e}^{-\mathrm{i}\mu t/\hbar}$ in
Equation~(\ref{eq:StationaryFieldExpansion}). In fact, the
time-independent GP equation~(\ref{eq:TIGP}) is precisely the condition
for $\Psi_0$ to be an extreme of the grand-canonical energy $K$, which
leads to the identification of the non-linear GP eigenvalue as the
chemical potential $\mu$, and to the absence of linear terms in the
field fluctuations in Equation~(\ref{eq:GrandCanonicalEnergy}). The
precise nature of the extreme is obtained by considering small
fluctuations of the stationary GP wavefunction: 
{\begin{eqnarray}\label{eq:energeticstability}
\delta K&\equiv& K[\Psi'_0]-K[\Psi_0]\simeq \tfrac{1}{2}(\Phi|M_0\Phi)=\tfrac{1}{2}\braket{\Phi|\Lambda|\Phi},\nonumber  \\
\Psi'_0(\mathbf{x})&=&\Psi_0(\mathbf{x})+\varphi(\mathbf{x}),\quad \Phi=\left[\begin{array}{@{}c@{}}\varphi\\ \varphi^*\end{array}\right].
\end{eqnarray}}\unskip
If $\Lambda$ is a positive-definite operator, then $\Psi_0$ is a
minimum and the system is energetically stable. In that case,
{\begin{equation}\label{eq:energeticdynamical}
    \braket{z_n|\Lambda|z_n}=(z_n|M_0z_n)=\epsilon_n(z_n|z_n)>0,
\end{equation}}\unskip
so all energies are positive, $\epsilon_n>0$, and the ground state is
the quasiparticle vacuum  $\hat{a}_{n}\ket{0}=0$. If $\Lambda$ is not
positive definite, we can have negative-energy modes, denoted as
anomalous, while positive-energy modes are denoted as normal. As a
result, the system is energetically unstable.

In general, the quantum state $\hat{\rho}$ of an ensemble of bosons at
thermal equilibrium at a temperature $T$ is 
{\begin{equation}
    \hat{\rho}=\frac{\mathrm{e}^{-\beta \hat{K}}}{Z},
\end{equation}}\unskip
with $\beta=1/k_B T$ and $Z=\mathrm{Tr}(\mathrm{e}^{-\beta \hat{K}})$ the
partition function. Within the Bogoliubov approximation, this leads in
an energetically stable condensate to a thermal Planckian distribution
for the quasiparticle occupation number:
{\begin{equation}\label{eq:PlanckianQuasiparticle}
\braket{\hat{a}^{\dagger}_{n}\hat{a}_{m}}=\mathrm{Tr}(\hat{a}^{\dagger}_{n}\hat{a}_{m}\hat{\rho})=\frac{\delta_{nm}}{\mathrm{e}^{\beta\epsilon_n}-1}.
\end{equation}}\unskip

\subsection{Gravitational analogy}

We now review how the original gravitational analogy~\cite{Unruh1981}
was established using a fluid flow. Specifically, we consider the Euler
equations for an ideal irrotational barotropic flow,
{\begin{eqnarray}\label{eq:EulerVelocity}
    0&=&\partial_t \rho+\boldsymbol{\nabla}\cdot\mathbf{J},\quad \mathbf{J}= \rho \mathbf{v},  \nonumber\\
  {}[\partial_t +\mathbf{v}\cdot\boldsymbol{\nabla}]\mathbf{v} &\equiv& D_{\mathrm{t}}\mathbf{v}=-\frac{1}{m}\boldsymbol{\nabla}V-\frac{1}{\rho}\boldsymbol{\nabla}P
\end{eqnarray}}\unskip
where $\rho$ is the mass density of the fluid, $\mathbf{J}$ is the
current, $V$ is some external potential (e.g., gravity), $D_{\mathrm{t}}$ is the
total derivative, and $P$ is the local pressure. The irrotationality
condition $\boldsymbol{\nabla}\times\mathbf{v}=0$ implies that the flow
is potential, i.e., $\mathbf{v}=\boldsymbol{\nabla}\phi$. By invoking
the barotropic condition, $P=P(\rho)$, we can arrive at a simplified
equation for the flow potential $\phi$,
{\begin{eqnarray}\label{eq:EulerPotential}
    0&=&\partial_t \rho+\boldsymbol{\nabla}\cdot(\rho \boldsymbol{\nabla}\phi), \nonumber\\
     \partial_t\phi &=& -\frac{{|\boldsymbol{\nabla}\phi|}^2}{2}-\frac{V}{m}-h(\rho),\quad \frac{\mathrm{d}h}{\mathrm{d}\rho}=\frac{1}{\rho}\frac{\mathrm{d}P}{\mathrm{d}\rho}.
\end{eqnarray}}\unskip
In the usual case of an isentropic flow, $h(\rho)$ is the specific
enthalpy. The gravitational analogy emerges when considering small
fluctuations $\delta\rho,\,\delta \phi$ of the density and the flow
potential around a certain background solution characterized by
$\rho,\,\phi$. Specifically, after expanding up to linear order, we get
{\begin{eqnarray}\label{eq:EulerPotentialLinear}
    D_{\mathrm{t}}\frac{\delta\rho}{\rho} &=&-\frac{1}{\rho}\boldsymbol{\nabla}\cdot(\rho\boldsymbol{\nabla}\delta\phi), \nonumber\\
     D_{\mathrm{t}}\delta \phi &=& -c^2\frac{\delta \rho}{\rho},\quad c^2\equiv \frac{\mathrm{d}P}{\mathrm{d}\rho},
\end{eqnarray}}\unskip
$c$ being the local speed of sound. By combining both equations, one
arrives at a single equation for the flow potential fluctuations that
can be rewritten as
{\begin{equation}
    \square \delta \phi=\nabla_\mu \nabla^\mu  \delta \phi=\frac{1}{\sqrt{-g}}\partial_\mu(\sqrt{-g}g^{\mu\nu}\partial_\nu\delta \phi)=0,
\end{equation}}\unskip
which is precisely the covariant form of the KG
equation~(\ref{eq:KleinGordon}) for a massless scalar field in a curved
spacetime described by the metric
{\begin{equation}\label{eq:RelativisticMetric}
g_{\mu\nu}(x)=\frac{\rho(x)}{c(x)}\left[\begin{array}{@{}cc@{}}-[c^2(x)-v^2(x)]& -\mathbf{v}^T(x) \\
-\mathbf{v}(x)& \delta_{ij}\\
\end{array}\right],\quad x\equiv(t,\mathbf{x}),
\end{equation}}\unskip
whose line element simply reads
{\begin{equation}\label{eq:LineElement}
\mathrm{d}s^2=\frac{\rho(x)}{c(x)}[-c^2(x)\mathrm{d}t^2+{|\mathrm{d}\mathbf{x}-\mathbf{v}(x)\mathrm{d}t|}^2].
\end{equation}}\unskip
The metric $g_{\mu\nu}$ is known as the hydrodynamic (or acoustic)
metric, and parametrizes a whole class of metrics. Thus, we can use
fluid flows, accessible to us in the laboratory, to study gravitational
phenomena that can be mimicked by an acoustic metric. Specifically, we
can address the physics of black holes since the acoustic metric
presents horizons at the subsonic/supersonic interfaces, where
$v(x)=c(x)$, denoted as acoustic horizons. In fact, the Schwarzschild
metric can be rewritten using the Gullstrand--Painlev\'e coordinates as
a stationary acoustic metric with
{\begin{equation}
    c(\mathbf{x})=c,\quad \mathbf{v}(\mathbf{x})=c \sqrt{\frac{r_{\mathrm{S}}}{r}}\frac{\mathbf{x}}{r},\quad r_{\mathrm{S}}=\frac{2GM}{c^2},
\end{equation}}\unskip
where $r_{\mathrm{S}}$ is the Schwarzschild radius. Thus, the exterior/interior of
a black hole is analogous to the subsonic/supersonic regions of a
flowing fluid. A more thorough discussion about ergoregions and
horizons in acoustic metrics is presented in
Ref.~\cite{Visser1998}.

In condensates, the gravitational analogy emerges within the Bogoliubov
formalism in the so-called hydrodynamic approximation. For that
purpose, we invoke the Madelung decomposition of the GP wavefunction,
{\begin{equation}
\Psi(\mathbf{x},t)=\sqrt{n(\mathbf{x},t)}\mathrm{e}^{\mathrm{i}\theta(\mathbf{x},t)},
\end{equation}}\unskip
which leads to a pair of hydrodynamic equations after rewriting the
time-dependent GP equation~(\ref{eq:TDGP}) in terms of the condensate
phase and density:
{\begin{eqnarray}\label{eq:EulerPhaseCondensate}
    0&=&\partial_t n+\boldsymbol{\nabla}\cdot(n \mathbf{v}),\quad \mathbf{v}=\frac{\hbar\boldsymbol{\nabla}\theta}{m},  \nonumber\\
  \frac{\hbar \partial_t \theta}{m} &=&\frac{\hbar^2}{2m^2\sqrt{n}}\nabla^2\sqrt{n}-\frac{v^2}{2}-\frac{V}{m}-\frac{gn}{m}.
\end{eqnarray}}\unskip
The first line is a continuity equation, from where we identify the
particle current  
{\begin{equation}\label{eq:SchrodingerCurrent}
    \mathbf{J}=n\mathbf{v}=-\frac{\mathrm{i}\hbar}{2m}\left[\Psi^*\boldsymbol{\nabla} \Psi-\Psi\boldsymbol{\nabla}\Psi^*\right]
\end{equation}}\unskip
and the flow velocity $\mathbf{v}$; notice that $n={|\Psi|}^2$ is the
particle density of the condensate, related to the mass density as
$\rho=n\cdot m$. Since $\mathbf{v}$ is proportional to the gradient of
the phase, the resulting flow is irrotational, with a flow potential
$\phi=\hbar \theta/m$. The second line provides the dynamics for the
potential flow, from where we can identify the local pressure as
{\begin{equation}
    h=\frac{gn}{m}\Longrightarrow P=\frac{gn^2}{2}.
\end{equation}}\unskip
The only genuine quantum-mechanical term involving $\hbar$ in
Equation~(\ref{eq:EulerPhaseCondensate}) is the so-called quantum
potential
{\begin{equation}
    Q\equiv -\frac{\hbar^2}{2m\sqrt{n}}\nabla^2\sqrt{n}.
\end{equation}}\unskip
Hence, in the hydrodynamic approximation, where $Q$ is negligible, the
GP equation reduces to the Euler equation for an ideal irrotational
barotropic flow, from where the gravitational analogy is retrieved, as
originally shown in Ref.~\cite{Garay2000}. 

Further insight on the quantum aspects of the gravitational analogy can
be obtained from the time-dependent BdG equations (\ref{eq:TDBdG}). By
using the relative quantum fluctuations,
$\hat{\varphi}(\mathbf{x},t)\equiv\Psi
(\mathbf{x},t)\hat{\chi}(\mathbf{x},t)$, we arrive at 
{\begin{equation}\label{eq:BdGfieldequationTIHydrodynamic}
\mathrm{i}\hbar D_{\mathrm{t}}\hat{\chi}=\left[T_n+mc^2\right]\hat{\chi}+mc^2\hat{\chi}^{\dagger},\quad T_n\equiv -\frac{\hbar^2}{2mn }\boldsymbol{\nabla}\cdot(n\boldsymbol{\nabla}),
\end{equation}}\unskip
where the speed of sound is simply found to be $c^2=h=gn/m$. The
relative quantum fluctuations can be in turn expressed in terms of the
more physical density and phase fluctuations\break from 
{\begin{eqnarray}
 \hat{\Psi}(\mathbf{x},t)=\Psi(\mathbf{x},t)+\hat{\varphi}(\mathbf{x},t)=\Psi(\mathbf{x},t)[1+\hat{\chi}(\mathbf{x},t)]   
   =\sqrt{n(\mathbf{x},t)+\delta\hat{n}(\mathbf{x},t)}\mathrm{e}^{\mathrm{i}[\theta(\mathbf{x},t)+\delta\hat{\theta}(\mathbf{x},t)]}.
\end{eqnarray}}\unskip
Expanding up to linear order in the density and phase fluctuations
yields
{\begin{eqnarray}\label{eq:BdGfieldequationTIHydrodynamic2}
\begin{array}{rcl}
\delta \hat{n}(\mathbf{x},t)&=&n(\mathbf{x},t)[\hat{\chi}(\mathbf{x},t)+\hat{\chi}^{\dagger}(\mathbf{x},t)],\vspace*{2.5pt}\\
 \delta \hat{\theta}(\mathbf{x},t)&=&\displaystyle -\frac{\mathrm{i}}{2}[\hat{\chi}(\mathbf{x},t)-\hat{\chi}^{\dagger}(\mathbf{x},t)].
\end{array}
\end{eqnarray}}\unskip
We can then rewrite the BdG equations
(\ref{eq:BdGfieldequationTIHydrodynamic}) as
{\begin{eqnarray}\label{eq:BdGfieldequationTIHydrodynamicPhase}
D_{\mathrm{t}}\frac{\delta \hat{n}}{n}&=&-\frac{1}{n}\boldsymbol{\nabla}\cdot\left[n\boldsymbol{\nabla}\frac{\hbar\delta \hat{\theta}}{m}\right],\nonumber\\
 D_{\mathrm{t}}\frac{\hbar \delta \hat{\theta}}{m}&=&-\left[\frac{T_n}{2m}+c^2\right]\frac{\delta \hat{n}}{n}.
\end{eqnarray}}\unskip
So far, within the Bogoliubov approximation, these equations are
\textit{exact}, resulting from a change of variables
$\{\hat{\varphi},\hat{\varphi}^\dagger\}\to\{\delta\hat{n},\delta\hat{\theta}\}$. 
Now, if we assume that we are in the Thomas--Fermi regime, where the
background condensate density smoothly varies on a sufficiently large
length scale, in the long-wavelength limit we can neglect the
contribution of $T_n$ at the r.h.s. of the second line. This precisely
amounts to work in the hydrodynamic approximation, where all the
contributions from the quantum potential are neglected, and we retrieve
the quantum version of Equation~(\ref{eq:EulerPotentialLinear}), from where
we find that $\square\delta\hat{\theta}=0$.

Therefore, in condensates, the gravitational analogy emerges in the
hydrodynamic approximation (equivalent to work in the long-wavelength
limit of the BdG equations above a condensate within the Thomas--Fermi
regime) as an equation of motion for the phase fluctuations which
mimics that of a massless scalar field in a curved spacetime described
by an acoustic\break metric (\ref{eq:RelativisticMetric}).

\subsection{Microscopic Hawking effect}

Due to the low temperature and genuine quantumness of Bose--Einstein
condensates, the gravitational analogy allows to study there the
Hawking effect, which is translated into the spontaneous emission of
phonon radiation by an acoustic horizon. In a pair of seminal works, it
was shown by Macher and Parentani~\cite{Macher2009a}, and by Recati,
Pavloff and Carusotto~\cite{Recati2009}, that the Hawking effect can be
studied within the full microscopic Bogoliubov framework without the
need of invoking the hydrodynamic approximation or even any metric at
all.

We now derive the Hawking effect from a microscopic approach along the
lines of Refs~\cite{Recati2009,Macher2009a}. For simplicity, hereafter
we focus on one-dimensional (1D) condensates, where $x$ will label the
1D spatial coordinate. Acoustic horizons are then reduced to discrete
points where the flow undergoes a subsonic/supersonic transition. 

We start by considering a stationary GP plane-wave solution
{\begin{equation}
\Psi_0(x)=\sqrt{n}\mathrm{e}^{\mathrm{i}(qx+\theta_0)}.
\end{equation}}\unskip
The associated BdG spectrum, resulting from Equation
(\ref{eq:BdGEigenmode}), is also described by plane waves with
wavevector $k$ and energy $\epsilon=\hbar\omega$, as given by the
dispersion relation
{\begin{equation}\label{eq:dispersionrelation}
\left(\omega-vk\right)^{2}=\Omega^2(k)=c^{2}k^{2}+\frac{\hbar^2k^{4}}{4m^2}=c^{2}k^{2}\left[1+\frac{(k\xi)^2}{4}\right],
\end{equation}}\unskip
where 
{\begin{equation}
    c=\sqrt{\frac{gn}{m}},\quad v=\frac{\hbar q}{m}
\end{equation}}\unskip
are the homogeneous sound and flow speeds. This is nothing else than
the usual Bogoliubov dispersion relation $\Omega(k)$ for a homogeneous
condensate at equilibrium plus a Doppler shift\break $\omega\to\omega-vk$,
resulting from the background condensate flow, which tilts the sound
cones. Remarkably, the Bogoliubov dispersion relation is superluminal,
with the healing length\break $\xi\equiv \hbar/mc$ playing the role of a
Planck length scale that controls the UV physics.

For given $\omega$, the dispersion relation
(\ref{eq:dispersionrelation}) provides $4$ wavevectors, labeled as
$k_a(\omega)$ and given by the roots of the fourth order polynomial,
which can be either real (describing propagating solutions) or complex
(describing exponentially growing/decaying solutions). The
corresponding BdG spinor for each wavevector $k_a$ reads
{\advance\jot by 8pt\begin{eqnarray} \label{eq:PlaneWaveSpinors}
s_{a,\omega}(x) & = & \frac{\mathrm{e}^{\mathrm{i}k_{a}\left(\omega\right)x}}{\sqrt{2\rmpi|w_{a}\left(\omega\right)|}}\left[\begin{array}{@{}c@{}}
\mathrm{e}^{\mathrm{i}(qx+\theta_0)}u_{a}(\omega)\vspace*{2pt}\\
\mathrm{e}^{-\mathrm{i}(qx+\theta_0)}v_{a}(\omega)
\end{array}\right] \nonumber\\
\left[\begin{array}{@{}c@{}}
\nonumber u_{a}(\omega)\\
v_{a}(\omega)
\end{array}\right]&=&N_a(\omega)\left[\begin{array}{@{}c@{}}
\frac{\hbar k_{a}^{2}\left(\omega\right)}{2m}+[\omega-vk_{a}\left(\omega\right)]\vspace*{5pt}\\
\frac{\hbar k_{a}^{2}\left(\omega\right)}{2m}-[\omega-vk_{a}\left(\omega\right)]
\end{array}\right]\\
N_a(\omega)&=&\left(\frac{m}
{2\hbar k_{a}^{2}\left(\omega\right)\left|\omega-vk_{a}\left(\omega\right)\right|}\right)^{\tfrac{1}{2}},
\end{eqnarray}}\unskip
with $u_a(\omega),v_a(\omega)$ the usual Bogoliubov components for a
homogeneous condensate, satisfying
${|u_{a}(\omega)|}^2-{|v_{a}(\omega)|}^2=\pm 1$, and
$w_{a}(\omega)\equiv[\mathrm{d}k_{a}(\omega)/\mathrm{d}\omega]^{-1}$ 
the group velocity, included here in order to normalize the propagating
modes in frequency domain,
$(s_{a,\omega}|s_{a,\omega'})=\pm\delta(\omega-\omega')$;
all normalization factors can be removed for complex wavevector
solutions, where they do not play any role. It is easy to check that
the dispersion relation (\ref{eq:dispersionrelation}) possess the
symmetry $k_{a}(\omega)=-k_a(-\omega)$, which implies
$\bar{s}_{a,\omega}=s_{a,-\omega}$. Hence, the $\pm$ branches of 
{\begin{equation}\label{eq:Dispersion}
    \omega(k)=vk\pm \Omega(k), 
\end{equation}}\unskip
depicted in blue (red) in Figure~\ref{fig:DispersionRelation},
respectively, are conjugate of each other, with the $\pm$ sign also
corresponding to the norm of the modes. 

\begin{figure}[t!]
\vspace*{2pt}
\includegraphics{fig01}
\vspace*{2pt}
\caption{\label{fig:DispersionRelation}Dispersion relation of a homogeneous flowing condensate. The
blue/red lines signal the $\pm$ branches of
Equation~(\ref{eq:Dispersion}). (a)~Subsonic regime with Mach number
$M=v/c=0.5$. (b) Supersonic regime with Mach number $M=v/c=2$. For a
certain frequency below the cutoff frequency $\omega_{{\max}}$, all
wavevectors are purely real (horizontal dashed line). The BCL mode has
zero frequency and finite wavevector $k_{\mathrm{BCL}}$ (black square). (c)
Schematic spacetime diagram of the modes in (a)--(b) when placed in an
analogue black hole, where the upstream/downstream region is
subsonic/supersonic. The position of the acoustic horizon is signaled
by a vertical solid black line. Incoming/outgoing modes are depicted as
solid/dashed lines. Horizontal arrows indicate the characteristic
correlations of both the Hawking and Andreev
effects.}
\vspace*{3pt}
\end{figure}

The dispersion relation displays two qualitatively different regimes,
depending on whether the flow is subsonic ($v<c$) or supersonic
($v>c$). In the subsonic regime, Figure~\ref{fig:DispersionRelation}a,
there are $2$ real wavevectors and $2$ complex ones. The propagating
solutions are labeled as $u-\mathrm{in}$ and $u-\mathrm{out}$, where ``in''
(``out'') indicates if the group velocity is positive (negative); the
motivation behind this notation will become clearer later. The flow is
energetically stable since all modes have positive energy. In the
supersonic regime, Figure~\ref{fig:DispersionRelation}b, for any
frequency $-\omega_{{\max}}<\omega<\omega_{{\max}}$, there are $4$
different propagating solutions, where the cutoff frequency
$\omega_{{\max}}$ is given by  
{\advance\jot by 8pt\begin{eqnarray}\label{eq:CutoffFrequency}
\begin{array}{rcl}
    \omega_{{\max}}&=&vk_{{\max}}-\Omega(k_{{\max}}),  \vspace*{2.5pt}\\
\displaystyle  (k_{{\max}}\xi)^2&=&\displaystyle \frac{M^2-4+M\sqrt{M^2+8}}{2},
\end{array}
\end{eqnarray}}\unskip
$M=v/c$ being the Mach number, with $M>1$ for a supersonic flow. The
``in'' (``out'') label is reverted here and now indicates if the group
velocity is negative (positive). The $d1$ modes are those normal, with
positive energy, while the $d2$ modes are those anomalous, with
negative energy. The presence of the anomalous $d2$ modes reveals the
energetic instability of a supersonic flow. This is a consequence of
the Landau criterion for superfluidity, which predicts the appearance
of a zero-frequency mode in a supersonic flow, the celebrated
Bogoliubov--Cherenkov--Landau (BCL) mode, with a finite wavevector $\hbar
k_{\mathrm{BCL}}=2m\sqrt{v^2-c^2}$ computed by
$\Omega(k_{\mathrm{BCL}})=vk_{\mathrm{BCL}}$ (black square in
Figure~\ref{fig:DispersionRelation}b). This results in the coherent
excitation of the BCL mode by the presence of any obstacle in a
supersonic flow, spoiling its superfluidity~\cite{Carusotto2006}; see
Equation~(\ref{eq:BCLStimulation}) and ensuing discussion for more
details. Nevertheless, supersonic flows are dynamically
stable~\cite{Mayoral2011}, since energetic instability is only a
necessary condition for dynamical instability, but not
sufficient~\cite{Wu2003}. This is directly seen from
Equation~(\ref{eq:energeticdynamical}): the presence of complex modes
(dynamical instability), which have zero norm by virtue of
Equation~(\ref{eq:EigenOrto}), necessarily implies that $\Lambda$
cannot be a positive-definite operator (energetic instability).

We are now in a position to study a black-hole (BH) solution, defined
here as a stationary 1D GP solution with two asymptotic homogeneous
regions, one subsonic and one supersonic, flowing from subsonic to
supersonic. When the flow travels from supersonic to subsonic, we have
a white hole (WH) solution, the time reversal of a BH (obtained simply
by conjugation of the wavefunction). Continuity of the GP wavefunction
implies that both BH and WH solutions always possess, at least, one
acoustic horizon, where $v(x)=c(x)$. We denote the region between the
two asymptotic regions, in which the acoustic horizon is located, as
the scattering region. The emergence of BH configurations is expected
on quite general grounds, since it was shown by Michel, Parentani and
Zegers~\cite{Michel2016} that stationary asymptotically uniform flows
are attractor solutions, providing an analogue version of the
celebrated no-hair\break theorem.

By convention, for BH solutions, we take the flow velocity always
positive, so the upstream subsonic region (labeled as ``u'') is located 
at $x\rightarrow-\infty$, while the downstream supersonic region
(labeled as ``d'') is located at $x\rightarrow\infty$. This convention
hence matches the notation previously introduced in
Figure~\ref{fig:DispersionRelation}, since ``in'' modes are incoming
(traveling towards the horizon), and ``out'' modes are outgoing
(traveling outwards). The asymptotic flow velocity, sound speed, and
healing length, are labeled as $v_{u,d},c_{u,d},\xi_{u,d}$. 

The BdG modes of a BH solution are asymptotically given in terms of
linear combinations of the plane-wave spinors
$s_{i-\mathrm{in/out},\omega}$ of Equation~(\ref{eq:PlaneWaveSpinors}),
$i=u,d1,d2$, representing the different incoming and outgoing
scattering channels. Throughout this work, we operate just with
positive frequencies $\omega>0$, and the remaining part of the spectrum
is obtained by conjugation. In particular, the scattering problem for
any positive frequency $0<\omega<\omega_{{\max}}$ involves the normal
$u,d1$ channels, and the conjugate of the anomalous $d2$ channel
(horizontal dashed line in Figure~\ref{fig:DispersionRelation}), so
$s_{d2-\mathrm{in/out},\omega}$ has negative norm. The retarded
(``in'') scattering states $z^{(+)}_{i,\omega}$ are global eigenmodes
of the stationary BdG equations (\ref{eq:BdGEigenmode}) with positive
frequency, presenting unit amplitude in the asymptotic incoming channel
$i$ and zero in the other incoming channels. The  amplitude of the
asymptotic ``out'' scattering channels is determined by the $S$-matrix,
as usual in scattering theory. For example, the scattering state
$z^{(+)}_{d2,\omega}$ asymptotically reads
{\begin{eqnarray}\label{eq:scatteringchannelstate}
\begin{array}{c}
\displaystyle z^{(+)}_{d2,\omega}\left(x\rightarrow-\infty\right)=S_{ud2}\left(\omega\right)s_{u-\mathrm{out},\omega}(x),\vspace*{2.5pt}\\
\displaystyle  z^{(+)}_{d2,\omega}\left(x\rightarrow\infty\right)=s_{d2-\mathrm{in},\omega}(x)+S_{d1d2}\left(\omega\right)s_{d1-\mathrm{out},\omega}(x)
+S_{d2d2}\left(\omega\right)s_{d2-\mathrm{out},\omega}(x).
\end{array}
\end{eqnarray}}\unskip
Similar expressions can be provided for the remaining ``in'' scattering
states. The advanced (``out'') scattering states
$z^{(-)}_{i,\omega}(x)$ are the outgoing analogues of the ``in''
states, having unit amplitude in the outgoing channel $i$ and zero in
the other outgoing channels. They are characterized by the inverse of
the scattering matrix $S(\omega)$,
{\begin{equation}
    z^{(-)}_{i,\omega}(x)=\sum_{j=u,d1,d2} S^{-1}_{ji}(\omega)z^{(+)}_{j,\omega}(x).
\end{equation}}\unskip
A schematic spacetime diagram of these modes is presented in
Figure~\ref{fig:DispersionRelation}c. Retarded scattering states are
characterized by one incoming channel (solid lines), which is
eventually scattered into the outgoing channels (dashed lines). In
contrast, advanced scattering states are labeled by just one outgoing
channel, whose time-reversed trajectory scatters at the horizon into
all incoming channels.

By invoking the conservation of the quasiparticle current
(\ref{eq:QuasiparticleCurrent}) for an arbitrary linear combination of
``in'' scattering states, it is shown that the $S$-matrix is
pseudo-unitary, i.e.,
{\begin{equation}\label{eq:pseudounitarity}
S^{\dagger}\eta S=\eta\equiv{\rm diag}(1,1,-1).
\end{equation}}\unskip
Thus, $S\in U(2,1)$, which implies $S^{-1}=\eta S^\dagger \eta$. A
direct consequence of the pseudounitarity of $S$ is that the scattering
states are orthonormal,
{\begin{equation}
(z^{(+)}_{i,\omega}|z^{(+)}_{j,\omega'})=(z^{(-)}_{i,\omega}|z^{(-)}_{j,\omega'})=\eta_{ij}\delta(\omega-\omega').
\end{equation}}\unskip

Since they form a complete orthonormal basis, the quantum fluctuations
of the field operator can be expanded in terms of the scattering states
as
{\begin{eqnarray}\label{eq:BHFieldOperator}
 \hat{\Phi}(x) = \displaystyle\sum_{I=u,d1}\int_{0}^{\infty}\mathrm{d}\omega\,
 [z^{(+)}_{I,\omega}(x)\hat{a}_{I}(\omega)+\bar{z}^{(+)}_{I,\omega}(x)
 \hat{a}_{I}^{\dag}(\omega)]
+\int_{0}^{\omega_{{\max}}}\mathrm{d}\omega\,[z^{(+)}_{d2,\omega}(x)\hat{a}_{d2}^{\dagger}(\omega)+\bar{z}^{(+)}_{d2,\omega}(x)\hat{a}_{d2}(\omega)].\nonumber\\
\end{eqnarray}}\unskip
A similar expression can be written using the ``out'' scattering states
after replacing $z^{(+)}_{i,\omega}(x)$ by $z^{(-)}_{i,\omega}(x)$, and
the ``in'' quantum amplitudes $\hat{a}_i(\omega)$ by the ``out'' ones
$\hat{b}_i(\omega)$, which are related through the scattering matrix as
{\begin{equation}\label{eq:inoutmodesrelation}
\left[\begin{array}{@{}c@{}}
\hat{b}_{u}\\
\hat{b}_{d1}\\
\hat{b}_{d2}^{\dagger}
\end{array}\right] = \left[\begin{array}{ccc}S_{uu}&S_{ud1}&S_{ud2}\\
S_{d1u}&S_{d1d1}&S_{d1d2}\\
S_{d2u}&S_{d2d1}&S_{d2d2}\end{array}\right]\left[\begin{array}{@{}c@{}}
\hat{a}_{u}\\
\hat{a}_{d1}\\
\hat{a}_{d2}^{\dagger}
\end{array}\right].
\end{equation}}\unskip
This is a Bogoliubov relation, mixing annihilation with creation
operators. It stems from the anomalous character of the
$z^{(\pm)}_{d2,\omega}$ scattering states, which have a negative norm
inherited from the corresponding anomalous scattering channels, and
hence their amplitudes behave as creation instead of annihilation
operators (see Equation~(\ref{eq:Aniquilacion}) and ensuing
discussion). In the following, we reserve lowercase Latin indices $i,j$
to label all channels, $i=u,d1,d2$, while uppercase Latin indices $I,J$
just label normal channels, $I=u,d1$. Lowercase Latin indices $a,b$
will label general BdG modes, either outgoing or incoming, either
propagating or not. In addition, the normal-normal and
anomalous-anomalous scattering processes (characterized by the
$S$-matrix elements $S_{IJ},S_{d2d2}$) are denoted as normal, while the
normal-anomalous and anomalous-normal scattering processes
(characterized by the $S$-matrix elements $S_{d2I},S_{Id2}$) are
denoted as anomalous.

The origin of the Hawking effect is the degeneracy of the vacuum of the
Bogoliubov theory, revealed by the anomalous sector of $\hat{K}$, 
{\begin{eqnarray}
 \hat{K}_{\mathrm{H}}\equiv \sum_{i,j}\int^{\omega_{{\max}}}_0\mathrm{d}\omega~\hbar\omega \hat{a}_{i}^{\dagger}(\omega)\eta_{ij} \hat{a}_{j}(\omega)   
 =\sum_{i,j}\int^{\omega_{{\max}}}_0\mathrm{d}\omega~\hbar\omega \hat{b}_{i}^{\dagger}(\omega)\eta_{ij} \hat{b}_{j}(\omega).
\end{eqnarray}}\unskip
Both the incoming vacuum $\hat{a}_{i}(\omega)\ket{0_{\mathrm{in}}}=0$
and the outgoing vacuum $\hat{b}_{i}(\omega)\ket{0_{\mathrm{out}}}=0$
satisfy $\hat{K}_{\mathrm{H}}\ket{0_{\mathrm{in}}}
=\hat{K}_{\mathrm{H}}\ket{0_{\mathrm{out}}}=0$. However, they do not
represent the same quantum state, as can be seen from the non-vanishing
population of normal outgoing modes in the incoming vacuum,
\def\bra#1{\mathinner{\mathop{\langle}{#1}|}}
{\begin{equation}\label{eq:hawkingef}
\bra{0_{\mathrm{in}}}
\hat{b}_{I}^{\dagger}(\omega)\hat{b}_{I}(\omega')
\ket{0_{\mathrm{in}}}
=\delta(\omega-\omega'){|S_{Id2}(\omega)|}^2\neq 0.
\end{equation}}\unskip
The Hawking effect is recovered for $I=u$, representing a spontaneous
outgoing flux of particles in the subsonic region (i.e., the exterior
of the black hole) in the absence of incoming radiation. This emission
is correlated with that of anomalous outgoing $d2$ modes into the
supersonic region (horizontal arrow in
Figure~\ref{fig:DispersionRelation}c), which are referred to as the
partner modes of the Hawking effect. 


The case $I=d1$ is characterized by the anomalous reflection
coefficient $S_{d1d2}$, representing the bosonic analogue of an
incident hole that is reflected as a particle in the normal side of a
normal/superconductor junction, the celebrated Andreev reflection. In
analogy with the Hawking effect, throughout this work we will refer to
the spontaneous emission of normal outgoing $d1$ modes into the
supersonic region as the Andreev effect. Once again, this emission is
correlated with that of anomalous outgoing $d2$ modes (horizontal arrow
in Figure~\ref{fig:DispersionRelation}c). A discussion of the
equivalence between the Andreev scattering picture and the spontaneous
emission of quasiparticle pairs in the context of normal/superconductor
interfaces was given in Refs~\cite{Samuelsson2003,Prada2004}. In
general, by invoking the Ginzburg--Landau (GL) order parameter, which
obeys a non-linear Schr\"odinger equation equivalent to the
time-independent GP equation~(\ref{eq:TIGP}), we can extend the analogy
with superconductivity and identify normal metals with supersonic
regions, and superconductors with subsonic regions. In fact, the
physics of real black holes and the physics of superconductors are in
close relationship~\cite{Manikadan2017,Manikadan2020}.


\subsection{Analytical BH solutions}

We present here canonical analytical BH solutions in condensates. In
the following, we set units and rescale the GP wavefunction as
{\begin{equation}
    \hbar=m=c_u=k_B=1,\quad \Psi_0\to \sqrt{n_u}\Psi_0.
\end{equation}}\unskip

The simplest BH solution is provided by the flat-profile model,
originally introduced in Ref.~\cite{Carusotto2008}, where the
plane wave $\Psi_0(x)=\mathrm{e}^{\mathrm{i}qx}$ is a solution of the GP
equation~(\ref{eq:TDGP}) at all times since the coupling constant
$g(x,t)$ and the external potential $V(x,t)$ are tuned in such a way
that
{\begin{equation}\label{eq:FlatProfileCondition}
    g(x,t)n_u+V(x,t)=E_b,
\end{equation}}\unskip
with $E_b$ some constant energy that can be subtracted from the
Hamiltonian. Nevertheless, even though the flow velocity is constant
and homogeneous, $v(x,t)=q$, the BdG modes do experience non-trivial
dynamics as the sound speed is $c^2(x,t)=g(x,t)n_u$; see
Equation~(\ref{eq:BdGfieldequationTIHydrodynamicPhase}). 

In particular, we can choose a time-independent piecewise homogeneous
dependence for the coupling constant, $g(x,t)=g(x)$, where the BdG
solutions within each homogeneous region are spanned by the plane-wave
spinors (\ref{eq:PlaneWaveSpinors}), thus allowing for simple
analytical calculations. Specifically, for the BH solution, we take
{\begin{equation}\label{eq:FlatProfile}
g(x)n_u=\left\{ \begin{array}{@{}lc}
1, & x< 0,\\
c_2^2, &  x \geq 0,
\end{array}\right.
\end{equation}}\unskip
with $c_2<q<1$, $c_2$ being the supersonic speed of sound and $q$
representing the subsonic Mach number. Hence, we reach a BH solution
whose acoustic horizon is placed at $x=0$, as depicted in
Figure~\ref{fig:BHModels}a.

\begin{figure}[t!]
\includegraphics{fig02}
\caption{\label{fig:BHModels}Upper row: Sound (solid blue) and flow
(dashed red) velocity profiles for different BH solutions with subsonic
Mach number $q=0.5$. (a) Flat profile. The shaded area indicates the
supersonic region, where $c_2=0.25$. (b) Waterfall potential. The shaded
area indicates the region where the step potential $V(x)=-V_0\Theta(x)$
is present. (c) Delta barrier. The arrow indicates the position of the
delta potential $V(x)=Z\delta(x)$. Lower row: (d)--(f) Hawking (solid
blue) and Andreev (solid red) spectra of the BH solutions above. Dashed
gold (cyan) line is a gray-body fit of the Hawking (Andreev) spectrum,
Equation~(\ref{eq:GrayBody}). Inset: Frequency-dependent Hawking (solid
black) and Andreev (solid magenta) temperatures $T_{u,d1}(\omega)$,
Equation~(\ref{eq:OmegaTemperature}). Horizontal dashed green line
marks the predicted Hawking temperature for a soliton,
Equation~(\ref{eq:SolitonHawking}).}
\end{figure}

However, in practice, the flat-profile condition
(\ref{eq:FlatProfileCondition}) is extremely challenging to implement
experimentally. More realistic models of BH configurations only involve
external potentials $V(x)$, which can be easily manipulated in the
laboratory, leaving the interaction strength $g$ constant.
Consequently, since the current $J$ is uniform for a stationary 1D
solution as dictated by the continuity equation (first line of
Equation~(\ref{eq:EulerPhaseCondensate})), we simply have that
{\begin{equation}
    c(x)=\sqrt{n(x)},\quad v(x)=\frac{J}{n(x)}=\frac{J}{c^2(x)}.
\end{equation}}\unskip
Typically, these BH solutions are described by a gray soliton in the
upstream region
{\begin{eqnarray}\label{eq:GraySoliton}
    \Psi_{0}(x)&=&\mathrm{e}^{\mathrm{i}(qx+\theta_0)}
    [q+\mathrm{i}\gamma_q(x-x_0)], \nonumber\\
    \gamma_q(x)&\equiv & \sqrt{1-q^2}\tanh(\sqrt{1-q^2}x).
\end{eqnarray}}\unskip
A gray soliton exponentially approaches a subsonic plane wave
$\Psi_0(x)\xrightarrow[x\to\pm\infty]{} \mathrm{e}^{\mathrm{i}qx}$
since $\gamma_q(\pm \infty)=\pm \sqrt{1-q^2}$, with $q< 1$ the
asympotic Mach number $M_u=q$. In our units, $q$ is also the minimum
soliton amplitude as well as the value of the conserved current, $J=q$.
In the downstream region, the BH solutions are given by the supersonic
plane wave 
{\begin{equation}
    \Psi_0(x)=c_d\mathrm{e}^{\mathrm{i}v_d x},\quad v_d=\frac{q}{c^2_d},
\end{equation}}\unskip
with $c_d<v_d$ the corresponding supersonic sound and flow velocities;
notice that $v_d$ is fixed by current conservation once $c_d$ is given.

For example, the waterfall configuration~\cite{Larre2012} accelerates
the atoms to supersonic speeds by means of an attractive step potential
of the form
{\begin{equation}\label{eq:waterfall}
V(x)=-V_0\Theta(x),\quad V_0=\frac{1}{2}\left(q^2+\frac{1}{q^2}\right)-1,
\end{equation}}\unskip
where $\Theta(x)$ is the Heaviside function. This gives rise to a
stationary GP wavefunction
{\begin{equation}\label{eq:CompactBHWF}
\Psi_0(x)=\left\{ \begin{array}{@{}cc}
\mathrm{e}^{\mathrm{i}qx}[q+\mathrm{i}\gamma_q(x)], & x< 0,\vspace*{2.5pt}\\
q\mathrm{e}^{\mathrm{i}\frac{x}{q}}, &  x \geq 0,
\end{array}\right.
\end{equation}}\unskip
which corresponds to half a gray soliton for $x<0$, and a homogeneous
supersonic flow with $c_d=q<v_d=1/q$ for $x>0$, where the attractive
potential is present. The resulting BH solution is represented in
Figure~\ref{fig:BHModels}b. The waterfall configuration is quite
relevant because it provides a simple theoretical model of the Technion
experiment~\cite{Lahav2010,Steinhauer2014,Steinhauer2016,deNova2019,Kolobov2021}.

Another possibility is to use a repulsive localized potential, which
can be modeled by a delta barrier of the form 
{\begin{equation}\label{eq:DeltaPotential}
    V(x)=Z\delta(x),\quad Z=\frac{(1-c^2_d)\sqrt{c^2_d-q^2}}{2c^2_d}.
\end{equation}}\unskip
This delta potential introduces a discontinuity in the derivative of
the GP wavefunction, $\Psi'_0(0^+)-\Psi'_0(0^-)=2Z\Psi_0(0)$. The
resulting BH solution is similar to that of the waterfall model:
{\begin{equation}\label{eq:CompactDelta}
\Psi_0(x)=\left\{ \begin{array}{@{}cc}
\mathrm{e}^{\mathrm{i}(qx+\theta_0)}[q+\mathrm{i}\gamma_q(x-x_0)] & x< 0,\vspace*{2.5pt}\\
c_d\mathrm{e}^{\mathrm{i}v_dx}, &  x \geq 0,
\end{array}\right.
\end{equation}}\unskip
where $x_0,\theta_0$ are such that the wavefunction is continuous and
{\begin{equation}\label{eq:DeltaConstraints}
    c_d=\frac{\sqrt{q^2+\sqrt{q^4+8q^2}}}{2},\quad v_d=\frac{q}{c^2_d}.
\end{equation}}\unskip
This BH solution is represented in Figure~\ref{fig:BHModels}c.

Remarkably, the scattering states associated to these BH solutions can
be also computed analytically because the stationary BdG solutions for
a gray soliton (\ref{eq:GraySoliton}) are known (see
Ref.~\cite{Zapata2011} for the technical details). They take a
similar form to the homogeneous solutions (\ref{eq:PlaneWaveSpinors}),
namely
{\begin{eqnarray}
 \zeta_{a,\omega}(x)&=&\frac{\mathrm{e}^{\mathrm{i}k_{a}(\omega)x}}{\sqrt{2\rmpi|w_{a}(\omega)|}}\left[\begin{array}{@{}c@{}}
\mathrm{e}^{\mathrm{i}(qx+\theta_0)}u_{a,\omega}(x)\\
\mathrm{e}^{-\mathrm{i}(qx+\theta_0)}v_{a,\omega}(x)
\end{array}\right], \nonumber
\\  \left[\begin{array}{@{}c@{}}
u_{a,\omega}(x)\\
v_{a,\omega}(x)
\end{array}\right]&=&N_a(\omega)\left[\begin{array}{r}
\left(1+\dfrac{k_a(\omega)}{\omega}\left[\dfrac{k_a(\omega)}{2}+\mathrm{i}\gamma_q(x-x_0)\right]\right)^2\vspace*{2.5pt}\\
-\left(1-\dfrac{k_a(\omega)}{\omega}\left[\dfrac{k_a(\omega)}{2}+\mathrm{i}\gamma_q(x-x_0)\right]\right)^2
\end{array}\right],\nonumber\\
N_a(\omega)&=&\frac{\omega}{\sqrt{8k^2_a(\omega)|\omega-vk_a(\omega)|}}.\label{eq:solitonspinors}
\end{eqnarray}}\unskip
In fact,
$\zeta_{a,\omega}(x)\xrightarrow[x\to\pm\infty]{}s_{a,\omega}(x)$ as
the soliton asymptotically approaches a subsonic plane wave
$\mathrm{e}^{\mathrm{i}qx}$, whose dispersion relation yields the
wavevectors $k_a(\omega)$. Thus, the scattering states
$z^{(\pm)}_{i,\omega}$ for the waterfall (\ref{eq:CompactBHWF}) and
delta (\ref{eq:CompactDelta}) BH solutions are obtained by matching at
$x=0$ the upstream soliton spinors $\zeta_{a,\omega}$ with the
corresponding downstream supersonic plane-wave spinors $s_{a,\omega}$;
notice that in the subsonic region one also needs to include the
evanescent solution $\zeta_{\mathrm{ev},\omega}$ with complex
wavevector $k_{\mathrm{ev}}(\omega)$ which exponentially decays at
$x\to-\infty$, $\mathrm{Im}\,k_{\mathrm{ev}}(\omega)<0$. The
flat-profile BH solution does not even involve the soliton spinors
$\zeta_{a,\omega}(x)$ because the upstream region is also homogeneous,
and the subsonic plane-wave spinors $s_{a,\omega}(x)$ are used instead.

The resulting anomalous scattering coefficients
${|S_{ud2}|}^2,{|S_{d1d2}|}^2$ characterizing the Hawking and Andreev
effects for each BH solution are depicted in
Figures~\ref{fig:BHModels}d--f as solid blue and red lines,
respectively. We can compare these results, fully derived within the
BdG microscopic framework, with the predictions from the gravitational
analogy in the hydrodynamic limit~\cite{Leonhardt2003a}
{\begin{equation}\label{eq:Hydrodynamic}
{|S_{ud2}(\omega)|}^2=\frac{1}{\mathrm{e}^{\frac{\omega}{T_{\mathrm{H}}}}-1},
\quad {|S_{d1d2}(\omega)|}^2=0,
\end{equation}}\unskip
where $T_{\mathrm{H}}$ is the Hawking temperature, given here by~\cite{Visser1998}
{\begin{equation}
    T_{\mathrm{H}}=\frac{1}{2\rmpi}\left|c'(x_{\mathrm{H}})-v'(x_{\mathrm{H}})\right|,
\end{equation}}\unskip
with $x_{\mathrm{H}}$ the position of the acoustic horizon. Notice that, for the
flat-profile configuration, this temperature is formally infinite due
to the discontinuity of the sound speed. On the other hand, for the
waterfall and delta configurations, it is easy to see that the acoustic
horizon of a gray soliton (\ref{eq:GraySoliton}) is placed where
$c(x)=v(x)=q^{\frac{1}{3}}$, yielding a predicted Hawking temperature
{\begin{equation}\label{eq:SolitonHawking}
 T_{\mathrm{H}}(q)=\frac{3}{2\rmpi}(1-q^{\frac{2}{3}})\sqrt{1-q^{\frac{4}{3}}}\leq T_{\mathrm{H}}(0)=\frac{3}{2\rmpi}<\frac{1}{2}.
\end{equation}}\unskip

A necessary condition for the validity of the hydrodynamic
approximation is that $T_{\mathrm{H}}\ll \omega_{{\max}}$. It is easily seen from
Equation~(\ref{eq:CutoffFrequency}) that $k_{{\max}}\leq v_d$ and that
$\omega_{{\max}}\leq v^2_d/2$; these inequalities are saturated in the
infinite supersonic Mach number limit, $c_d=0$. For the flat-profile
configuration, this implies a relatively small cutoff frequency as
$\omega_{{\max}}\leq q^2/2 < 1/2$. For the delta and waterfall
configurations, we respectively have that 
{\begin{equation}\label{eq:CutoffHawking}
    \omega_{{\max}}\leq \frac{8}{(q+\sqrt{q^2+8})^2}<1,\quad \omega_{{\max}}\leq \frac{1}{2q^2}.
\end{equation}}\unskip

As suggested in Ref.~\cite{Larre2012}, we can check the agreement
with the gravitational analogy by fitting the Hawking spectrum to a
gray-body distribution. We also extend this fit to the Andreev spectrum
as
{\begin{equation}\label{eq:GrayBody}
{|S_{Id2}(\omega)|}^2=\frac{\Gamma_I}{\mathrm{e}^{\frac{\omega}{T_I}}-1},
\end{equation}}\unskip
where $\Gamma_{I}$ is a gray-body factor and $T_{I}$ is the effective
temperature of the spectrum. In the hydrodynamic limit, we can expect
$\Gamma_u\simeq 1$ and $T_u\simeq T_{\mathrm{H}}$. The result of the fit for the
Hawking (Andreev) spectrum is depicted as a dashed gold (cyan) line in
lower row of Figure~\ref{fig:BHModels}. Another way to quantify the
Planckianity, proposed by Macher and Parentani~\cite{Macher2009a}, is
an $\omega$-dependent temperature, defined through
{\begin{equation}\label{eq:OmegaTemperature}
|S_{Id2}(\omega)|^2\equiv \frac{1}{\mathrm{e}^{\frac{\omega}{T_I(\omega)}}-1}.
\end{equation}}\unskip
$T_{u,d1}(\omega)$ is shown in the inset of lower
Figure~\ref{fig:BHModels} as a solid black (magenta) line. In the
Hawking case, we can compare the result with the predicted Hawking
temperature (\ref{eq:SolitonHawking}), indicated as a horizontal dashed
green line.

A gray-body distribution displays an excellent agreement with the
Hawking and Andreev spectra in all cases, even for the flat-profile
model (Figure~\ref{fig:BHModels}d), where the system is far from the
hydrodynamic regime and the predicted Hawking temperature is formally
infinite. The agreement is particularly good at low frequencies
because, in general, the scattering coefficients
$|S_{id1}(\omega)|^2,|S_{id2}(\omega)|^2$ display a universal scaling
${\sim}1/\omega$ at low frequencies~\cite{deNova2014} (the remaining
column $|S_{iu}(\omega)|^2$ approaches a finite value in this limit). A
gray-body distribution can reproduce exactly this behavior up to
corrections $O(\omega)$,
{\begin{equation}
    \frac{\Gamma}{\mathrm{e}^{\frac{\omega}{T}}-1}\simeq \frac{\Gamma T}{\omega}-\frac{\Gamma}{2}+O(\omega).
\end{equation}}\unskip

Regarding the effective temperatures $T_{I}(\omega)$, we observe that
$T_u(\omega)$ depends very mildly on the frequency, approaching the
predicted Hawking temperature in the low-frequency limit, while
$T_{d1}(\omega)$ depends more strongly on the frequency. This is quite
natural, since there is no analytical prediction of a Planckian
distribution for the Andreev spectrum. The effective temperature
$T_{d1}(\omega)$ also measures the strength of the Andreev effect,
which satisfies $|S_{d1d2}(\omega)|^2\ll|S_{ud2}(\omega)|^2$ for both
the waterfall and delta models, as expected from the gravitational
prediction (\ref{eq:Hydrodynamic}), and $|S_{d1d2}(\omega)|^2\sim
|S_{ud2}(\omega)|^2$ for the flat profile. This difference is due to
the sharpness of the flat-profile horizon, far from the hydrodynamic
limit. 

In the waterfall case, low subsonic Mach numbers imply large supersonic
Mach numbers, so $\omega_{{\max}}\simeq 1/2q^2$ (see
Equation~(\ref{eq:CutoffHawking})). This large cutoff frequency makes
dispersive effects important at intermediate frequencies $\omega
\lesssim \omega_{{\max}}$, explaining the strong deviations of
$T_u(\omega)$ from the predicted Hawking temperature as well as the
dominance of the Andreev effect, $T_{d1}(\omega)>T_u(\omega)$ (inset of
Figure~\ref{fig:BHModels}e). Nevertheless, the relevant part of the
Hawking spectrum is still Planckian (solid blue and dashed gold lines
in main Figure~\ref{fig:BHModels}e), with $T_u(\omega)\simeq T_{\mathrm{H}}$,
because it is restricted to the low-frequency dispersionless regime due
to the smallness of the Hawking temperature, $T_{\mathrm{H}}\ll \omega_{{\max}}$. 

Finally, close to $\omega_{{\max}}$, the Planckianity is necessarily
spoiled for both the Hawking and Andreev spectra since there the
anomalous coefficients $|S_{Id2}(\omega)|^2,|S_{d2I}(\omega)|^2$ vanish
as a power law\break ${\sim}\sqrt{\omega_{{\max}}-\omega}$~\cite{deNova2014}.
This is revealed by the departure from the gray-body fit, especially in
the flat-profile case, and by the sudden drop of $T_{I}(\omega)$ in all
configurations.


\section{Resonant Andreev--Hawking radiation}\label{sec:ResonantHawking}

The thermal character of the Andreev and Hawking spectra discussed
above makes quite difficult to isolate their signal in a real
experiment, as it can be quite easily misidentified or overshadowed by
another background thermal component. It was suggested by Zapata, Albert, 
Parentani and Sols~\cite{Zapata2011} that resonant BH
configurations could provide a strategic advantage due to their highly
non-thermal frequency dependence. We discuss in this section how
resonant BH configurations emerge in gravitational analogues and
propose possible experimental realizations based on the use of optical
lattices.

\subsection{Resonant BH configurations}

The first proposed model of resonant BH configuration~\cite{Zapata2011} consisted of a condensate flowing through a double-delta barrier
{\begin{equation}\label{eq:DoubleDeltaPotential}
    V(x)=Z\left[\delta\left(x+\frac{L}{2}\right)+\delta\left(x-\frac{L}{2}\right)\right].
\end{equation}}\unskip
The resulting BH solution, represented in Figure~\ref{fig:Resonant}a,
corresponds to a gray soliton for $x<-L/2$ and a supersonic homogeneous
plane wave for $x>L/2$, with the same relation between the parameters
$q,c_d,v_d$ as for the delta BH solution,
Equation~(\ref{eq:DeltaConstraints}). The difference is that now the
value of $Z$ is not fine-tuned as in
Equation~(\ref{eq:DeltaPotential}), and the stationary GP solution is
described by a cnoidal wave for $|x|<L/2$, given in terms of elliptic
functions (see Ref.~\cite{Zapata2011} for the technical details).
Specifically, for a certain value of the asymptotic subsonic flow
velocity $q$, there are only stationary GP solutions for $Z\in
[Z_{{\min}}(q),Z_{{\max}}(q)]$, and the number of available solutions
grows with the interbarrier distance $L$, displaying an increasing
number of cnoidal periods. These BH solutions may include several local
acoustic horizons, always an odd number of them to ensure that there is
a global transition from an asymptotic subsonic flow to a supersonic
one.

\begin{figure}[t!]
\includegraphics{fig03}
\caption{\label{fig:Resonant}Upper row: Sound (solid blue) and flow (dashed red) velocity
profiles for different resonant BH solutions. (a) Double delta. The
asymptotic subsonic flow velocity is $q=0.01$. The arrows indicate the
position of the delta barriers, whose amplitude is $Z=2.2$. The length
of the resonant cavity is $L\approx 3.62$. (b) Resonant flat-profile.
The global flow velocity is $q=0.75$. The shaded area indicates the
supersonic resonant cavity of length $L=20$, where the speed of sound
is $c_2=0.25$. The downstream supersonic sound speed is $c_3=0.5$.
Lower row: (c)--(d) Hawking (solid blue) and Andreev (solid red) spectra
of the BH solutions above.}
\end{figure}

The above configuration can be easily generalized to provide analytical
BH solutions by combining the three models of
Figure~\ref{fig:BHModels}, which yield piecewise homogeneous GP
equations, whose solutions are known. As a result, an analytical
resonant BH solution is typically given by either a gray soliton or a
homogeneous subsonic wave in the upstream region, by a homogeneous
supersonic plane wave in the downstream region, and by a cnoidal wave
in between; the GP wavefunction in each region is characterized by the
same global current $J$ and chemical\break potential $\mu$. 

A particularly simple model of resonant BH solution, represented in
Figure~\ref{fig:Resonant}b, was provided in Ref.~\cite{deNova2014}
by using a flat-profile piecewise configuration 
{\begin{equation}\label{eq:DoubleFlatProfile}
    g(x)n_u=\left\{ \begin{array}{cc}
1, & x< 0,\vspace*{2.5pt}\\
c_2^2, & 0\leq  x \leq L,\vspace*{3pt}\\
c_3^2, &  x > L,
\end{array}\right.
\end{equation}}\unskip
with $c_2,c_3<q$.


Regarding the Andreev and Hawking spectra, they are computed by solving
the scattering problem between the asymptotic upstream and downstream
channels, similar to that which was solved for the non-resonant BH
solutions. However, now there are two matchings, one with the upstream
region and one with the downstream region. Due to the periodic
character of a cnoidal wave, the corresponding BdG solutions take the
form of Bloch waves; analytical solutions can be obtained with the help
of mathematical tables (see for instance Ref.~\cite{Martone2021}). In
practice, a numerical integration of the time-independent BdG equations
for given frequency $\omega$ is quite efficient as these can be
recast as a simple $4\times 4$ linear system of first-order ordinary 
differential equations. In the case of the flat-profile resonant
configuration, this is not needed, as all the modes involved in each
region are the plane-wave spinors of
Equation~(\ref{eq:PlaneWaveSpinors}).

The Andreev and Hawking spectra for the BH solutions of
Figures~\ref{fig:Resonant}a,b are presented in
Figures~\ref{fig:Resonant}c,d, respectively. In the case of the double
delta barrier, Figure~\ref{fig:Resonant}c, we observe that apart from
the universal thermal $1/\omega$ peak of $|S_{Id2}(\omega)|^2$ at low
frequencies, there is a strongly non-thermal peak close to 
${\simeq}$0.6$\omega_{{\max}}$. This is because the two delta barriers behave as a
Fabry--Perot resonator for the anomalous scattering processes of the
Andreev and Hawking effects. We note that, as for the single delta
barrier, $|S_{d1d2}(\omega)|^2\ll|S_{ud2}(\omega)|^2$.

In the resonant flat-profile case, Figure~\ref{fig:Resonant}d, we
observe a similar trend, i.e., apart from the universal thermal peak at
low frequencies, there is a highly non-thermal peak at large
frequencies. However, unlike for the double delta barrier case, here
the Andreev signal is larger than the Hawking signal,
$|S_{d1d2}(\omega)|^2 > |S_{ud2}(\omega)|^2$. This highlights that
resonant structures may also be useful to enhance the Andreev effect,
which is suppressed with respect to the Hawking effect in typical
non-resonant configurations. It must be noted that, even though the
cavity is much longer than for the double delta barrier, the spectrum
only displays one peak. This is due to the smallness of the cutoff
frequency $\omega_{{\max}}$ for the flat-profile configuration,
resulting in low frequencies for the spectrum that are translated into
small wavevectors, whose inverse provides the typical length scale for
the occurrence of resonances.

\subsection{Black hole from an outcoupled condensate through an optical
lattice}

The resonant configurations discussed above displayed a single resonant
peak in the spectrum. Interestingly, the opposite limit of a long
cavity with many resonant peaks can be experimentally reproduced with
the help of an optical lattice, a major paradigm in AMO
physics~\cite{Orzel2001,Greiner2002,Bakr2009}. In particular,
Ref.~\cite{deNova2014a} provided a thorough numerical study of the
quasi-stationary BH resulting from the outcoupling of a condensate
through an optical lattice, which we proceed to discuss. 

A 1D optical lattice can be created from the interference of two
fixed-phase lasers of wavelength $\lambda$ whose wavevectors form an
angle $\theta$~\cite{Fabre2011}. The resulting potential can be written
as
{\begin{equation}\label{eq:OLPotential}
    V(x,t)=V(t)f(x-L)\cos^{2}\left[k_{\mathrm{L}}(x-L)\right]
\end{equation}}\unskip
with $k_{\mathrm{L}}=\rmpi/d$ and
$d=\lambda/\left[2\sin(\theta/2)\right]$ the lattice period. In the
above equation, $f(x)$ is a dimensionless function that characterizes
the global shape of the optical lattice, accounting for its finite size
in real experiments, while $V(t)$ represents its (possibly
time-dependent) amplitude. 

In our specific configuration, we assume that our condensate is
confined by a high-amplitude barrier placed at $x=0$, modeled by a
hard-wall boundary condition for the GP wavefunction $\Psi(0,t)=0$. The
position $L$ in Equation~(\ref{eq:OLPotential}) plays the role of the
approximate localization of the lattice, chosen here to be repulsive,
so the condensate is essentially confined between $0\leq x\lesssim L$.
The lattice amplitude is gradually lowered from $V_0$ to $V_{\infty}$
as
{\begin{equation}\label{eq:TDAmplitude}
    V(t)=\left\{ \begin{array}{cc}
V_0, & t\leq 0,\\
V_{\infty}+(V_{0}-V_{\infty})\mathrm{e}^{-t/\tau}, & t>0 ,
\end{array}\right.
\end{equation}}\unskip
where $\tau$ is the characteristic timescale of the process. This
causes the condensate to outcouple through the optical lattice from the
initial reservoir.

Quantitatively, the problem is described by the time-dependent GP equation
{\begin{equation}\label{eq:TDGPConfined}
\left[-\frac{\partial_{x}^{2}}{2}+V(x,t)+|\Psi(x,t)|^{2}\right]\Psi(x,t) = \mathrm{i}\partial_t\Psi(x,t),
\end{equation}}\unskip
where the initial condition $\Psi(x,0)=\Psi_0(x)$ is the wavefunction
describing the equilibrium condensate, which is solution of the
time-independent GP equation
{\begin{equation}\label{eq:TIGPConfined}
\left[-\frac{\partial_{x}^{2}}{2}+V(x,0)+|\Psi_0(x)|^{2}-\mu_{0}\right]\Psi_0(x) = 0,\quad \Psi_0(0)=0,
\end{equation}}\unskip
with the chemical potential $\mu_0$ determined by the normalization
condition (\ref{eq:Normalization}). This chemical potential also
defines some relevant physical scales: density $n_0\equiv \mu_{0}$,
length $\xi_0\equiv 1/\sqrt{\mu_{0}}$, velocity $c_{0}\equiv
\sqrt{\mu_{0}}$, time $t_0\equiv 1/\mu_{0}$, and temperature $T_0\equiv
\mu_0$. Typical orders of magnitude for $^{87}$Rb are $\xi_0\sim
0.1-1~\rmmu\mathrm{m}$, $t_0\sim 10^{-4}\ndash 10^{-3}~\mathrm{s}$,
$c_0\sim 0.1-1~\mathrm{mm}/\mathrm{s}$, and $T_0\sim 1\ndash 10
~\mathrm{nK}$. The lattice period satisfies $d>\lambda/2\gtrsim \xi_0$
while the lattice amplitude can essentially take any value; we choose
$V_0\gg \mu_0$ and $V_\infty\gtrsim \mu_0$, so the lattice goes from
providing a tight confinement to allow some leakage. On the other hand,
we require $\tau\gg t_0$, so the barrier lowering is adiabatic and does
not introduce further distortions in the condensate flow, and $L\gg
\xi_0$, so we are in the Thomas--Fermi regime where $\Psi_0(x)\simeq
\sqrt{n_0}$ is the bulk value (for $0\lesssim x \lesssim L$) of the
condensate amplitude. Therefore, we can regard $L$ as the approximate
size of the reservoir, containing $N\sim n_0 L$ particles.

The time evolution of $\Psi(x,t)$ is displayed in
Figure~\ref{fig:OLMF}, computed from numerical integration of
Equation~(\ref{eq:TDGPConfined}). In the first row, we consider an {\it
ideal} finite optical lattice, determined by an envelope
{\begin{equation}\label{eq:idealOL}
f(x)=\chi\left(\frac{x+\frac{d}{2}}{L_{\mathrm{lat}}}\right),
\end{equation}}\unskip
$\chi(x)$ being the characteristic function of the interval $[0,1]$.
Thus, a lattice with instantaneous uniform amplitude $V(t)$ extends
from $x_0=L-d/2$ to $x_1=x_0+L_{\mathrm{lat}}$, where the lattice
length is chosen such that it contains an integer number of periods
$n_{\rm osc}\sim~$10--50, $L_{\mathrm{lat}}\equiv n_{\rm osc}d$ and
$V(x_0)=V(x_1)=0$.


\begin{figure}[t!]
\includegraphics{fig04}
\caption{\label{fig:OLMF}Time evolution of an initially confined condensate which is
outcoupled through an optical lattice. Upper row: Ideal optical lattice
(\ref{eq:idealOL}) with $L\approx 125 \xi_0$, $d\approx 2.36\xi_0$, and
$n_{\rm osc}=30$. The lowering time is $\tau=500 t_0$. (a) Time-dependent 
profile of the sound speed, $c(x,t)/c_0$, with $c(x,t)=|\Psi(x,t)|$.
(b) Sound (solid blue) and flow (dashed red) velocity profiles of the
quasi-stationary regime, evaluated at the last snapshot of (a). The
initial sound speed profile is depicted in solid green (the initial
flow velocity is identically zero), while the optical lattice envelope
is shown as a black line. (c) Time evolution of the average chemical
potential $\bar{\mu}(t)$ (solid black) and its relative fluctuations
$\sigma(t)$ (solid red). The instantaneous conduction bands (energy
gaps) of the optical lattice are depicted as white (gray) bands. Lower
row: (d)--(f) Same as (a)--(c) but for a Gaussian optical lattice
(\ref{eq:actualpotential}) with $L\approx 1480 \xi_0$, $d\approx
1.73\xi_0$, and $\tilde{w}\approx 220.5\xi_0$.}
\end{figure}

In  Figure~\ref{fig:OLMF}a, we represent the time evolution of the
sound speed profile $c(x,t)$, proportional to the square root of the
local density, $c(x,t)=|\Psi(x,t)|=\sqrt{n(x,t)}$ (this choice improves
the visibility of the condensate outside the reservoir as compared to
using the density itself). After some transient times $t\sim 10^4 t_0$,
the condensate achieves a quasi-stationary regime in which it flows
through the lattice and eventually leaks outside. The sound and flow
velocity profiles in the quasi-stationary regime are shown in
Figure~\ref{fig:OLMF}b, where we observe that a BH configuration is
achieved, with the downstream supersonic region located outside the
lattice. 

We can quantify the degree of quasi-stationarity by defining a
\textit{local} chemical potential as
{\begin{equation}\label{eq:LocalChemicalPotential}
\mu(x,t)\equiv -\frac{1}{2}\frac{\partial_x^{2}\Psi(x,t)}{\Psi(x,t)}+V(x,t)+|\Psi(x,t)|^{2}.
\end{equation}}\unskip
For a stationary solution, $\mu(x,t)=\mu$ is real and constant. The
current is also constant and uniform for a 1D stationary solution.
However, this latter condition is impossible to fulfill strictly, since
the current is zero at $x=0$ due to the hard-wall boundary condition,
while the leaked downstream flow carries a non-zero flux. Hence, there
must be a current gradient, which, via the continuity equation, implies
a time-dependent density. This also implies a non-homogeneous
time-dependent complex chemical potential as
{\begin{equation}
    \partial_t \ln n=2\,\mathrm{Im}\,\mu .
\end{equation}}\unskip
Nevertheless, in practice, such dependence can become so weak that one
can neglect it, effectively achieving a quasi-stationary regime. This
regime should be characterized by a sufficiently uniform local chemical
potential $\mu(x,t)$, with small relative spatial fluctuations
$\sigma(t)$ around its instantaneous average value $\bar{\mu}(t)$,  
{\begin{eqnarray}
\bar{\mu}(t) & \equiv  & \frac{\int_{0}^{L_{g}}\mathrm{d}x~|\Psi(x,t)|^2\mu(x,t)}{\int_{0}^{L_{g}}\mathrm{d}x~|\Psi(x,t)|^2},\nonumber \\ 
\sigma(t) & \equiv  & \frac{1}{\bar{\mu}(t)}
\left[
\frac {\int_{0}^{L_{g}}\mathrm{d}x~|\Psi(x,t)|^2|\mu(x,t)-\bar{\mu}
(t)|^{2}} {\int_{0}^{L_{g}}\mathrm{d}x~|\Psi(x,t)|^2}
\right]^{\tfrac{1}{2}},
\label{eq:AverageChemicalPotential}
\end{eqnarray}}\unskip
where $L_g$ is the total length considered for the average. Typically,
$L_g$ is chosen well inside the downstream region, and the results are
quite insensitive to its specific value.


The time evolution of the real part of the average chemical potential
$\bar{\mu}(t)$ and its relative fluctuations $\sigma(t)$ is depicted in
Figure~\ref{fig:OLMF}c (imaginary values can be neglected as
$\mathrm{Im}\,\bar{\mu}\sim 10^{-6}\ndash 10^{-7}\mu_0$). In order to
understand their relation with $V(t)$, we also represent the
instantaneous band structure of the lattice, computed using the linear
Schr\"odinger equation since the non-linear interacting term is
negligible within the lattice due to the smallness of the density.
Specifically, the lowest conduction band is placed between $E_{0}(t)$
and $E_{1}(t)$, with a width $\Delta_c(t)=E_{1}(t)-E_{0}(t)$, where
dimensional arguments show that
{\begin{equation}\label{eq:BandStructure}
    E_{0,1}(t)=E_{\mathrm{R}}\,  F_{0,1}[\zeta(t)],\quad \zeta(t)\equiv \frac{V(t)}{16 E_{\mathrm{R}}},
\end{equation}}\unskip
with $F_{0,1}$ increasing functions of the dimensionless parameter $\zeta$ that can be computed perturbatively, 
{\begin{eqnarray}
\begin{array}{rcl}
F_0(\zeta)&=& 8\zeta-8\zeta^2+O(\zeta^4),\vspace*{2.5pt}\\
F_1(\zeta)&=&1+4\zeta-2\zeta^2+O(\zeta^4),
\end{array}
\end{eqnarray}}\unskip
and $E_{\mathrm{R}}\equiv k^2_{\mathrm{L}}/2$ the recoil energy of the lattice. In general,
the solutions to the 1D Schr\"odinger equation with a sinusoidal
potential can be obtained in terms of Mathieu functions. In practice,
the functions $F_{0,1}$ can be easily evaluated
numerically~\cite{deNova2014a}. The resulting time-dependent conduction
band (energy gap) is depicted as a white (gray) band. As we can see,
the condensate smoothly approaches the bottom of the asymptotic
conduction band, since tunneling is exponentially suppressed for
$\mu_0<E_{0}$. In this regime of small leaking, the fluctuations of the
chemical potential can become extremely small, $\sigma(t)\sim 10^{-4}$,
ensuring a high-degree of quasi-stationarity. 

The formation of a quasi-stationary BH can be then easily understood
from energetic arguments: in the upstream region, where the reservoir
is placed, the flow velocity is negligible and the chemical potential
is merely due to interactions, i.e., $\bar{\mu}\simeq n_u$ and
$v_u\simeq 0$. Due to the conservation of the chemical potential as
well as the small density there, in the downstream region the
condensate flows with a high velocity $v_d\sim \sqrt{2\bar{\mu}}\gg
c_d$, becoming supersonic. By continuity, this implies that there must
be an acoustic horizon somewhere within the lattice. In the bulk of the
lattice, the wavefunction is a Bloch wave, which is preferred to be
subsonic due to its energetic stability~\cite{Wu2003}. Thus, the
acoustic horizon must be placed at the right edge of the lattice, as
seen in Figure~\ref{fig:OLMF}b.

In the second row of Figure~\ref{fig:OLMF}, we analyze a more realistic
Gaussian envelope
{\begin{equation}\label{eq:actualpotential}
f(x)=\mathrm{e}^{-2\tfrac{x^2}{\tilde{w}^2}},
\end{equation}}\unskip
with $\tilde{w}$ the effective beam waist, which plays a similar role
to $L_{\mathrm{lat}}$ for the ideal optical lattice. We require the
length hierarchy 
{\begin{equation}\label{eq:LengthHierarchy}
    d\ll \tilde{w}\ll L,
\end{equation}}\unskip
where the second condition is imposed in order to have a sufficiently
large and homogeneous condensate reservoir. The first condition is
satisfied for typical waists, and implies that the overall Gaussian
amplitude behaves as a spatially adiabatic envelope, so the potential
can be regarded locally as an {\it ideal} optical lattice with an
amplitude $V_{\mathrm{A}}(x,t)$~\cite{Santos1998a,Santos1999,Carusotto2000}:
{\begin{eqnarray}
\begin{array}{rcl}
V(x,t)&=&V_{\mathrm{A}}(x,t)\cos^{2}[k_{\mathrm{L}}(x-L)],\vspace*{2.5pt}\\
\displaystyle V_{\mathrm{A}}(x,t)&\equiv&\displaystyle V(t)\exp\left[-2\left(\frac{x-L}{\tilde{w}}\right)^2\right].
\end{array}
\end{eqnarray}}\unskip
We observe in Figures~\ref{fig:OLMF}d,e the same trends as for the
ideal lattice case, namely, a quasi-stationary BH solution is achieved
for sufficiently long times. Due to the hierarchy
(\ref{eq:LengthHierarchy}), we only depict the vicinity of the lattice
peak $L-2.5\tilde{w}\leq x\leq L+2.5\tilde{w}$ to better observe the
structure of the horizon; in turn, at this scale, the lattice structure
cannot be resolved and the oscillations of the flow and sound
velocities appear as broadened lines. In Figure~\ref{fig:OLMF}f, the
average chemical potential also descends towards the bottom of the
conduction band, achieving a highly quasi-stationary regime where the
relative fluctuations $\sigma(t)$ are also insignificant,
$\sigma(t)\sim 10^{-4}$. The band structure is now evaluated at the
lattice peak $x=L$, where the local conduction band is determined by
the energies $E_{0,1}(x,t)$ obtained by taking
$\zeta(x,t)=V_A(x,t)/16E_{\mathrm{R}}$ in Equation~(\ref{eq:BandStructure}). The
fact that the lattice maximum is the transmission bottleneck can be
understood from the increasing character of the asymptotic energies
$E_{0,1}(x)\equiv E_{0,1}(x,\infty)$ with respect to the lattice
envelope $V_A(x)\equiv V_A(x,\infty)$. Thus, in order to place the
chemical potential within the local conduction band across the whole
lattice, it must be satisfied
{\begin{equation}
    E_{0}(x)\leq  E_{0}(L)< \mu_0 < E_{\mathrm{R}}<E_{1}(x),
\end{equation}}\unskip
where we have used that $F_1(0)=1$, so $E_1(0)=E_{\mathrm{R}}$. This
implies that the asymptotic lattice amplitude $V_\infty$ must be below
a certain critical value $V_{\mathrm{c}}$ so the condition
$E_{0}(L)<E_{\mathrm{R}}$ is met, which can be derived from
{\begin{equation}
  F_0\left(\frac{V_\infty}{16E_{\mathrm{R}}}\right)<1\Longrightarrow  V_\infty<V_{\mathrm{c}},\quad V_{\mathrm{c}}\approx 2.33 E_{\mathrm{R}}.
\end{equation}}\unskip

Another remarkable feature is that the acoustic horizon is now placed
exactly at the lattice peak (vertical dashed black line in
Figure~\ref{fig:OLMF}e). This is no coincidence and a detailed local
lattice calculation explicitly proves that the acoustic horizon must be
placed at the extremes of the lattice envelope~\cite{deNova2014a}. This
is in agreement with another result derived for a smooth potential
within the hydrodynamic approximation, stating that the acoustic
horizon must be located at the potential maximum~\cite{Giovanazzi2004}.

We summarize now the computation of the scattering matrix for the above
quasi-stationary BH solutions, where the interested reader can consult
Ref.~\cite{deNova2017b} for more details. From
Equation~(\ref{eq:LocalChemicalPotential}), and by invoking the
time-dependent GP equation~(\ref{eq:TDGPConfined}), it is easily shown
that $\mathrm{i}\partial_t \ln \Psi(x,t)=\mu(x,t)$. Hence,  
{\begin{equation}\label{eq:quasistationarywavefunction}
\Psi(x,t)=\Psi(x,t_s)\mathrm{e}^{-\mathrm{i}\int_{t_s}^{t}\mathrm{d}t'\mu(x,t')} \equiv\Psi_{\infty}(x,t)\mathrm{e}^{-\mathrm{i}{\rm Re}\,\bar{\mu}(t-t_s)},
\end{equation}}\unskip
where the rightmost term provides a definition for
$\Psi_{\infty}(x,t)$. If one chooses $t_s$ well inside the
quasi-stationary regime, when $\mu(x,t)\simeq \bar{\mu}\simeq {\rm
Re}\,\bar{\mu}$, then one can approximate $\Psi_{\infty}(x,t)\simeq
\Psi(x,t_s)=\Psi_{\infty}(x,t_s)\equiv \Psi_{\infty}(x)$. Thus, we can
work with a fully stationary BH solution whose chemical potential is
$\mu={\rm Re}\,\bar{\mu}$, and compute the Andreev and Hawking spectra
from the $S$-matrix by solving the associated BdG scattering problem,
where the external potential is the asymptotically stationary optical
lattice, $V(x)=V_\infty(x)\equiv V(x,t=\infty)$. Specifically, the
time-independent BdG equations for given $\omega$ are numerically
integrated and eventually matched with the corresponding scattering
channels at the asymptotic homogeneous subsonic (the upstream bulk of
the condensate reservoir) and supersonic (the downstream leaking flux)
regions. 

However, a major difficulty arises for the computation in the Gaussian
lattice, since its large size makes that exponentially growing modes,
corresponding to local Bloch waves with complex wavevector, explode
above the propagating modes, making the matching equations singular
within computer accuracy. This is known in general as the $\Omega d$
problem~\cite{Lowe1995}, emerging in a wide range of scenarios, ranging
from the propagation of ultrasonic and electromagnetic waves in
multilayered media~\cite{Lowe1995,Pernas2014} to Anderson
localization~\cite{Slevin2004}. Possible methods to deal with this
specific $\Omega d$ problem in the BdG context are the Global Matrix
method~\cite{Lowe1995} or QR decomposition~\cite{Slevin2004}.

As expected, the resulting Hawking and Andreev spectra display a highly
non-thermal structure. In particular, since the energy contribution 
from the speed of sound is negligible as compared to that from 
the lattice potential, one can analyze the
spectrum in terms of the underlying Schr\"odinger problem,
characterized by the band structure shown in Figures~\ref{fig:OLMF}c,f.
Thus, two main qualitative regimes arise: $\omega_{{\max}}<\Delta_c$,
where the spectrum is cut at its natural cutoff frequency, and
$\omega_{{\max}}>\Delta_c$, where the spectrum is abruptly cut by the
upper end of the conduction band. Hence, an optical lattice can behave
as a low-pass filter of Andreev--Hawking radiation, which may have
potential applications in quantum transport and
atomtronics~\cite{Amico2021}.

\section{Quantum Andreev--Hawking radiation}\label{sec:QuantumHR}

Although resonant configurations are highly non-thermal as shown in the
previous sections, this does not automatically imply that the
observation of a non-thermal Hawking spectrum is a signature of the
Hawking effect. This is because the BdG equations describe at the same
time the linear dynamics of both perturbations of the GP wavefunction
and quantum fluctuations of the field operator around the mean-field
expectation value. Moreover, the quantum state of the system may be a
highly thermal state. Thus, the observed non-resonant spectrum can
result from the coherent or thermal stimulation of Hawking radiation by
a classical source, instead of arising from a zero-point quantum
origin. The same applies to the Andreev effect. In this section, we
discuss how to unambiguously signal the genuine quantum character of
the Andreev and Hawking effects, as opposed to classical stimulation,
using different types of quantum correlations.

\subsection{Lessons from quantum optics}\label{subsec:QuantumOptics}

A major front of the quantum-classical frontier is present in the field
of quantum optics, from where we can borrow a number of concepts and
techniques; the interested reader is referred to
Ref.~\cite{Walls2008} for a pedagogical introduction to quantum
optics and a thorough discussion of a number of fundamental quantum
topics in that context. 

For simplicity, we begin by considering a single bosonic mode (as can
be that of a photon) whose amplitude is given by an annihilation
operator $\hat{a}$. The corresponding Hilbert space is the Fock space
spanned by the number states $\ket{n}$,
$\hat{a}^\dagger\hat{a}\ket{n}=n\ket{n}$, $n=0,1,2\ldots.$ A coherent
state is an eigenstate of the annihilation operator,
$\hat{a}\ket{\alpha}=\alpha\ket{\alpha},~\alpha\in\mathbb{C}$, and can
be expressed as
{\begin{equation}
\ket{\alpha}=\mathrm{e}^{-\frac{|\alpha|^2}{2}}\sum^{\infty}_{n=0}\frac{\alpha^n}{\sqrt{n!}}\ket{n}=D(\alpha)\ket{0},\quad D(\alpha)\equiv \mathrm{e}^{\alpha \hat{a}^\dagger-\alpha^*\hat{a}},
\end{equation}}\unskip
where $D(\alpha)$ is the displacement operator, which is unitary,
$D^\dagger(\alpha)=D^{-1}(\alpha)=D(-\alpha)$. The coherent states form
an overcomplete basis of the Hilbert space since
{\begin{equation}\label{eq:Overcomplete}
    |\braket{\alpha|\beta}|^2=\mathrm{e}^{-|\alpha-\beta|^2}\neq 0,\quad \int\frac{\mathrm{d}^2\alpha}{\rmpi}\ket{\alpha}\bra{\alpha}=\sum^{\infty}_{n=0}\ket{n}\bra{n}=1,
\end{equation}}\unskip
with
$\mathrm{d}^2\alpha\equiv\mathrm{d}\alpha_x\mathrm{d}\alpha_y,~\alpha=\alpha_x+\mathrm{i}\alpha_y$. 
Another relevant class of quantum states are the squeezed states
{\begin{equation}\label{eq:Squeezing}
    \ket{\varepsilon}=S(\varepsilon)\ket{0},\quad S(\varepsilon)\equiv \mathrm{e}^{\frac{\varepsilon^*\hat{a}^2-\varepsilon(\hat{a}^\dagger)^2}{2}},~ \varepsilon\equiv r \mathrm{e}^{\mathrm{i}2\theta},
\end{equation}}\unskip
where $S(\varepsilon)$ is the squeezing operator, also unitary as
$S^\dagger(\varepsilon)=S^{-1}(\varepsilon)=S(-\varepsilon)$. These
squeezed states are the vacuum of the annihilation operator $\hat{b}$,
$\hat{b}\ket{\varepsilon}=0$, arising from the Bogoliubov
transformation 
{\begin{equation}
\hat{b}=S(\varepsilon)\hat{a}S^\dagger(\varepsilon)=\hat{a}\cosh r+\hat{a}^\dagger \mathrm{e}^{\mathrm{i}2\theta}\sinh r.
\end{equation}}\unskip
By noticing that
$\hat{a}^2,(\hat{a}^\dagger)^2,\hat{a}^\dagger\hat{a}+\hat{a}\hat{a}^\dagger$ 
form a closed Lie algebra (actually, they form a representation of
$\mathfrak{su}(1,1)$), one can rewrite the squeezing operator in normal
order as
{\begin{equation}\label{eq:SqueezingOperatorNormal}
    S(\varepsilon)=\frac{1}{\sqrt{\cosh r}}\mathrm{e}^{-\frac{g}{2}(\hat{a}^\dagger)^2}\mathrm{e}^{f\hat{a}^\dagger\hat{a}}\mathrm{e}^{\frac{g^*}{2}\hat{a}^2},
\end{equation}}\unskip
with $g=\mathrm{e}^{\mathrm{i}2\theta}\tanh r$ and $f=-\ln \cosh r$. This allows to
readily express the squeezed states as
{\begin{equation}\label{eq:SqueezedFock}
    \ket{\varepsilon}=\frac{1}{\sqrt{\cosh r}}\sum^\infty_{n=0}(-1)^n\mathrm{e}^{\mathrm{i}2n\theta}\frac{\sqrt{2n!}\tanh^n r}{2^n\cdot n!}\ket{2n}.
\end{equation}}\unskip
The most general quantum state in this Hilbert space is described by a
density matrix $\hat{\rho}$ of the form
{\begin{equation}\label{eq:GeneralQuantumStateFock}
    \hat{\rho}=\sum^\infty_{n,m=0} \rho_{nm}\ket{n}\bra{m},
\end{equation}}\unskip
which is completely determined by its characteristic function
{\begin{equation}
    \chi(\eta)\equiv \braket{\mathrm{e}^{\eta\hat{a}^\dagger-\eta^*\hat{a}}}=\mathrm{Tr}[\mathrm{e}^{\eta\hat{a}^\dagger-\eta^*\hat{a}}\hat{\rho}].
\end{equation}}\unskip
One can also work with its normal and anti-normal versions
{\begin{eqnarray}\label{eq:CharacteristicParanormal}
\begin{array}{rcl}
      \chi_N(\eta)&\equiv& \braket{\mathrm{e}^{\eta\hat{a}^\dagger}\mathrm{e}^{-\eta^*\hat{a}}}, \vspace*{2.5pt}\\
      \chi_A(\eta)&\equiv& \braket{\mathrm{e}^{-\eta^*\hat{a}}\mathrm{e}^{\eta\hat{a}^\dagger}}.
\end{array}
\end{eqnarray}}\unskip

Alternative representations to the Fock expansion
(\ref{eq:GeneralQuantumStateFock}) are provided by the distributions
resulting from the Fourier transform of the characteristic functions:
{\begin{eqnarray}\label{eq:ProbabilitiesQuantumOptics}
\begin{array}{rcl}
P(\alpha)&\equiv&\displaystyle \int\frac{\mathrm{d}^2\eta}{\rmpi^2}\mathrm{e}^{(\eta^*\alpha-\eta\alpha^*)}\chi_N(\eta), \vspace*{2.5pt}\\
Q(\alpha)&\equiv&\displaystyle \int\frac{\mathrm{d}^2\eta}{\rmpi^2}\mathrm{e}^{(\eta^*\alpha-\eta\alpha^*)}\chi_A(\eta), \vspace*{2.5pt}\\
W(\alpha)&\equiv&\displaystyle \int\frac{\mathrm{d}^2\eta}{\rmpi^2}~\mathrm{e}^{(\eta^*\alpha-\eta\alpha^*)}\chi(\eta).
\end{array}
\end{eqnarray}}\unskip
All of them are quasi-probability distributions, properly normalized,
{\begin{equation}
\int\mathrm{d}^2\alpha~P(\alpha)=\int\mathrm{d}^2\alpha~W(\alpha)=\int\mathrm{d}^2\alpha~Q(\alpha)=1,
\end{equation}}\unskip
but do not describe disjoint events since coherent states are not
orthogonal, and can even take negative values, opening the door to
genuine non-classical behavior. 

The Glauber--Sudarshan $P$ function is equivalent to a diagonal
representation in the coherent basis,
{\begin{equation}
\hat{\rho}=\int\mathrm{d}^2\alpha~P(\alpha)\ket{\alpha}\bra{\alpha},
\end{equation}}\unskip
since its momenta provides the normal-ordered expectation values
{\begin{equation}
\braket{(\hat{a}^\dagger)^n\hat{a}^m}=\int\mathrm{d}^2\alpha~(\alpha^*)^n\alpha^m P(\alpha),
\end{equation}}\unskip
which are those typically characterizing correlation functions. This is
where the crucial role of the $P$ function in the understanding of the
classical-quantum frontier emerges: if it is non-negative,
$P(\alpha)\geq 0$, we can understand these correlations as statistical
averages over a continuous classical variable $\alpha$ with a
probability distribution given precisely by $P(\alpha)$. Thus, quantum
states with a non-negative $P$ function can be regarded as classical,
admitting a conventional probabilistic description in terms of a
stochastic complex amplitude. Examples of classical states are coherent
states, chaotic states, as well as quantum thermal states. On the other
hand, number states and squeezed states are intrinsically non-classical
as they do not even have a well-defined $P$-representation.

The $Q$-function yields the anti-normal expectation values
{\begin{equation}
\braket{\hat{a}^m(\hat{a}^\dagger)^n}=\int\mathrm{d}^2\alpha~(\alpha^*)^n\alpha^mQ(\alpha),
\end{equation}}\unskip
and it is easily evaluated by inserting the identity representation
(\ref{eq:Overcomplete}) in
Equation~(\ref{eq:CharacteristicParanormal}),
{\begin{equation}
Q(\alpha)=\frac{\braket{\alpha|\hat{\rho}|\alpha}}{\rmpi},
\end{equation}}\unskip
so it is non-negative and bounded, $0\leq Q(\alpha)\leq 1/\rmpi$.

Finally, the Wigner function $W(\alpha)$ characterizes symmetric
expectation values such as
{\begin{equation}
\frac{1}{2}\braket{\hat{a}^\dagger\hat{a}+\hat{a}\hat{a}^\dagger}=\int\mathrm{d}^2\alpha~|\alpha|^2\,W(\alpha).
\end{equation}}\unskip
Further insight on the physical meaning of the Wigner function is
obtained when working with the usual coordinate-momentum representation
{\begin{equation}
    \hat{a}=\frac{\hat{q}+\mathrm{i}\hat{p}}{\sqrt{2}},\quad \alpha=\frac{q+\mathrm{i}p}{\sqrt{2}},
\end{equation}}\unskip
whose eigenstates are labeled as $\ket{q}$, $\ket{p}$, respectively. After
proper normalization, the Wigner distribution in phase space reads
{\begin{equation}\label{eq:WignerPosition}
W(q,p)=\frac{1}{2\rmpi}\int\mathrm{d}q'~{\left\langle q-\frac{q'}{2}\right|\hat{\rho}\left|q+\frac{q'}{2}\right\rangle}\mathrm{e}^{\mathrm{i}pq'}.
\end{equation}}\unskip
Its marginal distributions
{\begin{eqnarray} 
\begin{array}{rcl}
W(q)&=&\displaystyle\int\mathrm{d}p~W(q,p)=\braket{q|\hat{\rho}|q}\geq 0,\vspace*{2.5pt}\\
W(p)&=&\displaystyle\int\mathrm{d}q~W(q,p)=\braket{p|\hat{\rho}|p}\geq 0
\end{array}
\end{eqnarray}}\unskip
are the spatial and momentum distributions of the quantum state. As a
result, the Wigner function represents a quantum version of the
classical Boltzmann distribution function. However, it is not
necessarily positive, and the presence of negative values $W(q,p)<0$ is
hence another quantum signature. Indeed, the negativity of the Wigner
function is a stronger condition than that of the Glauber--Sudarshan
function as they are related through the convolution
{\begin{equation}
    W(\alpha)=\frac{2}{\rmpi}\int\mathrm{d}^2\beta~\mathrm{e}^{-2|\alpha-\beta|^2}P(\beta),
\end{equation}}\unskip
derived by noting that
$\chi(\eta)=\mathrm{e}^{-\sfrac{|\eta|^2}{2}}\chi_N(\eta)$. For instance,
squeezed states do have a positive Wigner representation, while number
states do not. An interesting approach to quantum optics from phase
space can be found in Ref.~\cite{Schleich2001}. 

All the above concepts can be straightforwardly extended to
multipartite Hilbert spaces describing an ensemble of bosonic modes. Of
particular interest is the case of bipartite systems composed by two
modes, labeled as $i,j$, whose corresponding annihilation operators are
$\hat{a}_{i,j}$. Their Hilbert space is spanned by the Fock product
states $\ket{n_in_j}\equiv \ket{n_i}\otimes \ket{n_j}$, and accordingly
coherent states and quasi-probability distributions now have two
complex arguments $(\alpha_i,\alpha_j)$. For example, the
$P$-representation reads
{\begin{equation}\label{eq:PRepresentationBi}
\hat{\rho}=\int\mathrm{d}^2\alpha_i\mathrm{d}^2\alpha_j~P(\alpha_i,\alpha_j)\ket{\alpha_i\alpha_j}\bra{\alpha_i\alpha_j}.
\end{equation}}\unskip

Nevertheless, there are genuine bipartite quantum states which cannot
be expressed as a product of monomode states. One example is the
two-mode squeezed state
{\begin{equation}\label{eq:TwoModeSqueezed}
     \ket{\varepsilon}_{ij}\equiv U(\varepsilon)\ket{00},\quad U(\varepsilon)=\mathrm{e}^{\varepsilon^*\hat{a}_i\hat{a}_j-\varepsilon\hat{a}^\dagger_i\hat{a}^\dagger_j}.
\end{equation}}\unskip
In the same fashion of Equations (\ref{eq:SqueezingOperatorNormal}),
(\ref{eq:SqueezedFock}), it is shown that
{\begin{equation}
    \ket{\varepsilon}_{ij}=\frac{1}{\cosh r}\sum^\infty_{n=0}(-1)^n\mathrm{e}^{\mathrm{i}2n\theta}\tanh^n r\ket{nn}.
\end{equation}}\unskip
Remarkably, the reduced state
{\begin{equation}\label{eq:ThermalSqueezed}
\rho_i=\mathrm{Tr}_j(\ket{\varepsilon}_{ij} {\,}_{ij}{\bra{\varepsilon}})=\frac{1}{\cosh^2r}\sum^\infty_{n_i=0}(\tanh r)^{2n_i}\ket{n_i}  {\bra{n_i}}, 
\end{equation}}\unskip
is then a thermal state, whose equivalent temperature for a mode with
energy $\epsilon$ is obtained by $\mathrm{e}^{-\beta \epsilon}=\tanh^2r$. The
two-mode squeezed state is a non-classical state and it is of paramount
importance in quantum optics because it describes the time-evolution of
the so-called non-degenerate parametric amplifier, where one of the
modes is designed as the signal and the other as the
idler~\cite{Walls2008}.

\subsection{Cauchy--Schwarz violation}
Experimentally, the quantumness of a bipartite state, such as the
two-mode squeezed state (\ref{eq:TwoModeSqueezed}), can be
characterized through the measurement of the first 
{\begin{eqnarray}\label{eq:gdef}
g_{ij}\equiv \braket{\hat{a}_{i}^{\dagger}\hat{a}_{j}},\quad c_{ij}\equiv \braket{\hat{a}_{i}\hat{a}_{j}} ,
\end{eqnarray}}\unskip
and second-order correlation functions
{\begin{equation}\label{eq:GammaDef}
\Gamma_{ij}\equiv\braket{\hat{a}_{i}^{\dagger}
\hat{a}_{j}^{\dagger}\hat{a}_{j}\hat{a}_{i}}\geq 0\,.
\end{equation}}\unskip
Since they are normal-ordered expectation values, they can be computed
from the $P$-representation (\ref{eq:PRepresentationBi}).
Interestingly, for classical states, $P(\alpha_i,\alpha_j)\geq 0$, and
we can write the averages as a scalar product
{\begin{equation}
g_{ij}=\int\mathrm{d}^2\alpha~P(\alpha_i,\alpha_j)\alpha^*_i\alpha_j\equiv (\alpha_i,\alpha_j)_C,
\end{equation}}\unskip
with $c_{ij}=(\alpha^*_i,\alpha_j)_C$ and
$\Gamma_{ij}=(|\alpha_i|^2,|\alpha_j|^2)_C$. By invoking the
Cauchy--Schwarz (CS) inequality
{\begin{equation}
    |(\alpha_i,\alpha_j)_C|\leq \sqrt{(\alpha_i,\alpha_i)_C(\alpha_j,\alpha_j)_C},
\end{equation}}\unskip
one can prove that classical states obey the inequalities 
{\begin{eqnarray}\label{eq:ClassicalIneq}
 |g_{ij}|^2&\leq& g_{ii}g_{jj},\nonumber\\
    |c_{ij}|^2&\leq& g_{ii}g_{jj}, \nonumber\\
    \Gamma_{ij}&\leq& \sqrt{\Gamma_{ii}\Gamma_{jj}}.
\end{eqnarray}}\unskip
The violation of any of the above \textit{classical} CS inequalities
requires a negative-valued $P$ function, a genuine signature of
quantumness. The first violation of a CS inequality in photons was
observed in 1974~\cite{Clauser1974}. In condensates, CS violation has
been also observed in 2012~\cite{Kheruntsyan2012}.

Actual mathematical CS inequalities that are never violated can be
proven for quantum operators, which we now review along the lines of
the enlightening discussion from Adamek,\break 
Busch and Parentani~\cite{Adamek2013}. For two operators $\hat{A},\hat{B}$, one
can associate a scalar product to a quantum state $\hat{\rho}$ as
{\begin{equation}\label{eq:ScalarProductOperator}
(\hat{A},\hat{B})_Q\equiv\braket{\hat{A}^{\dagger}\hat{B}}=\mathrm{Tr}[\hat{A}^{\dagger}\hat{B}\hat{\rho}],
\end{equation}}\unskip
which satisfies the usual properties of a scalar product, including the
\textit{quantum} CS inequality
{\begin{equation}\label{eq:CSmathematical}
|\braket{\hat{A}^{\dagger}\hat{B}}|^2\leq \braket{\hat{A}^{\dagger}\hat{A}}\braket{\hat{B}^{\dagger}\hat{B}}.
\end{equation}}\unskip
Notice that the only assumption here is that $\hat{\rho}$ is a
\textit{physical} quantum state, specifically, a non-negative operator,
and this is satisfied by definition. We can now derive the quantum
versions of the CS inequalities (\ref{eq:ClassicalIneq}), which are
then strict mathematical inequalities. For instance, by substituting
$\hat{A}=\hat{a}_{i}$ and $\hat{B}=\hat{a}_{j}$ in
Equation~(\ref{eq:CSmathematical}), we find
{\begin{equation}\label{eq:CSgijimpossible}
|g_{ij}|^2=|\braket{\hat{a}^{\dagger}_{i}\hat{a}_{j}}|^2
\leq\braket{\hat{a}^{\dagger}_{i}\hat{a}_{i}}\braket{\hat{a}^{\dagger}_{j}\hat{a}_{j}}=g_{ii}g_{jj}.
\end{equation}}\unskip
This is the same CS inequality as in the classical case, so it is
always verified. However, taking $\hat{A}=\hat{a}^{\dagger}_{i}$ and
$\hat{B}=\hat{a}_{j}$ yields
{\begin{equation}\label{eq:CScijpos}
|c_{ij}|^2=|\braket{\hat{a}_{i}\hat{a}_{j}}|^2
\leq\braket{\hat{a}_{i}\hat{a}^{\dagger}_{i}}
\braket{\hat{a}^{\dagger}_{j}\hat{a}_{j}}=(g_{ii}+1)g_{jj} .
\end{equation}}\unskip
Interestingly, the \textit{quantum} CS inequality leaves the
possibility of violating the \textit{classical} CS inequality
$|c_{ij}|^2\leq g_{ii}g_{jj}$. In order to quantify this violation, we
define the CS witness
{\begin{eqnarray}\label{eq:CSviolation2}
\Delta_{ij}\equiv |c_{ij}|^2-g_{ii}g_{jj},
\end{eqnarray}}\unskip
and denote the condition $\Delta_{ij}>0$ as {\it quadratic} CS
violation. For the second-order correlation function, the associated
quantum CS inequality is obtained by setting
$\hat{A}=\hat{a}^{\dagger}_{i}\hat{a}_{i}$ and
$\hat{B}=\hat{a}^{\dagger}_{j}\hat{a}_{j}$,
{\vspace*{2pt}}
{\begin{eqnarray}\label{eq:CSGammaijpos}
 |\Gamma_{ij}|^2=|\braket{\hat{a}^{\dagger}_{i}\hat{a}_{i}\hat{a}^{\dagger}_{j}\hat{a}_{j}}|^2
\leq\braket{\hat{a}^{\dagger}_{i}\hat{a}_{i}\hat{a}^{\dagger}_{i}\hat{a}_{i}}\braket{\hat{a}^{\dagger}_{j}\hat{a}_{j}\hat{a}^{\dagger}_{j}\hat{a}_{j}}
=(\Gamma_{ii}+g_{ii})(\Gamma_{jj}+g_{jj}),
\end{eqnarray}}\unskip
which also leaves the possibility of violating the classical CS
inequality $|\Gamma_{ij}|^2\leq\Gamma_{ii}\Gamma_{jj}$, as quantified
by the CS witness
{\begin{equation}\label{eq:CSviolation4}
\Theta_{ij}\equiv \Gamma_{ij}-\sqrt{\Gamma_{ii}\Gamma_{jj}}.
\end{equation}}\unskip
In analogy to the quadratic CS violation, we refer to the condition
$\Theta_{ij}>0$ as {\it quartic} CS violation. Remarkably, the origin
of the violation of classical CS inequalities can be pin-pointed to the
non-commutativity of quantum operators, a property present at the very
core of quantum\break mechanics.

{\vspace*{2pt}}

\subsection{Entanglement}

{\vspace*{2pt}}

Entanglement is perhaps the most genuine quantum feature. It has been
observed in a wide variety of systems as different as
photons~\cite{Aspect1982}, neutrinos~\cite{Formaggio2016},
quarks~\cite{ATLAS2024,CMS2024}, mesons~\cite{Go2007},
atoms~\cite{Hagley1997}, molecules~\cite{Bao2023,Holland2023},
superconductors~\cite{Steffen2006}, nitrogen-vacancy centers in
diamond~\cite{Pfaff2013}, and even macroscopic diamond
itself~\cite{Lee2011}. In general, entanglement is defined as the
non-separability of the quantum state of a system~\cite{Werner1989}. In
turn, a quantum state in a bipartite Hilbert space is said to be
separable \textit{iff} it can be written as a convex sum of product
states,
{\vspace*{2pt}}
{\begin{equation}\label{eq:separability}
\hat{\rho}=\sum_n p_n\hat{\rho}^{(i)}_n\otimes\hat{\rho}^{(j)}_n,\quad  \sum_n p_n=1,\quad p_n\geq 0,
\end{equation}}\unskip
where $\hat{\rho}^{(i),(j)}_n$ are states within the Hilbert subspaces
associated to the $i,j$ modes, respectively. Classical states are
separable, as directly seen from Equation~(\ref{eq:PRepresentationBi}).

In order to characterize entanglement, we make use of the generalized
Peres--Horodecki (GPH) criterion~\cite{Simon2000}, which extends the
celebrated Peres--Horodecki (PH)
criterion~\cite{Peres1996,Horodecki1997} to continuous systems. The PH
criterion results from the fact that, if $\hat{\rho}$ is separable, its
partial transpose $\hat{\rho}_{t}$ with respect to one of the
subsystems is also a physical density matrix and, in particular, a
non-negative operator. Thus, the PH criterion states that, if
$\hat{\rho}_{t}$ is not non-negative, then $\hat{\rho}$ is necessarily
entangled. For $2\times 2$ and $2\times 3$ systems, the PH criterion is
a necessary and sufficient condition for entanglement; in general, it
is only a sufficient condition.

In order to obtain the partial transpose $\hat{\rho}_t$ of a density
matrix $\hat{\rho}$, we make use of its Wigner function in phase space,
$W(X)$, computed through the analogous version for bipartite systems of
Equation~(\ref{eq:WignerPosition}), where we gather the phase-space
variables in a single vector $X\equiv [q_i,p_i,q_j,p_j]^{T}$. Without
loss of generality, we take the partial transpose with respect to the
subsystem $j$, which amounts to transpose the matrix elements of
$\hat{\rho}$ with respect to the Hilbert subspace of the mode $j$. It
is straightforward to show then that the Wigner distribution $W_t(X)$
associated to $\hat{\rho}_t$ is simply given by 
{\vspace*{2pt}}
{\begin{equation}\label{eq:Wignertransposed}
W_t(X)=W(\Lambda X),\quad \Lambda=\mathrm{diag}[1,1,1,-1].
\end{equation}}\unskip
Another way to put it is that the transposition operation amounts to a
time reversal transformation in the Wigner function. 

The effects of this seemingly innocuous transformation can become
critical, as revealed when evaluating the uncertainties of the
phase-space operators
{\vspace*{2pt}}
{\begin{equation}
\hat{Y}=\sum_{\alpha=1}^4u^\alpha\Delta\hat{X}_\alpha,
\end{equation}}\unskip\noindent
where  $\hat{X}\equiv [\hat{q}_i,\hat{p}_i,\hat{q}_j,\hat{p}_j]^{T}$ is
the quantum version of the phase-space vector $X$, $\Delta
\hat{X}\equiv \hat{X}-\braket{\hat{X}}$, and $u^\alpha$ are the
components of an arbitrary four-dimensional complex vector $u$. Due to
the positiveness of the scalar product
(\ref{eq:ScalarProductOperator}), 
$\braket{\hat{Y}^{\dagger}\hat{Y}}\geq 0$, which in compact vector\break
notation reads
{\begin{equation}\label{eq:Uncer}
u^{\dagger}Mu\geq0,\quad M_{\alpha\beta}=\braket{\Delta\hat{X}_{\alpha}\Delta\hat{X}_{\beta}}.
\end{equation}}\unskip
This is an alternative expression of the uncertainty principle, holding
for any complex vector $u$, which implies that $M$ must be a
non-negative matrix, $M\geq 0$. 

For its computation, we separate the matrix $M$ into its symmetric and
antisymmetric part as
{\begin{eqnarray}\label{eq:Uncerx}
\begin{array}{rcl}
 M&=&V+\mathrm{i}\displaystyle\frac{L}{2},\vspace*{2.5pt}\\
V_{\alpha\beta}&=&\displaystyle\frac{\braket{\{\Delta\hat{X}_{\alpha},\Delta\hat{X}_{\beta}\}}}{2},\vspace*{2.5pt}\\
  \mathrm{i}L_{\alpha\beta}&=&\displaystyle\braket{[\Delta\hat{X}_{\alpha},\Delta\hat{X}_{\beta}]}
  =\braket{[\hat{X}_{\alpha},\hat{X}_{\beta}]}.
\end{array}
\end{eqnarray}}\unskip
where $\{\ldots\}$ is the anticommutator. The matrix $V$ is the
symmetric covariance matrix and, since it contains symmetric
expectation values, it is readily evaluated with the help of the Wigner
distribution,
{\begin{eqnarray}\label{eq:SymmetricCovariances}
\begin{array}{rcl}
\displaystyle\braket{X_\alpha}&=&\displaystyle\int\mathrm{d}^4X~X_\alpha W(X),\vspace*{2.5pt}\\
\displaystyle \frac{\braket{\{\Delta\hat{X}_{\alpha},\Delta\hat{X}_{\beta}\}}}{2}&=&\displaystyle\int\mathrm{d}^4X~\Delta X_\alpha \Delta X_\beta W(X).
\end{array}
\end{eqnarray}}\unskip
On the other hand, the commutators between phase-space operators are
proportional to the identity, and thus $L$ is a $4\times 4$ matrix
independent of the state $\hat{\rho}$, 
{\begin{eqnarray}\label{eq:XCommutators}
L=\left[\begin{array}{@{}cc@{}} J & 0\\
0 & J \end{array}\right],\quad J=\left[\begin{array}{@{}cc@{}}
0 & 1\\
-1 & 0
\end{array}\right],
\end{eqnarray}}\unskip
$J$ being the symplectic matrix in two dimensions.

The above decomposition not only allows to evaluate $M$, but it also
provides a straightforward way to derive the corresponding uncertainty
principle for $\hat{\rho}_t$. Indeed, since its Wigner function
$W_t(X)$ satisfies Equation~(\ref{eq:Wignertransposed}), it is
immediate to see that the condition
$\braket{\hat{Y}^\dagger\hat{Y}}_t\equiv
\mathrm{Tr}[\hat{Y}^\dagger\hat{Y}\hat{\rho}_t]\geq 0$ is equivalent to
$M_t\geq 0$, with
{\begin{eqnarray}\label{eq:Uncertrans}
M_t&\equiv&V_t+\mathrm{i}\frac{L}{2},\quad V_t=\Lambda V \Lambda.
\end{eqnarray}}\unskip
We can finally formulate quantitatively the GPH criterion: if
$\hat{\rho}$ is separable, $\hat{\rho}_t$ must be a physical state
satisfying the uncertainty principle, which implies $M_t\geq 0$.
Therefore, if $M_t$ is not non-negative, the state is entangled. Notice
that, since $M$ is always non-negative as the original $\hat{\rho}$ is
a physical density matrix, by the Sylvester--Jacobi criterion, $M_t$ is
non-negative \textit{iff} $\det M_t\geq 0$. The conditions $\det M,\det
M_t\geq 0$ are respectively equivalent to
$\mathcal{P}^{\pm}_{ij}\geq0$,\break where
{\begin{eqnarray}\label{eq:GPHpm}
\mathcal{P}^{\pm}_{ij}\equiv\det A_i\det A_j+(\tfrac{1}{4}\mp \det C_{ij})^2
-\mathrm{tr}(JA_iJC_{ij}JA_jJC_{ij}^{T})-\tfrac{1}{4}(\det A_i+\det A_j)
\end{eqnarray}}\unskip
and the matrices $A_i,A_j,C_{ij}$ are the $2\times 2$ blocks forming
the covariance matrix $V$, 
{\begin{equation}\label{eq:WBlocks}
V=\left[\begin{array}{@{}cc@{}} A_{i} & C_{ij}\\
C^{T}_{ij} & A_{j} \end{array}\right].
\end{equation}}\unskip
We can put together both conditions by defining the GPH function
$\mathcal{P}_{ij}$ as
{\begin{eqnarray}\label{eq:GPH}
\mathcal{P}_{ij}\equiv\det A_i\det A_j+(\tfrac{1}{4}-|\det C_{ij}|)^2
- \mathrm{tr}(JA_iJC_{ij}JA_jJC_{ij}^{T})-\tfrac{1}{4}(\det A_i+\det A_j),
\end{eqnarray}}\unskip
where $\mathcal{P}_{ij}<0$ is a sufficient condition for entanglement.
This is the entanglement witness that results from the GPH criterion.
Notice that, whenever $\det C_{ij}\geq 0$, $\hat{\rho}$ is separable,
because then $\mathcal{P}_{ij}=\mathcal{P}^{+}_{ij}\geq0$, so only
states with $\det C_{ij}< 0$ can be entangled. 

In the usual case where the expectation values of the operators
$\braket{\hat{X}_\alpha}=0$ vanish, the matrices $A_{k},C_{ij}$ are
expressed in terms of the first-order correlation functions as
{\begin{eqnarray}\label{eq:WBlocksx}
\begin{array}{rcl}
 A_{k}&=&\displaystyle \left(g_{kk}+\tfrac{1}{2}\right)\mathbb{I}_2+\left[\begin{array}{@{}cc@{}}
\mathrm{Re}~c_{kk} & \mathrm{Im}~c_{kk}\\
\mathrm{Im}~c_{kk} & -\mathrm{Re}~c_{kk}
\end{array}\right],\quad k=i,j,\vspace*{2.5pt}\\
C_{ij}&=&\displaystyle 
\left[\begin{array}{@{}cc@{}}
\mathrm{Re}(g_{ij}+c_{ij}) & \mathrm{Im}(g_{ij}+c_{ij})\\
\mathrm{Im}(-g_{ij}+c_{ij}) & \mathrm{Re}(g_{ij}-c_{ij})
\end{array}\right].
\end{array}
\end{eqnarray}}\unskip
Alternatively, we can work directly with the operators $\hat{X}_\alpha$
instead of their fluctuations $\Delta \hat{X}_\alpha$. This allows to
prove that quadratic CS violation is a sufficient condition for the
fulfillment of the GPH criterion. Indeed, suppose that $M_t\geq 0$.
Then, we can define an associated scalar product as $(u,v)_t\equiv
u^{\dagger}M_tv$, satisfying the CS inequality $|(u,v)_t|^2\leq
(u,u)_t(v,v)_t$. By choosing 
{\begin{equation}\label{eq:CSvectors4}
u=\frac{1}{\sqrt{2}}\left[\begin{array}{@{}c@{}}
0\\ 0\\ 1\\ \mathrm{i}
\end{array}\right],\quad v=\frac{1}{\sqrt{2}}\left[\begin{array}{@{}c@{}}
1\\ \mathrm{i}\\ 0\\ 0
\end{array}\right],
\end{equation}}\unskip
we obtain the quadratic CS inequality $|c_{ij}|^2\leq g_{ii}g_{jj}$.
Thus, quadratic CS violation implies that the matrix $M_t$ is not
non-negative, and hence the GPH criterion is satisfied. 



More generally, from the definition of separability,
Equation~(\ref{eq:separability}), we can apply the following chain of
CS inequalities if the state is separable,  
{\begin{eqnarray}\label{eq:separabilitycs}
|c_{ij}|&=&|\braket{\hat{a}_i\hat{a}_j}|=\left|\mathrm{Tr}[\hat{a}_i\hat{a}_j\hat{\rho}]\right|=\left|\sum_np_n\braket{\hat{a}_i}_n\braket{\hat{a}_j}_n\right|\nonumber\\
 &\leq&\sum_n p_n|\braket{\hat{a}_i}_n\braket{\hat{a}_j}_n|\leq\sum_n p_n\sqrt{\braket{\hat{a}^\dagger_i\hat{a}_i}_n\braket{\hat{a}^\dagger_j\hat{a}_j}_n}\nonumber\\
 &\leq&  \sqrt{\sum_n p_n \braket{\hat{a}^\dagger_i\hat{a}_i}_n}\sqrt{\sum_n p_n \braket{\hat{a}^\dagger_j\hat{a}_j}_n}=\sqrt{g_{ii}g_{jj}},
\end{eqnarray}}\unskip
where $\braket{\hat{A}}_n\equiv
\mathrm{Tr}[\hat{A}(\hat{\rho}^{(i)}_n\otimes\hat{\rho}^{(j)}_n)]$.
Thus, quadratic CS violation is a sufficient condition for
entanglement. However, the above derivation does not work for quartic
CS violation. As a counterexample, the product of two number states
$\hat{\rho}=\ket{nm}\bra{nm}$, with $n,m>0$, is clearly a separable
state that nevertheless violates the quartic CS inequality.

Technically, the GPH criterion here explained is only a sufficient
condition for the negativity of $\hat{\rho}_t$; in general, an infinite
set of sufficient and necessary conditions for the negativity of
$\hat{\rho}_t$, based on higher-order CS violations, can be
derived~\cite{Shchukin2005}. For multipartite entanglement, the role of
the CS inequality is played by the H\"older inequality~\cite{Wolk2014},
which is a generalization of the CS inequality. However, in practice,
the GPH criterion provides a very powerful  and simple tool to signal
entanglement. In particular, for bipartite Gaussian states, the GPH
criterion is a sufficient and necessary condition for
entanglement~\cite{Simon2000}.


\subsection{Cauchy--Schwarz violation and entanglement in Andreev--Hawking radiation}

We finally switch back to analogue gravity and, in particular, we apply
the above techniques to the study of Andreev--Hawking radiation in
condensates. We first note that the mean-field formalism can be
alternatively described in the Schr\"odinger picture by a coherent
ansatz of the form
{\begin{equation}
\ket{\Psi}=\mathrm{e}^{\int\mathrm{d}\mathbf{x}\left[\Psi(\mathbf{x},t)\hat{\Psi}^\dagger(\mathbf{x})-\Psi^*(\mathbf{x},t)\hat{\Psi}(\mathbf{x})\right]}\ket{0}_{\mathrm{MB}},
\end{equation}}\unskip
where $\ket{0}_{\mathrm{MB}}$ is the many-body vacuum containing no bosons,
$\hat{\Psi}(\mathbf{x})\ket{0}_{\mathrm{MB}}=0$. When this ansatz is inserted
into a variational principle, such as the Dirac-Frenkel one, the
time-dependent GP equation~(\ref{eq:TDGP}) is retrieved.


The genuine quantum character of the Andreev and Hawking effects is
revealed by reexamining the relation of
Equation~(\ref{eq:inoutmodesrelation}). After some
algebra~\cite{deNova2015}, it is shown that any matrix $S\in U(2,1)$
can be written as $S=B^{\dagger}S_{PA}A^{\dagger}$, with
{\begin{eqnarray}\label{eq:param}
\begin{array}{c}
\displaystyle S_{PA}=\left[\begin{array}{ccc} \mathrm{e}^{\mathrm{i}\phi} &0&0\\
0&\mathrm{e}^{-\mathrm{i}2\theta}\cosh r&\sinh r\\
0&\sinh r &\mathrm{e}^{\mathrm{i}2\theta} \cosh r\end{array}\right],\vspace*{3.5pt}\\
\displaystyle  A=\left[\begin{array}{@{}cc|c@{}}
{U_A} && 0  \\
  && 0 \\
  \hline
    0 & 0 & 1 
\end{array}\right],\quad\displaystyle  U_A\equiv \frac{1}{\sinh r}\left[\begin{array}{@{}cc@{}} -S_{d2d1}&S^*_{d2u}\\
\displaystyle S_{d2u}&S^*_{d2d1}\end{array}\right],\vspace*{3.5pt}\\
 B=\left[\begin{array}{@{}cc|c@{}}
{U_B} && 0  \\
  && 0 \\
  \hline
    0 & 0 & 1 
\end{array}\right],\quad\displaystyle  U_B\equiv\frac{1}{\sinh r}\left[\begin{array}{@{}cc@{}} -S_{d1d2}&S_{ud2}\\
\displaystyle S^*_{ud2}&S^*_{d1d2}\end{array}\right],
\end{array}
\end{eqnarray}}\unskip
where
{\begin{eqnarray}
\begin{array}{rcl}
\mathrm{e}^{\mathrm{i}\phi}&=&\det S, \vspace*{2.5pt}\\
\mathrm{e}^{\mathrm{i}2\theta}\cosh r&=&S_{d2d2},\vspace*{2.5pt}\\
\sinh r&=&\sqrt{|S_{ud2}|^2+|S_{d1d2}|^2}=\sqrt{|S_{d2u}|^2+|S_{d2d1}|^2}.
\end{array}
\end{eqnarray}}\unskip
Since the matrices $U_A,U_B$ are unitary, the above transformation
amounts to a change of basis in the normal $u-d1$ sector,
{\begin{equation}\label{eq:AHTransform}
\left[\begin{array}{@{}c@{}}
\hat{a}_{s}\\
\hat{a}_{AH}
\end{array}\right]=U^\dagger_A 
\left[\begin{array}{@{}c@{}}
\hat{a}_{u}\\
\hat{a}_{d1}
\end{array}\right],\quad \left[\begin{array}{@{}c@{}}
\hat{b}_{s}\\
\hat{b}_{AH}
\end{array}\right]=U^\dagger_B 
\left[\begin{array}{@{}c@{}}
\hat{b}_{u}\\
\hat{b}_{d1}
\end{array}\right].
\end{equation}}\unskip
In this new basis, the scattering of the spectator (s) channel is a
fully normal, unitary process. On the other hand, the normal hybrid
Andreev--Hawking (AH) channel couples to the anomalous $d2$ channel as
{\begin{equation}\label{eq:ParametricInOut}
\left[\begin{array}{@{}c@{}}
\hat{b}_{AH}\\
\hat{b}^{\dagger}_{d2}
\end{array}\right] =\left[\begin{array}{@{}cc@{}}
\mathrm{e}^{-\mathrm{i}2\theta}\cosh r&\sinh r\\
\sinh r &\mathrm{e}^{\mathrm{i}2\theta} \cosh r\end{array}\right]\left[\begin{array}{@{}c@{}}
\hat{a}_{AH}\\
\hat{a}_{d2}^{\dagger}
\end{array}\right] .
\end{equation}}\unskip
Apart from a trivial phase, this is precisely the same Bogoliubov
relation arising from a two-mode squeezed operator $U(\varepsilon)$,
Equation~(\ref{eq:TwoModeSqueezed}). Moreover, the unitary
transformation (\ref{eq:AHTransform}) does not change the incoming or
outgoing vacua. As a result, the ``in'' vacuum can be regarded as a
two-mode squeezed state for the ``out'' modes, and the spontaneous
production of outgoing AH modes simultaneously describes both the
Andreev and Hawking effects,
{\begin{eqnarray}\label{eq:AHEffect}
\bra{0_{\mathrm{in}}}\hat{b}_{AH}^{\dagger}(\omega)\hat{b}_{AH}(\omega')\ket{0_{\mathrm{in}}}=\delta(\omega-\omega')\sinh^2r(\omega)
=\delta(\omega-\omega')\left[|S_{ud2}(\omega)|^2+|S_{d1d2}(\omega)|^2\right].
\end{eqnarray}}\unskip
Thus, in the quantum optics jargon, the joint AH effect is nothing else
than a non-degenerate parametric amplifier, where the outgoing AH mode
is the signal and the outgoing partner $d2$ mode is the idler.
Furthermore, as discussed in Equation~(\ref{eq:ThermalSqueezed}), the
reduced state for the AH mode is a thermal state, in analogy with the
original prediction of Hawking. Notice, however, that this
transformation applies independently to each frequency, which results
in a $\omega$-dependent effective temperature.


Although the AH mode is very appealing from the theoretical point of
view, in practice, the channels $i,j=u,d1,d2$ are still more convenient
since (i) they are located in separated physical regions and (ii) quantum
states are typically defined in this basis. Specifically, a most
relevant class of quantum states is the family of incoherent Gaussian
states characterized by the following first and second-order momenta in
the incoming basis:
{\begin{eqnarray} \label{eq:incoherent}
\begin{array}{rcl}
 \braket{\hat{a}_{i}(\omega)\hat{a}_{j}(\omega')}&=& \braket{\hat{a}_{i}(\omega)}=0\vspace*{2.5pt}\\
\braket{\hat{a}_{i}^{\dagger}(\omega)\hat{a}_{j}(\omega')}
&=&n_i(\omega)\delta_{ij}\delta(\omega-\omega') .
\end{array}
\end{eqnarray}}\unskip
Since the state is Gaussian, any higher-order momentum can be put in
terms of these expectation values via Wick theorem. The above class of
states includes thermal states of incoming modes. A typical
choice~\cite{Macher2009a,Busch2014} is that where the incoming modes
have thermalized in the comoving frame of the condensate,
{\begin{equation}\label{eq:thermalin}
    n_i(\omega)=\frac{1}{\mathrm{e}^{\beta\Omega_{i-in}(\omega)}-1},
\end{equation}}\unskip
with $\Omega_{i-in}(\omega)=\Omega(k_{i-in}(\omega))$ the comoving Bogoliubov frequency in the corresponding asymptotic region; see Equation~(\ref{eq:dispersionrelation}). In order to assess the quantumness of the AH effect for incoherent Gaussian states (\ref{eq:incoherent}), we compute the first- and second-order correlation functions of Equations (\ref{eq:gdef}) and (\ref{eq:GammaDef}) for the ``out'' modes at given $\omega$,
{\begin{eqnarray}\label{eq:CSGamma}
\begin{array}{rcl}
 g_{ij}(\omega)&\equiv&\braket{\hat{b}_{i}^{\dagger}(\omega)\hat{b}_{j}(\omega)},\vspace*{2.5pt}\\
c_{ij}(\omega)&\equiv&\braket{\hat{b}_{i}(\omega)\hat{b}_{j}(\omega)},\vspace*{2.5pt}\\
\Gamma_{ij}(\omega)&\equiv&\braket{\hat{b}_{i}^{{{\dagger}}}(\omega)\hat{b}_{j}^{\dagger}(\omega)\hat{b}_{j}(\omega)\hat{b}_{i}(\omega)},
\end{array}
\end{eqnarray}}\unskip
where we obviate all Dirac delta factors; a discussion on how to
regularize the infinities resulting from $\delta(\omega=0)$ is present
at the beginning of Section~\ref{subsec:ExperimentalCS}. 

From the previous definitions, it is easy to show that the only
non-zero first-order correlations are
{\begin{eqnarray}\label{eq:quadraticorrs}
\begin{array}{rcl}
g_{IJ}&=&\alpha_{I}^{\dagger}\cdot\alpha_{J}, \vspace*{2.5pt}\\
g_{d2d2}&=&|\alpha_{d2}|^{2}-1, \vspace*{2.5pt}\\
c_{Id2}&=&\alpha_{d2}^{\dagger}\cdot\alpha_{I}.
 \end{array}
\end{eqnarray}}\unskip
On the other hand, since we are working with Gaussian states, the
second-order correlation functions are expressed in terms of the
first-order ones as
{\begin{eqnarray}\label{eq:CSWick}
\begin{array}{rcl}
\Gamma_{IJ} & = & g_{II}g_{JJ}+|g_{IJ}|^2=|\alpha_{I}|^{2}|\alpha_{J}|^{2}+|\alpha_{I}^{\dagger}\cdot\alpha_{J}|^2, \vspace*{3pt}\\
 \Gamma_{d2d2} & = & 2g^2_{d2d2}=2(|\alpha_{d2}|^{2}-1)^{2},\vspace*{3pt}\\
 \Gamma_{Id2} & = & |c_{Id2}|^2+g_{II}g_{d2d2}
 =|\alpha_{d2}^{\dagger}\cdot\alpha_{I}|^{2}+|\alpha_{I}|^{2}(|\alpha_{d2}|^{2}-1).
\end{array}
\end{eqnarray}}\unskip
In the above equations, we have compacted the notation by defining the complex vector
{\begin{eqnarray}\label{eq:defcorrs}
\alpha_i(\omega)&\equiv&\left[\begin{array}{@{}c@{}}
S_{iu}(\omega)\sqrt{n_u(\omega)}\vspace*{2pt}\\
S_{id1}(\omega)\sqrt{n_{d1}(\omega)}\vspace*{2pt}\\
S_{id2}(\omega)\sqrt{n_{d2}(\omega)+1}
\end{array}\right].
\end{eqnarray}}\unskip
This is equivalent to working with complex vectors $\alpha_i$ with
components $(\alpha_i)_j=S_{ij}$, and whose scalar product is given by
the metric $g=\mathrm{diag}[n_u,n_{d1},n_{d2}+1]$.

With the help of the above results, we quantify the quartic and
quadratic CS violations using $\Delta_{ij}(\omega)$ and
$\Theta_{ij}(\omega)$, Equations (\ref{eq:CSviolation2}),
(\ref{eq:CSviolation4}), respectively. Notice that, as explained after
Equation~(\ref{eq:CSgijimpossible}), no CS violation is possible for
$g_{ud1}$. Moreover, the condition for quartic CS violation in the
normal $u-d1$ correlations is equivalent to that for quadratic CS
violation, $\Theta_{ud1}=\Delta_{ud1}=|g_{ud1}|^2-g_{uu}g_{d1d1}\leq
0$. Thus, only the anomalous correlations characterizing the Andreev
and Hawking effects can give rise to genuine quantum correlations,
quantified by
{\begin{eqnarray}\label{eq:CSnonseparable}
\begin{array}{rcl}
\Delta_{Id2}&=&|\alpha_{d2}^{\dagger}\cdot\alpha_{I}|^{2}-|\alpha_{I}|^{2}(|\alpha_{d2}|^{2}-1),\vspace*{2.5pt}\\
 \Theta_{Id2}&=&\Delta_{Id2}.
\end{array}
\end{eqnarray}}\unskip
Similarly, only anomalous processes can be entangled. Since our state
is Gaussian, entanglement is equivalent to the condition
$\mathcal{P}_{Id2}(\omega)<0$, Equation~(\ref{eq:GPH}), where the GPH
function reads
{\begin{equation}\label{eq:GPHanomalous1}
\mathcal{P}_{Id2}=-\Delta_{Id2}[(g_{II}+1)(g_{d2d2}+1)-|c_{Id2}|^2].
\end{equation}}\unskip
Since $(g_{II}+1)(g_{d2d2}+1)\geq |c_{Id2}|^2$,
Equation~(\ref{eq:CScijpos}), we observe that, for the class of states
(\ref{eq:incoherent}), entanglement, quadratic and quartic CS
violations are all equivalent conditions, $\Delta_{Id2}(\omega)>0$.
Therefore, we can use $\Delta_{Id2}(\omega)$ simultaneously as both
entanglement and CS witness. By invoking the pseudounitarity of $S$,
Equation~(\ref{eq:pseudounitarity}), a simple explicit expression for
$\Delta_{Id2}$ can be\break derived:
{\begin{eqnarray}\label{eq:WitnessExplicit}
    \Delta_{Id2}&=&|S_{Id2}|^2(1+n_u+n_{d1}+n_{d2})
    -|S_{d2d1}|^2n_u-|S_{d2u}|^2n_{d1}\nonumber\\
    &&-\,|S_{I'u}|^2n_{d1}n_{d2}-|S_{I'd1}|^2n_{u}n_{d2}-|S_{I'd2}|^2n_{u}n_{d1},
\end{eqnarray}}\unskip
where here $I'=d1,u$ is the complementary normal channel to $I=u,d1$.

We compute $\Delta_{Id2}(\omega)$ for different BH solutions assuming a
thermal quantum state in the incoming channels, Equations
(\ref{eq:incoherent}), (\ref{eq:thermalin}). In particular, at zero
temperature, there is entanglement across the whole Andreev--Hawking
spectrum,
{\begin{equation}
    \Delta_{Id2}(\omega,T=0)=|S_{Id2}(\omega)|^2>0.
\end{equation}}\unskip
At finite temperature, $n_i(\omega)\neq 0$; in particular,
$n_{u}(\omega)$ is the only divergent occupation number since
$\Omega_{u-\mathrm{in}}(\omega)\sim \omega$ at low frequencies, while
the other comoving frequencies are finite at zero frequency since their
wavevector is $k_{i-\mathrm{in}}(\omega=0)=k_{\mathrm{BCL}}$,
$i=d1,d2$. Moreover, by invoking pseudounitarity and noting that
$|S_{ij}(\omega)|^2\sim 1/\omega$ only for $j=d1,d2$, we have that
$|S_{Id2}|^2<|S_{d2d1}|^2$ in this regime, which implies that there is
no entanglement close to $\omega\simeq 0$ at finite temperature. On the
other side of the spectrum, $\omega\simeq \omega_{{\max}}$,
$|S_{Id2}(\omega)|^2\sim \sqrt{\omega_{{\max}}-\omega}$, and
entanglement is lost again.

The results for the non-resonant BH solutions of
Figure~\ref{fig:BHModels} are depicted in Figure~\ref{fig:Witness}. We
observe that, for the flat-profile model (left column), entanglement
only survives at low temperatures and high frequencies for both Andreev
and Hawking radiation. However, for the delta and waterfall models
(center and right columns, respectively), entanglement is present even
at relatively high temperatures of the order of the chemical potential
for the Hawking effect (upper row), with a significant reduction of the
entanglement signal for the Andreev effect (lower row). These features
can be understood from (i) the  smallness of $\omega_{\max}$ for the
flat-profile model results in large occupation numbers $n_i(\omega)$
even at low temperatures, and (ii) for the delta and waterfall models,
$|S_{d1d2}(\omega)|^2\ll |S_{ud2}(\omega)|^2$ in the relevant part of
the spectrum, explaining the reduction of the Andreev entanglement.

\begin{figure}[t!]
\includegraphics{fig05}
\caption{\label{fig:Witness}
Entanglement witness $\Delta_{Id2}$ as function of $\omega$
for an incoherent thermal Gaussian state, Equations (\ref{eq:incoherent}),
(\ref{eq:thermalin}). Different colors label different temperatures, as
indicated in the legend of each panel. (a)--(c) $\Delta_{ud2}$ for the BH
solutions of Figures~\ref{fig:BHModels}a--c. (d)--(f) Same as (a)--(c) but for
$\Delta_{d1d2}$.}
\end{figure}

The above picture is modified for the resonant structures of
Figure~\ref{fig:Resonant}, as shown in
Figure~\ref{fig:ResonantWitness}. For the double-delta configuration
(left column), there is a strong entanglement signal near the resonant
peak, close to the zero-temperature value even at high temperatures for
both the Andreev and Hawking effects. For the resonant flat-profile
configuration (right column), although entanglement is lost for low
temperatures because of the smallness of $\omega_{\max}$, the Andreev
entanglement signal is now larger than the Hawking one close to the
resonant peak since there $|S_{d1d2}(\omega)|^2> |S_{ud2}(\omega)|^2$.
These results clearly demonstrate the potential of resonant structures
for studying the quantumness of the AH effect, in particular that of
the Andreev effect, greatly attenuated in non-resonant structures.

\begin{figure}[t!]
{\vspace*{5pt}}
\includegraphics{fig06}
\caption{\label{fig:ResonantWitness}Entanglement witness $\Delta_{Id2}$
as function of $\omega$ for an incoherent thermal Gaussian state, Equations
(\ref{eq:incoherent}), (\ref{eq:thermalin}). Different colors label
different temperatures, as indicated in the legend of each panel.
(a)--(b) $\Delta_{ud2}$ for the resonant BH solutions of
Figures~\ref{fig:Resonant}a,b. (c)--(d) Same as (a)--(b) but for
$\Delta_{d1d2}$.}
\end{figure}

To conclude the discussion, we examine the physical implications of
both incoherent and Gaussianity conditions. Gaussianity results from
the BdG approximation, see Equation~(\ref{eq:GrandCanonicalEnergy}),
and it is only expected to be broken in a strongly interacting regime
beyond Bogoliubov. In 1D, this regime is only achieved at low densities
(see discussion after Equation~(\ref{eq:RelativeCorrelationsBdG})). On
the other hand, incoherence results from a stationary description in
which the incoming modes that eventually scatter at the horizon are
populated in the asymptotic regions, for instance by a thermal bath. 

If we remove the incoherence in the incoming basis but not Gaussianity,
the GPH criterion is still equivalent to entanglement. However, the
quadratic CS violation becomes only a sufficient condition for the GPH
criterion and the quartic CS violation is then independent of the GPH
criterion (see Equation~(\ref{eq:CSvectors4}) and subsequent
discussion). On the other hand, if we remove Gaussianity but not
incoherence, the GPH criterion, now only a sufficient entanglement
condition, is independent from the quartic CS violation but equivalent
to the quadratic CS violation. For a general state which is neither
Gaussian nor incoherent in the incoming basis, the quadratic CS
violation is only a sufficient condition for the GPH criterion, which
in turn is a sufficient condition for the presence of entanglement; all
of them are independent from quartic CS violation. The above logical
relations are summarized in Table~\ref{TCS}.

\begin{table}
\caption{\label{TCS}Logical relations between the different quantum criteria
considered here, where CS2, CS4 stand for quadratic and quartic CS
violations, and EI means ``Entanglement independent''}
\begin{tabular}{cccrc}
\thead
Incoherent &  Gaussian & $\text{CS4}$ & $\text{CS2}$ & GPH\\
\endthead
${\checkmark}$ & ${\checkmark}$& $\Leftrightarrow$ GPH &$\Leftrightarrow$ GPH &$\Leftrightarrow$ Entanglement\\
${\checkmark}$ & ${\times}$ &~EI & $\Leftrightarrow$  GPH & $\Rightarrow$ Entanglement\\
${\times}$& ${\checkmark}$ &~EI & $\Rightarrow$  GPH &$\Leftrightarrow$ Entanglement\\
${\times}$ &${\times}$ &~EI & $\Rightarrow$ GPH &$\Rightarrow$ Entanglement
\botline
\end{tabular}
\end{table}

\subsection{Experimental considerations and conceptual exports}\label{subsec:ExperimentalCS}

Experimental proposals for the detection of CS violation and
entanglement in an analogue context were performed using time-of-flight
(TOF) techniques in Ref.~\cite{deNova2014} (see also
Ref.~\cite{Boiron2015}), and density--density correlations in
Ref.~\cite{Steinhauer2015}. A regularization procedure for the
infinities arising from the Dirac delta factors ignored in
Equation~(\ref{eq:CSGamma}), based on the use of windowed Fourier
transforms, was presented in Refs~\cite{deNova2014,deNova2015} for
each experimental scheme, respectively. The first claimed observation
of the Hawking effect~\cite{Steinhauer2016} was indeed based on the
detection of the quadratic CS violation $\Delta_{ud2}>0$ from the
measurement of density--density correlations. Later observations of the
Hawking effect~\cite{deNova2019,Kolobov2021}, although exhibiting a
more accurate agreement with the theoretical prediction for the Hawking
correlations $c_{ud2}(\omega)$, did not address the question of
entanglement or CS violation. Entanglement in the Andreev--Hawking
effect still represents an active topic of research, and it can be of
interest for quantum technologies, since then an analogue horizon
behaves as a source of entangled phonons.

Remarkably, the concepts and techniques discussed here to signal
quantum correlations can be exported to qudits through the
$P$-representation developed in Ref.~\cite{Giraud2008}. For
simplicity, we consider the particular case of a qubit, i.e., a
two-level quantum system consisting of two states $\ket{+},\ket{-}$.
Furthermore, for the sake of definiteness, we identify these states as
spin projections along the $z$-axis of a spin-$1/2$ particle,
$\sigma_z\ket{\pm}=\pm\ket{\pm}$, although the discussion can be
trivially adapted to any type of qubit by using the pseudospin formalism. 

The quantum state of a spin-$1/2$ particle is described by a $2\times 2$
density matrix of the form
{\vspace*{-1.5pt}}
{\begin{equation}\label{eq:QubitState}
    \rho=\sum_{n,m=\pm} \rho_{nm}\ket{n}\bra{m}=\frac{I_2+\mathbf{B}\cdot \boldsymbol{\sigma}}{2}
\end{equation}}\unskip
with $I_n$ the $n\times n$ identity matrix and $\boldsymbol{\sigma}$ a
vector containing the Pauli matrices. The Bloch vector $\mathbf{B}$
fully determines the quantum state, representing the spin polarization
of the particle,
$\mathbf{B}=\braket{\boldsymbol{\sigma}}=\mathrm{Tr}[\boldsymbol{\sigma}\rho]$, 
where we have invoked the trace orthogonality of the Pauli matrices,
$\mathrm{Tr}[\sigma_i]=0$,
$\mathrm{Tr}[\sigma_i\sigma_j]=2\delta_{ij}$. The above expression can
be regarded as the qubit version of the Fock expansion
(\ref{eq:GeneralQuantumStateFock}), where the spin states $\ket{\pm}$
play the role of the number states $\ket{n}$. Actually, the Fock space
of a fermion mode is also a qubit, spanned by the two number states
$\ket{n}$, $n=0,1$.

The analogue of the coherent states $\ket{\alpha}$ are the
spin-coherent states $\ket{\hat{\mathbf{n}}}$, 
{\vspace*{-1.5pt}}
{\begin{equation}\label{eq:CoherentState1/2}
    \ket{\hat{\mathbf{n}}}=\cos\frac{\theta}{2}\mathrm{e}^{-\mathrm{i}\frac{\phi}{2}}\ket{+}+\sin\frac{\theta}{2}\mathrm{e}^{\mathrm{i}\frac{\phi}{2}}\ket{-},
\end{equation}}\unskip
which are spin eigenstates with maximum projection along the direction
of the unit vector $\hat{\mathbf{n}}=[\sin \theta \cos\phi,\sin \theta
\sin\phi,\cos \theta]$, $(\hat{\mathbf{n}}\cdot
\boldsymbol{\sigma})\ket{\hat{\mathbf{n}}}=\ket{\hat{\mathbf{n}}}$.
These states also form an overcomplete basis as
{\vspace*{-1.5pt}}
{\begin{equation}\label{eq:Identity}     
|{\braket{\hat{\mathbf{n}}|\hat{\mathbf{n}}'}}|^2
=\frac{1+\hat{\mathbf{n}}\cdot \hat{\mathbf{n}}'}{2},\quad \frac{1}{2\rmpi}\int \mathrm{d}\Omega\ket{\hat{\mathbf{n}}}\bra{\hat{\mathbf{n}}}=\sum_{n=\pm}\ket{n}\bra{n}=1,
\end{equation}}\unskip
where $\Omega$ is the solid angle associated to the unit vector
$\hat{\mathbf{n}}$. Moreover, there are spin-squeezed states for larger
total spin~\cite{Kitagawa1993}, which play a central role in metrology.


The spin-coherent basis allows for a $P$-representation, 
{\vspace*{-1.5pt}}
{\begin{equation}\label{eq:PRepresentationQubit}
    \rho=\int \mathrm{d}\Omega~P(\hat{\mathbf{n}})\ket{\hat{\mathbf{n}}}\bra{\hat{\mathbf{n}}},\quad \int \mathrm{d}\Omega~P(\hat{\mathbf{n}})=1.
\end{equation}}\unskip
The function $P(\hat{\mathbf{n}})$ is also a quasi-probability
distribution, representing the spin analogue of the Glauber--Sudarshan
$P$-function. Since the density matrix (\ref{eq:QubitState}) is
diagonalized in the $\ket{\pm \hat{\mathbf{n}}}$ basis, with
$\hat{\mathbf{n}}\parallel \mathbf{B}$, $P(\mathbf{n})$ can be always
chosen as non-negative for one qubit.

Genuine quantum signatures arise when considering two spin-$1/2$
particles, labeled as $i,j$, whose total quantum state is described by
a $4\times 4$ density matrix 
{\vspace*{-1.5pt}}
{\begin{equation}\label{eq:TwoQubitState}
    \rho=\frac{I_4+\mathbf{B}_i\cdot \boldsymbol{\sigma}_i+\mathbf{B}_j\cdot \boldsymbol{\sigma}_j+\boldsymbol{\sigma}_i\cdot \mathbf{C}\cdot \boldsymbol{\sigma}_j}{4},
\end{equation}}\unskip
with $\boldsymbol{\sigma}_{i,j}$ a vector containing the Pauli matrices
in each subspace, $\mathbf{B}_{i,j}$ the individual spin polarizations,
and $\mathbf{C}$ the spin-correlation matrix. From the definition of
separability, Equation~(\ref{eq:separability}), and by using that any
one-qubit state is diagonalized in the spin-coherent basis, we find
that separability for two-qubit states is equivalent to the existence
of a non-negative $P$-representation, 
{\vspace*{-1.5pt}}
{\begin{equation}\label{eq:PRepresentation2Qubit}
\rho=\int\mathrm{d}\Omega_i\mathrm{d}\Omega_j~P(\mathbf{n}_i,\mathbf{n}_j)\ket{\mathbf{n}_i\mathbf{n}_i}\bra{\mathbf{n}_i\mathbf{n}_j},\quad P(\mathbf{n}_i,\mathbf{n}_j)\geq 0.
\end{equation}}\unskip
Hence, the absence of a non-negative $P$ representation is
automatically a signature of entanglement. This implies that any CS
violation is then a sufficient condition for entanglement. For
instance, 
{\vspace*{-3pt}}
{\advance\jot by -2pt\begin{eqnarray}\label{eq:CSViolationSpin}
 |\mathrm{Tr}[\mathbf{C}]|&=&|\braket{\boldsymbol{\sigma}_i\cdot \boldsymbol{\sigma}_j}|=\left|\int\mathrm{d}\Omega_i\mathrm{d}\Omega_j~P(\mathbf{n}_i,\mathbf{n}_j)\mathbf{n}_i\cdot\mathbf{n}_j\right|\nonumber\\
 &\leq &\int\mathrm{d}\Omega_i\mathrm{d}\Omega_j~P(\mathbf{n}_i,\mathbf{n}_j)|\mathbf{n}_i\cdot\mathbf{n}_j|\nonumber\\
&\leq& \int\mathrm{d}\Omega_i\mathrm{d}\Omega_j~P(\mathbf{n}_i,\mathbf{n}_j)=1
\end{eqnarray}}\unskip
is a classical CS inequality, based on the non-negativity of the
$P$-function. Qualitatively, we can understand this CS inequality as
the fact that the classical average of the scalar product 
of two unit vectors (such as the spin orientations
$\mathbf{n}_i,\mathbf{n}_j$) is never larger than one. Thus, 
{\vspace*{-1.5pt}}
{\begin{equation}\label{eq:EntanglementWitness}
    \Delta\equiv -\mathrm{Tr}[\mathbf{C}]-1>0
\end{equation}}\unskip
represents a CS violation that provides an entanglement witness.

Using the analogies above as a pipeline, and inspired by the techniques
discussed here for the study of quantum Andreev--Hawking radiation, as
well as by the fact that the Standard Model is based on a relativistic
quantum field theory in a flat spacetime, it was recently shown that
quantum correlations can be also studied at the 
LHC~\cite{Afik2021}. In particular, it was proven that the spin
quantum state of a pair of top-antitop quarks, the most massive
fundamental particles known to exist, can be fully reconstructed from
their decay products, implementing the so-called quantum tomography in
quantum information jargon. This is possible because the large top mass
is translated into a short lifetime that avoids any other process,
including hadronization, to affect its spin before the decay. 

Another remarkable source of inspiration was the study of quantum
steering in Hawking radiation by Robertson, Michel and 
Parentani~\cite{Robertson2017}, which directly motivated the analysis
of steering in top quarks~\cite{Afik2023}. In general, the study of
quantum information in high-energy physics is becoming an active topic
of research (see for instance
Refs~\cite{Fabbrichesi2021,Severi2022,Aguilar2022,Aoude2022,Barr2022,Bernal2023,Morales2023,Cheng2023,Sakurai2024,Dong2024,Afik2024});
the interested reader is referred to Ref.~\cite{Afik2022} for a
pedagogical introduction to the topic, aimed at an audience outside
particle physics. Moreover, the experimental proposal of
Ref.~\cite{Afik2021} has been implemented by both the ATLAS and
CMS collaborations~\cite{ATLAS2024,CMS2024}, leading to the first
observation of entanglement in quarks and to the highest-energy
entanglement detection ever achieved. Specifically, the entanglement
witness (\ref{eq:EntanglementWitness}) was directly measured from the
angular distribution of the separation between the leptons arising from
the top-antitop decay, obtaining $\Delta>0$ with more than $5\sigma$
(the standard candle for discovery in particle physics), which also
represents the violation of a CS\break inequality.

\section{Black-hole lasers}\label{sec:BHL}

Another main topic of research involving resonant analogue
configurations is the so-called black-hole laser
(BHL)~\cite{Corley1999}. The BHL emerges in a configuration similar to
that of resonant BH solutions, but now the asymptotic downstream region
is again subsonic. As a result, a BHL displays a pair of BH/WH
horizons, and the resulting finite-size supersonic cavity becomes
unstable due to the successive bouncing of Andreev--Hawking radiation
between the\break horizons. 


Qualitatively, we can understand the BHL instability as the partner
$d2$ modes from the Andreev--Hawking effect being reflected at the WH
horizon as $d2$-in modes that bounce back towards the BH, further
stimulating the production of Andreev--Hawking radiation and thus
leading to a process of self-amplification, similar to that occurring
in a lasing cavity. Quantitatively, the BHL effect is characterized by
a discrete BdG spectrum of dynamical instabilities, computed by
extending the usual scattering problem (see
Equation~(\ref{eq:solitonspinors}) and subsequent discussion) to
complex frequencies and retaining only the asymptotically bounded modes
outside the cavity. This procedure bears some resemblance to the
computation of the discrete spectrum of bound states for an attractive
potential in the Schr\"odinger equation, where the usual scattering
problem for positive energies is extended to negative energies, keeping
only the exponentially decaying solutions at infinity.

A systematic procedure for the quantization of the unstable lasing
modes was provided by Finazzi and Parentani~\cite{Finazzi2010}. For
simplicity, we discuss the case of a single unstable BdG mode $z_I$
with complex frequency $\omega=\gamma+\mathrm{i}\Gamma$, where $\gamma$ is the
real, oscillatory part of the frequency, and $\Gamma$ is the imaginary
part of the frequency, determining the growth rate of the instability.
We also assume that the mode is non-degenerate, which means that
$\gamma\neq0$ so $\omega\neq -\omega^*$. In general, any dynamically
unstable mode $z_I$ has associated a stable mode $z_S$ with frequency
$\omega^*$~\cite{Leonhardt2003}. Their eigenvalue equation reads 
{\begin{equation}
M_0z_I=\omega z_I,\quad M_0z_S=\omega^*z_S,
\end{equation}}\unskip
with $M_0$ the BdG matrix operator, Equation~(\ref{eq:BdGEigenmode}).
Because of their complex frequency, both modes have zero norm
$(z_S|z_S)=(z_I|z_I)=0$; see Equation~(\ref{eq:EigenOrto}). However, we
can choose their normalization such that
{\begin{equation}\label{eq:NormalizationLasing}
    (z_S|z_I)=-(\bar{z}_I|\bar{z}_S)=1.
\end{equation}}\unskip
Properly normalized states are defined through
{\begin{equation}\label{eq:UnstableToNormalModes}
Z_{+}\equiv \tfrac{1}{\sqrt{2}}(z_I+z_S),\quad Z_{-}\equiv\tfrac{1}{\sqrt{2}}(\bar{z}_I-\bar{z}_S),
\end{equation}}\unskip
satisfying 
{\begin{equation}\label{eq:BHLOrto}
    (Z_{+}|Z_{+})=(Z_{-}|Z_{-})=1,\quad (\bar{Z}_{-}|Z_{+})=(Z_{-}|Z_{+})=0.
\end{equation}}\unskip

As a result, their quantum amplitudes 
{\begin{equation}
\hat{a}_{\pm}=(Z_{\pm}|\hat{\Phi})
\end{equation}}\unskip
do behave as proper annihilation operators (see
Equation~(\ref{eq:Aniquilacion})). Their time evolution is easily
derived from the quantum amplitudes of the original complex eigenmodes,
{\begin{equation}\label{eq:UnstableToNormalAmplitudes}
\hat{a}_{I}=(z_S|\hat{\Phi})=\frac{1}{\sqrt{2}}(\hat{a}_{+}+\hat{a}^{\dagger}_{-}),\quad \hat{a}_{S}=(z_I|\hat{\Phi})=\frac{1}{\sqrt{2}}(\hat{a}_{+}-\hat{a}^{\dagger}_{-})
\end{equation}}\unskip
which are not annihilation operators as
$[\hat{a}_{I},\hat{a}^\dagger_{I}]=[\hat{a}_{S},\hat{a}^\dagger_{S}]=0$, $[\hat{a}_{I},\hat{a}^\dagger_{S}]=[\hat{a}_{S},\hat{a}^\dagger_{I}]=1$. 
Nevertheless, they evolve as expected from
Equation~(\ref{eq:BogoliubovQuantum}),
{\begin{eqnarray}
\begin{array}{rcl}
\displaystyle
 \mathrm{i}\partial_t\hat{a}_{I}&=&(z_S|M_0\hat{\Phi})=\omega \hat{a}_{I}\Longrightarrow
\hat{a}_{I}(t)=\hat{a}_{I}\mathrm{e}^{-\mathrm{i}\omega t}, \vspace*{2.5pt}\\
\displaystyle  \mathrm{i}\partial_t\hat{a}_{S}&=&(z_I|M_0\hat{\Phi})=\omega^*\hat{a}_{S}\Longrightarrow
\hat{a}_{S}(t)=\hat{a}_{S}\mathrm{e}^{-\mathrm{i}\omega^* t}.
\end{array}
\end{eqnarray}}\unskip
To invert the relation, it is quite convenient to employ matrix
notation. First, Equation~(\ref{eq:UnstableToNormalAmplitudes}) can be
rewritten as 
{\begin{equation}\label{eq:LasingNormalRelation}
   \left[\begin{array}{@{}c@{}} \hat{a}_{I}\\
\hat{a}_{S}\end{array}\right]=U \left[\begin{array}{@{}c@{}} \hat{a}_{+}\\
\hat{a}^\dagger_{-}\end{array}\right],\quad  U\equiv \mathrm{e}^{-\mathrm{i}\frac{\rmpi}{4}\sigma_y}\sigma_z=\frac{1}{\sqrt{2}}\left[\begin{array}{rr} 1 & 1\\
1 &-1\end{array}\right].
\end{equation}}\unskip
The matrix $U=U^\dagger=U^{-1}$ describes a spin inversion in the
$x$--$y$ plane plus a rotation of $\rmpi/2$ around the $y$-axis,
satisfying $U\sigma_{x,z} U^\dagger=U^\dagger\sigma_{x,z}
U=\sigma_{z,x}$ and $U\sigma_y U^\dagger=U^\dagger\sigma_y
U=-\sigma_y$. In this notation, the time evolution of the  amplitudes
of the complex modes simply reads
{\begin{equation}
    \left[\begin{array}{@{}c@{}} \hat{a}_{I}(t)\\
\hat{a}_{S}(t)\end{array}\right]=\left[\begin{array}{@{}c@{}} \mathrm{e}^{-\mathrm{i}\omega t}\hat{a}_{I}\\
\mathrm{e}^{-\mathrm{i}\omega^* t}\hat{a}_{S}\end{array}\right]=\mathrm{e}^{-\mathrm{i}\gamma t}\mathrm{e}^{\Gamma t \sigma_z}\left[\begin{array}{@{}c@{}} \hat{a}_{I}\\
\hat{a}_{S}\end{array}\right].
\end{equation}}\unskip
As a result, we trivially find
{\begin{eqnarray}
\left[\begin{array}{@{}c@{}} \hat{a}_{+}(t)\\
\hat{a}^\dagger_{-}(t)\end{array}\right]=U^\dagger \left[\begin{array}{@{}c@{}} \hat{a}_{I}(t)\\
\hat{a}_{S}(t)\end{array}\right]=\mathrm{e}^{-\mathrm{i}\gamma t}\mathrm{e}^{\Gamma t \sigma_x}\left[\begin{array}{@{}c@{}} \hat{a}_{+}\\
\hat{a}^\dagger_{-}\end{array}\right]
 =\mathrm{e}^{-\mathrm{i}\gamma t}\left[\begin{array}{@{}cc@{}} \cosh \Gamma t & \sinh \Gamma t\\
\sinh \Gamma t &\cosh \Gamma t\end{array}\right]\left[\begin{array}{@{}c@{}} \hat{a}_{+}\\
\hat{a}^\dagger_{-}\end{array}\right].
\end{eqnarray}}\unskip
This is a similar evolution to that of a non-degenerate parametric
amplifier. This can be better seen by examining the contribution from
these modes to the field spinor $\hat{\Phi}$ of
Equation~(\ref{eq:QuantumFieldFluctuations}),
{\begin{eqnarray}
 \hat{\Phi}_{\mathrm{L}}=Z_{+}\hat{a}_{+}+Z_{-}\hat{a}_{-}+\bar{Z}_{+}\hat{a}^{\dagger}_{+}+\bar{Z}_{-}\hat{a}^{\dagger}_{-}
=z_{I}\hat{a}_{I}+z_{S}\hat{a}_{S}+\bar{z}_{I}\hat{a}^{\dagger}_{I}+\bar{z}_{S}\hat{a}^{\dagger}_{S}.
\end{eqnarray}}\unskip
When inserted into the Bogoliubov expansion for the grand-canonical
Hamiltonian (\ref{eq:GrandCanonicalEnergy}), which governs the
dynamics, we obtain an orthogonal contribution to that of the regular
Bogoliubov sector with real frequencies,
{\begin{eqnarray}
  \hat{K}_{\mathrm{L}}&=&\frac{1}{2}(\hat{\Phi}_{\mathrm{L}}|M_0\hat{\Phi}_{\mathrm{L}})=[\hat{a}^\dagger_{I}~\hat{a}^\dagger_{S}]\left[\begin{array}{@{}cc@{}}0& \Omega^*  \nonumber\\
\Omega & 0\end{array}\right]\left[\begin{array}{@{}c@{}} \hat{a}_I\\
\hat{a}_S\end{array}\right]  \nonumber\\
&=&[\hat{a}^\dagger_{+}~\hat{a}_{-}](\gamma \sigma_z-\Gamma \sigma_y)\left[\begin{array}{@{}c@{}} \hat{a}_{+}\\
\hat{a}^\dagger_{-}\end{array}\right]  \nonumber\\
&=&\gamma(\hat{a}^{\dagger}_{+}\hat{a}_{+}-\hat{a}^{\dagger}_{-}\hat{a}_{-})
+\mathrm{i}\Gamma(\hat{a}^{\dagger}_{+}\hat{a}^{\dagger}_{-}-\hat{a}_{+}\hat{a}_{-}),
\end{eqnarray}}\unskip
where we neglect zero-point contributions. This is precisely the same
Hamiltonian of a non-degenerate parametric amplifier (see
Equation~(\ref{eq:TwoModeSqueezed}) and ensuing discussion). Another
remarkable feature of dynamical instability is that there is no
well-defined vacuum for the unstable modes~\cite{Ribeiro2022}. This is
immediately seen by noticing that any Bogoliubov transformation
{\vspace*{2pt}}
{\begin{equation}
\left[\begin{array}{@{}c@{}} \hat{b}_{+}\\
\hat{b}^\dagger_{-}\end{array}\right]=\mathrm{e}^{u \sigma_x}\left[\begin{array}{@{}c@{}} \hat{a}_{+}\\
\hat{a}^\dagger_{-}\end{array}\right]=\left[\begin{array}{@{}cc@{}} \cosh u& \sinh u\\
\sinh u &\cosh u\end{array}\right]\left[\begin{array}{@{}c@{}} \hat{a}_{+}\\
\hat{a}^\dagger_{-}\end{array}\right],
\end{equation}}\unskip
leaves invariant $\hat{K}_{\mathrm{L}}$, with each $ \hat{b}_{\pm}$ giving rise to
a different vacuum.

The above derivations can be easily adapted to the case of a degenerate
unstable mode, $\gamma=0$. This implies that $\omega=-\omega^*$, so
$z_I,z_S$ are proportional to $\bar{z}_{I},\bar{z}_{S}$. In particular,
one can set $\bar{z}_{I}=z_I$ and $\bar{z}_{S}=-z_S$, and thus
Equation~(\ref{eq:NormalizationLasing}) is again satisfied. However,
now there is only one independent normalizable mode
{\vspace*{2pt}}
{\begin{equation}\label{eq:UnstableToNormalModesDegenerate}
Z\equiv \tfrac{1}{\sqrt{2}}(z_I+z_S),\quad \bar{Z}\equiv\tfrac{1}{\sqrt{2}}(z_I-z_S),
\end{equation}}\unskip
whose amplitude $\hat{a}=(Z|\hat{\Phi})$ is a proper annihilation
operator. Consequently, Equation~(\ref{eq:LasingNormalRelation}) now
reads
{\vspace*{2pt}}
{\begin{equation}\label{eq:LasingNormalRelationDegenerate}
\left[\begin{array}{@{}c@{}} \hat{a}_{I}\\
\hat{a}_{S}\end{array}\right]=U \left[\begin{array}{@{}c@{}} \hat{a}\\
\hat{a}^\dagger\end{array}\right],
\end{equation}}\unskip
giving rise to the time evolution
{\vspace*{2pt}}
{\begin{equation}
\left[\begin{array}{@{}c@{}} \hat{a}\\
\hat{a}^\dagger\end{array}\right]=\mathrm{e}^{\Gamma t \sigma_x}\left[\begin{array}{@{}c@{}} \hat{a}\\
\hat{a}^\dagger\end{array}\right]=\left[\begin{array}{@{}cc@{}} \cosh \Gamma t & \sinh \Gamma t\\
\sinh \Gamma t &\cosh \Gamma t\end{array}\right]\left[\begin{array}{@{}c@{}} \hat{a}\\
\hat{a}^\dagger\end{array}\right],
\end{equation}}\unskip
which results from the grand-canonical Hamiltonian
{\vspace*{2pt}}
{\begin{equation}
    \hat{K}=\mathrm{i}\Gamma\frac{(\hat{a}^\dagger)^2-\hat{a}^2}{2}.
\end{equation}}\unskip
This is indeed the Hamiltonian of a degenerate parametric amplifier,
where the same mode is both the signal and the idler. Interestingly,
this Hamiltonian is, after a trivial phase transformation, that of an
unstable harmonic oscillator,
{\vspace*{2pt}}
{\begin{equation}\label{eq:UnstablePendulum}
    \hat{H}=\omega\frac{\hat{p}^2-\hat{q}^2}{2}=-\omega\frac{(\hat{a}^\dagger)^2+\hat{a}^2}{2}.
\end{equation}}\unskip
The resulting time-evolution operator $\mathrm{e}^{-\mathrm{i}\hat{H}
t}$ is just the squeezing operator (\ref{eq:Squeezing}) with a linearly
increasing amplitude $\varepsilon=-\mathrm{i}\omega t$. Therefore,
dynamically unstable modes can be understood as unstable harmonic
oscillators, in contrast to the regular Bogoliubov modes of purely real
frequency, which behave as normal harmonic oscillators.



Returning to the BHL, the above quantization procedure is applied
separately to each unstable lasing mode. Due to the exponential
parametric amplification of the lasing modes, at some point the
Bogoliubov approximation ceases to be valid and one needs to take into
account higher-order interacting terms to describe the dynamics.
Consequently, we separate our discussion of the BHL effect following
the three different stages of its time evolution: (i) short times, when
the dynamics is still governed by the linear BdG equations; (ii)
intermediate times, when the evolution  is driven solely by the dominant
unstable mode up to the saturation regime, where the full interacting
Hamiltonian is required again; and (iii) long times, when the system
reaches its final state after the collapse of the metastable state
achieved in the saturation\break regime.


\subsection{Short times: linear and non-linear spectra}

The microscopic derivation of Ref.~\cite{Finazzi2010} was extended
in another seminal work by Michel and Parentani~\cite{Michel2013} using
a simple analytical model based on the flat-profile configuration,
Figure~\ref{fig:BHModels}a, where the flow velocity is homogeneous,
$v(x)=q$, and the speed of sound is changed to $c(x)=c_2<q$ for
$|x|<L/2$, Figure~\ref{fig:BHLModels}a (we still set the asymptotic
subsonic sound speed to $c_u=1$). We label this stationary BHL solution
as $\Psi_{\mathrm{BHL}}$. 


\begin{figure}[t!]
\includegraphics{fig07}
\caption{\label{fig:BHLModels}(a)--(c) BHL solutions resulting from mirroring the BH
solutions of Figures~\ref{fig:BHModels}a--c. (d)--(f) Linear BdG
spectrum of dynamical instabilities for the BHL solutions of (a)--(c).
Solid (dashed) lines represent the imaginary (real) part $\Gamma_n$
($\gamma_n$) of the complex frequency $\omega_n$. The dash-dotted blue
line is the inverse of the round-trip time. Solid (dashed) vertical
lines mark the critical lengths $L_n$ ($L_{n+1/2}$). (g)--(i)
Non-linear spectrum of stationary GP solutions $\Psi_n$ for the
background configurations of (a)--(c), where the color code is chosen
to match the associated lasing modes in (d)--(f). The density profile
of the initial BHL solutions \mbox{(a)--(c)} is depicted as solid
blue.}
\end{figure}

The discrete BdG spectrum of complex frequencies arising from
$\Psi_{\mathrm{BHL}}$ is depicted in Figure~\ref{fig:BHLModels}d as a
function of $L$. The critical lengths of the cavity $L=L_n$ at which
the $n$th dynamically unstable mode emerges (vertical solid lines) are
given by
{\begin{equation}\label{eq:CriticalLengths}
    L_n=\frac{\varphi_0+2\rmpi n}{k_{\mathrm{BCL}}}=L_0+n\lambda_{_{\mathrm{BCL}}},\quad n=0,1\ldots
\end{equation}}\unskip
where
{\begin{equation}\label{eq:CriticalFP}
    \varphi_0=2\arctan\sqrt{\dfrac{1-q^2}{q^2-c_2^2}},\quad k_{\mathrm{BCL}}=2\sqrt{q^2-c_2^2}.
\end{equation}}\unskip
This equation can be simply understood as that, after some threshold
length $L_0=\varphi_0/k_{\mathrm{BCL}}$ at which the first unstable
lasing mode appears, the cavity gives birth to a new unstable mode each
BCL wavelength $\lambda_{{\mathrm{BCL}}}=2\rmpi/k_{\mathrm{BCL}}$. The
lasing modes are initially degenerate, i.e., they have purely imaginary
frequency, $\omega_n=\mathrm{i}\Gamma_n$. For lengths $L>L_{n+1/2}$, with
$L_{n+1/2}$ obtained by inserting half-integer values $n+1/2$ in the
above equation (vertical dashed lines), the $n$th unstable mode
becomes non-degenerate, developing a non-vanishing real part of the
frequency\break $\gamma_n\neq 0$. 

The dominant mode is that with the largest growth rate $\Gamma_n$, and
determines the overall growth rate $\Gamma$ of the lasing instability,
$\Gamma=\max_{n}\Gamma_n$. For short cavities, this is typically the
mode with the largest $n$. However, as the cavity becomes longer and
longer, the competition between the different unstable modes becomes
stronger and stronger. We can compare these exact results with an
estimation for the growth rate resulting from the qualitative picture
of bouncing Hawking radiation, $\Gamma\sim 1/\tau_{\mathrm{RT}}$, with
$\tau_{\mathrm{RT}}$ the roundtrip time for a zero-frequency $d2$ mode
to travel back and forth between the horizons; the zero-frequency
choice is motivated by the small value of $\gamma_n$ for the dominant
mode observed in the plot. The result is depicted in dashed-dotted blue
line, finding that it provides a decent estimation for long cavities.
An elaborated WKB calculation shows a much more accurate agreement with
the exact BdG results~\cite{Michel2013}; however, it completely misses
the existence of degenerate unstable modes, which are the dominant ones
in short cavities. Thus, WKB prescriptions can only be used reliably in
the long-cavity\break limit.

Interestingly, the work by Michel and Parentani~\cite{Michel2013}
further established a perfect correspondence between the emergence of
dynamical instabilities in the BdG spectrum and the emergence of
non-linear stationary solutions in the GP equation
$\Psi_n(x),~n=0,1\ldots.$ Specifically, these are stationary GP
solutions for the same underlying Hamiltonian that are smoothly
connected (as a function of $L$) to $\Psi_{\mathrm{BHL}}$, sharing the
same conserved current $J$ and chemical potential $\mu$. The
correspondence is shown in Figure~\ref{fig:BHLModels}g, where the
spectrum of stationary non-linear GP solutions for the largest value of
$L$ in Figure~\ref{fig:BHLModels}d is represented using the same color
code of the associated lasing modes. The GP solution $\Psi_n(x)$ first
emerges at $L=L_n$ as a sinusoidal oscillation around the supersonic
cavity with the BCL wavelength (see for instance $\Psi_3(x)$ in
dashed-dotted magenta), eventually becoming a non-linear cnoidal wave
as the cavity enlarges. These solutions have lower grand-canonical
energy (\ref{eq:energeticstability}) than the original BHL solution
$\Psi_{\mathrm{BHL}}$, following the\break hierarchy
{\begin{equation}
 K[\Psi_0]<K[\Psi_1]<K[\Psi_2]<\cdots<K[\Psi_{\mathrm{BHL}}].
\end{equation}}\unskip
Michel and Parentani~\cite{Michel2015} conjectured that all of them are
also dynamically unstable except for the ground state solution
$\Psi_0(x)$, which accumulates particles in the cavity in order to
become fully subsonic and evaporate the horizons. It can be proven that
the degeneracy breaking at $L=L_{n+1/2}$ of the $n$th unstable mode
can be also attributed to the emergence of a new non-linear GP
solution, but this one is asymmetric and contains one soliton minimum
outside the cavity, thus being energetically unfavored with respect to
the symmetric solutions $\Psi_n(x)$ (i.e., their density
$|\Psi_n(x)|^2$ has even parity with respect to the center of the
lasing cavity). Remarkably, this non-linear spectrum of solutions is
similar to that arising in a superconductor/normal/superconductor
junction~\cite{Sols1994} due to the GP--GL\break correspondence.


More realistic BHL models were developed in
Ref.~\cite{deNova2017a}, where it was shown that any BH solution
leading to a homogeneous supersonic region, as those of
Figure~\ref{fig:BHModels}, can be mirrored to produce a WH solution by
parity inversion of the Hamiltonian and time-reversal symmetry of the
GP wavefunction. By matching the two BH/WH solutions in the homogeneous
supersonic region, one obtains a symmetric BHL solution
$\Psi_{\mathrm{BHL}}(x)$ with a homogeneous lasing cavity of arbitrary
length $L$. Indeed, the flat-profile BHL solution of
Figure~\ref{fig:BHLModels}a is a particular example of this general
result. Two more examples are the attractive well and double-delta BHL
solutions of Figures~\ref{fig:BHLModels}b,c, obtained from the
corresponding waterfall and delta BH solutions of
Figures~\ref{fig:BHModels}b,c. Their spectrum of dynamical
instabilities is computed in Figures~\ref{fig:BHLModels}e,f, which
exhibits the same trends as the flat-profile case. In particular, the
critical lengths $L_n$ at which a new dynamical instability emerges are
also given by Equation~(\ref{eq:CriticalLengths}), where
{\begin{equation}
    \varphi_0=\rmpi,\quad k_{\mathrm{BCL}}=2\sqrt{v^2_d-c^2_d}=2\sqrt{\frac{1}{q^2}-q^2},
\end{equation}}\unskip
for the attractive-well BHL solution, while for the double-delta BHL solution 
{\begin{eqnarray}
    \varphi_0=4\arcsin\sqrt{\frac{1-r}{2}},\quad k_{\mathrm{BCL}}=2\sqrt{v^2_d-c^2_d},\quad
    r=c^2_d\sqrt{\frac{2(M^2_d-1)}{2Z^2c^2_d+\sqrt{q^4+8q^2}(1-c^2_d)}},
\end{eqnarray}}\unskip
with $M_d=v_d/c_d=q/c^3_d$ the supersonic Mach number and $Z$ the
amplitude of the delta barrier; see Equations 
(\ref{eq:CompactBHWF})--(\ref{eq:CompactDelta}) and 
ensuing discussion for the details of the
waterfall and delta configurations. Each lasing mode has again
associated a non-linear symmetric GP solution $\Psi_n(x)$ smoothly
connected to $\Psi_{\mathrm{BHL}}(x)$ as a function of the cavity
length, Figures~\ref{fig:BHLModels}h,i. Moreover, the $n$th unstable
mode also becomes non-degenerate at $L=L_{n+1/2}$, coinciding with the
emergence of an asymmetric GP solution. The inverse of the roundtrip
time still provides an estimation for the growth rate that improves for
long cavities. In summary, all the trends predicted by Michel and
Parentani in Ref.~\cite{Michel2013} are further confirmed by these
alternative BHL solutions. The only exception is the appearance of a
dynamical instability in the attractive square well for $0<L<L_0$,
labeled as the short-length (SL) mode in Figure~\ref{fig:BHLModels}e.
This instability is not related to the BHL effect itself but its origin
lies in the fact that  $\Psi_{\mathrm{BHL}}(x)$ is here smoothly
connected to the soliton solution (\ref{eq:GraySoliton}) when no
potential is present ($L=0$). However, the soliton has larger energy
than the homogeneous plane wave $\Psi_0(x)=\mathrm{e}^{\mathrm{i}qx}$, which in turn is
smoothly connected to the actual ground state, labeled as GS in
Figure~\ref{fig:BHLModels}h. Thus, the SL mode is simply a consequence
of the energetic instability of the BHL solution for any $L>0$. This
provides further numerical evidence for the conjecture of
Ref.~\cite{Michel2015}: in flowing scattering configurations,
energetic and dynamical instability are equivalent conditions.

\subsection{Intermediate times: quantum amplification in the BHL--BCL
crossover}

The exponential amplification of the dominant lasing mode drives the
linear Bogoliubov dynamics for times $\Gamma t\gtrsim 1$ until the
saturation regime, when it typically reaches the non-linear stationary
GP solution with the largest $n$. However, due to the energetic
instability of the supersonic cavity, the exponential growth of the
dominant lasing mode can be overshadowed by the coherent stimulation of
the BCL wave resulting from the presence of an obstacle in the flow.
For instance, as originally shown in
Refs~\cite{Wang2016,Wang2017}, this can be the case of the WH
horizon itself in highly time-dependent configurations, far from the
fine-tuned stationary BHL solutions of Figure~\ref{fig:BHLModels}.
Moreover, since the real part of the frequency of the dominant mode is
very small, it contains wavevectors close to that of the BCL mode,
something that greatly complicates their clear distinction in real
setups~\cite{Steinhauer2014,Kolobov2021,Steinhauer2022,Tettamanti2016,Steinhauer2017,Wang2016,Wang2017,Llorente2019,Tettamanti2021,deNova2023}. 
All of this leads to a strong competition between the BCL and BHL
mechanisms.

We can understand coherent BCL stimulation from a simple model based on
linear response theory~\cite{Carusotto2006}. We consider a general
stationary condensate, solution of the time-independent GP
equation~(\ref{eq:TIGP}). At $t=0$, a small external perturbation
described by a potential $W(x,t)$ is introduced. Expanding the GP
wavefunction as
{\begin{equation}\label{eq:ClassicalExpansion}
\Psi(\mathbf{x},t)=\left[\Psi_0(\mathbf{x})+\varphi(\mathbf{x},t)\right]\mathrm{e}^{-\mathrm{i}\mu t}
\end{equation}}\unskip
in the time-dependent GP equation~(\ref{eq:TDGP}) leads, at linear
order in the external perturbation, to the classical BdG equations with
a source,
{\begin{eqnarray}\label{eq:TIBdGSource} 
[\mathrm{i}\hbar\partial_t-M_0]\Phi(\mathbf{x},t)&=&\mathrm{i}W(\mathbf{x},t)z_\theta(\mathbf{x}),\nonumber \\
\Phi(\mathbf{x},t)&=&\left[\begin{array}{r}\varphi(\mathbf{x},t)\\ \varphi^*(\mathbf{x},t)\end{array}\right],\quad z_\theta(\mathbf{x})=\left[\begin{array}{r}-\mathrm{i}\Psi_0(\mathbf{x})\\ \mathrm{i}\Psi^*_0(\mathbf{x})\end{array}\right],
\end{eqnarray}}\unskip
where $z_\theta$ is the zero-frequency Nambu--Goldstone mode associated
to the spontaneous $U(1)$-symmetry breaking by the coherent GP
wavefunction~\cite{Lewenstein1996}; see
Equation~(\ref{eq:NambuGoldstone}) and ensuing discussion. This
equation can be solved by performing a classical expansion in terms of
the complete set of BdG eigenmodes, analogous to that of
Equation~(\ref{eq:QuantumFieldFluctuations}),
{\begin{equation}
    \Phi(\mathbf{x},t)=\sum_{n}a_{n}(t)z_{n}(\mathbf{x})\mathrm{e}^{-\mathrm{i}\omega_nt}+a^*_{n}(t)\bar{z}_{n}(\mathbf{x})\mathrm{e}^{\mathrm{i}\omega_nt},
\end{equation}}\unskip
where we subtract here the intrinsic time evolution of each mode. This
leads to
{\begin{eqnarray}\label{eq:BCLStimulation}
\partial_t a_n=\mathrm{e}^{\mathrm{i}\omega_n t}(z_n|W(t) z_\theta)\Longrightarrow
a_n(t)= \int^t_0\mathrm{d}t'~\mathrm{e}^{\mathrm{i}\omega_n t'}(z_n|W(t') z_\theta).
\end{eqnarray}}\unskip
Thus, for sufficiently long times, the amplitude of the collective
modes, corresponding here to the BdG modes, is determined by the
Fourier spectrum of the external perturbation, as expected from the
usual theory of linear response. In the specific case of a supersonic
condensate, since the BCL mode has zero-frequency, it is resonantly
stimulated by any static obstacle in the supersonic flow: this is
precisely the origin of Landau criterion for superfluidity. Notice that
this stimulation is coherent, imprinted on the GP wavefunction, and
thus it has a completely classical\break nature.


In order to study the BHL--BCL crossover, we first model the initial
background condensate (before the BHL and/or BCL onsets) within the
bulk of the lasing cavity by a supersonic plane wave of the form
$\Psi_0(x)\simeq \sqrt{n_0}\mathrm{e}^{\mathrm{i}q_0x}$, 
with an associated healing
length $\xi_0$. For the analysis, we focus on the expectation values of
the density and its correlations, measurable in the laboratory through
{in situ} imaging after averaging over ensembles of repetitions
of the
experiment~\cite{Shammass2012,Steinhauer2014,Steinhauer2016,deNova2019,Kolobov2021}. 
After expanding the density in terms of the quantum fluctuations of the
background condensate, one obtains
{\begin{eqnarray}\label{eq:Density expansion}
\begin{array}{l}
      \hat{n}(x,t)=\hat{\Psi}^\dagger(x,t)\hat{\Psi}(x,t)= n_0+\delta\hat{n}^{(1)}(x,t)+\delta\hat{n}^{(2)}(x,t),\vspace*{2.5pt}\\
    \delta\hat{n}^{(1)}(x,t)=\Psi_0^*(x)\hat{\varphi}(x,t)+\Psi_0(x)\hat{\varphi}^\dagger(x,t)=
    -\mathrm{i}z^\dagger_\theta \sigma_z \hat{\Phi},\vspace*{2.5pt}\\
   \delta\hat{n}^{(2)}(x,t)=\hat{\varphi}^\dagger(x,t)\hat{\varphi}(x,t),
\end{array}
\end{eqnarray}}\unskip
where we separate the linear contribution in the field fluctuations
$\delta\hat{n}^{(1)}(x,t)$ (the same of
Equation~(\ref{eq:BdGfieldequationTIHydrodynamic2}), obtained within
the BdG approximation), from the quadratic contribution
$\delta\hat{n}^{(2)}(x,t)$. Since dimensional analysis implies that the
quantum fluctuations around $\Psi_0$ scale as $\hat{\varphi}\sim
1/\sqrt{\xi_0}$, we have the scalings 
{\begin{equation}\label{eq:DensityScalings} 
    \delta\hat{n}^{(1)}\sim \sqrt{\frac{n_0}{\xi_0}},\quad \delta\hat{n}^{(2)}\sim \frac{1}{\xi_0}.
\end{equation}}\unskip
Using these results, we characterize the BHL and BCL mechanisms through
the first-order correlation function
{\begin{equation}
     G^{(1)}(x,t)\equiv \frac{\braket{\hat{n}(x,t)}}{n_0}-1,
\end{equation}}\unskip
which measures the amplitude of the developing density modulation above
the background condensate, and the \textit{normalized} density--density
correlation function 
{\begin{equation}\label{eq:NormalizedRelativeCorrelationsBdG}
    G^{(2)}(x,x',t)\equiv n_0\xi_0 g^{(2)}(x,x',t),
\end{equation}}\unskip
which in turn measures the quantum fluctuations around the density
modulation, $g^{(2)}(x,x',t)$ being the \textit{relative}
density--density correlation function 
{\begin{eqnarray}\label{eq:RelativeCorrelationsBdG}
 g^{(2)}(x,x',t)\equiv \frac{\braket{\hat{n}(x,t)\hat{n}(x',t)}-\braket{\hat{n}(x,t)}\braket{\hat{n}(x',t)}}{n^2_0}  
    \simeq \frac{\braket{\delta\hat{n}^{(1)}(x,t)\delta\hat{n}^{(1)}(x',t)}}{n^2_0}\sim \frac{1}{n_0\xi_0},
\end{eqnarray}}\unskip
where we take the leading contribution in the Bogoliubov approximation.
The relative density--density correlation function $g^{(2)}$ provides
the relative amplitude of the quantum fluctuations, so we can regard
$(n_0\xi_0)^{-1}$ as the initial strength of the quantum fluctuations,
which must be small for the Bogoliubov approximation to be valid,
$n_0\xi_0\gg 1$. On the other hand, the normalization of $G^{(2)}$
ensures that it is a dimensionless function that does not depend
explicitly on $n_0$ in the Bogoliubov approximation, only implicitly
through the healing length $\xi_0$. Since both BHL and BCL mechanisms
involve modes with well-defined wavevectors within the lasing cavity,
we will use the Fourier transforms in the supersonic region of the
above observables as figures of merit,
{\begin{eqnarray}
\begin{array}{rcl}
G^{(1)}_{\mathrm{peak}}(t)&\equiv& \max_{k} |G^{(1)}(k,t)|,\vspace*{3.5pt}\\
G^{(2)}_{\mathrm{peak}}(t)&\equiv& \max_{k,k'} |G^{(2)}(k,k',t)|.
\end{array}
\end{eqnarray}}\unskip
In real space, this peaked Fourier structure is translated into a
ripple in the ensemble-averaged density profile and into a checkerboard
pattern in the density--density correlations, respectively.

By borrowing the analogy with an unstable pendulum from
Equation~(\ref{eq:UnstablePendulum}), we can distinguish three main
regimes in the BHL--BCL crossover, represented in
Figure~\ref{fig:Pendula}, depending on the interplay between quantum
fluctuations and classical BCL stimulation, where the former are
controlled by the dimensionless amplitude 
{\begin{equation}\label{eq:QFAmplitude}
     A_{\mathrm{QF}}\sim \frac{\hat{\varphi}}{\Psi_0}\sim \frac{\delta\hat{n}^{(1)}}{n_0} \sim \frac{1}{\sqrt{n_0\xi_0}}\ll 1,
\end{equation}}\unskip
while the latter is controlled by the relative amplitude of the
coherent BCL wave with respect to the background condensate,  
{\begin{equation}
    A_{\mathrm{BCL}}\sim \frac{\varphi}{\Psi_0}.
\end{equation}}\unskip


\begin{figure}[t!]
\includegraphics{fig08}
\caption{\label{fig:Pendula}Schematic depiction of the three different
regimes of the BHL--BCL crossover using an analogy with an unstable
pendulum. (a) Quantum BHL: Quantum fluctuations cause the unstable
equilibrium position to collapse due to the Heisenberg uncertainty
principle. (b) Classical BHL: A small kick on the pendulum displaces it
some angle $\theta$ from the unstable equilibrium position, falling
down with a well-defined classical trajectory as a result. (c) BCL: An
external force (horizontal arrows) pushes the pendulum out of
equilibrium, governing the dynamics instead of gravity.} 
\end{figure}

\subsubsection{Quantum BHL}

When $A_{\mathrm{BCL}}\ll A_{\mathrm{QF}}\ll 1$, the BHL instability is
purely triggered by quantum fluctuations (e.g., there is no BCL
stimulation, $A_{\mathrm{BCL}}=0$), and the dynamics is driven by the
parametric amplification of the dominant lasing mode $z_I$, whose
frequency and quantum amplitude are $\omega=\gamma+\mathrm{i}\Gamma$ and
$\hat{a}_I$, respectively. The contribution of the dominant mode to the
quantum field fluctuations is
{\begin{equation}\label{eq:QuantumDominant}
    \hat{\Phi}(x,t)\simeq \mathrm{e}^{\Gamma t}\left[z_I(x)\mathrm{e}^{-\mathrm{i}\gamma t}\hat{a}_I+\bar{z}_I(x)\mathrm{e}^{\mathrm{i}\gamma t}\hat{a}^\dagger_I\right].
\end{equation}}\unskip
In a quantum BHL, the phase of the amplitude of the dominant mode is
expected to be random and hence we can take
$\braket{\hat{a}_I\hat{a}_I}\simeq 0$. This assumption yields that
$G^{(1)}_{\mathrm{peak}},G^{(2)}_{\mathrm{peak}}$ behave as
{\begin{eqnarray}\label{eq:QBHLGrowth}
 G^{(1)}_{\mathrm{peak}}(t)&\sim& \frac{\braket{\delta\hat{n}^{(2)}}}{n_0}\sim \frac{\mathrm{e}^{2\Gamma t}}{n_0\xi_0}, \nonumber\\
    G^{(2)}_{\mathrm{peak}}(t)&\sim& n_0\xi_0 \frac{\braket{\delta\hat{n}^{(1)}\delta\hat{n}^{(1)}}}{n^2_0}\sim \mathrm{e}^{2\Gamma t}.
\end{eqnarray}}\unskip
Hence, for a quantum BHL, the correlation functions
$G^{(1)}_{\mathrm{peak}}(t)$, $G^{(2)}_{\mathrm{peak}}(t)$ scale
quadratically in the field fluctuations. In the case of
$G^{(1)}_{\mathrm{peak}}$, this is because the $\mathbb{Z}_2$ symmetry
of a purely quantum BHL sets $\braket{\hat{\varphi}(x,t)}=0$ and thus,
$\braket{\delta\hat{n}^{(1)}(x,t)}=0$, as originally discussed by
Michel and Parentani~\cite{Michel2015}. 

In the pendulum analogy, the unstable equilibrium position is the
initial BHL solution and gravity is the lasing instability. Due to the
Heisenberg uncertainty principle, the unstable equilibrium position
collapses at the quantum level, and then the pendulum falls,
Figure~\ref{fig:Pendula}a. This is akin to the parametric amplification
of the lasing instability, where the $\mathbb{Z}_2$ symmetry can be
understood as that of the unstable equilibrium position of the
pendulum. 

The exponential growth of the dominant mode ceases when the system
reaches the saturation regime, corresponding to one of the stationary
GP solutions of the spectrum (typically, that with the largest $n$),
where the density modulation becomes of the order of the background
density itself, $G^{(1)}_{\mathrm{peak}}\sim 1$; see lower row of
Figure~\ref{fig:BHLModels}. Since this saturation stems from the
amplification of the quantum fluctuations of the dominant mode, we also
have that $g_{\mathrm{sat}}^{(2)}\sim 1$. As a result, the saturation
values of both correlation functions are roughly
{\begin{eqnarray}
\begin{array}{rcl}
 G^{(1)}_{\mathrm{sat}}&\sim&\displaystyle 1\sim \frac{\mathrm{e}^{2\Gamma t_{\mathrm{sat}}}}{n_0\xi_0 }, \vspace*{3.5pt}\\
    G^{(2)}_{\mathrm{sat}}&\sim&\displaystyle n_0\xi_0 \sim \mathrm{e}^{2\Gamma t_{\mathrm{sat}}},
\end{array}
\end{eqnarray}}\unskip
where $t_{\mathrm{sat}}$ is the time needed to reach saturation,
{\begin{equation}
    t_{\mathrm{sat}}\sim \frac{\ln n_0\xi_0}{2\Gamma}.
\end{equation}}\unskip

\subsubsection{Classical BHL} 

Here, $A_{\mathrm{QF}}\ll A_{\mathrm{BCL}}\ll 1$, so BHL amplification
still dominates the dynamics but the seed of the instability is now the
classical amplitude of the BCL wave in the condensate, leading to a
well-defined mean-field trajectory. Specifically, the perturbation of
the flow gives a classical coherent amplitude to the dominant lasing
mode through the stimulation of the BCL wave (see
Equation~(\ref{eq:ClassicalExpansion}) and subsequent results), which
is then exponentially amplified as
{\begin{equation}\label{eq:ClassicalDominant}
    \Phi(x,t)\simeq \mathrm{e}^{\Gamma t}\left[z_I(x)\mathrm{e}^{-\mathrm{i}\gamma t}a_I+\bar{z}_I(x)\mathrm{e}^{\mathrm{i}\gamma t}a^*_I\right].
\end{equation}}\unskip
Hence, $\braket{\delta\hat{n}^{(1)}}\neq 0$, and the $\mathbb{Z}_2$ symmetry is broken, 
{\begin{equation}\label{eq:Z2BreakingLinear}
   G^{(1)}_{\mathrm{peak}}(t)\sim \frac{\braket{\delta\hat{n}^{(1)}}}{n_0}\sim A_{\mathrm{BCL}}\mathrm{e}^{\Gamma t}.
\end{equation}}\unskip
Precisely because of its classical deterministic character, at the
linear level the BCL amplitude does not show up in the density--density
correlation function, and $G^{(2)}_{\mathrm{peak}}(t)$ still follows
Equation~(\ref{eq:QBHLGrowth}) in this regime. Therefore, the
$\mathbb{Z}_2$ symmetry-breaking implies now
{\begin{equation}\label{eq:Z2Breaking}
    G^{(1)}_{\mathrm{peak}}(t)\sim \mathrm{e}^{\Gamma t},\quad  G^{(2)}_{\mathrm{peak}}(t)\sim \mathrm{e}^{2\Gamma t}.
\end{equation}}\unskip

In the pendulum analogy, a classical BHL is akin to separate the
pendulum some small angle $\theta$ from its equilibrium position, which
consequently falls following a well-defined classical trajectory,
Figure~\ref{fig:Pendula}b. Here, the angle $\theta$ plays the role of
the Cherenkov amplitude $A_{\mathrm{BCL}}$ that seeds the BHL
instability, breaking the original $\mathbb{Z}_2$ symmetry of the
problem. 

In the saturation regime, $G^{(1)}_{\mathrm{peak}}\sim 1$, which now implies that
{\begin{equation}\label{eq:SaturationTimeCBHL}
    G^{(1)}_{\mathrm{sat}}\sim 1 \sim A_{\mathrm{BCL}}\mathrm{e}^{\Gamma t_{\mathrm{sat}}}.
\end{equation}}\unskip
Therefore, the saturation time is predicted to behave as
{\begin{equation}
    t_{\mathrm{sat}}\sim - \frac{\ln A_{\mathrm{BCL}}}{\Gamma},
\end{equation}}\unskip
and then
{\begin{equation}\label{eq:SaturationAmplitudeCBHL}
    G^{(2)}_{\mathrm{sat}}\sim \mathrm{e}^{2\Gamma t_{\mathrm{sat}}}\sim A^{-2}_{\mathrm{BCL}}.
\end{equation}}\unskip

\subsubsection{BCL}  

When the BCL amplitude is highly non-linear, $A_{\mathrm{QF}}\ll
A_{\mathrm{BCL}}\sim 1$, it dominates the mean-field dynamics towards
the saturation regime, overshadowing the lasing mechanism. Since the
BCL stimulation depends on the specific perturbation of the flow (see
Equation~(\ref{eq:BCLStimulation})), in general, no analytical formula
is available for the evolution of
$G^{(1)}_{\mathrm{peak}}(t),G^{(2)}_{\mathrm{peak}}(t)$. In the
pendulum analogy, the BCL regime is akin to applying a strong external
force to the pendulum that overcomes gravity,
Figure~\ref{fig:Pendula}c, so its evolution depends on the specific
force.



Regarding saturation, a large BCL amplitude is no longer described by
linear response theory, but it rather requires the full GP equation.
The saturation amplitude then scales as $G^{(1)}_{\mathrm{sat}}\sim
A^2_{\mathrm{BCL}}$ by definition of BCL amplitude. Regarding the
density fluctuations $G^{(2)}$, the sharp peaked structure of the BCL
wave now acts as a new mean-field background over which fluctuations
evolve. This gives rise to a checkerboard pattern in the correlation
function whose origin is completely different to that from a BHL, which
there stems from the exponential amplification of the quantum
fluctuations of the lasing modes, with wavevectors close to the BCL
one. Therefore, we can expect
$G^{(2)}_{\mathrm{sat}}=F(A_{\mathrm{BCL}})$, where $F$ is in general a
monotonically increasing function of $A_{\mathrm{BCL}}$ that also
depends on the other parameters of the background flow. Due to the
strong BCL stimulation, the saturation time is essentially limited by
$t_{\mathrm{sat}}\gtrsim \tau_{\mathrm{BCL}}$, where
$\tau_{\mathrm{BCL}}$ is the time that it takes the BCL wave to expand
along the whole lasing cavity.


\subsubsection{Upshot}

All the above scalings were originally derived in
Ref.~\cite{deNova2023} within a model based on a flat-profile BHL
solution (Figure~\ref{fig:BHLModels}a) with a delta barrier placed
exactly at the WH horizon to stimulate BCL radiation, which allows to
isolate the contribution of both mechanisms to the dynamics, finding a
good agreement with numerical results. Interestingly, it was also
stressed there that each regime of the BHL--BCL crossover can be
characterized according to its efficiency as a quantum amplifier.
Specifically, we measure the quantum amplification using the relative
density--density correlations, following
Equation~(\ref{eq:RelativeCorrelationsBdG}). In the initial state,
{\begin{equation}
    g_{\mathrm{peak}}^{(2)}(t=0)\sim \frac{1}{n_0\xi_0}\ll 1, 
\end{equation}}\unskip
which is the input of the quantum amplifier, while the output is the
saturation value $g_{\mathrm{sat}}^{(2)}$. This implies that the gain
of the quantum amplifier $\mathcal{G}$ is directly proportional to the
saturation value $G_{\mathrm{sat}}^{(2)}$ since
{\begin{equation}
    \mathcal{G}\equiv \frac{g_{\mathrm{sat}}^{(2)}}{g_{\mathrm{peak}}^{(2)}(t=0)}\propto n_0\xi_0 g_{\mathrm{sat}}^{(2)}= G_{\mathrm{sat}}^{(2)}. 
\end{equation}}\unskip
For each regime, we find
{\begin{eqnarray}
\begin{array}{rcl}
   \mathcal{G}_{\mathrm{QBHL}}&\sim& n_0\xi_0,  \vspace*{2.5pt}\\
        \mathcal{G}_{\mathrm{CBHL}}&\sim& \mathrm{e}^{2\Gamma t_{\mathrm{sat}}}, \vspace*{2.5pt}\\
    \mathcal{G}_{\mathrm{BCL}}&\sim &F(A_{\mathrm{BCL}}).
\end{array}
\end{eqnarray}}\unskip
This means that a quantum BHL behaves as a non-linear quantum
amplifier, since it amplifies the initial quantum fluctuations up to
the same saturation amplitude $g_{\mathrm{sat}}^{(2)}\sim 1$, so the
gain depends on the input amplitude $1/n_0\xi_0$. On the other hand,
classical BHL and BCL do behave as linear quantum amplifiers (i.e.,
their gain does not depend on the initial quantum strength
$1/n_0\xi_0$). In a classical BHL, the gain is exponentially large in
the saturation time $t_{\mathrm{sat}}$, which in turn decreases with
the BCL amplitude $A_{\mathrm{BCL}}$. This is because
$t_{\mathrm{sat}}$ is the lasing time during which the exponential
amplification of quantum fluctuations takes place. Hence, for
increasing $A_{\mathrm{BCL}}$, the system starts closer to the
saturation regime and less amplification is needed to reach it. In the
BCL regime, the gain is exponentially smaller as compared to that of a
classical BHL since there is no microscopic mechanism of exponential
amplification, and the enhancement of quantum fluctuations stems just
from the large BCL modulation of the background mean-field density.
This also implies that the function $F(A_{\mathrm{BCL}})$ determining
the gain increases with $A_{\mathrm{BCL}}$, in stark contrast with the
decrease expected for lasing amplification. Thus, the dependence of the
gain with respect to the BCL amplitude can be used to distinguish
between classical BHL and BCL in experiments. 

It is also remarkable that, in the quantum BHL and BCL regimes, the
behavior of $G_{\mathrm{peak}}^{(1)},G_{\mathrm{peak}}^{(2)}$ is
correlated since they are dominated by the same mechanism (either
exponential amplification of quantum fluctuations or classical BCL
stimulation), while classical BHL is a hybrid regime where
$G_{\mathrm{peak}}^{(1)}$ has a classical nature while
$G_{\mathrm{peak}}^{(2)}$ has a quantum one. Another qualitative
criterion of distinguishability is the non-monotonic behavior of the
growth rate of the density ripple and the checkerboard pattern of the
density--density correlations with respect to the parameters determining
the background flow (for instance, the cavity length $L$, as shown in
central row of Figure~\ref{fig:BHLModels}), in contrast with the smooth
behavior expected for BCL stimulation. This non-monotonicity is a
typical feature of resonant structures, and it is also observed in the
peak structure of the Andreev--Hawking spectrum discussed in
Section~\ref{sec:ResonantHawking}~\cite{Zapata2011}. A summary of the
main results for each regime is presented in Table~\ref{TBCL}.

\begin{table}
\caption{\label{TBCL}Summary of the different scalings for the three
regimes of the BHL--BCL crossover: quantum BHL, classical BHL, and BCL
\vspace*{-2pt}}
\tabcolsep 4pt
\begin{tabular}{ccccccc}
\thead
 &\raisebox{-.8em}{}\raisebox{1.25em}{} $G^{(1)}_{\mathrm{peak}}(t)$ & $G^{(2)}_{\mathrm{peak}}(t)$ & $G^{(1)}_{\mathrm{sat}}$ & $G^{(2)}_{\mathrm{sat}}$ & $t_{\mathrm{sat}}$ & Monotonic \\
\endthead
Quantum BHL &\raisebox{1.2em}{} ${\sim} \mathrm{e}^{2\Gamma t}/n_0\xi_0 $ & ${\sim} \mathrm{e}^{2\Gamma t}$ & ${\sim}$1 & $\sim n_0\xi_0 $ &  ${\sim}\ln n_0\xi_0 /2\Gamma$  & No\vspace*{2pt}\\
Classical BHL  & ${\sim} A_{\mathrm{BCL}}\mathrm{e}^{\Gamma t}$ & ${\sim} \mathrm{e}^{2\Gamma t}$  & ${\sim}$1 & ${\sim}\mathrm{e}^{2\Gamma t_{\mathrm{sat}}} \sim A_{\mathrm{BCL}}^{-2}$ &  ${\sim}-\ln A_{\mathrm{BCL}}/\Gamma$ & No\vspace*{2pt}\\
BCL  & -- & -- & ${\sim} A^2_{\mathrm{BCL}}$ & $F(A_{\mathrm{BCL}}) $ &  ${\gtrsim}\tau_{\mathrm{BCL}}$ & Yes\vspace*{2pt}
\botline
\end{tabular}
\tabnote{There is no analytical prediction for
$G^{(1)}_{\mathrm{peak}}(t),G^{(2)}_{\mathrm{peak}}(t)$ in the BCL
regime since they depend on the particular stimulation source.
$F(A_{\mathrm{BCL}})$ is an increasing function of  $A_{\mathrm{BCL}}$
and $\tau_{\mathrm{BCL}}$ is the time that it takes the BCL wave to
reach the BH horizon. The column ``Monotonic'' indicates a monotonic
dependence on the background parameters of the flow.}
\vspace*{-3pt}
\end{table}


\subsection{Long times: spontaneous Floquet states}

As explained, the lasing instability grows up to the saturation regime,
where it reaches a certain stationary GP solution. However, this
solution is also dynamically unstable, with a lifetime much longer than
that of the initial BHL solution, and eventually collapses. After some
(possibly very long) transient, where the system may intercept a number
of intermediate metastable GP solutions, it was numerically
observed~\cite{deNova2016} that the flat-profile BHL of
Figure~\ref{fig:BHLModels}a only has two possible fates: it either
reaches the true ground state (solid black line in
Figure~\ref{fig:BHLModels}g) or the so-called Continuous Emission of
Solitons (CES) state, where the system self-oscillates periodically
while continuously emitting solitons into the downstream region. Due to
its periodic nature, the CES state has been argued to be the
\textit{bona-fide} black-hole laser~\cite{deNova2016}.

Both trajectories are represented in Figures~\ref{fig:CES}a,b, where we
show the time evolution of the density $|\Psi(x,t)|^2$ for two
different flat-profile BHL configurations where some initial noise has
been added to trigger the BHL instability. We take a short cavity $L=2$
in both cases since that implies a large growth rate for the lasing
instability as well as only one stationary GP solution in the
non-linear spectrum, the ground state $\Psi_0(x)$, considerably
shortening the transient towards the final state. In
Figure~\ref{fig:CES}a, for $q=0.7,c_2=0.2$, the system directly reaches
the stable ground state $\Psi_0(x)$ in the saturation regime (vertical
red stripe centered at $x=0$). After further increasing the initial
flow velocity to $q=0.9$, Figure~\ref{fig:CES}b, the system also
approaches the ground state by increasing the density in the cavity
while expelling a gray soliton upstream to conserve particle number;
however, now the flow velocity is high enough to drag the soliton back
to the cavity (blue half-rings at the left of the vertical red stripe),
which then passes to the downstream region and travels along the flow
(diagonal blue lines downstream). The process is accompanied by the
emission of waves (diagonal yellow lines upstream) to ensure
conservation of particle number and energy. The passage of the dragged
soliton through the cavity restarts the cycle, giving rise to a
periodic behavior; this is the physical mechanism behind the CES\break state.

\begin{figure}[t!]
\includegraphics{fig09}
\caption{\label{fig:CES}(a) Spatio-temporal density profile $|\Psi(x,t)|^2$ for an
initial flat-profile BHL solution with $c_2=0.2$, $L=2$ and $q=0.7$.
Some random noise is initially added to trigger the BHL instability. (b)
Same as (a) but now $q=0.9$. (c) Dynamical phase diagram as a function of
$(L,q)$ with fixed $c_2=0.2$ for the final state of a flat-profile BHL,
where SB denotes that the cavity is subsonic and DS denotes the
dynamically stable region. (d) Fourier spectrum $|\Psi(x,\omega)|^2$ of
the simulation in (b) once in the CES state.}
\vspace*{5pt}
\end{figure}

The final state only depends on the background parameters of the flow
$(L,q,c_2)$, being quite insensitive to the particular details of the
transient or the initial noise. This gives rise to a dynamical phase
diagram, shown in Figure~\ref{fig:CES}c as a function of $(L,q)$ for
fixed $c_2=0.2$. Above the green region of dynamical stability (denoted
as DS), whose upper boundary is given by the condition $L=L_0(q,c_2)$
(see Equations (\ref{eq:CriticalLengths}), (\ref{eq:CriticalFP})), the
initial BHL solution first asymptotically reaches the ground state,
denoted as GS (blue region). Above some critical flow velocity\break
$q=q_c(L,c_2)$, numerically obtained, the final state is the CES state
(red region), which can be then regarded as a non-linear extension of
the Landau criterion. Later work~\cite{deNova2021} numerically showed
that the final fate of BHL solutions from an attractive well
(Figure~\ref{fig:BHLModels}h) is also either the ground state or the
CES state, suggesting the universality of the long-time behavior of a
BHL. Moreover, it was observed that the same CES state can be reached
even without starting from a BHL solution, indicating that the CES
state is an intrinsic state of the system, and not some fine-tuned\break
trajectory. 

This observation, along with the periodicity of the CES state, led to
identify a novel type of quantum state~\cite{deNova2022}: the so-called
spontaneous Floquet state, which is a state of a time-independent
Hamiltonian that oscillates like a Floquet state due to many-body
interactions, without the need of external periodic driving. The
emergence of a spontaneous Floquet state can be easily understood from
the time-dependent GP equation for a time-independent\break Hamiltonian,
which is a non-linear Schr\"odinger equation of the form
{\begin{eqnarray}\label{eq:HamiltonianEffective}
   \mathrm{i}\partial_t\Psi(x,t)=\displaystyle H_{\mathrm{GP}}(x,t)\Psi(x,t),\quad
   H_{\mathrm{GP}}(x,t)=\displaystyle-\frac{\partial_x^2}{2}+V(x)+g(x)|\Psi(x,t)|^2,
\end{eqnarray}}\unskip
where we allow for a possible inhomogeneous coupling constant (as it is
the case of a flat-profile BHL configuration). Notice that the only
possible time dependence of $H_{\mathrm{GP}}(x,t)$ results from that of
the GP wavefunction itself. A periodic density $|\Psi(x,t)|^2$, as that
of the CES state, implies an effective periodic Hamiltonian,
$H_{\mathrm{GP}}(x,t+T)=H_{\mathrm{GP}}(x,t)$, where $T$ is the
oscillation period.\break Self-consistently, $\Psi(x,t)$ becomes a Floquet
state of its own periodic Hamiltonian,
{\begin{equation}\label{eq:FloquetWaveMu}
\Psi(x,t)=\Psi_0(x,t)\mathrm{e}^{-\mathrm{i}\tilde{\mu}t},\quad \Psi_0(x,t)=\sum^{\infty}_{n=-\infty}u_n(x)\mathrm{e}^{-\mathrm{i}n\omega_0t},
\end{equation}}\unskip
with $\Psi_0(x,t+T)=\Psi_0(x,t)$, $\omega_0=2\rmpi/T$, and $\tilde{\mu}$
the quasi-chemical potential, defined modulo $\omega_0$. By inserting
this expansion into the GP equation, we get self-consistent equations
for the Floquet components $u_n(x)$,
{\begin{eqnarray}\label{eq:FloquetEquation}
    n\omega_0u_n(\mathbf{x})= \left[-\frac{\partial_x^2}{2}+V(x)-\tilde{\mu}\right]u_n(x)
    + \sum^{\infty}_{m=-\infty}\sum^{\infty}_{k=-\infty} g(x)u^*_{k}(x)u_{k+n-m}(x)u_{m}(x).
\end{eqnarray}}\unskip
The CES state provides a particular example of spontaneous Floquet
state, as can be seen from its Fourier transform $\Psi(x,\omega)$,
represented in Figure~\ref{fig:CES}d. We observe that the Fourier
spectrum consists of a series of equispaced lines at frequencies
$\omega=\tilde{\mu}+n\omega_0$, which can be identified as the Floquet
components $u_n(x)$. The dominant Floquet component has a frequency of
the order of the initial chemical potential $\omega\simeq q^2/2$
(notice that for the flat-profile BHL solution we subtract the
interacting contribution to the chemical potential, as discussed after
Equation~(\ref{eq:FlatProfileCondition})), which we can identify as a
non-trivial quasi-chemical potential $\tilde{\mu}\neq 0$. 

Interestingly, the concept of spontaneous Floquet state can be extended
to any many-body system whose dynamics can be described by a
variational ansatz that leads to an effective self-consistent
Hamiltonian like the GP equation, as proven in
Ref.~\cite{deNova2022}. This includes several canonical many-body
descriptions such as the MultiConfiguration Time-Dependent Hartree
method for bosons and fermions~\cite{Caillat2005,Alon2008}, the
Hartree-Fock equations for fermions~\cite{Thouless2014}, or the
Gutzwiller ansatz in Bose--Hubbard models~\cite{Jaksch1998}.

Furthermore, since a spontaneous Floquet state breaks the
time-translation symmetry of the underlying Hamiltonian, the CES state
represents a realization of continuous time
crystal~\cite{Wilczek2012,Kongkhambut2022}. This in stark contrast with
discrete time crystals~\cite{Sacha2015,Else2016,Zhang2017,Choi2017},
arising in conventional, driven Floquet systems, where the periodicity
is imposed by a subharmonic response to the external driving, and the
resulting symmetry breaking is discrete, not continuous. Actually, the
CES state was shown~\cite{deNova2022} to satisfy typical time-crystal
criteria of robustness (against the presence of time-dependent disorder
or variations of the parameters of the Hamiltonian), independence from
the initial state, and universality (i.e., it is a feature of a wide
class of\break Hamiltonians). 

For instance, a simple realization of the CES state is achieved by
quenching an attractive delta barrier $V(x)=-Z\delta(x)$ at $t=0$ in a
homogeneous flowing condensate with velocity $q$, described by an
initial GP wavefunction $\Psi(x,0)=\mathrm{e}^{\mathrm{i}qx}$. The resulting dynamics is
deterministic, and the final state of the system is described by a
similar dynamical phase diagram, solely function of $(Z,q)$, which only
displays the ground state at low velocities and the CES state at high
velocities, Figure~\ref{fig:Critical}a. The GS/CES phase diagram is an
example of dynamical phase
transition~\cite{Moeckel2008,Sciolla2010,Lang2018}, where the ground
state is the symmetry-unbroken phase, with continuous time translation
symmetry, while the CES state is the time-crystalline phase, with
discrete time translation symmetry. Indeed, as predicted by Renaud
Parentani himself during our visit to Orsay in 2015, the oscillation
frequency $\omega_0$ of a CES state exhibits a critical behavior close
to the phase transition, where the critical exponents
$\alpha_q,\alpha_Z$ for $q,Z$ (obtained from a fit in
Figure~\ref{fig:Critical}b) are both approximately $\alpha_q\simeq
\alpha_Z\simeq 0.50$, strongly suggesting a possible analytical\break
derivation.

\begin{figure}
\includegraphics{fig10}
\caption{\label{fig:Critical}(a) Dynamical phase diagram as a function of $(Z,q)$ for the
final state of an initially subsonic flowing condensate
$\Psi(x,0)=\mathrm{e}^{\mathrm{i}qx}$ in which an attractive delta barrier
$V(x)=-Z\delta(x)$ is quenched at $t=0$. (b) Critical behavior of the
CES frequency $\omega_0$ close to the phase transition along the green
lines in (a), where the red line represents a fit to a power law. Main
panel: Velocity dependence. Inset: Delta-strength
dependence.}
\end{figure}


Another remarkable feature of a spontaneous Floquet state is that it
conserves energy due to the time-independence of the underlying
Hamiltonian, unlike conventional Floquet systems. After developing the
so-called $(t,\phi)$ formalism within the generalized Gibbs
ensemble~\cite{Cazalilla2006,Rigol2007,Langen2015},
Ref.~\cite{deNova2024} showed that   spontaneous Floquet states
have a unique conserved magnitude, the Floquet charge $F$, whose
conjugate thermodynamic variable is the frequency $\omega$. This allows
to identify spontaneous Floquet states as isofloquetic, conserving the
total energy, and driven Floquet states as isoperiodic, conserving the
Floquet enthalpy, $I=E-\omega F$, in analogy with isochoric and
isobaric systems, respectively. This characterization gives rise to the
so-called Floquet thermodynamics, which describes Floquet systems with
the same thermodynamical tools as in stationary states.

The quantum fluctuations of a spontaneous Floquet state in an atomic
condensate are described by the BdG equations resulting from the
expansion $\hat{\Psi}=[\Psi_0(x,t)+\hat{\varphi}(x,t)]
\mathrm{e}^{-\mathrm{i}\tilde{\mu}t}$,\break namely
{\begin{equation}\label{eq:BdGfieldequation}
\mathrm{i}\partial_t\hat{\Phi}=M_0(t)\hat{\Phi},\quad 
M_0(t)=\left[\begin{array}{@{}cc@{}} N_0(t) & A_0(t)\\
-A_0^*(t) &-N_0^*(t)\end{array}\right],
\end{equation}}\unskip
where
{\begin{eqnarray}
 N_0(t)&=&-\frac{\partial_x^2}{2}+V(x)+2g(x)|\Psi_0(x,t)|^2-\tilde{\mu},\nonumber\\
A_0(t)&=&g(x)\Psi_0^{2}(x,t)
\end{eqnarray}}\unskip
are periodic operators. Consequently, the BdG matrix $M_0(t)$ is a
periodic linear operator, described by conventional Floquet theory, so
its spectrum is given in terms of quasi-energy bands,
{\begin{equation}\label{eq:BogoliubovFloquetSectorTime}
    [M_0(t)-\mathrm{i}\partial_t]z_{\varepsilon,\nu}(t)=\varepsilon z_{\varepsilon,\nu}(t),
\end{equation}}\unskip
with $z_{\varepsilon,\nu}(\mathbf{x},t+T)=z_{\varepsilon,\nu}(\mathbf{x},t)$, 
$\varepsilon$ the quasi-energy (defined again modulo $\omega_0$), and
$\nu$ a discrete index labeling the solution. This BdG expansion also
results from a conventional Floquet state, which takes the same form of
Equation~(\ref{eq:FloquetWaveMu}) but the period $T$ there is imposed
by the external driving, while for a spontaneous Floquet state $T$ is
spontaneously chosen by the system.

Of particular interest is the presence of Nambu--Goldstone (NG) modes.
In a stationary context, if the GP wavefunction spontaneously breaks
one of the continuous symmetries of the Hamiltonian in such a way that,
if $\Psi_0(x)$ is a stationary GP solution, then
$\Psi_\alpha(x)=\mathrm{e}^{-\mathrm{i}\alpha G}\Psi_0(x)$ 
is another stationary GP
solution, a zero-energy NG mode emerges in the BdG spectrum:
{\begin{equation}\label{eq:NambuGoldstone}
    M_0z_\alpha=0,\quad z_\alpha=\left[\begin{array}{@{}c@{}} \partial_\alpha \Psi_0 \\
\partial_\alpha \Psi^*_0\end{array}\right]=\left[\begin{array}{@{}c@{}} -\mathrm{i}G \Psi_0 \\
\mathrm{i} (G\Psi_0)^*\end{array}\right].
\end{equation}}\unskip
This can be proven by expanding the symmetry transformation to linear
order in $\alpha$, where $G$ is the infinitesimal generator of the
transformation. Since NG modes have zero norm, their amplitude does not
behave as an annihilation operator but instead as a coordinate
operator, with a conjugate momentum that describes the fluctuations of
the conserved charge $Q$ associated to the broken continuous
symmetry~\cite{Lewenstein1996,Dziarmaga2004}. A major example is the
Goldstone mode $z_\theta$ corresponding to the $U(1)$-symmetry
breaking, see Equation~(\ref{eq:TIBdGSource}), resulting from the fact
that $\mathrm{e}^{-\mathrm{i}\theta}\Psi_0$ is also a stationary GP solution for
arbitrary $\theta$. When several symmetries are spontaneously broken,
the quantization procedure is elegantly described by a geometric
formalism involving the so-called Berry-Gibbs connection, which is the
Berry connection associated to the GP wavefunction, the continuous
parameters of the manifold being the conserved charges associated to
the broken symmetries~\cite{deNova2024}.

For Floquet states, both spontaneous and conventional, spontaneous
symmetry breaking is translated into the emergence of
Floquet--Nambu--Goldstone (FNG) modes with zero quasi-energy. In a
condensate, this means $[M_0(t)-\mathrm{i}\partial_t]z_\alpha(t)=0$, where
$z_\alpha(t)$ is now periodic. In the specific case of a spontaneous
Floquet state, a genuine temporal FNG mode arises from the spontaneous
symmetry breaking of time-translational invariance~\cite{deNova2024}.
This is the characteristic hallmark of a \textit{bona-fide} spontaneous
Floquet state, reflecting also its time-crystalline nature, which
distinguishes it from trivial periodic behavior such as traveling
soliton waves in a ring, which do not possess a proper temporal FNG
mode. The CES state indeed has a genuine temporal FNG mode, stemming
from the fact that $\Psi_0(x,t+t_0)\mathrm{e}^{-\mathrm{i}\tilde{\mu}t}$ is also a
solution of the time-dependent GP
equation~(\ref{eq:HamiltonianEffective}) for arbitrary $t_0$. 

Interestingly, the quantum amplitude of the temporal FNG mode provides
a unique realization of a time operator in quantum mechanics, which
commutes with the linear fluctuations of the grand-canonical energy, no
longer vanishing since $\Psi_0$ is not here an extreme of the
grand-canonical energy as it is time-dependent~\cite{deNova2024}. In
general, quantum mechanics forbids the existence of a time operator
$\hat{T}$ since the Hamiltonian should be its canonical conjugate,\break
$[\hat{T},\hat{H}]=-\mathrm{i}$. This implies that
{\begin{equation}
\mathrm{e}^{-\mathrm{i}\hat{T}E_0}\hat{H}\mathrm{e}^{\mathrm{i}\hat{T}E_0}=\hat{H}-E_0.
\end{equation}}\unskip
for
arbitrary value of $E_0$. Hence, if $E$ is an eigenenergy, then $E-E_0$
is also an eigenergy, which contradicts the fact that Hamiltonians are
bounded from below by the ground-state energy.

Nevertheless, this problem is similar to that of the phase
operator~\cite{Mora2003}. If a phase operator $\hat{\theta}$ exists,
then it satisfies $[\hat{N},\hat{\theta}]=\mathrm{i}$, which implies
{\begin{equation}
\mathrm{e}^{-\mathrm{i}\hat{\theta}N_0}\hat{N}\mathrm{e}^{\mathrm{i}\hat{\theta}N_0}=\hat{N}-N_0
\end{equation}}\unskip
for arbitrary value of $N_0$. Hence, if $N$ is an
eigenvalue of $\hat{N}$, then $N-N_0$ is also, violating two
fundamental properties of the spectrum of $\hat{N}$: its positive
definiteness and its discreteness.  However, it is well known that one
can define phase fluctuations in condensates involving large particle
numbers, where one can neglect the discreteness of $N$ and states with
low occupation number. In analogy, by identifying energy with particle
number and time with phase, one can define time fluctuations for large
energies well above the ground state. This correspondence indicates
that spontaneous Floquet states must be highly excited states, shifted
by a macroscopic energy from the ground state. This is indeed the case
of the CES state, whose energy is\break shifted $\Delta E\simeq Nq^2/2$ 
above that of a condensate at rest. 

\section{Discussion}\label{sec:conclu}

Resonant configurations represent a rich paradigm in analogue gravity.
In addition to the universal thermal behavior at low frequencies, they
display a highly non-thermal peaked structure in the Andreev and
Hawking spectra, since they act as a Fabry--Perot resonator for the
negative-energy partners of the Andreev--Hawking effect. This contrasts
with standard analogue configurations, whose spectra have a marked
thermal character which is easy to misidentify with any other
background thermal component, as in the real astrophysical scenario.
Another interesting feature of resonant configurations is that they
highly enhance the Andreev signal, which can even overcome the Hawking
one, while in standard analogue configurations the former is typically
highly suppressed.

A feasible experimental scheme to implement a resonant configuration
using atomic condensates relies on outcoupling a large boson reservoir
through an optical lattice, eventually achieving a quasi-stationary
black-hole configuration. Remarkably, the optical lattice acts as a
low-pass filter of Andreev--Hawking radiation~\cite{deNova2017b},
something that could have potential applications in quantum transport
and atomtronics~\cite{Amico2021}, further motivating the interest for
its experimental implementation. An alternative approach to achieve a
non-thermal spectrum  is provided by dipolar condensates, where the
departure from thermality is induced by the presence of a roton minimum
in the dispersion relation~\cite{Ribeiro2023}. 


By applying concepts and techniques originally derived in quantum
optics, we can understand the Andreev and Hawking effects as a joint
phenomenon, resulting from the non-degenerate parametric amplification
of the outgoing modes, where the hybrid Andreev--Hawking channel is the
signal, and the anomalous channel is the idler. Regarding quantum
correlations, we have analyzed the occurrence of violations of
Cauchy--Schwarz inequalities and entanglement, which are equivalent
conditions for a broad and relevant class of quantum states. We have
observed that resonant configurations highly enhance entanglement near
the resonant peaks of the spectrum for both the Andreev and Hawking
effects, allowing for its detection even at high temperatures
comparable to the chemical potential. Thus, they improve the
performance of standard analogue configurations, where entanglement is
highly attenuated with temperature, fading away at low temperatures in
the Andreev case. 

The characterization of quantum correlations in the Andreev--Hawking
effect, including tripartite entanglement~\cite{Isoard2021} and Bell
non-locality~\cite{Ciliberto2024}, is still an active topic of
research, which could lead to potential applications in quantum
technologies, since an analogue horizon behaves as a spontaneous source
of entangled phonons. 

An interesting spin-off of the study of quantum correlations in
analogue gravity is the research on quantum information in high-energy
colliders~\cite{Afik2021}, which is rapidly becoming a whole topic of
research by itself. It has already led to the first observation of
entanglement in quarks, in turn the highest-energy entanglement
detection ever, by the ATLAS and CMS
collaborations~\cite{ATLAS2024,CMS2024}. This observation paves the way
to use high-energy colliders for the study of foundational quantum
problems, something of great interest due to their genuine relativistic
nature and fundamental character, operating at the current frontier of
knowledge in Physics. In fact, a number of experimental analyses
searching for genuine quantum signatures at the LHC are currently
ongoing.  

A most important phenomenon arising in resonant configurations is the
black-hole laser effect. For its discussion, we have separated the
three main stages of its time evolution. At short times, the dynamics
is governed by the spectrum of dynamical instabilities in the linear
BdG equations. By analyzing several BHL models, we have observed the
generality of the original predictions by Michel and
Parentani~\cite{Michel2013}, namely: (i) the lasing modes emerge as
degenerate at critical lengths equispaced by the BCL wavelength,
becoming non-degenerate at halfway between the critical lengths; and
(ii) there is a perfect correspondence between the emergence and later
degeneracy breaking of the lasing modes, and the emergence of
stationary GP solutions in the non-linear spectrum. Our results also
confirm the conjecture of Michel and Parentani~\cite{Michel2015}: in
flowing scattering configurations, energetic and dynamical instability
are equivalent conditions, and the only stable solution is the ground
state, which evaporates all the acoustic horizons to become fully
subsonic.

At intermediate times, we have theoretically studied the BHL--BCL
crossover, originally characterized in Ref.~\cite{deNova2023}. By
invoking the analogy with an unstable
pendulum~\cite{Leonhardt2003,Finazzi2010,Burkle2018,deNova2023}, three
regimes can be identified: quantum BHL, classical BHL, and BCL. Their
most characteristic trait is their efficiency as quantum amplifiers: a
quantum BHL is a non-linear quantum amplifier, increasing quantum
fluctuations up to the same saturation amplitude regardless of their
initial strength, while classical BHL and BCL are linear quantum
amplifiers, where the output is proportional to the input. In
particular, quantum amplification in a classical BHL is exponentially
large in the lasing time, and much larger than in the BCL regime, where
there is no microscopic amplification mechanism, and the amplification
just stems from the strong background modulation induced by the BCL
wave. The characterization of each regime as a quantum amplifier can be
complemented with a qualitative analysis of the monotonicity of their
growth rate and their quantum gain, providing practical experimental
criteria for the unambiguous detection of the BHL effect, a major
remaining challenge in the analogue field. Furthermore, our analysis
identifies the classical BHL regime as the most reachable target, where
the BCL wave fuels the BHL effect by providing it with a classical
seed, instead of\break undermining it.

Remarkably, the study of the BHL--BCL crossover in
Ref.~\cite{deNova2023} has also allowed to identify novel analogue
phenomena such as Hawking-stimulated white-hole (HSWH) radiation at the
start of the BHL process (when the partner modes of the Andreev--Hawking
effect stimulate the continuous spectrum of white-hole
radiation~\cite{Mayoral2011}), or quantum BCL-stimulated Hawking
radiation (the spontaneous resonant Hawking radiation above the
non-linear saturated BCL wave). The analysis of the BHL--BCL 
crossover can be of interest for
other analogue setups in which low-frequency undulations similar to the
BCL wave hinder the BHL
effect~\cite{Coutant2012,Coutant2014,Bossard2023}. In general, a
quantum/classical BHL provides an ideal testing ground for the study of
quantum/classical
backreaction~\cite{Balbinot2005a,Patrick2021,Baak2022,Butera2023}.
Apart from its intrinsic interest for the analogue field, a BHL behaves
as a quantum amplifier, which could have potential applications in
atomtronics.

At long times, a black-hole laser exhibits a dynamical phase diagram
with two states: the ground state, with continuous time-translation
symmetry, and the CES state, with discrete time-translation symmetry,
resulting from its periodic nature. Indeed, the CES state is a
universal feature of a flowing condensate, representing a particular
example of the much more general concept of spontaneous Floquet
state~\cite{deNova2022}: a Floquet state arising from a
time-independent Hamiltonian, whose periodicity is spontaneously set by
many-body interactions.

Spontaneous Floquet states are by themselves a novel non-equilibrium
paradigm. For instance, they conserve energy, in contrast to
conventional Floquet states, which in turn conserve the so-called
Floquet enthalpy as they arise from periodically driven Hamiltonians,
operating at fixed frequency. These conserved magnitudes allow for a
thermodynamic description of Floquet states completely analogous to
that of stationary states, which has been labeled Floquet
thermodynamics~\cite{deNova2024}. 

In addition, spontaneous Floquet states spontaneously break continuous
time-translation symmetry, representing a specific realization of a
continuous time crystal. This results in the emergence of a genuine
temporal Floquet--Nambu--Goldstone mode with zero quasi-frequency, whose
amplitude provides a unique realization of a time operator in a
tangible condensed-matter setup~\cite{deNova2024}. We note that the
construction of a time operator is a fundamental subject in quantum
mechanics~\cite{Aharonov1961,Susskind1964,Unruh1989,Honh2021}.
Therefore, the identification of a time operator in an analogue gravity
setup provides a rich scenario that could lead to fundamental research
on the quantum foundations of spacetime.

Although our discussion is restricted to atomic condensates, the
results of this work can be easily translated to optical and
polaritonic analogues due to the similarity of the equations of motion.
For instance, Andreev reflection has been studied in polaritonic
condensates~\cite{Septembre2021,Septembre2023}. The excitation of a
quasi-normal mode from vacuum fluctuations has been recently
numerically observed in a resonant polaritonic
configuration~\cite{Jacquet2023}. Quantum correlations in optical and
polaritonic analogues have been also
studied~\cite{Busch2014a,Agullo2022,Brady2022,Delhom2024}. Similar
periodic states to the CES state are predicted for
polaritons~\cite{Opala2018}.

As a global remark, it must be noted that the low-pass filter of
Andreev--Hawking radiation provided by an optical lattice, the
stationary source of entangled phonons provided by the\break Hawking
effect~\cite{Kolobov2021}, the behavior of a black-hole laser as a
quantum amplifier, and the spontaneous Floquet state represented by the
CES state, demonstrate the potential of interdisciplinary applications
of analogue gravity concepts. This is further supported by the
establishment of a line of research on quantum information in
high-energy colliders, conceptually based on the study of quantum
correlations in Andreev--Hawking radiation.

We would like to conclude by emphasizing the profound influence
provided by the intellectual leadership of Renaud Parentani in the
shaping of the field of analogue gravity. His ideas and inspiration
permeate any coherent narration of the evolution of this research
field, and in particular are ubiquitous in all the results discussed in
the present article.

\section*{Declaration of interests}

The authors do not work for, advise, own shares in, or receive funds
from any organization that could benefit from this article, and have
declared no affiliations other than their research organizations.

\section*{Acknowledgments}
We devote this article to the memory of Renaud Parentani. This work has
received funding from European Union's Horizon 2020 research and
innovation programme under the Marie Sk\l{}odowska-Curie Grant
Agreement No. 847635, from Spain's Agencia Estatal de Investigaci\'on
through Grant No. PID2022-139288NB-I00, and from Universidad
Complutense de Madrid through Grant No. FEI-EU-19-12.

\CDRGrant[EU]{847635}
\CDRGrant[AEI]{PID2022-139288NB-I00}
\CDRGrant[UCM]{FEI-EU-19-12}

\back{}

\bibliographystyle{crunsrt}
\bibliography{crphys20240461}
\refinput{crphys20240461-reference.tex}

\end{document}
