Black-hole lasing in Bose-Einstein condensates: analysis of the role of the dynamical instabilities in a nonstationary setup

TL;DR

This paper investigates how dynamical instabilities influence black-hole lasing phenomena in Bose-Einstein condensates, emphasizing nonstationary effects and distinguishing between classical noise and quantum fluctuations in experimental setups.

Contribution

It introduces a toy model to analyze nonstationary effects on black-hole lasing, revealing the roles of dynamical instabilities and differentiating classical noise from quantum fluctuations.

Findings

Dynamical instabilities significantly affect black-hole lasing behavior.

Nonstationary effects alter the characterization of stimulated Hawking radiation.

Classical noise impacts the amplification process differently than quantum fluctuations.

Abstract

We present a theoretical study on the origin of some findings of recent experiments on sonic analogs of gravitational black holes. We focus on the realization of a black-hole lasing configuration, where the conclusive identification of stimulated Hawking radiation requires dealing with the implications of the nonstationary character of the setup. To isolate the basic mechanisms responsible for the observed behavior, we use a toy model where nonstationarity can be described in terms of departures from adiabaticity. Our approach allows studying which aspects of the characterization of black-hole lasing in static models are still present in a dynamical scenario. In particular, variations in the role of the dynamical instabilities can be traced. Arguments to conjecture the twofold origin of the detected amplification of sound are given: the differential effect of the instabilities on the mean field and on the quantum fluctuations gives some clues to separate a deterministic component from self-amplified Hawking radiation. The role of classical noise, present in the experimental setup, is also tackled: we discuss the emergence of differences with the effect of quantum fluctuations when various unstable modes are relevant to the dynamics.

Equations

$$\hat{\Psi }={\Psi }_0+\widehat{\delta \Psi },$$

$$i\hslash \frac{\partial {\Psi }_0(\mathbf{r},t)}{\partial t}=\left[-\frac{{\hslash }^2}{2m}{\nabla }_{\mathbf{r}}^2+V_{ex}(\mathbf{r},t)+g{\left|{\Psi }_0(\mathbf{r},t)\right|}^2\right]{\Psi }_0(\mathbf{r},t),$$

$$\hat{\Psi }(x,t)$$

