Entanglement and violation of classical inequalities in the Hawking radiation of flowing atom condensates

The detection of spontaneous Hawking radiation (HR) remains a major challenge in modern physics. Originally predicted for cosmological black-holes (BH) [1], it was soon noticed [2, 3] that the emission of HR is a kinematic effect that can be observed in an ordinary laboratory. For a quantum fluid traversing a subsonic–supersonic interface (which amounts to a sonic event horizon), the spontaneous correlated emission of phonons into the subsonic and the supersonic regions has been predicted [4–11]. Sonic event horizons have been realized by accelerating a Bose–Einstein condensate [12]. In a similar setup, the self-amplifying stimulated HR resulting from the BH laser effect [10] has been observed [13]. However, the emission of spontaneous HR still remains unobserved.

In the context of quasi-stationary flow scenarios, it has been proposed that a large leaking condensate can provide a suitable subsonic–supersonic interface, where HR production could be conveniently observed [14, 15]. Realistic protocols to produce such quasi-stationary flow regimes have been investigated recently [16].

It has also been noticed [17] that the violation of classical Cauchy–Schwarz (CS) inequalities [18, 19] by outgoing quasiparticles can provide unambiguous evidence of deep quantum behavior and ultimately of the existence of a spontaneous contribution to HR. The distinction between spontaneous and thermal (stimulated) contributions is a fundamental requirement in the search for spontaneous HR. Proposals based on density–density correlations in real space do not meet that requirement [7, 20]. The mere measurement of phonon [21] or atom [14] intensity spectra would not permit that distinction either.

An alternative scheme to identify deep quantum behavior relies on the detection of entanglement between the various outgoing radiation channels from a sonic horizon [22, 23]. A recent work has addressed the possibility of detecting entanglement through density–density correlations in Fourier space [24]. Approaches based on entanglement detection have been also proposed in analogous contexts such as inflationary cosmology [25], astrophysical BH [26], general relativistic quantum fields [27], or other BH analogs [28, 29].

Our present work aims at clearly establishing the theoretical relation between CS violation and entanglement in the outgoing quasiparticle modes of analog HR, unifying the existent work of [17, 22, 23]. We also study their potential measurability in specific detection schemes [24]. Importantly, by introducing envelope-modulated Fourier transforms we explicitly take into account the role played by the spatial location of the asymptotic regions. This allows us to show that the violation of only quadratic CS inequalities can be measured. In particular, we prove that the complete implementation of the generalized Peres–Horodecki (GPH) criterion [22, 23] or the quartic CS violation [17] requires the knowledge of parameters that are impossible or very difficult to measure, at least within the class of detection schemes here considered. However, we also show that this does not represent an important limitation in practice, since, under rather broad conditions, all the previous criteria become equivalent.

Although the present work is motivated by the quest for the observation of spontaneous HR, the results here obtained are also of relevance in neighboring fields such as quantum optics [30, 31] and quantum information physics [32, 33], as well as in the broader topic of bosonic condensates [34, 35]. Atom flow through a sonic black hole may also be viewed as one of the paradigms of atom quantum transport through a barrier [36]. In general, a subsonic/supersonic interface can behave as a basic element that provides novel functionalities in atom circuits within the emerging field of atomtronics [37, 38]. For instance, the BH configuration of [16] can be used to produce a quasi-stationary supersonic atom current with well-defined velocity. Another recent numerical work [39] has shown that a BH laser configuration is able to reach a regime of continuous emission of solitons, providing a hydrodynamic analog of an optical laser.

This paper is arranged as follows. In section 2, we present the physical model that we use in this work. Section 3 contains a detailed comparison among the existent works on the detection of the spontaneous HR and unifies the criteria under certain conditions. In section 4, we discuss the possibility of experimentally verifying the discussed criteria. Finally, in section 5, we present numerical results on the detection of the CS violation in typical HR setups.

We start by considering the stationary flow of a one-dimensional BEC in a regime known as 1D-mean field regime [36, 40]. The condensate is described by the stationary GP wave function ${{\rm{\Psi }}}_{0}(x)=\sqrt{{n}_{0}(x)}{{\rm{e}}}^{{\rm{i}}\theta (x)}.$ We can define the local sound speed and flow velocity as $c(x)={[{gn}(x)/m]}^{1/2}$ and $v(x)={\hslash }{\theta }^{\prime }(x)/m,$ respectively, with g the coupling constant and m the mass of the atoms. In our convention, $v(x)\gt 0.$ A sonic BH configuration is that involving two asymptotically homogeneous regions such that one is subsonic (the upstream region, $x\to -\infty$), with ${v}_{{\rm{u}}}\lt {c}_{{\rm{u}}},$ while the other one is supersonic (downstream region, $x\to \infty$), with ${v}_{{\rm{d}}}\gt {c}_{{\rm{d}}}.$ The magnitudes ${c}_{{\rm{u}},{\rm{d}}},{v}_{{\rm{u}},{\rm{d}}}$ represent the asymptotic homogeneous values of the sound and flow velocities. Hereafter, we use units such that ${\hslash }=m={c}_{{\rm{u}}}={k}_{{\rm{B}}}=1.$ Throughout this work, we follow the notation of [14, 15, 20].

Within the Bogoliubov–de Gennes (BdG) approximation, the field operator is decomposed as:

$$\begin{eqnarray}&&\hat{{\rm{\Psi }}}(x)={{\rm{\Psi }}}_{0}(x)+\hat{\varphi }(x){{\rm{e}}}^{{\rm{i}}\theta (x)}=[\sqrt{{n}_{0}(x)}+\hat{\varphi }(x)]{{\rm{e}}}^{{\rm{i}}\theta (x)},\end{eqnarray} \tag{ 1 }$$

where we have removed the mean-field condensate phase, ${{\rm{e}}}^{{\rm{i}}\theta (x)},$ from the definition of the field fluctuations $\hat{\varphi }(x)$ for computational purposes. The expression of $\hat{\varphi }(x)$ in a BH configuration is given by:

$$\begin{eqnarray}\hat{\varphi }(x) & = & {\displaystyle \int }_{0}^{\infty }{\rm{d}}\omega \displaystyle \sum _{I}\left[{u}_{I-{\rm{in}},\omega }(x){\hat{\gamma }}_{I-{\rm{in}}}(\omega )+{v}_{I-{\rm{in}},\omega }^{*}(x){\hat{\gamma }}_{I-{\rm{in}}}^{\dagger }(\omega )\right]\\ & & +{\displaystyle \int }_{0}^{{\omega }_{{\rm{max}}}}{\rm{d}}\omega \left[{u}_{d2-{\rm{in}},\omega }(x){\hat{\gamma }}_{d2-{\rm{in}}}^{\dagger }(\omega )+{v}_{d2-{\rm{in}},\omega }^{*}(x){\hat{\gamma }}_{d2-{\rm{in}}}(\omega )\right].\end{eqnarray} \tag{ 2 }$$

The components of the spinor ${z}_{a,\omega }(x)\equiv {[{u}_{a,\omega }(x),{v}_{a,\omega }(x)]}^{\top }$ are solutions to the BdG equations at a given frequency ω. Here we have adopted the convention of working with positive frequencies at the expense of introducing the negative-normalization channel d2. The index I is summed over the conventional, positive-normalization channels (u and d1); see figure 1 for the mode notation. In the asymptotically homogeneous regions, the macroscopic wave function has the form

$$\begin{eqnarray}&&{{\rm{\Psi }}}_{0}(x)=\sqrt{{n}_{r}(x)}{{\rm{e}}}^{{\rm{i}}{v}_{r}x+{\theta }_{r}},\quad \quad r=u,d\end{eqnarray} \tag{ 3 }$$

where ${\theta }_{r}$ is a constant phase. Thus, the scattering states ${z}_{i,\omega }(x)$ are combinations of plane waves (scattering channels) in the asymptotic regions u and d, while ${\hat{\gamma }}_{i}(\omega )$ is their corresponding a nnihilation operator, with the index i taking values $u-,d1-,d2-{\rm{in}}.$ For instance, the scattering state $d2-{\rm{in}}$ has the asymptotic form

$$\begin{eqnarray}{z}_{d2-{\rm{in}},\omega }\left(x\to -\infty \right) & = & {S}_{{ud}2}\left(\omega \right){s}_{u-{\rm{out}},\omega }\left(x\right),\\ {z}_{d2-{\rm{in}},\omega }\left(x\to \infty \right) & = & {s}_{d2-{\rm{in}},\omega }\left(x\right)+{S}_{d1d2}\left(\omega \right){s}_{d1-{\rm{out}},\omega }\left(x\right)\\ & & +{S}_{d2d2}\left(\omega \right){s}_{d2-{\rm{out}},\omega }\left(x\right).\end{eqnarray} \tag{ 4 }$$

In this expression, ${s}_{a,\omega }\left(x\right)$ is the spinor (plane) wave function of the corresponding scattering channel a:

$$\begin{eqnarray}{s}_{a,\omega }\left(x\right) & = & \displaystyle \frac{{{\rm{e}}}^{{\rm{i}}{k}_{a}\left(\omega \right)x}}{\sqrt{2\pi | {w}_{a}\left(\omega \right)| }}\left[\begin{array}{c}{u}_{a}(\omega )\\ {v}_{a}(\omega )\end{array}\right],\\ &&\left[\begin{array}{c}{u}_{a}(\omega )\\ {v}_{a}(\omega )\end{array}\right] & = & {N}_{a}\left[\begin{array}{c}\displaystyle \frac{{k}_{a}^{2}\left(\omega \right)}{2}+[\omega -{v}_{r}{k}_{a}\left(\omega \right)]\\ \displaystyle \frac{{k}_{a}^{2}\left(\omega \right)}{2}-[\omega -{v}_{r}{k}_{a}\left(\omega \right)]\end{array}\right],\\ {N}_{a} & = & {[2{k}_{a}^{2}\left(\omega \right)\left|\omega -{v}_{r}{k}_{a}\left(\omega \right)\right|]}^{-1/2},\end{eqnarray} \tag{ 5 }$$

