\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 ($vc$). 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_20$, 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, \end{array}\right. \end{equation}}\unskip with $c_2,c_3 |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\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_2L_{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]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}