$$\left[-\frac{{\hslash }^2}{2m}\frac{{\partial }^2}{\partial x^2}+V_{ex}(x,\mathbf{\Lambda }(t))+g{\left|{\Phi }_{0,ad}\right|}^2-\mu (\mathbf{\Lambda }(t)\right]{\Phi }_0^{ad}(x,\mathbf{\Lambda }(t))=0.$$

$$\widehat{\mathbf{\delta }\mathbf{\Phi }}\equiv \left(\begin{array}{c}\widehat{\delta \Phi }\\ {\widehat{\delta \Phi }}^{†}\end{array}\right),$$

$$\mathbf{\delta }\mathbf{\Phi }\equiv \left(\begin{array}{c}\delta \Phi \\ \delta {\Phi }^{\ast }\end{array}\right).$$

$$\frac{\partial \mathbf{\delta }\mathbf{\Phi }}{\partial t}=-\frac{i}{\hslash }{\mathcal{L}}_{\mathbf{BdG}}^{\mathbf{ad}}\mathbf{\delta }\mathbf{\Phi }+\mathbf{S},$$

$${\mathcal{L}}_{\mathbf{BdG}}^{\mathbf{ad}}\equiv \left(\begin{array}{cc}{H}_0-\mu +2g{\left|{\Phi }_0^{ad}\right|}^2& g{\left({\Phi }_0^{ad}\right)}^2\\ -g{\left({\Phi }_0^{ad\ast }\right)}^2& -\left(H_0-\mu +2g{\left|{\Phi }_0^{ad}\right|}^2\right)\end{array}\right),$$

$$\mathbf{S}\equiv \left(\begin{array}{c}S\\ -S^{\ast }\end{array}\right),$$

$$S=-i\hslash \dot{\mathbf{\Lambda }}\frac{d{\Phi }_{0,ad}}{d\mathbf{\Lambda }},$$

$$\widehat{\delta \Phi }={\widehat{\delta \Phi }}^{ad}+\delta {\Phi }_0^{nad}.$$

$${\widehat{\delta \Psi }}^{ad}(x,t)$$

$${\mathbf{V}}_{\mathbf{a}}(x,\mathbf{\Lambda }(t))$$

$${\mathbf{Z}}_{\mathbf{a}}(x,\mathbf{\Lambda }(t))$$

$${\vartheta }_a(t)={\int }_0^t{\lambda }_a(\mathbf{\Lambda }(t^{\prime }))d{t}^{\prime }.$$

$$\frac{d\mathbf{U}}{dt}=-\frac{i}{\hslash }{\mathcal{L}}_{\mathbf{BdG}}^{\mathbf{ad}}\mathbf{U}.$$

$$\mathbf{\delta }{\mathbf{\Phi }}_0^{\mathbf{nad}}=\mathbf{U}(t){\int }_0^t{\mathbf{U}}^{-1}(t^{\prime })\mathbf{S}(t^{\prime })d{t}^{\prime }.$$

$${\mathbf{U}}_{jk}=\exp \left(-i{\int }_0^t{\lambda }_j(t^{\prime })d{t}^{\prime }\right){\delta }_{jk}=\exp \left(-i{\vartheta }_j(t)\right){\delta }_{jk}$$

$$({\mathbf{U}}^{-1}{)}_{jk}=\exp \left(i{\int }_0^t{\lambda }_j(t^{\prime })d{t}^{\prime }\right){\delta }_{jk}=\exp \left(i{\vartheta }_j(t)\right){\delta }_{jk}.$$

$$\delta {\Psi }_0^{nad}(x,t)$$

$$B_a(t)$$

$$C_a(t)$$

$$\hat{\rho }={\rho }_0^{ad}+{\hat{\rho }}_1^{ad}+{\rho }_1^{nad},$$

$${\hat{\rho }}_1^{ad}=({\Phi }_0^{ad}{\widehat{\delta \Phi }}^{ad†}+\text{h.c.}),$$

$${\rho }_1^{nad}=\left[{\Phi }_0^{ad}\delta {\Phi }_0^{nad\ast }+\text{c.c.}\right].$$

$$⟨0|{\rho }_0^{ad}+{\hat{\rho }}_1^{ad}+{\rho }_1^{nad}|0⟩$$

$$⟨0|{\hat{\rho }}_1^{ad}|0⟩=0.$$

$$⟨0|\hat{\rho }|0⟩\sim {\rho }_0^{ad}(x,\mathbf{\Lambda }(t))\left[1+2\text{Re}\left\{B_a(t){\sigma }_a(x,\mathbf{\Lambda }(t))\right\}\right],$$

$$⟨0|{\hat{\rho }}_1^{ad}(x,t){\hat{\rho }}_1^{ad}(x^{\prime },t)|0⟩$$

$${\Delta }_a(t)={\int }_0^t{\Gamma }_a(\mathbf{\Lambda }(t^{\prime }))d{t}^{\prime }.$$

$$\mathbf{\delta }{\mathbf{\Psi }}^{\mathbf{cn}}\equiv \left(\begin{array}{c}\delta {\Psi }^{cn}\\ \delta {\Psi }^{cn\ast }\end{array}\right),$$

Full text

I Introduction

The predictions on Hawking radiation (HR), i.e., on spontaneous emission of thermal radiation from a black hole (BH)key-Hawking1 ; key-Hawking2 , have not had observational verification. The low temperatures of emission constitute the main handicap to the observation. Yet, the fundamental character of the involved physics allows searching for alternative strategies to test the theory via the detection of similar effects in parallel systems key-Unruh . Specially promising are the proposals for building sonic analogs of gravitational BHs with atomic Bose-Einstein condensates (BECs) key-Garay1 ; key-Garay2 ; key-Barcelo ; key-Gardiner . The basis of the proposals is the creation of a sonic BH horizon, where the character of a flowing condensate changes from subsonic to supersonic. As in the gravitational context, the quantum treatment of the field, in contrast with the classical description, predicts the emission of radiation from the supersonic region key-Carusotto2 ; key-Carusotto3 ; key-MacherParentani . Given the low temperatures and the access to measure in any region in the condensates, the proposed schemes can be expected to have practical applicability. Indeed, a first realization was reported in key-Steinhauer3 . Advances in this line have been achieved recently with the implementation of a black-hole lasing (BHL) configuration key-Steinhauer1 . In this variation of the basic scheme, a condensate flow is made to cross twice the speed of sound. Its characteristics in a stationary regime have been analyzed in previous theoretical work. The presence of the second (white-hole) (WH) horizon has been predicted to lead to amplification of spontaneous HR through a resonance mechanism similar to that of a lasing cavity key-Corley . The analyses are based on the use of the Bogoliubov–de Gennes (BdG) approach to describe the excitations key-Michel1 ; key-Carusotto4 ; key-Finazzi1 . The appearance of complex eigenvalues of the BdG operator marks the emergence of lasing: the instabilities, associated with the imaginary parts of the eigenfrequencies, induce the amplification of the perturbations, in particular, the self-amplification of the quantum fluctuations. As the amplification is frequency dependent, the thermal character of the spontaneous HR is lost in BHL. In the experiments key-Steinhauer1 , the BH horizon was generated by displacing a step-like potential along an atomic BEC initially at rest. In this form, the fluid was accelerated from a subsonic to a supersonic velocity in the reference frame of the step. The subsequent deceleration, due to the effect of the trap, led to the second (WH) horizon. The observed magnitudes were the density and the nonlocal density correlation function (DCF) in the supersonic region. The evolution and spatial dependence of these observables have been compared with classical-field simulations and predictions of studies on static models. Interestingly, despite the nonstationary character of the implementation, the static picture seems to provide valuable insight into the underlying mechanisms. For instance, the concept of dynamical instability has been used to interpret the measured increase in the DCF. From those analyses, there is broad agreement on the detection of lasing. Still, there is some debate on whether the amplified radiation is actually rooted in quantum noise. (Alternative interpretations tracing the origin to classical noise or to a deterministic seed have been formulated from numerical simulations key-Carusotto1 ; key-Clark ). A related question which requires a detailed explanation is the growth of the mean density, absent in the predictions for a static model. (In the analysis of this issue, the role of the nonlinear backreaction of the excited modes on the condensate has been considered key-Michel2 ). Also, the spectral structure of the DCF must be clarified. Here, we aim at establishing a link between the (static) theoretical approaches and the (dynamical) experimental realization. As the questions that have been the subject of debate have fundamental character, we have opted for tackling them in a simple scenario, where the basic physical mechanisms can be isolated. In particular, we deal with a toy model where the corrections to an adiabatic approximation can already give the clues to the effects of nonstationarity, and, in turn, to tracing the origin of some of the observed features. The predictions of this simplified description will be tested through a numerical simulation of the experiment using a more complete model.

The outline of the paper is as follows. In Sec. II, we present our model system. A nonadiabatic approach with a generic set of variable parameters is developed in Sec. III. In Sec. IV, we evaluate the effects of nonadiabaticity on the density and on the DCF. Differential effects of the instabilities on quantum fluctuations and diverse types of classical noise are tackled in Sec. V. In Sec. VI, we go beyond the postadiabatic scenario: numerical results are presented for the simulation of the experimental setup with a model with no restrictions on the time-variation regime. Finally, some general conclusions are summarized in Sec. VI. (To give a self-contained presentation, we summarize some results of previous descriptions of BHL in static systems in Appendix A. Those results are taken as starting point in our study).

II The model system

We consider an atomic BEC in a confining potential $V_{ex}(\mathbf{r},t)$. We assume that the field operator of the condensate $\hat{\Psi}$ can be separated as

[TABLE]

where $\Psi_{0}$ is a classical field and $\hat{\delta\Psi}$ is a perturbative quantum contribution. The mean-field description of the classical component is given by the time-dependent Gross-Pitaevskii (GP) equation key-Stringari

[TABLE]

where $m$ is the mass of a condensate atom, and $g$ is the strength that characterizes the atom-atom interaction. It is assumed that the evolution of the quantum term $\hat{\delta\Psi}$ can be analyzed via a linearization procedure.

In the experimental realization, $V_{ex}(\mathbf{r},t)$ was conveniently varied in order to implement the black-hole laser configuration, i.e., to achieve two crossing points between the flow velocity and the speed of sound. Specifically, a step potential was displaced along the condensate in the harmonic trap. The consequent acceleration (in the frame of the step) of the atomic flow led to the emergence of the BH horizon. Subsequently, as the trapping effect set in, the flow was decelerated and the inner (WH) horizon was formed. Appropriate for implementing an adiabatic approximation is to formally incorporate the time variation in the external potential through a set of parameters $\boldsymbol{\Lambda}(t)$. Therefore, we rewrite the confining term as $V_{ex}(\mathbf{r},\boldsymbol{\Lambda}(t))$. Additionally, a dimensional reduction, allowed by the characteristics of the harmonic trap, is applied: the strong confinement in the directions transversal to the step displacement can be incorporated through effective parameters in the equation for the reduced dynamics in the longitudinal $x$-direction. [The effective monodimensional version of the trapping potential will be denoted as $V_{ex}(x,\boldsymbol{\Lambda}(t))$].

A crucial aspect of the practical arrangement is the nonstationary character of the flow generated by the changing potential. Previous work on the characterization of instabilities has been carried out on models where stationary flows and permanent horizon-schemes are assumed. In the present case, the whole process of horizon formation presents nontrivial dynamical aspects, which constitute a handicap to the interpretation of the findings. Actually, there is an open debate on the relevance of different mechanisms to the observed behavior key-Carusotto1 ; key-Clark ; key-SteinhauerNova . The discussion has incorporated arguments developed from various numerical simulations. Those studies, based on solving the GP equation, have included different components of the setup. Early results for a completely classical field were argued to reproduce the detected amplification. The appearance and growth of undulations in the mean density were reported. Additionally, the DCF, obtained with added classical noise, seemed to evolve in agreement with the observations. The similar spatial patterns of the density and of the DCF pointed to a common origin for the growth of the mean values of both observables. A global conclusion was that the emergence of amplification in a purely classical context ruled out the necessary identification of the observed radiation with self-amplified HR key-Carusotto1 . In contrast with this picture, subsequent numerical work, which included elements simulating quantum fluctuations, revealed qualitative differences between the responses of the system to quantum and classical noise key-SteinhauerNova . Quantum fluctuations were found to be necessary to account for salient spectral features. Also, it was argued that, since the density undulation (the ripple) appeared before the WH formation and evolved steadily with no qualitative changes, its origin must be different from that of the DCF evolution, linked initially to spontaneous HR, and, developing subsequently (different) characteristics associated with BHL. The questions raised by this debate show the convenience of presenting additional arguments to conform a conclusive interpretation of the experiments. Here, we work with a simplified version of the setup where the roles of different components of the dynamics can be singled out. To deal with implications intrinsic to nonstationarity, we consider a working regime where the adiabatic approximation gives the zero-order description of the dynamics and the departures from adiabaticity are taken as corrections. Additionally, the quantum fluctuations from the adiabatic stationary solutions are assumed to be well described by the BdG approach. The presence of classical fluctuations is also considered. The experimental results, satisfactorily reproduced by some of the numerical studies key-Carusotto1 ; key-Clark ; key-SteinhauerNova , seem to indicate a more regular and slower evolution of the flow once the horizons have been formed. From this feature we extract some clues to simplify the approach: instead of attempting to describe the whole process of the emergence of the BHL, we will focus on the system once the BHL regime has been reached. Moreover, as it will be specified further on, in order to avoid the breakdown of the adiabatic approximation key-Anglin , we will concentrate on a temporal range where no changes in the unstable or stable character of the eigenstates of the BdG operator take place. Although the evolution in the practical setup cannot be assumed to correspond to a perturbed adiabatic regime, the study of nonadiabaticity as a correction in our model will be shown to uncover some basic mechanisms potentially relevant to general time variations in the BHL configuration. Moreover, even though our approach cannot account for the potential transition of the emitted radiation from spontaneous to self-amplified HR, it can trace the differences with the characteristics of the density growth.

III Non-adiabatic approach to a non-stationary black-hole lasing configuration

Following the above considerations, we use for the field operator in Eq. 1 the ansatz

[TABLE]

where $\Phi_{0}^{ad}(x,\boldsymbol{\Lambda}(t))$ is the adiabatic wavefunction, i.e., the solution to the time-independent Gross-Pitaevskii equation for a frozen set of parameters $\boldsymbol{\Lambda}(t)$

[TABLE]

$\mu(\boldsymbol{\Lambda}(t)$ is the related chemical potential. Moreover, $\hat{\delta\Phi}(x,t)$ is the perturbative quantum term. From it, we define the two-component field $\mathbf{\boldsymbol{\hat{\delta\Phi}}}$ as

[TABLE]

and, following the standard procedure to solve for it (see the Appendix), we work initially with its classical counterpart

[TABLE]

It is straightforwardly shown that, to first order, $\boldsymbol{\delta\Phi}$, as does $\boldsymbol{\hat{\delta\Phi}}$, obeys the equation

[TABLE]

where two contributions to the evolution can be differentiated. First, the homogeneous part, characterized by the operator

[TABLE]

where $H_{0}=-\frac{\hbar^{2}}{2m}\frac{\partial^{2}}{\partial x^{2}}+V_{ex}(x,\boldsymbol{\Lambda}(t))$, corresponds to the evolution that would experience the perturbative term in a strictly adiabatic regime. Note that it parallels the equation for the perturbation in a static scenario.

Second, the source matrix

[TABLE]

with

[TABLE]

specifically incorporates departures from adiabaticity. This term affects the dynamics irrespective of the system preparation. Actually, it constitutes a seed for perturbations.

In the following, we will build up the general solution for the perturbative field operator $\hat{\delta\Phi}$ as the sum of a general solution to the homogeneous equation, which will be denoted as $\hat{\delta\Phi}^{ad}$, and a particular solution to the complete equation, represented by $\delta\Phi_{0}^{nad}$. Hence, we will write

[TABLE]

III.1 The homogeneous equation

The parallelism existent between the adiabatic BdG operator $\boldsymbol{\mathcal{L}_{BdG}^{ad}}$ and its counterpart $\boldsymbol{\mathcal{L}_{BdG}}$ for a strictly static system allows applying the methods developed for the case of a stationary flow key-Leonhardt ; key-Finazzi1 , summarized in the Appendix, to solve the homogeneous equation. Accordingly, we will expand the field in the eigenmodes of $\boldsymbol{\mathcal{L}_{BdG}^{ad}}$. As we are interested primarily in the effect of the instabilities on the quantum fluctuations, we work with the quantum form of the expansion and retain in it only the contribution of the unstable modes. Therefore, we write

[TABLE]

where the operators $\hat{b}_{a}$ ($\hat{b}_{a}^{\dagger}$) and $\hat{c}_{a}$($\hat{c}_{a}^{\dagger}$) respectively correspond to the unstable modes

[TABLE]

with respective eigenfrequencies $\lambda_{a}(\boldsymbol{\Lambda}(t))$ and $\lambda_{a}^{*}(\boldsymbol{\Lambda}(t))$ [$\lambda_{a}(\boldsymbol{\Lambda}(t))=\omega_{a}(\boldsymbol{\Lambda}(t))+i\Gamma_{a}(\boldsymbol{\Lambda}(t))$, see the Appendix]. Moreover, we have used

[TABLE]

It is assumed that an adiabatic approximation for the wave functions of the modes is feasible. Indeed, we can regard the nonadiabatic corrections to the wavefunctions of the unstable modes as having a secondary effect on the dynamics compared with the role played by the complex character of their eigenvalues. As the applicability of this approximation can be jeopardized by changes in the stability of the modes, we will focus on a temporal range where no changes in the unstable or stable character of the eigenmodes of $\boldsymbol{\mathcal{L}_{BdG}^{ad}}$ take place. Despite its apparent restrictive character, this sound time regime will be shown to be appropriate to uncover general characteristics of the mechanisms underlying the observed features.

To account, in the proposed adiabatic framework, for the evolution of a classical perturbation, e.g., of classical noise in the initial preparation, we can use the classical version of Eq. (12), where the operators are replaced by $\mathbb{C}$-functions corresponding to the projections of the perturbation on the unstable modes. This analysis is postponed to Sec. V.

III.2 The source term

Convenient to deal with the inhomogeneous term in Eq. (7) is the use of the time-evolution operator $\boldsymbol{U}(t)$ associated with $\boldsymbol{\mathcal{L}_{BdG}^{ad}}$, i.e., of the operator defined through the equation

[TABLE]

Using it, one can readily show that a particular solution to Eq. (7) is given by

[TABLE]

The characterization of $\boldsymbol{U}(t)$, which, in a general regime, can be considerably involved because of the time dependence of $\boldsymbol{\mathcal{L}_{BdG}^{ad}}$, is trivial in the considered adiabatic approximation for the eigenmodes of $\boldsymbol{\mathcal{L}_{BdG}^{ad}}$. Indeed, in the representation of (instantaneous) adiabatic eigenstates, $\boldsymbol{U}(t)$, and, in turn, $\boldsymbol{U}^{-1}(t)$, are diagonal. Accordingly, using Eq. (15), we write the matrix elements of the operators $\boldsymbol{U}(t)$ and $\boldsymbol{U}^{-1}(t)$ in terms of the eigenvalues $\lambda_{j}(t)$ of the instantaneous $\boldsymbol{\mathcal{L}_{BdG}^{ad}}$ as

[TABLE]

[TABLE]

Consequently, the projection of Eq. (17) on the eigenmodes of $\boldsymbol{\mathcal{L}_{BdG}^{ad}}$, and, in turn, the expansion of $\delta\Phi_{0}^{nad}$ in them, is straightforward. As we deal here with corrections to the classical field, it is the classical version of the expansion that is applicable. Accordingly, we write

[TABLE]

where instead of the operators $\hat{b}_{a}$ ($\hat{b}_{a}^{\dagger}$) and $\hat{c}_{a}$ ($\hat{c}_{a}^{\dagger}$), present in the expansion of the quantum contribution to the field, we have here the (c-number) time-dependent functions $B_{a}$ ($B_{a}^{*}$) and $C_{a}$ ($C_{a}^{*}$), which give the projections of Eq. (17) on the eigenmodes of $\boldsymbol{\mathcal{L}_{BdG}^{ad}}$, specifically, on the (discrete) modes $\boldsymbol{V_{a}}$ and $\boldsymbol{Z_{a}}$. Again, as we are interested in the effect of the instabilities, we have retained only the contribution of the unstable (discrete) modes in the expansion. According to Eq. (17), those functions are given by

[TABLE]

[TABLE]

where the eigenfunctions $\xi_{a}(x,\boldsymbol{\Lambda}(t))$, $\eta_{a}(x,\boldsymbol{\Lambda}(t))$, $\psi_{a}(x,\boldsymbol{\Lambda}(t))$, and $\zeta_{a}(x,\boldsymbol{\Lambda}(t))$ have the same characteristics as those corresponding to $\boldsymbol{\mathcal{L}_{BdG}}$, introduced in the Appendix. Note that it is the projection of the eigenmodes with the source term $\boldsymbol{S}$ that determines the relevance of the nonadiabatic corrections.

IV The effects of nonadiabaticity on the density and on the nonlocal

density correlation function

To first order in $\hat{\delta\Phi}$, the density, $\hat{\rho}=\hat{\Psi}(x,t)\hat{\Psi}^{\dagger}(x,t)$, can be written as

[TABLE]

where $\rho_{0}^{ad}(x,\boldsymbol{\Lambda}(t))=\left|\Phi_{0}^{ad}(x,\boldsymbol{\Lambda}(t))\right|^{2}$ and

[TABLE]

[TABLE]

One of the most conspicuous experimental features is the growth of the density. The identification of its origin, in particular, the characterization of its differential aspects with the mechanism responsible for the DCF evolution, is one of the open questions. Let us see that some clues to its understanding can be obtained by using our approach. Specifically, we consider the system in the vacuum state $\left|0\right\rangle$ of the annihilation operators $\hat{d}_{a+}$ and $\hat{d}_{a-}$, defined from $\hat{b}_{a}$($\hat{b}_{a}^{\dagger}$) and $\hat{c}_{a}$($\hat{c}_{a}^{\dagger}$) as $\hat{d}_{a+}=\frac{\hat{b}_{a}+i\hat{c}_{a}}{\sqrt{2}}$, and $\hat{d}_{a-}=\frac{\hat{b}_{a}^{\dagger}+i\hat{c}_{a}^{\dagger}}{\sqrt{2}}$. (This is the reference vacuum considered in the studies of stationary models). Although it cannot be assumed that, in the real scenario, the system adiabatically follows the vacuum of excitations from the initial preparation, this simplification can serve to isolate one of the basic mechanisms that can be responsible for the increase in the density. Therefore, we calculate the mean value of $\hat{\rho}$ in the adiabatically evolved state $\left|0\right\rangle$. Using Eqs. (12) and (20), we find

[TABLE]

where we have taken into account that, in parallel with the result for a strictly static setup, the mean value of the quantum term is zero in the vacuum state, i.e.,

[TABLE]

The functions $\sigma_{a}(x,\boldsymbol{\Lambda}(t))$ and $\nu_{a}(x,\boldsymbol{\Lambda}(t))$ present in Eq. (26) are the counterparts for $\boldsymbol{\mathcal{L}_{BdG}^{ad}}$ of the functions respectively given in Eq. (58) for the eigenmodes of $\boldsymbol{\mathcal{L}_{BdG}}$.

It is important to emphasize that, as can be shown from the analysis of Eqs. (21) and (22), the functional forms of $B_{a}(t)$ and $C_{a}(t)$ depart from pure exponentials. Namely, in Eq. (21), the time dependence is contained not only in the factor $e^{-i\vartheta_{a}(t)}$, but also in the integral of $e^{i\vartheta_{a}(t^{\prime})}$, in the adiabatic wavefunctions, and, importantly, in the source term. Yet, assuming that the increasing (decreasing) character of $B_{a}(t)$ $[C_{a}(t)]$ is robust against nonadiabatic corrections, and, therefore, that, at sufficiently large times, the terms that include $B_{a}(t)$ dominate, we can approximate the total density by

[TABLE]

where only the contribution of the most unstable growing term with nonzero projection with the source has been kept. Hence, in agreement with the experimental results, and, in contrast with the picture corresponding to a stationary flow (see the Appendix), changes in the mean value of the density are observed here. (It is worth recalling that the variation in the density is a salient feature of the experimental realization). Two components can be singled out in the spatial pattern of those changes. The first one is the adiabatic wave function of the background field, which enters Eq. (28) through $\rho_{0}^{ad}(x,\boldsymbol{\Lambda}(t))$. The second component is the wave function of the most unstable mode that projects with the nonadiabatic source term, which is incorporated in Eq. (28) via $\sigma_{a}(x,\boldsymbol{\Lambda}(t))$. It is the second component that gives the peculiar character to the density evolution: because of the dynamical instability, the effect of the nonadiabatic variation of the system parameters is amplified; then, it can lead to drastic changes in the density. Extrapolating these conclusions, derived in a post-adiabatic picture, to a general time-varying scenario, we conjecture that the mere dynamical implementation of the black-hole lasing configuration can provide the seed for an instability-determined variation in the system density.

Note that, at the considered order of approximation, the nonadiabatic correction has a purely deterministic effect. Therefore, it is not seen in the DCF if the mean value of the density is extracted in the evaluation of the correlation, as it is done in the standard procedure. Consequently, the DCF is merely determined by the solutions to the homogeneous equation; namely, it reads

[TABLE]

where

[TABLE]

Therefore, in parallel with the results observed in the experiments, we find that, in our toy model, the form of the time dependence of the nonlocal correlation function, given by the exponential $e^{2\Delta_{a}(t)}$, differs from that of the density, incorporated by the function $B_{a}(t)$ in Eq. (28). (We recall that an approximate exponential behavior of experimental DCF is reported in Fig. 5 in key-Steinhauer1 ). In contrast, the presence of $\sigma_{a}(x,\boldsymbol{\Lambda}(t))$ in the obtained expressions for both observables, implies a similar spatial form, also in agreement with the experimental findings. (The analysis of the experimental results has actually uncovered a similar pattern in the spatial dependence of the DCF and of the density key-Clark ). In previous work on static models key-Coutant ; key-Michel1 ; key-Carusotto4 , the wavefunctions of the unstable modes have been evaluated using different approximate methods. Despite the departures of the practical implementation from those simplified setups, a certain similarity still exists between the model-proposed profiles of the velocity of the flow and the speed of sound and those measured in the experiments. Hence, the undulations observed in practice can be linked to the spatial oscillations of the calculated theoretical wavefunctions.

V The effect of the instabilities on the classical noise

Classical fluctuations are present in the practical arrangements for BHL. The study of their implications is required: the identification of differential effects of the instabilities on classical and quantum fluctuations can be crucial to trace the presence of self-amplified HR. A variety of forms of classical noise can be considered. Depending on their characteristics, their effect on the analyzed observables can vary. Let us exemplify this diversity of noisy responses by dealing with two types of fluctuations potentially relevant to the experimental realization key-Steinhauer1 .

V.1 Noise in the system preparation

The limited precision in the preparation of the initial state in the experiments can be accounted for in our model by including classical fluctuations $\delta\Phi^{cn}$ in the initial mean-field wavefunction. For instance, shot-to-shot variations in the number of atoms can be simulated in this form. This type of classical noise is incorporated into our scheme via the two-component spinor

[TABLE]

which is modified in each noise realization. Using the classical version of Eq. (12), $\boldsymbol{\delta\Psi^{cn}}$ is expanded in terms of the eigenmodes of $\boldsymbol{\mathcal{L}_{BdG}^{ad}}$. Then, retaining only the contribution of the unstable modes, we write for the evolved noisy perturbation

[TABLE]

where the coefficients $b_{a}$ and $c_{a}$ respectively denote the projections of the noise spinor on the modes $\boldsymbol{V_{a}}$ and $\boldsymbol{Z_{a}}$, i.e., they are given by

[TABLE]

[TABLE]

Let us see how this kind of classical fluctuations can affect the observables measured in the experiments. Eq. (23) for the density is replaced by $\hat{\rho}=\rho_{0}^{ad}+\hat{\rho}_{1}^{ad}+\rho_{1}^{nad}+\rho_{1}^{cn},$ where to first order in $\delta\Psi^{cn}(x,t)$, the classical-noise contribution to the density is $\rho_{1}^{cn}=(\Psi_{0}^{ad}\delta\Psi^{cn*}+\textrm{c.c.})$. Consequently, the variation in the mean value of $\hat{\rho}$ due to the fluctuations is obtained as

[TABLE]

which, averaged over stochastic realizations $\left\langle\>\right\rangle_{nr}$, gives

[TABLE]

Then, if the noise has zero mean value, i.e., if $\left\langle\delta\Psi^{cn}\right\rangle_{nr}=0$, which, in turn, implies $\left\langle b_{a}\right\rangle_{nr}=\left\langle c_{a}\right\rangle_{nr}=0$, it follows that

[TABLE]

which parallels the result obtained for quantum fluctuations (see Eq. (27)). Additionally, the contribution of noise to the DCF is straightforwardly obtained if the existence of a clearly dominating unstable mode, i.e., of a mode with a growing rate $\Gamma_{a}$ much larger than the others, is assumed. In that case, we find

[TABLE]

where have retained only the secular (nonoscillating) terms. Hence, the pattern is given by the form of the wavefunction of the most unstable mode. Again, this is an analog of the result corresponding to quantum noise, as can be shown by retaining only the dominant term in Eq. (29). Therefore, in the considered temporal range, namely, for times sufficiently large for having the dynamics determined by a clearly dominant unstable mode, there are not qualitative differences between the effects of quantum fluctuations and (classical) noise in the system preparation.

V.2 Technical noise

We consider now fluctuations in the elements that constitute the practical arrangement. For instance, let us deal with stochastic variations in the realization of the optical step-potential used to implement the two-horizon configuration. In this case, the different noise realizations can be regarded as implementations of diverse setups. Therefore, a different set of eigenmodes is applicable to each experimental run. Additionally, there is effective noise in the preparation: the difference between the prepared wave function and the stationary solution of the GP equation for the actually realized set of parameters is taken as a stochastic perturbation. The expansions given by Eqs. (32), (33), and (34) are still applicable. Furthermore, the density is given by Eq. (35) if only the contribution of the unstable modes is retained. The differences with the effect of shot-to-shot noise become evident when the averages over noise realizations are carried out. Namely, the average of density outputs, keeping only the increasing terms, is now given by

[TABLE]

where $b_{a}$ is the projection of the effective noise in the preparation on the mode $\boldsymbol{V_{a}}$. It is important to realize that, even when the fluctuations have zero average over realizations, they can change the density mean value. Actually, as each realization implies not only a different value of $b_{a}$, but also of $\sigma_{a}(x,\boldsymbol{\Lambda}(t))$, the average of density outputs does not have to be zero, i.e., in general,

[TABLE]

Moreover, again, for times sufficiently large for having the evolution determined by only a dominant growing rate $\Gamma_{a}$, we have for the DCF

[TABLE]

where $\Delta_{a}$, and $\sigma_{a}(x,\boldsymbol{\Lambda}(t))$ are characteristics of the most unstable mode for each realization of system parameters. Note that, in this case, as the pattern incorporates an average over wave-functions corresponding to the differently realized set of parameters, the form of the density correlation function differs from that obtained in the previously analyzed cases, i.e., for quantum fluctuations and for noise in the system preparation. In particular, in the present case, since slightly different wave-functions are being averaged, the contrast of the pattern lines can be expected to weaken.

V.3 Interference effects in the system response to classical noise

In the above analysis of the DCF, we have considered the case where there is one unstable mode whose frequency has an imaginary part $\Gamma_{a}$ much larger than those of the rest of modes. Then, at times sufficiently large, the DCF is determined by the characteristics of that dominant mode. No spectral structure is then observed. Let us evaluate now the appearance of qualitative differences in the situation corresponding to having, at least, two unstable modes whose frequencies have imaginary parts, $\Gamma_{a1}$ and $\Gamma_{a2}$, of the same magnitude. The analysis is not specific to the adiabatic scenario. In fact, it also applies to a stationary setup. To emphasize its generality, we suppress the references to adiabaticity in the notation.

For quantum noise, when there are two predominant unstable modes, we obtain from Eq. (59)

[TABLE]

Hence, no interferences effects emerge in the mean value of the DCF if the system can be assumed to be prepared in the vacuum state. (Here, it is worth pointing out that the experimental conditions correspond to a more complex situation. Indeed, since, in the practical realization, the system departs from a stable regime and enters the instability region through the variation in the external potential, the adiabatic following of the vacuum state, assumed in our toy model, does not hold. Therefore, one cannot discard the emergence of spectral structure rooted in the nontrivial evolution of the initial state).

In contrast with the prediction for quantum fluctuations, for (classical) noise in the system preparation, we find

[TABLE]

where an oscillating term with a frequency given by the difference between the real parts of the two considered eigenvalues is apparent. The counterpart expression for technical noise reads

[TABLE]

where, again, oscillations are observed. It is then concluded that, for the two types of classical noise previously considered, interferences appear in the DCF. In the evaluation of the practical importance of these results, one must keep in mind that the approximation of retaining only the most unstable mode in the description becomes worse as the time of observation is reduced: at shorter times, more unstable modes have a non negligible role in the dynamics. Given the time limitations to maintain the two-horizon scheme in practice, the interference effects cannot be neglected. Then, in the analysis of the spectral characteristics of the DCF observed in the experiments, the relevance of more than one unstable mode to the dynamics, and, as a consequence, the resulting interference effects on classical fluctuations, must be taken into account. As previously indicated, one can also contemplate the potential emergence of spectral structure in the quantum scenario when the adiabatic following of the vacuum state does not hold: the breakdown of the adiabatic approximation implies the population of higher modes, and, consequently, the appearance of nontrivial spectral features. An open question is how the quantum fluctuations in the preparation are affected by the system entering the instability region. In this context of discriminating the classical and quantum origin of the findings, it is worth recalling the applicability of the measurement of entanglement as an unambiguous signature for quantum behavior key-Steinhauer2 .

VI Numerical results

In order to confirm the predictions of the above analytical study, we present in this section a numerical simulation of the referred experiments key-Steinhauer1 . We depart from the three-dimensional GP equation where we incorporate the time variation of the trapping potential implemented in the practical setup. The dimensional reduction allowed by the strong confinement in the directions transversal to the step displacement leads to the Non-Polynomial Schrödinger Equation (NPSE) for the longitudinal coordinate. The system is considered to be prepared in the ground state corresponding to the initial confining potential. To emulate the time dependent potential we use the functional form and parameters, in particular the step velocity $v_{s}$, given in Ref. key-Carusotto1 . The NPSE is solved using split-operator and fast Fourier-transform techniques. Results for the sound speed and for the fluid velocity at two different times are presented in Fig. 1.

The similarity with the experimental findings is evident. In particular, a conspicuous feature of the practical realization of Ref. key-Steinhauer1 is reproduced: although an adiabatic approximation is not applicable, a relatively slow variation of the observables is apparent once the white horizon is formed. This aspect of the dynamics has been a key element in the design of our toy model. Indeed, one can think of using a postadiabatic scenario, with an appropriately chosen set of slowly varying parameters, to obtain insight into the observed behavior. Now, focusing on that time regime, we will show that, as analytically predicted, any departure from adiabaticity can become a seed for instability. It is worth recalling that the general objective of our approach has been to establish a link between previous theoretical studies of instability, carried out in static models, and the nonstationary experimental realization. Here, in order to apply the theory to the practical setup, we must turn to a static parallel system where the characterization of the modes can be feasible. Moreover, nonstationarity must be introduced in that analogue system in a form that emulates the detected behavior. Accordingly, we proceed as follows:

i) We employ a monodimensional stationary model used in former work on static BHL settings key-Finazzi1 . In it, the sound speed $c(x)$ and the fluid velocity $v(x)$ are given by the equations

[TABLE]

[TABLE]

[TABLE]

In this framework, the white horizon is located at $-L$ with surface gravity $-c_{H}\kappa_{W}$, and the black horizon is placed at $L$ with surface gravity $c_{H}\kappa_{B}$. [$c_{H}$ denotes the sound speed at the horizons, $D_{c}/c_{H}$ gives the extension of the range (close to the horizons) where both velocities are not flat, and $q$ specifies which part of the sum $c(x)+v(x)$ corresponds to each velocity]. The connection of this model with the practical setup is evident in Fig. 2 (upper panel), where we depict $c(x)$ and $v(x)$ as given by Eqs. (45), (46), and (47). Actually, the characterization of the modes in this system provided interesting clues (key-Finazzi1 ) to modify an early implementation of the BHL setup (key-Steinhauer3 ) in order to facilitate the detection of the amplification. Namely, a change of parameters was proposed to enhance the instability, i.e., to increase the imaginary part of the relevant eigenfrequencies.

ii) By time varying some of the parameters of the model, we incorporate nonstationarity in the setup. Since, in the experimental curves, both, $c$ and $v$, are seen to increase with time, it is appropriate to introduce nonstationarity through the time variation of $c_{H}$. (The implications of changes in other parameters will be discussed further on). Specifically, we consider that $c_{H}$ is modified according to