and the corresponding wave vector, ${k}_{a}(\omega )$ is given implicitly by the dispersion relation

$$\begin{eqnarray}&&{{\rm{\Omega }}}_{a}^{2}={\left(\omega -{v}_{r}{k}_{a}\right)}^{2}={c}_{r}^{2}{k}_{a}^{2}+{k}_{a}^{4}/4,\end{eqnarray} \tag{ 6 }$$

where ${{\rm{\Omega }}}_{a}\gt 0$ is the comoving frequency, defined positive, and the index $r=u,d$ characterizes the asymptotic region corresponding to the scattering channel a. The group velocity of the scattering channel a is ${w}_{a}\left(\omega \right)\equiv {[{\rm{d}}{k}_{a}\left(\omega \right)/{\rm{d}}\omega ]}^{-1}$ and it is included in the definition of the scattering channels in order to properly normalize them to a given quasiparticle flux. The dispersion relation is depicted in figure 1 for the subsonic and supersonic asymptotic regions.

Zoom In Zoom Out Reset image size

The expressions for the other scattering states are similar to equation (4). They all are characterized by an incoming channel carrying unit flux and several scattering amplitudes (elements of the S-matrix) describing the transition to the outgoing modes, which are also assumed to carry unit flux. The 'out' (advanced) scattering states are related to the 'in' (retarded) scattering states through the scattering matrix S for a given value of the frequency ω:

$$\begin{eqnarray}&&\left[\begin{array}{c}{\hat{\gamma }}_{u-{\rm{out}}}\\ {\hat{\gamma }}_{d1-{\rm{out}}}\\ {\hat{\gamma }}_{d2-{\rm{out}}}^{\dagger }\end{array}\right]=\left[\begin{array}{ccc}{S}_{{uu}} & {S}_{{ud}1} & {S}_{{ud}2}\\ {S}_{d1u} & {S}_{d1d1} & {S}_{d1d2}\\ {S}_{d2u} & {S}_{d2d1} & {S}_{d2d2}\end{array}\right]\left[\begin{array}{c}{\hat{\gamma }}_{u-{\rm{in}}}\\ {\hat{\gamma }}_{d1-{\rm{in}}}\\ {\hat{\gamma }}_{d2-{\rm{in}}}^{\dagger }\end{array}\right].\end{eqnarray} \tag{ 7 }$$

Here, as often in the rest of the paper, the frequency dependence is understood.

We note that equation (7) is a Bogoliubov-type relation linearly combining destruction and creation operators. This is due to the negative normalization of the d2 modes, see figure 1. Thus, the vacuum of the 'in' modes can be regarded as a squeezed state in the representation of the 'out' modes. This non-uniqueness of the vacuum is at the root of the Hawking effect itself. In order to detect the intrinsic quantum behavior corresponding to the squeezed character of the zero-point HR, we focus on the correlations between the 'out' modes.

Several criteria have been proposed in order to distinguish the spontaneous from the stimulated (thermal) HR or the coherent collective modes. First, it was argued in [17] that the violation of CS inequalities can be regarded as an unequivocal signature of the presence of zero-point dynamics. Specifically, the following second-order correlation function between two given modes $i,j$ was considered:

$$\begin{eqnarray}&&{{\rm{\Gamma }}}_{{ij}}\equiv \langle {\hat{\gamma }}_{i}^{\dagger }{\hat{\gamma }}_{j}^{\dagger }{\hat{\gamma }}_{j}{\hat{\gamma }}_{i}\rangle \gt 0,\end{eqnarray} \tag{ 8 }$$

where ${\hat{\gamma }}_{i}$ is the annihilation operator of the mode i. The expectation value of an operator $\hat{O}$ is taken over the state of the system, described by the density matrix $\hat{\rho },$ $\langle \hat{O}\rangle =\mathrm{Tr}(\hat{\rho }\hat{O}).$ The relation

$$\begin{eqnarray}&&{{\rm{\Gamma }}}_{{ij}}\leqslant \sqrt{{{\rm{\Gamma }}}_{{ii}}{{\rm{\Gamma }}}_{{jj}}},\end{eqnarray} \tag{ 9 }$$

is a CS-type inequality which is always satisfied in a classical context. The violation of (9) is characterized by the positiveness of the difference

$$\begin{eqnarray}&&{{\rm{\Theta }}}_{{ij}}\equiv {{\rm{\Gamma }}}_{{ij}}-\sqrt{{{\rm{\Gamma }}}_{{ii}}{{\rm{\Gamma }}}_{{jj}}}\gt 0.\end{eqnarray} \tag{ 10 }$$

We will refer to equation (10) as the quartic CS violation. We can also define the corresponding first-order correlation functions,

$$\begin{eqnarray}&&{g}_{{ij}}\equiv \langle {\hat{\gamma }}_{i}^{\dagger }{\hat{\gamma }}_{j}\rangle ,\;{c}_{{ij}}\equiv \langle {\hat{\gamma }}_{i}{\hat{\gamma }}_{j}\rangle ,\end{eqnarray} \tag{ 11 }$$

whose associated quadratic CS violations, in analogy to equation (10), are given by:

$$\begin{eqnarray}&&{{\rm{\Delta }}}_{{ij}}\equiv | {c}_{{ij}}{| }^{2}-{g}_{{ii}}{g}_{{jj}}\gt 0.\end{eqnarray} \tag{ 12 }$$

The inequality $| {g}_{{ij}}{| }^{2}\leqslant {g}_{{ii}}{g}_{{jj}}$ is always satisfied; see appendix A.

Another possible signature of the quantum character of the system is the presence of entanglement between two modes $i,j.$ In this work, we follow the convention of defining entanglement as the non-separability of the state of the system. We say that a two-mode state $\hat{\rho }$ is separable when it can be decomposed as:

$$\begin{eqnarray}&&\hat{\rho }=\displaystyle \sum _{n}{p}_{n}{\hat{\rho }}_{n}^{(i)}\otimes {\hat{\rho }}_{n}^{(j)},\end{eqnarray} \tag{ 13 }$$

where ${\hat{\rho }}_{n}^{(i)}$ is a state of the Hilbert subspace spanned by the i mode. The use of the GPH criterion [41, 42] was proposed in [22, 23] to identify the entanglement between two modes within the context of analog HR emission. The GPH criterion asserts that the state is entangled if ${{\mathcal{P}}}_{{ij}}\lt 0,$ where ${{\mathcal{P}}}_{{ij}}$ is the GPH function defined in equation (A.11). In the particular case of Gaussian states, the reverse implication is also true, i.e., the GPH criterion is also a necessary condition for entanglement [42].

Some general remarks about CS violations and entanglement are in order. The violation of a quadratic CS inequality, equation (12), is a sufficient condition for the non-separability of the system [43]. In particular, we show in appendix A that the quadratic CS violation is a sufficient condition for the fulfillment of the GPH criterion, which in turn is known to imply the presence of entanglement. It can also be shown that the quartic CS violation is independent of the entanglement of the system (see appendix A), so separable states can violate quartic CS inequalities and viceversa.

Now we focus on the specific case of analog HR in a BEC. We evaluate the correlation functions of equations (8) and (11) for the 'out' modes at given ω:

$$\begin{eqnarray}{{\rm{\Gamma }}}_{{ij}}(\omega ) & = & \langle {\hat{\gamma }}_{i-{\rm{out}}}^{\dagger }(\omega ){\hat{\gamma }}_{j-{\rm{out}}}^{\dagger }(\omega ){\hat{\gamma }}_{j-{\rm{out}}}(\omega ){\hat{\gamma }}_{i-{\rm{out}}}(\omega )\rangle ,\\ {g}_{{ij}}(\omega ) & \equiv & \langle {\hat{\gamma }}_{i-{\rm{out}}}^{\dagger }(\omega ){\hat{\gamma }}_{j-{\rm{out}}}(\omega )\rangle ,\\ {c}_{{ij}}(\omega ) & \equiv & \langle {\hat{\gamma }}_{i-{\rm{out}}}(\omega ){\hat{\gamma }}_{j-{\rm{out}}}(\omega )\rangle .\end{eqnarray} \tag{ 14 }$$

Hereafter, $i,j=u,d1,d2.$ The associated quartic and quadratic CS violations are characterized by the positivity of ${{\rm{\Theta }}}_{{ij}}(\omega )$ and ${{\rm{\Delta }}}_{{ij}}(\omega ),$ respectively, as defined in equations (10), (12). In a similar fashion, we study the GPH function ${{\mathcal{P}}}_{{ij}}(\omega )$ for the $i,j$ 'out' modes at the same frequency ω. For simplicity, in this section we will obviate the Dirac delta factors appearing in all equal-frequency correlation functions (see equation (15) for a complete expression).

One may wonder to what extent the previous criteria differ when applied to spontaneous HR. We devote the rest of this section to prove that, when the state of the system $\hat{\rho }$ is Gaussian and incoherent in the 'in' modes, the quartic and quadratic CS violations become equivalent to the GPH criterion. Specifically, the requirement of incoherent incoming modes can be expressed as:

$$\begin{eqnarray}\langle {\hat{\gamma }}_{I-{\rm{in}}}(\omega ){\hat{\gamma }}_{j-{\rm{in}}}({\omega }^{\prime })\rangle & = & 0,\\ \langle {\hat{\gamma }}_{I-{\rm{in}}}^{\dagger }(\omega ){\hat{\gamma }}_{j-{\rm{in}}}({\omega }^{\prime })\rangle & = & {n}_{i}(\omega ){\delta }_{{ij}}\delta (\omega -{\omega }^{\prime })\ne 0.\end{eqnarray} \tag{ 15 }$$

These correlation functions fully characterize the state provided it is Gaussian. It is shown in appendix B that, for that class of states, we only need seven parameters for the computation of the CS violations and the GPH function. The case of Gaussian states is important because, within the BdG approximation, the dynamics is Gaussian.