[TABLE]

Our procedure to describe the dynamics resulting from the variation of $c_{H}$ consists in obtaining first the confining potential $V(x)$ and the effective interaction strength $g(x)$ which lead to a ground-state solution of the GP equation with the characteristics given by Eqs. (45), (46), (47) with $c_{H}=c_{H_{0}}$ key-Finazzi1 . Additionally, the (thus derived) potential and interaction strength are made to incorporate the time-varying $c_{H}(t)$. Then, taking as initial preparation the wave function determined by Eqs. (45), (46), (47) with $c_{H}=c_{H_{0}}$, we solve the GP equation with the time-dependent potential and interaction strength, $V(x,t)$ and $g(x,t)$. We have chosen system parameters appropriate to emulate the practical setup. They correspond to the existence of unstable modes. We intend to confirm that, as shown in Sec. III, [see Eq. (33)], the perturbation of the system, given by the projection of the source term on the modes, is amplified provided that any of the involved modes is unstable, and, that, in turn, it leads to a significant variation in the density. The magnitude of nonadiabaticity is controlled with the parameter $\alpha$. Results for the sound speed and fluid velocity at two different times are presented in Fig. 2 (middle and lower panels). (The thiner lines stand for the strictly adiabatic regime. The thicker lines correspond to $\alpha=0.02$). The differences between the results for the two regimes of time variation show that, as predicted by our analytical study, a departure from adiabaticity can become a seed for activating the instability of the system. It is worth emphasizing the time growing character of those differences.

The properties of the modes in the static system, i.e., in the model with frozen parameters, give further insight into the observed behavior. The use in a former study of a semiclassical approximation to analyze them has provided useful information on the real and imaginary parts of the eigenfrequencies and on the associated eigenfunctions. That study focused on the dependence of the mode properties on $L$. The application of that approach, adapted to the present scenario of fixed $L$ and slightly varying $c_{H}$, allows us to test our predictions on the role of the modes in the emergent dynamics. We have produced two different proofs of consistency. First, from the number of oscillations in the density, we have estimated the wavelength of the leading unstable mode. Then, from the approximately known dependence of the modes on $L$, we have evaluated the separation between horizons that corresponds to having that mode as dominant. (From the application of the quantization rule of Bohr-Sommerfeld, we can assume that no changes take place in the number of unstable modes as $c_{H}$ is varied in the considered time interval). We have found that the required length has the same magnitude as that corresponding to the experiment. Second, from our characterization of the modes, we have been able to alter the system parameters to reduce in a controlled way the number of oscillations in the density. The results are shown in Fig. 3. (It must be noticed that to obtain the curves in Fig. 3 we have kept the former model parameters and have altered the reduced units to work with $\hbar=2$).

It is apparent that despite the restrictions on its applicability, the study offers a framework where a certain degree of control over the system can be achieved. Advances in the design of optimized setups can also be expected from the study of the effect of the variation of other parameters. That analysis requires the characterization of the dependence of the modes on those parameters. Indeed, in order to predict the emergence of instability, the overlapping of the source term with the mode eigenfunctions must be evaluated. The same argument is applicable to the evaluation of the effect of other, deterministic or noisy, perturbation.