For convenience, we define the complex vector

$$\begin{eqnarray}&&{\alpha }_{i}(\omega )\equiv \left[\begin{array}{c}{S}_{{iu}}(\omega )\sqrt{{n}_{u}(\omega )}\\ {S}_{{id}1}(\omega )\sqrt{{n}_{d1}(\omega )}\\ {S}_{{id}2}(\omega )\sqrt{{n}_{d2}(\omega )+1}\end{array}\right].\end{eqnarray} \tag{ 16 }$$

From the previous definitions, it is easy to show that the only non-zero quadratic correlations for 'out' modes are:

$$\begin{eqnarray}{g}_{{{II}}^{\prime }} & = & {\alpha }_{I}^{\dagger }\cdot {\alpha }_{{I}^{\prime }},\\ {g}_{d2d2} & = & | {\alpha }_{d2}{| }^{2}-1,\\ {c}_{{Id}2} & = & {\alpha }_{d2}^{\dagger }\cdot {\alpha }_{I},\end{eqnarray} \tag{ 17 }$$

where the index I stands for a normal, positive-normalization (u or d1) mode. First, we study the correlation between a normal and an anomalous (negative-normalization) mode. The case I = u corresponds to the proper Hawking effect [14, 20] and the case $I=d1$ corresponds to the bosonic equivalent of the Andreev reflection [44].

As we are working with Gaussian states, we can apply Wick's theorem to compute the quartic correlations as a function of the quadratic correlations

$$\begin{eqnarray}{{\rm{\Gamma }}}_{{II}} & = & 2\;{g}_{{II}}^{2}=2| {\alpha }_{I}{| }^{4},\\ {{\rm{\Gamma }}}_{d2d2} & = & 2\;{g}_{d2d2}^{2}=2{({| {\alpha }_{d2}| }^{2}-1)}^{2},\\ {{\rm{\Gamma }}}_{{Id}2} & = & | {c}_{{Id}2}{| }^{2}+{g}_{{II}}{g}_{d2d2},\\ & = & | {\alpha }_{d2}^{\dagger }\cdot {\alpha }_{I}{| }^{2}+{| {\alpha }_{I}| }^{2}(| {\alpha }_{d2}{| }^{2}-1).\end{eqnarray} \tag{ 18 }$$

Therefore, the condition for the quartic CS violation ${{\rm{\Theta }}}_{{Id}2}(\omega )\gt 0$ reduces to the simpler quadratic CS violation

$$\begin{eqnarray}&&{{\rm{\Theta }}}_{{Id}2}(\omega )={{\rm{\Delta }}}_{{Id}2}(\omega )\gt 0\end{eqnarray} \tag{ 19 }$$

which, using equations (17) and (18), can be rewritten as:

$$\begin{eqnarray}&&| {\alpha }_{d2}^{\dagger }\cdot {\alpha }_{I}{| }^{2}\gt {| {\alpha }_{I}| }^{2}(| {\alpha }_{d2}{| }^{2}-1).\end{eqnarray} \tag{ 20 }$$

The $-1$ within the second bracket results from the anomalous character of the d2 mode and is responsible for making the violation of the CS inequality possible. Thus, we have proven that the quartic and the quadratic violations are equivalent for Id2, normal-anomalous correlations.

Now we turn our attention to the GPH criterion. From equation (A.18) we compute the GPH function for the pair of modes Id2:

$$\begin{eqnarray}&&{{\mathcal{P}}}_{{Id}2}=({g}_{{II}}{g}_{d2d2}-{| {c}_{{Id}2}| }^{2})[({g}_{{II}}+1)({g}_{d2d2}+1)-{| {c}_{{Id}2}| }^{2}].\end{eqnarray} \tag{ 21 }$$

As the second (square) bracket in the r.h.s. of equation (21) is always positive (see equation (A.5)), we conclude that ${{\mathcal{P}}}_{{Id}2}\lt 0$ if and only if equation (19) is satisfied. As the state is Gaussian, the GPH criterion is equivalent to the entanglement of the state of the system $\hat{\rho }.$ We conclude that the quadratic and quartic CS violations are equivalent to entanglement. We note that the equivalence between condition (19) and entanglement was already pointed out in [22] but the connection with the CS violation criterion of [17] was not made explicit. A similar result appeared in [45] on the equivalence between non-separability and CS violation when studying the correlation function at two different times in fluids of light.

In regard to the correlation between the two normal modes, we obtain, following similar arguments, that the quartic CS inequality (9) reduces to:

$$\begin{eqnarray}&&| {g}_{{ud}1}{| }^{2}\leqslant {g}_{{uu}}{g}_{d1d1},\end{eqnarray} \tag{ 22 }$$

which can never be violated, as previously explained (see equation (12) and accompanying discussion). In particular, the inequality (22) can be rewritten as:

$$\begin{eqnarray}&&| {\alpha }_{u}^{\dagger }\cdot {\alpha }_{d1}{| }^{2}\leqslant | {\alpha }_{u}{| }^{2}| {\alpha }_{d1}{| }^{2},\end{eqnarray} \tag{ 23 }$$

which is always satisfied for two complex vectors. Thus, there is no CS violation in the correlation between normal modes.

When considering the entanglement between the two normal modes, we have $| {g}_{{ud}1}| \gt | {c}_{{ud}1}| =0$ and thus there is no entanglement (see equation (A.18) and the ensuing discussion).

From the previous discussion we conclude that, when considering Gaussian and incoming incoherent states (i.e., states satisfying equation (15)), all the here considered criteria for characterizing the quantum character of the system become equivalent. This result is important because it unifies the work of [17, 22, 23].

A most important particular case is that involving a thermal distribution of incoming quasiparticles that have thermalized in the comoving reference frame [22, 46] so that their occupation factor is:

$$\begin{eqnarray}&&{n}_{i}(\omega )=\displaystyle \frac{1}{\mathrm{exp}({{\rm{\Omega }}}_{i}(\omega )/T)-1}\end{eqnarray} \tag{ 24 }$$

with ${{\rm{\Omega }}}_{i}(\omega )$ the commoving frequency of the mode i-in at laboratory frequency ω, as given by equation (6). This state satisfies the mentioned conditions, so there is no difference between using the GPH criterion or the CS violation to characterize the quantum behavior of the system.

Finally, we discuss the differences that appear when removing some restrictions. If the state is Gaussian but not incoherent in the incoming channels, the GPH criterion is still equivalent to the non-separability of the system. On the other hand, in the same case (Gaussian and not incoherent), the quadratic CS violation is no longer equivalent to the GPH criterion; rather, it is only a sufficient condition for it. As a consequence, the GPH criterion is a more powerful criterion than the quadratic CS violation to detect non-separability (it can identify non-separability in cases where the quadratic CS violation would fail).

However, when relaxing the requirement of incoherence, the quartic CS violation becomes independent of the GPH criterion (see appendix A). In this context, we wish to note, using quantum optics terminology, that the presence of quartic CS violation or entanglement requires that the system is described by a non-classical Glauber–Sudarshan function, i.e., a Glauber–Sudarshan function that takes negative values.

On the other hand, if the state of the system is not Gaussian but is incoherent in the incoming channels, the quartic CS violation is also independent of the entanglement. The quadratic CS violation equation (19) and the GPH criterion are still equivalent between them, but no longer equivalent to the entanglement of the system.

For a general state which does not satisfy any of the previous conditions, 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. We summarize all these logical relations in table 1.

Table 1.  Logical relations between the different criteria studied throughout this work in the various physical cases. The three rightmost entries in the upper row indicate the different criteria for the identification of quantum behavior: $\mathrm{CS}2,$ $\mathrm{CS}4$ stand for quadratic and quartic CS violations, respectively, while GPH stands for generalized Peres–Horodecki criterion; NS means non-separability. The two leftmost columns define the various physical cases here considered. By 'incoherent' state, we understand a density matrix which is diagonal in the representation of retarded quasiparticle scattering modes, each characterized by a single incoming channel; see equation (15). The symbol • stands for the uppermost entry in the corresponding column. The abbreviation 'indep.' means that, in the three lower cases, the quartic CS violation is independent of the non-separability of the system.

Gaussian	Incoherent	$\mathrm{CS}4$	$\mathrm{CS}2$	GPH
yes	yes	•⇔ GPH	$\bullet \;\iff$ GPH	$\bullet \iff$ NS
yes	no	•indep. NS	$\bullet \;\Rightarrow$ GPH	$\bullet \;\iff$ NS
no	yes	•indep. NS	$\bullet \;\iff$ GPH	$\bullet \;\Rightarrow$ NS
no	no	•indep. NS	$\bullet \;\Rightarrow$ GPH	$\bullet \;\Rightarrow$ NS

Remarkably, the quadratic CS violation equation (19) reveals at the same time two different aspects of quantum behavior: the violation of a classical inequality and the entanglement of the system.

In this section, we analyze possible detection schemes of the criteria discussed in the previous section. A particular type of CS violation between two colliding condensates was measured using time-of-flight detection [34]. The possible detection of the quartic CS violation in a TOF experiment for analog HR was discussed in [17]. A recent work has studied the measurement of quartic CS violations [47] using phonon evaporation, following a related work on the dynamical Casimir effect [48]. On the other hand, the detection of entanglement using the GPH criterion was analyzed in [23] using density fluctuations or the optomechanical detection of phonons. Recently, it has been proposed that the entanglement of HR in a BEC can be detected experimentally by measuring the density–density correlations through in situ imaging [24]. All these setups have in common that they involve collective atom flow and thus can be viewed as elementary components of a larger atomtronic circuit.

For illustrative purposes, we focus on the comparison between the detection of the CS violation and the GPH criterion for the particular scenario of density–density correlations, but the analysis which we present in this work can be generalized to other schemes. The role of the spatial density correlation function in analog models has been extensively studied in several works [7, 15, 20]. In the BEC context, the measurement of the spatial density correlations has been experimentally used to characterize the BH Laser [13]. The density–density correlation function also plays an important role in the polariton analog [49, 50].