In summary, the results of our numerical study confirm the validity of the picture given by our analytical approach. We emphasize that the mere use of terms like instability or unstable modes, which is frequent in the analysis of the experimental findings and which has been applied in the design of the experimental setup, implies assuming that an adiabatic scenario can give a useful ground for describing the essential physics of the system. Our model has given support to that use.

VII Concluding remarks

The non-stationary character of the implementation of BHL has been shown to imply departures from the behavior predicted for a static configuration. Specifically, the study of a post-adiabatic model has uncovered differential aspects in the role of the dynamical instabilities. Apart from leading to the amplification of quantum fluctuations, as in a static setup, the instabilities induce nontrivial variations in the density. The mechanism responsible for the growth of the density mean-value incorporates different elements. First, the variation in the adiabatic wave function of the substrate leads to changes in the density, as expected also in a stable configuration. Second, specific to the presence of instabilities is that the non-adiabatic corrections get amplified provided that they have nonzero projection on the unstable modes. Hence, the density grows because the nonstationary deterministic seed is amplified via the unstable modes. It is this combination of factors that makes the form of the time dependence of the density to differ from that of the nonlocal DCF, which simply corresponds to the exponential growth associated with the instabilities. On the other hand, the similarity between the spatial patterns of both, the density and the DCH, is basically due to the dependence of both observables on the wave-functions of the dominant unstable modes. The generality of the identified mechanism provides a certain predictive power on the implications of the different system components. For instance, although, in our model, the set of variable parameters $\boldsymbol{\Lambda}(t)$ has been assumed to enter the system through the external potential, the consequences of changes in the interaction strength could equally be evaluated.

We have not found qualitative differences between the effects of quantum noise and (classical) fluctuations in the preparation when only one unstable mode dominates the dynamics. This analogy does not persist in the case of technical noise: changes in the density and in the DCF can take place. Furthermore, differential effects of the instabilities on quantum and classical noise are apparent when more than one unstable mode intervene. Namely, for quantum fluctuations from the (adiabatically evolved) vacuum state, no interference effects appear in the DCF. On the contrary, for the two types of classical noise previously considered, oscillations, and, therefore, spectral structure, are apparent in the DCF.

A comment on the limitations of our model is pertinent. Given the presence of instabilities in the system, the use of the adiabatic approximation has required assuming restrictive conditions key-Anglin . In this sense, we have opted for working in a temporal range where the applicability of our approach is guaranteed, and, still, where valuable information on the physical mechanisms specific to non-stationarity could be extracted. Whereas we have dealt with non-adiabatic corrections to the classical field $\Phi_{0}^{ad}(x,\boldsymbol{\Lambda}(t))$, the adiabatic following of, both, the vacuum state and the eigenstates of the operator $\boldsymbol{\mathcal{L}_{BdG}^{ad}}$, has been assumed. As our approach put the focus on a temporal range subsequent to the formation of the two-horizon scheme, it does do not allow us to assess the potential transition of the HR from spontaneous to self-amplifying. A more complete description, which incorporates the breakdown of the adiabatic approximation, and, in particular, the nontrivial evolution of the quantum fluctuations from the initial state, can be expected to account for more complex spectral characteristics. In spite of being a drastic simplification of the experimental setup, our picture is sufficiently accurate to provide arguments useful in the discussion of the findings.