We start by considering the density in the same BdG approximation of section 2 and expand up to first order in the fluctuations of the field operator:

$$\begin{eqnarray}&&\hat{n}(x)\simeq {n}_{0}(x)+\sqrt{{n}_{0}(x)}\hat{\phi }(x),\;\hat{\phi }(x)=\hat{\varphi }(x)+{\hat{\varphi }}^{\dagger }(x).\end{eqnarray} \tag{ 25 }$$

The expression of $\hat{\phi }(x)$ in terms of the BdG scattering states is formally similar to that of $\hat{\varphi }(x),$ equation (2):

$$\begin{eqnarray}\hat{\phi }(x) & = & {\displaystyle \int }_{0}^{\infty }{\rm{d}}\omega \displaystyle \sum _{I}{r}_{I-{\rm{in}},\omega }(x){\hat{\gamma }}_{I-{\rm{in}}}(\omega )\\ & & +{\displaystyle \int }_{0}^{{\omega }_{{\rm{max}}}}{\rm{d}}\omega \;{r}_{d2-{\rm{in}},\omega }(x){\hat{\gamma }}_{d2-{\rm{in}}}^{\dagger }(\omega )+{\rm{h.c.}}.\end{eqnarray} \tag{ 26 }$$

with ${r}_{a,\omega }(x)={u}_{a,\omega }(x)+{v}_{a,\omega }(x).$ Importantly, we note that the density in the BdG approximation is linear in the destruction operators, rather than quadratic. Thus, we can extract the quadratic correlations of section 3 by measuring density–density correlations [23, 24]. For that purpose, we restrict ourselves to asymptotic (subsonic and supersonic) regions, where the GP wave function adopts the form of equation (3). We refer to the region between the two asymptotic regions as the scattering region, where the sonic black hole is placed; i.e., within the scattering region the flow velocity crosses the sound velocity at least once. In the following, we assume that the scattering region is placed near x = 0 and its size is much smaller than the size of the asymptotic, homogeneous regions.

By taking the Fourier transform of the density in the asymptotic regions at $k\ne 0,$ we can get rid of the condensate signal and extract the phonon signal. For definiteness, we focus on the correlations $u-d2$ but the procedure for the other cases is similar. We define:

$$\begin{eqnarray}&&{\hat{n}}_{r}\left(k\right)\equiv \int {\rm{d}}x\;\hat{n}(x){f}_{r}^{*}\left(x\right){{\rm{e}}}^{-{\rm{i}}{kx}},\end{eqnarray} \tag{ 27 }$$

where ${f}_{r}\left(x\right)$ are normalized functions ($\int {\rm{d}}x| {f}_{r}(x){| }^{2}=1$) localized in the corresponding asymptotic homogeneous regions, sufficiently far from the scattering region. They represent the envelope of the Fourier transform, which has to be introduced to explicitly take into account the fact that the asymptotic subsonic and supersonic regions are placed in different spatial regions in a realistic situation. We choose them such that

$$\begin{eqnarray}&&{f}_{r}(x)=\displaystyle \frac{1}{\sqrt{{L}_{r}}}f\left(\displaystyle \frac{x-{x}_{r}}{{L}_{r}}\right),\end{eqnarray} \tag{ 28 }$$

with f a symmetric and real dimensionless function, x_r the point where the envelope is centered, and L_r the size of the wave packet, which is taken sufficiently large for the Fourier transform of the envelope function,

$$\begin{eqnarray}&&{f}_{r}(k)=\displaystyle \frac{1}{\sqrt{2\pi }}\int {\rm{d}}x\;{{\rm{e}}}^{-{\rm{i}}{kx}}{f}_{r}\left(x\right),\end{eqnarray} \tag{ 29 }$$

to be sufficiently well peaked at zero momentum. Since the subsonic (supersonic) region is placed at the left (right) of the BH, we have ${x}_{u}\lt 0$ (${x}_{d}\gt 0$).

Taking into account all the previous considerations, we obtain:

$$\begin{eqnarray}{\hat{n}}_{u}\left({k}_{u-{\rm{out}}}(\omega )\right) & \simeq & \sqrt{{n}_{u}}\left(\displaystyle \int \right.{\rm{d}}{\omega }^{\prime }\;{f}_{u}^{*}({k}_{u-{\rm{out}}}({\omega }^{\prime })-{k}_{u-{\rm{out}}}(\omega ))\displaystyle \frac{{r}_{u-{\rm{out}}}({\omega }^{\prime })}{{\left|{w}_{u-{\rm{out}}}\left(\omega ^{\prime} \right)\right|}^{1/2}}{\hat{\gamma }}_{u-{\rm{out}}}({\omega }^{\prime })\\ & & +{f}_{u}^{*}(-{k}_{u-{\rm{in}}}({\omega }^{\prime })-{k}_{u-{\rm{out}}}(\omega ))\displaystyle \frac{{r}_{u-{\rm{in}}}({\omega }^{\prime })}{{| {w}_{u-{\rm{in}}}\left(\omega ^{\prime} \right)| }^{1/2}}{\hat{\gamma }}_{u-{\rm{in}}}^{\dagger }\left.({\omega }^{\prime })\right),\\ {\hat{n}}_{d}\left(-{k}_{d2-{\rm{out}}}(\omega )\right) & \simeq & \sqrt{{n}_{d}}\left({\displaystyle \int }_{0}^{{\omega }_{{\rm{max}}}}\right.{\rm{d}}{\omega }^{\prime }\;{f}_{d}^{*}({k}_{d2-{\rm{out}}}(\omega )-{k}_{d2-{\rm{out}}}({\omega }^{\prime }))\displaystyle \frac{{r}_{d2-{\rm{out}}}({\omega }^{\prime })}{{| {w}_{d2-{\rm{out}}}\left(\omega ^{\prime} \right)| }^{1/2}}{\hat{\gamma }}_{d2-{\rm{out}}}({\omega }^{\prime })\\ & & +\displaystyle \int \;{\rm{d}}{\omega }^{\prime }\;{f}_{d}^{*}({k}_{d2-{\rm{out}}}(\omega )-{k}_{d1-{\rm{out}}}({\omega }^{\prime }))\displaystyle \frac{{r}_{d1-{\rm{out}}}({\omega }^{\prime })}{{| {w}_{d1-{\rm{out}}}\left(\omega ^{\prime} \right)| }^{1/2}}{\hat{\gamma }}_{d1-{\rm{out}}}^{\dagger }\left.({\omega }^{\prime })\phantom{\rule[-0000]{0ex}{1.28em}}\right),\end{eqnarray} \tag{ 30 }$$

where ${r}_{a}(\omega )$ is (see equation (5))

$$\begin{eqnarray}&&{r}_{a}(\omega )={u}_{a}(\omega )+{v}_{a}(\omega )=\sqrt{\displaystyle \frac{{k}_{a}^{2}}{2{{\rm{\Omega }}}_{a}(\omega )}},\end{eqnarray} \tag{ 31 }$$

and we have used $| \omega -{v}_{r}{k}_{a}(\omega )| ={{\rm{\Omega }}}_{a}(\omega )=\sqrt{{c}_{r}^{2}{k}^{2}+\frac{{k}^{4}}{4}}.$ Now, we can connect the correlations studied in section 3 with the density correlations by taking into account that ${\hat{n}}_{r}^{\dagger }(k)={\hat{n}}_{r}(-k).$ We obtain:

$$\begin{eqnarray}{G}_{{uu}}(\omega ) & \equiv & \langle {\hat{n}}_{u}\left(-{k}_{u-{\rm{out}}}\left(\omega \right)\right){\hat{n}}_{u}\left({k}_{u-{\rm{out}}}\left(\omega \right)\right)\rangle ={n}_{u}{r}_{u-{\rm{out}}}^{2}(\omega )[{g}_{{uu}}(\omega )+{n}_{u}({\omega }_{{uu}}(\omega ))+1],\\ {G}_{d2d2}(\omega ) & \equiv & \langle {\hat{n}}_{d}\left({k}_{d2-{\rm{out}}}(\omega )\right){\hat{n}}_{d}\left(-{k}_{d2-{\rm{out}}}(\omega )\right)\rangle ={n}_{d}{r}_{d2-{\rm{out}}}^{2}(\omega )[{g}_{d2d2}(\omega )+{g}_{d1d1}({\omega }_{d1d2}(\omega ))+1],\\ {G}_{{ud}2}(\omega ) & \equiv & \langle {\hat{n}}_{d}\left(-{k}_{d2-{\rm{out}}}(\omega )\right){\hat{n}}_{u}\left({k}_{u-{\rm{out}}}\left(\omega \right)\right)\rangle =F(\omega )\sqrt{{n}_{u}{n}_{d}}{r}_{u-{\rm{out}}}{r}_{d2-{\rm{out}}}{c}_{{ud}2}(\omega ),\end{eqnarray} \tag{ 32 }$$

with ${\omega }_{d1d2}(\omega )={\omega }_{d1}({k}_{d2-{\rm{out}}}(\omega ))$ and ${\omega }_{{uu}}(\omega )={\omega }_{u}(-{k}_{u-{\rm{out}}}(\omega )),$ ${\omega }_{i}(k)$ being the dispersion relation of the mode i. Their corresponding values are obtained graphically in figure 2. On the other hand, the overlap function between the subsonic and supersonic regions, $F(\omega ),$ is given by:

$$\begin{eqnarray}F(\omega ) & = & \displaystyle \int {\rm{d}}k\;{f}_{u}^{*}\left(\displaystyle \frac{k}{\zeta (\omega )}\right){f}_{d}^{*}\left(k\zeta (\omega )\right),\\ & = & \displaystyle \frac{1}{\sqrt{{L}_{u}{L}_{d}}}\displaystyle \int {\rm{d}}x\;{f}^{*}\left(\displaystyle \frac{\zeta \left(\omega \right)x+{x}_{u}}{{L}_{u}}\right){f}^{*}\left(\displaystyle \frac{\zeta {\left(\omega \right)}^{-1}x-{x}_{d}}{{L}_{d}}\right),\end{eqnarray} \tag{ 33 }$$

with

$$\begin{eqnarray}&&\zeta (\omega )\equiv {\left|\displaystyle \frac{{w}_{u-{\rm{out}}}(\omega )}{{w}_{d2-{\rm{out}}}(\omega )}\right|}^{1/2}.\end{eqnarray} \tag{ 34 }$$

Zoom In Zoom Out Reset image size

The previous integral can be interpreted as the scalar product of two normalized functions, which satisfies $F(\omega )\leqslant 1.$ It is easy to prove that $F(\omega )=1$ when

$$\begin{eqnarray}&&\zeta (\omega )=\sqrt{\displaystyle \frac{{L}_{u}}{{L}_{d}}}=\sqrt{\left|\displaystyle \frac{{x}_{u}}{{x}_{d}}\right|}.\end{eqnarray} \tag{ 35 }$$

This condition has been interpreted as the condition for the wave-packets to have equal width in frequency space [17]. However, resulting from the well-known properties of the spatial density correlations, equation (35) can be also interpreted as expressing that the envelopes in the subsonic and supersonic regions have to be placed along the correlation lines that maximize the spatial density correlation function, see [7, 15, 20].

In the following, we suppose that condition (35) is satisfied and $F(\omega )=1.$ It is worth noting that the correlation functions ${G}_{{ij}}(\omega )$ can never violate the CS inequality by themselves since ${G}_{{ud}2}\lt \sqrt{{G}_{{uu}}{G}_{d2d2}}.$ Thus, we try to relate these functions, which can be measured, to the phonon correlations considered in the previous section. For that purpose, we normalize the density–density correlations of equation (32):

$$\begin{eqnarray}&&{{\mathcal{G}}}_{{ij}}(\omega )\equiv \displaystyle \frac{{G}_{{ij}}(\omega )}{\sqrt{{n}_{i}{n}_{j}}{r}_{i-{\rm{out}}}(\omega ){r}_{j-{\rm{out}}}},\end{eqnarray} \tag{ 36 }$$

where ${n}_{i,j}$ must be interpreted as ${n}_{i,j}={n}_{d}$ when the indices $i,j$ take values $d1,d2.$ From equation (32), we see that ${{\mathcal{G}}}_{{ud}2}(\omega )={c}_{{ud}2}(\omega ).$ For the extraction of the other correlations, in a similar way to [24], we define the magnitudes:

$$\begin{eqnarray}{\tilde{g}}_{{uu}}(\omega ) & \equiv & {g}_{{uu}}(\omega )+{n}_{u}({\omega }_{{uu}}(\omega ))={{\mathcal{G}}}_{{uu}}(\omega )-1,\\ {\tilde{g}}_{d2d2}(\omega ) & \equiv & {g}_{d2d2}(\omega )+{g}_{d1d1}({\omega }_{d1d2}(\omega ))={{\mathcal{G}}}_{d2d2}(\omega )-1.\end{eqnarray} \tag{ 37 }$$

The $-1$ appearing in the definitions of ${\tilde{g}}_{{uu}}(\omega ),{\tilde{g}}_{d2d2}(\omega )$ reflects the subtraction of the atomic depletion contribution.

We see that ${\tilde{g}}_{{uu}},{\tilde{g}}_{d2d2}$ are over-estimations of the correlation functions ${g}_{{uu}},{g}_{d2d2},$ since ${n}_{u},{g}_{d1d1}$ are always positive. From this, we can define

$$\begin{eqnarray}&&{\tilde{{\rm{\Delta }}}}_{{ud}2}\equiv | {c}_{{ud}2}{| }^{2}-{\tilde{g}}_{{uu}}{\tilde{g}}_{d2d2},\end{eqnarray} \tag{ 38 }$$

noting that ${{\rm{\Delta }}}_{{ud}2}\gt {\tilde{{\rm{\Delta }}}}_{{ud}2}.$ Thus, measuring ${\tilde{{\rm{\Delta }}}}_{{ud}2}\gt 0$ amounts to experimentally observing the quadratic CS violation. We note that, for these calculations, the sole assumption has been made that the state of the system is stationary.

The upshot of this discussion is that we can observe the violation of a quadratic CS inequality through the measurement of the function ${\tilde{{\rm{\Delta }}}}_{{ud}2}.$ Similar claims can be made about the measurement of the quadratic CS inequality involving the $d1,d2$ modes. We can repeat the same strategy as before to obtain the correlation functions ${\tilde{g}}_{d1d1},{c}_{d1d2}.$ However, it is not possible to obtain the correlation functions ${c}_{{uu}},{c}_{d1d1},{c}_{d2d2},{g}_{{ud}2},{g}_{d1d2}$ from this scheme. This is due to the vanishing overlap integrals that appear when trying to obtain the corresponding correlations. For instance, if we compute $\langle {\hat{n}}_{u}\left({k}_{u-{\rm{out}}}\left(\omega \right)\right){\hat{n}}_{u}\left({k}_{u-{\rm{out}}}\left(\omega \right)\right)\rangle ,$ we face an overlap integral of the type

$$\begin{eqnarray}&&\int {\rm{d}}k\;{f}_{u}^{*}\left(k\right){f}_{u}^{*}\left(k\right)=\int {\rm{d}}x\;{f}_{u}^{*}\left(x\right){f}_{u}^{*}\left(-x\right)=0,\end{eqnarray} \tag{ 39 }$$

because the f_u function is, by construction, well localized in the subsonic region, far from the scattering region around x = 0, see equation (27) and the discussion below. Thus, we cannot obtain c_uu by this method. A similar reasoning applies to other correlations.

It is important to remark on the crucial role played by the spatial location of the asymptotic regions. If we did not introduce the envelopes in equation (27) and rather did take the Fourier transforms as ideally infinite, we would obtain:

$$\begin{eqnarray}{\hat{n}}_{u}\left({k}_{u-{\rm{out}}}(\omega )\right) & \simeq & \sqrt{{n}_{u}}\left(\displaystyle \frac{{r}_{u-{\rm{out}}}(\omega )}{{\left|{w}_{u-{\rm{out}}}\left(\omega \right)\right|}^{1/2}}{\hat{\gamma }}_{u-{\rm{out}}}(\omega )+\displaystyle \frac{{r}_{u-{\rm{in}}}({\omega }_{{uu}}(\omega ))}{{| {w}_{u-{\rm{in}}}\left({\omega }_{{uu}}(\omega )\right)| }^{1/2}}{\hat{\gamma }}_{u-{\rm{in}}}^{\dagger }({\omega }_{{uu}}(\omega ))\right),\\ {\hat{n}}_{d}\left(-{k}_{d2-{\rm{out}}}(\omega )\right) & \simeq & \sqrt{{n}_{d}}\left(\displaystyle \frac{{r}_{d2-{\rm{out}}}(\omega )}{| {w}_{d2-{\rm{out}}}\left(\omega \right){| }^{1/2}}{\hat{\gamma }}_{d2-{\rm{out}}}(\omega )+\displaystyle \frac{{r}_{d1-{\rm{out}}}({\omega }_{d1d2}(\omega ))}{| {w}_{d1-{\rm{out}}}\left({\omega }_{d1d2}(\omega )\right){| }^{1/2}}{\hat{\gamma }}_{d1-{\rm{out}}}^{\dagger }({\omega }_{d1d2}(\omega ))\right).\end{eqnarray} \tag{ 40 }$$

We can infer that such an approach would lead to contradictory results. For instance, the commutator

$$\begin{eqnarray}&&[{\hat{n}}_{u}\left({k}_{u-{\rm{out}}}(\omega )\right),{\hat{n}}_{u}\left(-{k}_{u-{\rm{in}}}(\omega )\right)]={n}_{u}\displaystyle \frac{{r}_{u-{\rm{out}}}(\omega )}{{\left|{w}_{u-{\rm{out}}}\left(\omega \right)\right|}^{1/2}}\displaystyle \frac{{r}_{u-{\rm{in}}}(\omega )}{{| {w}_{u-{\rm{in}}}\left(\omega \right)| }^{1/2}}{S}_{{uu}}(\omega )\ne 0\end{eqnarray} \tag{ 41 }$$

would be non-vanishing, in contradiction with the well-known fact that it has to be zero since $[\hat{n}(x),\hat{n}(x^{\prime} )]=0$ for every $x,x^{\prime} .$ By way of contrast, when using the full expressions of equation (30) we arrive at the correct result

$$\begin{eqnarray}[{\hat{n}}_{u}\left({k}_{u-{\rm{out}}}(\omega )\right),{\hat{n}}_{u}\left(-{k}_{u-{\rm{in}}}(\omega )\right)] & = & {n}_{u}\displaystyle \frac{{r}_{u-{\rm{out}}}(\omega )}{{\left|{w}_{u-{\rm{out}}}\left(\omega \right)\right|}^{1/2}}\displaystyle \frac{{r}_{u-{\rm{in}}}(\omega )}{{| {w}_{u-{\rm{in}}}\left(\omega \right)| }^{1/2}}{S}_{{uu}}(\omega )\\ & & \times \displaystyle \int {\rm{d}}{\omega }^{\prime }\;{f}_{u}^{*}\left(-\displaystyle \frac{{\omega }^{\prime }-\omega }{| {w}_{u-{\rm{in}}}\left(\omega \right)| }\right){f}_{u}^{*}\left(-\displaystyle \frac{{\omega }^{\prime }-\omega }{| {w}_{u-{\rm{out}}}\left(\omega \right)| }\right).\end{eqnarray} \tag{ 42 }$$

The integral of the rhs, after returning to real space, is similar to that considered in equation (39) and gives zero.