Appendix A DESCRIPTION OF THE STATIC SCENARIO

A.1 Characterization of the eigenmodes of the Bogoliubov-De Gennes operator

$\boldsymbol{\mathcal{L}_{BdG}}$

To describe the elementary excitations of the condensate in a static scenario, i.e., of the considered system for fixed values of the parameters $\boldsymbol{\Lambda}(t)$, we simply rewrite the ansatz in Eq. (11) as

[TABLE]

where $\mu$ denotes the chemical potential. Then, it is found that

[TABLE]

with $H_{GP}$ representing the GP Hamiltonian, i.e.,

[TABLE]

[The (fixed) confining potential is denoted as $V_{ex}(x)$]. Assuming that $\Phi_{0}$ and $\mu$ are known, the evolution of the quantum contribution $\hat{\delta\Psi}$ can be formally analyzed key-Finazzi1 ; key-Leonhardt . In a linear approximation, the dynamics of the excitations are given by the BdG approach. Following the standard procedure, we define the two-component field $\hat{\boldsymbol{W}}\equiv\begin{pmatrix}\hat{\delta\Phi}\\ \hat{\delta\Phi}^{\dagger}\end{pmatrix}$, and work initially with its classical version $\boldsymbol{W}\equiv\begin{pmatrix}\delta\Phi\\ \delta\Phi^{*}\end{pmatrix}$. Using Eqs. (49), (50), and (51), we find