We arrive at the same conclusion when considering the density–density correlations in real space. As noted in [20], only the terms with a stationary phase should be kept when computing $\langle \hat{n}(x)\hat{n}(x^{\prime} )\rangle ,$ which yields the condition (35). In particular, only the 'proper' correlations between out-out modes can be extracted ${c}_{{ud}2},{c}_{d1d2},{g}_{{ud}1}.$ The other correlations ${c}_{{uu}},{c}_{d1d1},{c}_{d2d2},{g}_{{ud}2},{g}_{d1d2},{c}_{{ud}1},$ even when they are non-zero, cannot be obtained because the associated exponential terms do not present a stationary phase when integrating over frequencies.

We conclude that an important consequence of taking into account the spatial location of the subsonic and supersonic regions is that only a limited number of correlations can be experimentally observed. In particular, this implies that we cannot measure all the correlation functions appearing in the GPH criterion (A.11) by this procedure. Moreover, and for the same reason, if the state is Gaussian, we cannot obtain all the correlation functions appearing in the quartic CS violation, see equation (A.17).

The results of the present section show that, within this kind of detection schemes, we can only aim at observing a quadratic CS violation. One can expect this limitation to be not exclusive of the specific procedure here considered. The reason is that any realistic attempt to obtain the correlation functions between phonons must necessarily take into account the spatial location of the asymptotic regions, which implies that a similar reasoning does apply.

Nevertheless, the previous considerations do not pose a problem for the detection of the quantum–Hawking effect for two reasons: (a) if the state of the system is Gaussian and incoherent over the incoming channels, we have proven that the GPH criterion and the quartic CS violation are equivalent to the quadratic CS violation and (b) even if the state of the system does not belong to that class, the quadratic CS violation is still a signature of the presence of the entanglement; in particular, it is a sufficient condition for the fulfillment of the GPH criterion. Even more than that, as remarked at the end of section 3, the quadratic CS violation is also by itself a clear indication of the quantum nature of the system, as CS inequalities are always satisfied in a classical system.

We investigate here the possibility of experimental detection of CS violation in different BH setups. For that purpose, we compute the measurable quantity ${\tilde{{\rm{\Delta }}}}_{{ud}2},$ defined in equation (38). We suppose that the state of the system is given by a thermal distribution in the incoming modes, as expressed in equation (24). Since we only focus on the identification of some physical trends and not on the study of the whole parameter space, characterized by seven variables (see appendix B), we will consider for simplicity some typical configurations studied in the literature in order to compute ${\tilde{{\rm{\Delta }}}}_{{ud}2}.$ For this analysis we distinguish between non-resonant and resonant structures. For the non-resonant case, we consider in this work two scenarios: the single delta barrier configuration and the waterfall configuration, schematically depicted in figure 3.

Zoom In Zoom Out Reset image size

In the delta barrier configuration, the black hole forms near a localized potential of the form $V(x)=Z\delta (x)$ [14, 15, 51]; see left panel of figure 3. This configuration permits us to study theoretically the flow of a condensate through a potential barrier of finite size. On the other hand, the waterfall configuration creates the sonic horizon by using a negative step potential $V(x)=-{V}_{0}\theta (x),$ where $\theta (x)$ is the step function [15]. This kind of scenario, where the black hole is produced by introducing a negative potential, has been already experimentally realized by the Technion group [12].

In figure 4, we represent the Hawking spectrum, $| {S}_{{ud}2}{| }^{2}$ and the functions ${{\rm{\Delta }}}_{{ud}2},{\tilde{{\rm{\Delta }}}}_{{ud}2}$ for the delta barrier and the waterfall setups, both at temperature $T=0.2,$ which is of the order of magnitude of typical experimental setups, where it can be as low as $T\sim 0.1$ [13]. We see that, for some range of frequencies, the measurable function ${\tilde{{\rm{\Delta }}}}_{{ud}2}$ satisfies ${\tilde{{\rm{\Delta }}}}_{{ud}2}\gt 0,$ which implies the presence of CS violation. As noted previously, ${\tilde{{\rm{\Delta }}}}_{{ud}2}\gt 0$ is a more restrictive condition than the bare quadratic CS violation ${{\rm{\Delta }}}_{{ud}2}\gt 0.$ However, we observe that for sufficiently large ω, both ${{\rm{\Delta }}}_{{ud}2},{\tilde{{\rm{\Delta }}}}_{{ud}2}$ converge to $| {S}_{{ud}2}{| }^{2},$ since for high ω we have that the occupation numbers are negligible and thus, we recover the zero-temperature limit, for which ${{\rm{\Delta }}}_{{ud}2}(T=0)=| {S}_{{ud}2}{| }^{2}.$ Moreover, we see that, at high ω, $| {S}_{d1d2}({\omega }_{d1d2}(\omega )){| }^{2}=0,$ because for ω such that ${\omega }_{d1d2}(\omega )\gt {\omega }_{{\rm{max}}}$ the anomalous scattering channel disappears, and then we have ${\tilde{{\rm{\Delta }}}}_{{ud}2}\simeq {{\rm{\Delta }}}_{{ud}2}(T=0)=| {S}_{{ud}2}{| }^{2}.$ Finally, we note that there is no violation near $\omega =0,$ as argued quite generally in [17].

Zoom In Zoom Out Reset image size

We now focus on resonant configurations. In particular, we focus on the case of a double delta barrier, which was first considered in [14]. In this situation, the black hole is formed by two single delta barriers separated by a distance d, so the potential is given by $V(x)=Z[\delta (x)+\delta (x-d)].$ This setup is schematically depicted in figure 5. It was shown in [17] that resonant spectra can be expected to present a strong signal of CS violation. This trend can be observed in figure 6, where we represent the same magnitudes as in figure 4 but now for a higher temperature T = 0.7. We see that, even for this relatively high temperature, the experimental signal ${\tilde{{\rm{\Delta }}}}_{{ud}2}$ is, at the resonance frequency, substantially larger than in the non-resonant case.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We have analyzed the range of validity of different existent criteria for the identification of deep quantum behavior as applied to the outgoing channels involved in the spontaneous HR. Specifically, we have compared the GPH criterion with the violation of quadratic and quartic CS inequalities. We have shown that, under certain physical conditions (Gaussian processes and, simultaneously, incoherent incoming channels), all the considered criteria are equivalent. When such conditions are not fulfilled, we have shown that the quartic CS violation and the non-separability of the state represent independent mathematical conditions.

We have also investigated the possible measurement of the different criteria in realistic scenarios. By taking into account the different spatial location of the subsonic and supersonic regions, we have shown that only certain correlation functions can be obtained. For simplicity we have focused only on one detection scheme, but we expect similar problems to arise in other kind of measurements, since in any realistic situation the supersonic and subsonic regions are necessarily placed in different regions. However, our work also shows that this limitation is not a major problem, as we can often measure the quadratic CS violation, which is also a sufficient condition for the entanglement of the system. Finally, our numerical results of section 5 show that, in typical analog configurations, the CS violation can be detected in an achievable range of temperatures.

Sonic black holes, as well as related setups such as BH lasers, have the potential to become an important element of atomtronic circuits, for instance, as a source of entangled quasiparticles in the broader context of quantum communication.

We thank A Amo, F Michel, R Parentani, D Guéry-Odelin and J Steinhauer for valuable discussions. This work has been supported by MINECO (Spain) through grants FIS2010-21372 and FIS2013-41716-P, and by Comunidad de Madrid through grant MICROSERES-CM (S2009/TIC-1476).

A.1. Cauchy–Schwarz inequality for operators

We briefly review the mathematical CS inequalities for operators in a Hilbert space, which are always fulfilled, see [43] for example. We start by defining the scalar product associated to a state $\hat{\rho }$ for two operator $\hat{A},\hat{B}$ as:

$$\begin{eqnarray}&&(\hat{A},\hat{B})\equiv \langle {\hat{A}}^{\dagger }\hat{B}\rangle =\mathrm{Tr}(\hat{\rho }{\hat{A}}^{\dagger }\hat{B}).\end{eqnarray} \tag{ A.1 }$$

It is easy to check that the previous product satisfies the usual properties of a scalar product. Thus, it satisfies the mathematical CS inequality:

$$\begin{eqnarray}&&{| \langle {\hat{A}}^{\dagger }\hat{B}\rangle | }^{2}\leqslant \langle {\hat{A}}^{\dagger }\hat{A}\rangle \langle {\hat{B}}^{\dagger }\hat{B}\rangle .\end{eqnarray} \tag{ A.2 }$$

In particular, the previous inequality implies for B = 1:

$$\begin{eqnarray}&&{| \langle \hat{A}\rangle | }^{2}={| \langle {\hat{A}}^{\dagger }\rangle | }^{2}\leqslant \langle {\hat{A}}^{\dagger }\hat{A}\rangle ,\end{eqnarray} \tag{ A.3 }$$

and, for $A={\hat{\gamma }}_{i}$ and $B={\hat{\gamma }}_{j},$

$$\begin{eqnarray}&&| {g}_{{ij}}{| }^{2}={| \langle {\hat{\gamma }}_{i}^{\dagger }{\hat{\gamma }}_{j}\rangle | }^{2}\leqslant \langle {\hat{\gamma }}_{i}^{\dagger }{\hat{\gamma }}_{i}\rangle \langle {\hat{\gamma }}_{j}^{\dagger }{\hat{\gamma }}_{j}\rangle ={g}_{{ii}}{g}_{{jj}}.\end{eqnarray} \tag{ A.4 }$$

Thus, the CS inequality $| {g}_{{ij}}{| }^{2}\leqslant {g}_{{ii}}{g}_{{jj}}$ is always satisfied and can never be violated. On the other hand, taking $A={\hat{\gamma }}_{i}^{\dagger }$ and $B={\hat{\gamma }}_{j},$ we arrive at:

$$\begin{eqnarray}&&| {c}_{{ij}}{| }^{2}={| \langle {\hat{\gamma }}_{i}{\hat{\gamma }}_{j}\rangle | }^{2}\leqslant \langle {\hat{\gamma }}_{i}{\hat{\gamma }}_{i}^{\dagger }\rangle \langle {\hat{\gamma }}_{j}^{\dagger }{\hat{\gamma }}_{j}\rangle =({g}_{{ii}}+1){g}_{{jj}}.\end{eqnarray} \tag{ A.5 }$$

Thus, the mathematical CS inequality states $| {c}_{{ij}}{| }^{2}\leqslant ({g}_{{ii}}+1){g}_{{jj}},$ leaving the possibility of quadratic CS violation, equation (12). The same type of argument holds for quartic CS violations, where an analogous reasoning leads to:

$$\begin{eqnarray}| {{\rm{\Gamma }}}_{{ij}}{| }^{2} & = & | \langle {\hat{\gamma }}_{i}^{\dagger }{\hat{\gamma }}_{i}{\hat{\gamma }}_{j}^{\dagger }{\hat{\gamma }}_{j}\rangle {| }^{2}\leqslant \langle {\hat{\gamma }}_{i}^{\dagger }{\hat{\gamma }}_{i}{\hat{\gamma }}_{i}^{\dagger }{\hat{\gamma }}_{i}\rangle \langle {\hat{\gamma }}_{j}^{\dagger }{\hat{\gamma }}_{j}{\hat{\gamma }}_{j}^{\dagger }{\hat{\gamma }}_{j}\rangle ,\\ & = & ({{\rm{\Gamma }}}_{{ii}}+{g}_{{ii}})({{\rm{\Gamma }}}_{{jj}}+{g}_{{jj}})\end{eqnarray} \tag{ A.6 }$$

which also leaves the possibility of quartic CS violations; see equation (10).

A.2. Quadratic CS violation and the GPH criterion

We now consider bipartite states of two modes $i,j.$ We define the vector $\hat{X}\equiv {[{\hat{q}}_{i},{\hat{p}}_{i},{\hat{q}}_{j},{\hat{p}}_{j}]}^{T},$ where the phase space operators ${\hat{q}}_{k},{\hat{p}}_{k}$ are related to the annihilation operator of the modes $k=i,j$ through ${\hat{\gamma }}_{k}=({\hat{q}}_{k}+i{\hat{p}}_{k})/\sqrt{2}.$ Following [42], we define the 4 × 4 matrix:

$$\begin{eqnarray}M & \equiv & V+{\rm{i}}\displaystyle \frac{W}{2}\\ W & = & \left[\begin{array}{cc}J & 0\\ 0 & J\end{array}\right],J=\left[\begin{array}{cc}0 & 1\\ -1 & 0\end{array}\right]\end{eqnarray} \tag{ A.7 }$$

where V is the covariance matrix with matrix elements ${V}_{\alpha \beta }=\frac{1}{2}\{{\rm{\Delta }}{\hat{X}}_{\alpha },{\rm{\Delta }}{\hat{X}}_{\beta }\},$ with ${\rm{\Delta }}\hat{X}\equiv \hat{X}-\langle \hat{X}\rangle$ and $\{{\rm{\Delta }}{\hat{X}}_{\alpha },{\rm{\Delta }}{\hat{X}}_{\beta }\}$ the anticommutator. The matrix J is the simplectic matrix in two-dimensions. The matrix M is non-negative, $M\geqslant 0,$ as this fact represents an alternative expression of the uncertainty principle [42].

If the state $\hat{\rho }$ is separable, it has the form (13), and by taking partial transpose of the state with respect to the subsystem j, we obtain an operator ${\hat{\rho }}_{t}$ which is also necessarily a physical state. The GPH criterion is based on this fact. It can be shown that if ${\hat{\rho }}_{t}$ is a physical state, the matrix M_t is also non-negative, where M_t is given by [42]:

$$\begin{eqnarray}{M}_{t} & \equiv & {V}_{t}+{\rm{i}}\displaystyle \frac{W}{2},\\ {V}_{t} & = & {\rm{\Lambda }}V{\rm{\Lambda }},\;{\rm{\Lambda }}\equiv {\rm{diag}}[1,1,1,-1].\end{eqnarray} \tag{ A.8 }$$

The condition ${M}_{t}\geqslant 0$ is the uncertainty principle for the state ${\hat{\rho }}_{t}.$ As M is non-negative, M_t is non-negative if and only if $\mathrm{det}{M}_{t}\geqslant 0.$ The conditions $\mathrm{det}M,{M}_{t}\geqslant 0$ are equivalent to ${{\mathcal{P}}}_{{ij}}^{\pm }\geqslant 0,$ respectively, where:

$$\begin{eqnarray}{{\mathcal{P}}}_{{ij}}^{\pm } & \equiv & \mathrm{det}{A}_{i}\mathrm{det}{A}_{j}+{\left(\displaystyle \frac{1}{4}\mp \mathrm{det}{C}_{{ij}}\right)}^{2}\\ & & -\mathrm{tr}({{JA}}_{i}{{JC}}_{{ij}}{{JA}}_{j}{{JC}}_{{ij}}^{T})-\displaystyle \frac{1}{4}(\mathrm{det}{A}_{i}+\mathrm{det}{A}_{j}).\end{eqnarray} \tag{ A.9 }$$

The matrices ${A}_{i},{A}_{j},{C}_{{ij}}$ are 2 × 2 submatrices of the covariance matrix V:

$$\begin{eqnarray}V=\left[\begin{array}{cc}{A}_{i} & {C}_{{ij}}\\ {C}_{{ij}}^{T} & {A}_{j}\end{array}\right].\end{eqnarray} \tag{ A.10 }$$

We note that ${{\mathcal{P}}}_{{ij}}^{+}\geqslant 0$ is always satisfied since $\hat{\rho }$ is a physical state. However, when ${{\mathcal{P}}}_{{ij}}^{-}\lt 0,$ the state ${\hat{\rho }}_{t}$ is not a physical state, which implies that the original state $\hat{\rho }$ is not separable. We can put together both conditions by defining the GPH function ${{\mathcal{P}}}_{{ij}}$ as:

$$\begin{eqnarray}{{\mathcal{P}}}_{{ij}} & \equiv & \mathrm{det}{A}_{i}\mathrm{det}{A}_{j}+{\left(\displaystyle \frac{1}{4}-| \mathrm{det}{C}_{{ij}}| \right)}^{2}\\ & & -\mathrm{tr}({{JA}}_{i}{{JC}}_{{ij}}{{JA}}_{j}{{JC}}_{{ij}}^{T})-\displaystyle \frac{1}{4}(\mathrm{det}{A}_{i}+\mathrm{det}{A}_{j}),\end{eqnarray} \tag{ A.11 }$$

Thus, if ${{\mathcal{P}}}_{{ij}}\lt 0,$ the state is entangled. This result is the GPH criterion. Note that when $\mathrm{det}{C}_{{ij}}\geqslant 0,$ the state is separable, as ${{\mathcal{P}}}_{{ij}}={{\mathcal{P}}}_{{ij}}^{+}\geqslant 0,$ so only states with $\mathrm{det}{C}_{{ij}}\lt 0$ can be entangled.

On the basis of the previous results, we now prove that the quadratic CS violation implies the fulfillment of the GPH criterion. Suppose that ${{\mathcal{P}}}_{{ij}}\geqslant 0.$ In that case, the matrix M_t is non-negative and we can define an associated scalar product for vectors $u,v\in {{\mathbb{C}}}^{4}$ as:

$$\begin{eqnarray}&&{(u,v)}_{t}\equiv {u}^{\dagger }{M}_{t}v\end{eqnarray} \tag{ A.12 }$$

which satisfies the associated CS inequality:

$$\begin{eqnarray}&&| {(u,v)}_{t}{| }^{2}\leqslant {(u,u)}_{t}{(v,v)}_{t}.\end{eqnarray} \tag{ A.13 }$$

Assuming the usual case where the annihilation and destruction operators have zero expectation values, the matrices ${A}_{i,j}$ and C_ij can be written in terms of the correlation functions of equation (11) as:

$$\begin{eqnarray}{A}_{k}=\left({g}_{{kk}}+\displaystyle \frac{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],\end{eqnarray} \tag{ A.14 }$$

$$\begin{eqnarray}{C}_{{ij}}=\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{eqnarray} \tag{ A.15 }$$

with the index $k=i,j$ and ${{\mathbb{I}}}_{2}$ the 2 × 2 identity matrix. By inserting

$$\begin{eqnarray}&&u=\displaystyle \frac{1}{\sqrt{2}}\left[\begin{array}{c}0\\ 0\\ 1\\ i\end{array}\right],\;v=\displaystyle \frac{1}{\sqrt{2}}\left[\begin{array}{c}1\\ i\\ 0\\ 0\end{array}\right]\end{eqnarray} \tag{ A.16 }$$

into equation (A.13), we obtain the quadratic CS inequality $| {c}_{{ij}}{| }^{2}\leqslant {g}_{{ii}}{g}_{{jj}}.$ Thus, if there is quadratic CS violation, the matrix M_t cannot be non-negative, which implies that ${{\mathcal{P}}}_{{ij}}\lt 0$ and the GPH criterion is satisfied. We conclude that the quadratic CS violation is a sufficient condition for the fulfillment of the GPH criterion.

A.3. Quartic CS violation and entanglement

The previous arguments cannot be applied to the quartic CS violation equation (10). As a counterexample, the direct product of two pure number states of the modes $i,j,$ $\hat{\rho }=| n\rangle \langle n| \otimes | n^{\prime} \rangle \langle n^{\prime} | ,$ which is a manifestly separable state, violates the quartic CS inequality. Even for Gaussian states, the two conditions, quartic CS violation and entanglement, are still independent. For instance, for a general Gaussian state, the quartic CS violation of equation (10) can be expressed, via Wick's theorem, in terms of the quadratic correlations:

$$\begin{eqnarray}&&{{\rm{\Theta }}}_{{ij}}=| {c}_{{ij}}{| }^{2}+| {g}_{{ij}}{| }^{2}+{g}_{{ii}}{g}_{{jj}}-\sqrt{2{g}_{{ii}}^{2}+{| {c}_{{ii}}| }^{2}}\sqrt{2{g}_{{jj}}^{2}+{| {c}_{{jj}}| }^{2}}.\end{eqnarray} \tag{ A.17 }$$

For Gaussian states the GPH condition is equivalent to the non-separability of the system. Using equations (A.11), (A.14), we can compute explicitly the GPH function:

$$\begin{eqnarray}{{\mathcal{P}}}_{{ij}} & = & [{g}_{{ii}}{g}_{{jj}}-{| {c}_{{ij}}| }^{2}][({g}_{{ii}}+1)({g}_{{jj}}+1)-{| {c}_{{ij}}| }^{2}]+4\left({g}_{{ii}}+\displaystyle \frac{1}{2}\right)\mathrm{Re}({g}_{{ij}}{c}_{{jj}}^{*}{c}_{{ij}})+4\left({g}_{{jj}}+\displaystyle \frac{1}{2}\right)\mathrm{Re}({g}_{{ij}}{c}_{{ii}}{c}_{{ij}}^{*})\\ & & -2\mathrm{Re}({c}_{{ij}}^{2}{c}_{{ii}}^{*}{c}_{{jj}}^{*})-2\mathrm{Re}({g}_{{ij}}^{2}{c}_{{ii}}{c}_{{jj}}^{*})+| {c}_{{ii}}{| }^{2}| {c}_{{jj}}{| }^{2}+| {g}_{{ij}}{| }^{4}-| {g}_{{ij}}{| }^{2}({g}_{{ii}}+{g}_{{jj}}+2{g}_{{ii}}{g}_{{jj}}+2{| {c}_{{ij}}| }^{2})\\ & & -{g}_{{ii}}({g}_{{ii}}+1)| {c}_{{jj}}{| }^{2}-{g}_{{jj}}({g}_{{jj}}+1)| {c}_{{ii}}{| }^{2}.\end{eqnarray} \tag{ A.18 }$$

This expression applies whenever $| {c}_{{ij}}| \geqslant | {g}_{{ij}}| .$ For $| {c}_{{ij}}| \lt | {g}_{{ij}}| ,$ we have $\mathrm{det}{C}_{{ij}}\gt 0$ and the state is separable, as previously explained.

As we can see, the expressions for the quartic CS violation and the GPH criterion represent different conditions for arbitrary Gaussian states. In particular, it is easy to find states violating the quartic CS inequality with $\mathrm{det}{C}_{{ij}}\gt 0,$ which means that they are separable. Also, we can find entangled states that satisfy the quartic CS inequality.

We discuss in this appendix the parameters needed to compute the correlation functions defined in the main text. A related discussion appeared in [22]. The scattering matrix S relates the 'out' scattering states with the 'in' scattering states through equation (7). The conservation of the commutation relations implies the relation:

$$\begin{eqnarray}&&{S}^{\dagger }\eta S=\eta \equiv {\rm{diag}}(1,1,-1),\end{eqnarray} \tag{ B.1 }$$

which means that $S\in U(2,1).$ It follows that ${S}^{-1}=\eta {S}^{\dagger }\eta$ and, using standard linear algebra, we can prove:

$$\begin{eqnarray}&&{S}_{{ij}}^{*}=\displaystyle \frac{{\eta }_{{ii}}{m}_{{ij}}{\eta }_{{jj}}}{{\rm{det}}\;S},\end{eqnarray} \tag{ B.2 }$$

where m_ij is the minor associated to the S matrix element S_ij. In order to simplify the notation, in this section we relabel the indices $u,d1,d2$ as $1,2,3$ in order to match the matrix indices ordering. With the help of the previous results, it can be proven that any matrix $S\in U(2,1)$ satisfies:

$$\begin{eqnarray}{BSA} & = & {S}_{{\rm{PA}}} & = & \left[\begin{array}{ccc}\mathrm{det}S & 0 & 0\\ 0 & {S}_{33}^{*} & N\\ 0 & N & {S}_{33}\end{array}\right],\\ & B & \equiv & \displaystyle \frac{1}{N} & \left[\begin{array}{ccc}-{S}_{23} & {S}_{13} & 0\\ {S}_{13}^{*} & {S}_{23}^{*} & 0\\ 0 & 0 & N\end{array}\right],\\ & A & \equiv & \displaystyle \frac{1}{N} & \left[\begin{array}{ccc}-{S}_{32} & {S}_{31}^{*} & 0\\ {S}_{31} & {S}_{32}^{*} & 0\\ 0 & 0 & N\end{array}\right],\end{eqnarray} \tag{ B.3 }$$

where $N=\sqrt{| {S}_{13}{| }^{2}+| {S}_{23}{| }^{2}}=\sqrt{| {S}_{31}{| }^{2}+| {S}_{32}{| }^{2}}=\sqrt{| {S}_{33}{| }^{2}-1}.$ We see that the matrices $A,B$ are unitary. Inverting the previous relation, we obtain:

$$\begin{eqnarray}&&S={B}^{\dagger }{S}_{{\rm{PA}}}{A}^{\dagger }.\end{eqnarray} \tag{ B.4 }$$

Equations (B.3) and (B.4) show explicitly that we only need nine real parameters to parametrize the S matrix: $| {S}_{31}| ,| {S}_{13}| ,| {S}_{33}| ,\chi ,{\phi }_{13},{\phi }_{23},{\phi }_{31},{\phi }_{32},{\phi }_{33},$ with ${\phi }_{{ij}}$ the phase of S_ij, ${S}_{{ij}}=| {S}_{{ij}}| {{\rm{e}}}^{{\rm{i}}{\phi }_{{ij}}}$ and χ that of $\mathrm{det}S,$ $\mathrm{det}S={{\rm{e}}}^{{\rm{i}}\chi }.$ That is, we use three amplitudes and six phases to characterize the complete scattering matrix. Interestingly, the matrix ${S}_{{\rm{PA}}}$ shows clearly that the HR acts as a non-degenerate parametric amplifier [19].

A convenient way to write the S matrix for our problem is:

$$\begin{eqnarray}S & = & {U}_{{\rm{out}}}^{\dagger }\tilde{S}{U}_{{\rm{in}}}\\ {U}_{k} & = & {\rm{diag}}[{{\rm{e}}}^{{\rm{i}}{\delta }_{1-k}},{{\rm{e}}}^{{\rm{i}}{\delta }_{2-k}},{{\rm{e}}}^{{\rm{i}}{\delta }_{3-k}}]\end{eqnarray} \tag{ B.5 }$$

(with $k={\rm{in,}}\;{\rm{out}}$), since we can fix the elements of ${U}_{{\rm{in,}}\;{\rm{out}}}$ in such a way that they absorb the phases ${\phi }_{31},{\phi }_{32},{\phi }_{13},{\phi }_{23},{\phi }_{33},$ and so ${\tilde{S}}_{31},{\tilde{S}}_{32},{\tilde{S}}_{13},{\tilde{S}}_{23},{\tilde{S}}_{33}$ become purely real. This means that the matrix $\tilde{S}$ is characterized by only four parameters, $| {S}_{31}| ,| {S}_{13}| ,| {S}_{33}| ,\tilde{\chi },$ with

$$\begin{eqnarray}&&\mathrm{det}\tilde{S}={{\rm{e}}}^{{\rm{i}}\tilde{\chi }}=\mathrm{det}S\;{{\rm{e}}}^{-{\rm{i}}({\phi }_{13}+{\phi }_{23}+{\phi }_{13}+{\phi }_{23}-{\phi }_{33})}.\end{eqnarray} \tag{ B.6 }$$

Nevertheless, we do not need all nine parameters for computing the physical quantities of the main text. For instance, the CS violations associated to the correlations of equation (14) are invariant under phase transformations of the 'out' states. It can be also shown that the GPH function is invariant under these transformations, see equation (A.18). As we are always considering incoherent incoming modes, all the quantities are also invariant under phase transformations of the incoming modes. From equation (B.5), and taking into account the previous observations, we conclude that all the requested quantities depend only on the matrix $\tilde{S},$ which is characterized by just 4 parameters, as noted in the previous paragraph. Using equations (B.1), (B.4) we can write $\tilde{S}$ as:

$$\begin{eqnarray}\tilde{S}=\left[\begin{array}{ccc}-\mathrm{sin}\alpha & \mathrm{cos}\alpha & 0\\ \mathrm{cos}\alpha & \mathrm{sin}\alpha & 0\\ 0 & 0 & 1\end{array}\right]\left[\begin{array}{ccc}{{\rm{e}}}^{{\rm{i}}\tilde{\chi }} & 0 & 0\\ 0 & \mathrm{cosh}\gamma & \mathrm{sinh}\gamma \\ 0 & \mathrm{sinh}\gamma & \mathrm{cosh}\gamma \end{array}\right]\left[\begin{array}{ccc}-\mathrm{sin}\beta & \mathrm{cos}\beta & 0\\ \mathrm{cos}\beta & \mathrm{sin}\beta & 0\\ 0 & 0 & 1\end{array}\right]\end{eqnarray} \tag{ B.7 }$$

with $| {S}_{33}| =\mathrm{cosh}\gamma ,| {S}_{13}| =\mathrm{sinh}\gamma \mathrm{cos}\alpha$, and $| {S}_{31}| =\mathrm{sinh}\gamma \mathrm{cos}\beta$. Thus, we only need four parameters in order to parametrize the scattering matrix for the calculations of this work: $\alpha ,\beta ,\gamma ,\tilde{\chi }$. We note that this parametrization differs slightly from that used in [22], where four amplitudes were used to characterize the state of the system. Nevertheless, that election is equivalent to that used here.

Finally, taking into account that the state of the system $\hat{\rho }$ is characterized by only three numbers, n_i (with $i=1,2,3$, see equation (15)), we have that the whole problem is completely determined by seven parameters, four arising from the S matrix and three arising from the specification of the incoherent incoming Gaussian state.