[TABLE]

where the BdG operator $\boldsymbol{\mathcal{L}_{BdG}}$ is given by

[TABLE]

The spectrum of $\boldsymbol{\mathcal{L}_{BdG}}$ has been analyzed in previous research. From the non-Hermitian character of $\boldsymbol{\mathcal{L}_{BdG}}$, it follows that its eigenvalues can be real or complex. The complex eigenfrequencies constitute a discrete spectrum, and, given the characteristics of $\boldsymbol{\mathcal{L}_{BdG}}$, appear in pairs, which will be denoted as $\lambda_{a}$ and $\lambda_{a}^{*}$ ($\lambda_{a}=\omega_{a}+i\Gamma_{a}$). Without loss of generality, we assume $\Gamma_{a}>0$. The respective associated eigenmodes can be written as

[TABLE]

The imaginary parts of the eigenvalues lead to instabilities which are responsible for the lasing effect. The growth rate of the perturbations is given by $\Gamma_{a}$. Additionally, the real eigenfrequencies are shown to form a continuous spectrum. They will be denoted as $\omega$, and the associated eigenmodes can be expressed as

[TABLE]

The variable $\alpha$ accounts for possible degeneracies. The field $\hat{\boldsymbol{W}}$ can be expanded as $\hat{\boldsymbol{W}}=\sum_{n}\boldsymbol{W_{n}}\hat{a}_{n}+\bar{\boldsymbol{W}}_{\boldsymbol{n}}\hat{a}_{n}^{\dagger}$, where $\boldsymbol{W_{n}}$ denotes the different eigenmodes (both of the continuum and discrete spectrum) and $\bar{\boldsymbol{W}}_{\boldsymbol{n}}$ stands for the conjugate spinor of $\boldsymbol{W_{n}}$, which, using the Pauli matrix $\sigma_{x}$, is defined as $\bar{\boldsymbol{W}}_{\boldsymbol{n}}=\sigma_{x}\boldsymbol{W_{n}^{*}}$. Moreover, $\hat{a}_{n}$ is a compact notation for the two kinds of operators potentially present in the description: it incorporates the operators $\hat{a}_{\omega}^{\alpha}$ corresponding to the continuous set of modes, and it also stands for the operators $\hat{b}_{a}$ and $\hat{c}_{a}$ linked to the discrete set. The commutators are $\left[\hat{a}_{\omega}^{\alpha},\hat{a}_{\omega^{\prime}}^{\alpha^{\prime}\dagger}\right]=\delta_{\alpha,\alpha^{\prime}}\delta(\omega-\omega^{\prime})$, and $\left[\hat{b}_{a},\hat{c}_{a^{\prime}}^{\dagger}\right]=i\delta_{a,a^{\prime}}$, vanishing all the others. As opposed to $\hat{a}_{\omega}^{\alpha}$, $\hat{b}_{a}$ and $\hat{c}_{a}$ are not annihilation operators key-Finazzi1 .

It is shown that, in a general case [i.e., when $\boldsymbol{\mathcal{L}_{BdG}}$ has a continuous set of modes $\boldsymbol{W_{\omega}^{\alpha}}$ and a discrete set of (unstable) modes $\boldsymbol{V_{a}}$ and $\boldsymbol{Z_{a}}$], the expansion of the field operator $\hat{\delta\Psi}(x,t)=e^{-i\mu t/\hbar}\hat{\delta\Phi}(x,t)$ is given by

[TABLE]

A.2 The density and the density-density correlation function

To first order in $\hat{\delta\Psi}(x,t)$, we have for the density $\hat{\rho}=\rho_{0}+(\Phi_{0}\hat{\delta\Phi}^{\dagger}+\textrm{h.c.})\equiv\rho_{0}+\hat{\rho}_{1}$, ($\rho_{0}(x)=\left|\Phi_{0}(x)\right|^{2}$). As we intend to account for Hawking radiation, we consider the system prepared in the vacuum state $\left|0\right\rangle$. It is worth clarifying that $\left|0\right\rangle$ refers to the vacuum of the real annihilation operators $\hat{d}_{a+}$ and $\hat{d}_{a-}$ defined in the standard form $\hat{d}_{a+}=\frac{\hat{b}_{a}+i\hat{c}_{a}}{\sqrt{2}}$, and $\hat{d}_{a-}=\frac{\hat{b}_{a}^{\dagger}+i\hat{c}_{a}^{\dagger}}{\sqrt{2}}.$ Introducing those changes into the expansion in Eq. (56), we find $\left\langle 0\right|\hat{\rho}_{1}\left|0\right\rangle=0.$ Additionally, for the nonlocal DCF, we have

[TABLE]

where we have retained only the contribution of the unstable modes. The functions

[TABLE]

respectively correspond to the increasing and decreasing modes. Keeping only the increasing terms (recall that $\Gamma_{a}>0$), we obtain key-Finazzi1

[TABLE]

Furthermore, if there is one unstable mode whose frequency has an imaginary part $\Gamma_{a}$ much larger than that of any other mode, we can, for sufficiently large times, keep only the contribution of that mode in Eq. (59). An exponential growth of the nonlocal density correlation is then apparent. Moreover, the spatial pattern is determined by $\sigma_{a}(x)$, i.e., by a combination of the wave-functions of the dominant mode.

A conclusion particularly relevant to the comparison with the non-stationary configuration must be singled out. Namely, in the stationary regime, the existence of instabilities does not lead the quantum fluctuations to contribute to the mean density. In contrast, a self-amplifying effect of quantum noise is apparent in the nonlocal density correlation function.

Acknowledgments

One of us (JMGL) acknowledges the support of the Spanish Ministerio de Economía y Competitividad and the European Regional Development Fund (Grant No. FIS 2016-79596-P).