Violation of Bell inequalities in an analogue black hole
Giorgio Ciliberto Stephanie Emig Nicolas Pavloff Mathieu Isoard




Université Paris-Saclay, CNRS, LPTMS, 91405, Orsay, France
Physikalisches Institut, Albert-Ludwigs-Universität Freiburg,
Hermann-Herder-Straße 3, 79104 Freiburg, Germany

Physikalisches Institut, Albert-Ludwigs-Universität Freiburg,
Hermann-Herder-Straße 3, 79104 Freiburg, Germany

Université Paris-Saclay, CNRS, LPTMS, 91405, Orsay, France
Institut Universitaire de France

Laboratoire Kastler Brossel, Sorbonne Université, CNRS, ENS-PSL Research University,
Collège de France, 4 place Jussieu, F-75252 Paris, France

  Signals of entanglement and nonlocality are quantitatively evaluated
  at zero and finite temperature in an analogue black hole realized in
  the flow of a quasi one-dimensional Bose-Einstein condensate. The
  violation of Lorentz invariance inherent to this analog system opens
  the prospect to observe 3-mode quantum correlations and we study the
  corresponding violation of bipartite and tripartite Bell
  inequalities. It is shown that the long wavelength modes of the
  system are maximally entangled, in the sense that they realize a
  superposition of continuous variable versions of
  Greenberger-Horne-Zeilinger states whose entanglement resists
  partial tracing.


Violation of Bell inequalities in an analogue black hole
    Mathieu Isoard
    
========================================================





§ INTRODUCTION



The domain of analog gravity aims at providing laboratory models for
gaining insight on general relativity phenomena that cannot be
directly observed in the usual gravitational context or for which
there exists no complete theoretical framework. Two such phenomena are
black hole superradiance
<cit.> and Hawking radiation
<cit.>. In that context, the most successful
analogous experimental platforms have been surface gravity
<cit.> and
acoustic <cit.> waves, nonlinear optical systems
<cit.>, cavity polaritons
<cit.> and Bose-Einstein
condendates (BEC) of atomic vapors
<cit.>. Because of their low
temperature and coherence properties BECs appear particularly well
suited for demonstrating quantum features, for instance for deciding
whether the phenomenon of analogous Hawking radiation is mostly
stimulated or does actually correspond to spontaneous emission of
pairs of entangled (quasi-)particles as in Hawking’s original
scenario.  We do not directly address this point in the present work
(see, e.g., Ref. <cit.> for a recent discussion) but rather
take the theoretical analysis a little further by asking: which
general insight can we reach by studying quantum correlations of the
Hawking signal emitted by an analogue black hole? A natural approach
for such an investigation is a test of nonlocality via violation of
Bell inequality. The epistemological query of refutation of local
hidden variable theories has already received an unambiguous answer in
many contexts (see, e.g., Ref. <cit.> and references
therein) and important progresses have also been achieved in the field
of BEC matter waves
<cit.> we consider
here. Such a test is nonetheless a nontrivial extension of the scope
of analogue gravity and would constitute a primer for continuous
variables entanglement in a matter wave environment (see also Refs.
<cit.> for related proposals).  In view
of future experimental studies it is relevant to quantitatively
evaluate in realistic configurations to what extent Bell inequalities
can be violated in BEC analogues. This is a natural question to ask,
all the more so as we argue in the following that in some (exotic)
limits the analog black hole we consider exactly realizes a
Einstein-Podolsky-Rosen (EPR) pair.  Furthermore, as will be shown,
the specifics of the system provide a natural testing ground for
genuine tripartite nonlocality. Our theoretical investigation of
the matter reveals an unexpected generic feature of black hole
analogues: in the long wave-length limit, the state of the system
realizes an infinite sum of degenerate Greenberger-Horne-Zeilinger
(GHZ) states. Interestingly, thanks to the continuous nature of its
degrees of freedom, and despite the clear GHZ nature of its long
wavelength modes, the analog system remains entangled after partial
tracing.

The paper is organized as follows. Section <ref> presents the BEC analogue system we consider and the theoretical tools we employ for its study. Section <ref> is devoted to the study of measures of bipartite nonlocality and Section <ref> to tripartite nonlocality. We present our conclusions in Sec. <ref>. Technical aspects are summarized in the appendices. Appendices <ref> and <ref> are devoted to a brief presentation of previous results. 
The black hole analogue we consider can be modeled by an equivalent optical system 
 <cit.> whose relevance for the aspects we consider in the main text is discussed in Appendix <ref>. Appendix <ref> present elementary properties of the pseudo-spin we use in the main text. Appendices <ref>, 
<ref> and <ref> treat specific technical aspects of the procedures used in Secs. <ref>  and <ref> for maximizing the violation of bipartite and tripartite Bell inequalities. Appendix <ref> also concerns technical aspects, but these are of crucial importance for the computation of all the quantities we present and plot in the main text, such as the Clauser-Horne-Shimony-Holt (CHSH), Svetlichny and Mermin parameters. The last appendix (<ref>) presents results in analogue black hole configurations different from the one on which we focus in the main text.



§ THE MODEL


We consider a one-dimensional (1D) BEC described by a quantum field Ψ̂(x,t) solution of

    iħ∂_t Ψ̂= -ħ^2/2m∂_x^2Ψ̂+[U(x) + g  n̂ - μ]Ψ̂.

In this equation U(x) is an external potential whose precise form depends
on the black hole configuration considered, μ is the chemical potential, n̂(x,t)=Ψ̂^†Ψ̂ is the density operator and g is the nonlinear constant describing the effective interaction within the system. We consider an effective repulsion, with g>0. In the so-called 1D mean-field regime <cit.> the quantum field can be written as

    Ψ̂(x,t)=Φ(x)+ψ̂(x,t),

where Φ(x) describes the background flow while ψ̂
accounts for small quantum fluctuations. Whereas the decomposition
(<ref>) is legitimate in 3D, it has a finite range of
validity in the 1D configurations we consider; however, its conditions
of applicability are commonly met in standard experimental situations,
see, e.g., the discussion in <cit.>. In decomposition
(<ref>) the classical field Φ is a stationary solution
of the Gross-Pitaevskii equation, see Eq. (<ref>).  The
associated flow realizes an analogue black hole horizon if it is
subsonic in the asymptotic upstream region and supersonic in the
asymptotic downstream region. Several such flows have been considered
in the past. We will here focus on the “waterfall configuration”
which is close to the experimental realization of
Refs. <cit.>. In order to assess the
generality of our results, we present in Appendix <ref> the
results obtained for two other configurations, we denote as “delta
peak” and “flat profile”.  All these configurations are described in
Appendix <ref>.



 §.§ Propagation channels and quantum modes


The decomposition (<ref>) is meaningful in a regime of small quantum fluctuations where the operator ψ̂ can be treated within a Bogoliubov approach. In this case  ψ̂ is naturally expanded along the asymptotic ingoing and outgoing channels of the flow. The dispersion relation of elementary excitations in the asymptotic upstream subsonic and downstream supersonic regions 
(x→ -∞ and +∞, respectively) 
is of Bogoliubov type, with a Doppler shift accounting for the finite velocity of the background:

    (ω - V_α q)^2=ω_ B,α^2(q),

where α=u far upstream and α=d downstream. In this expression  V_α is the asymptotic velocity of the flow in region α and

    ω_ B,α(q)=c_α q √(1+q^2ξ^2_α/4)

is the Bogoliubov dispersion relation <cit.>. In Eq. (<ref>) c_α  and ξ_α are the asymptotic speed of sound and healing length in region α, respectively, see Appendix <ref>.
These dispersion relations are represented in Fig. <ref>.


From the identification of the relevant channels and of their
direction of propagation it is possible to define quantum modes
forming a basis enabling to describe all the elementary excitations of
the background flow. To each such mode is associated a quantum
operator: b̂_i^†(ω) and b̂_i(ω) (i=0,
1 or 2) are the creation and annihilation operators of an excitation
of energy ħω which is ingoing in channel i| in and
scattered by the horizon onto the three outgoing channels
0| out, 1| out and 2| out. Since each b̂
mode is associated with a single ingoing channel, it is denoted as an
“ingoing mode”. It is also relevant to define “outgoing modes”
associated with a single outgoing channel. The corresponding operators
are denoted as ĉ_i(ω) and
ĉ_i^†(ω). For instance ĉ^ †_0 is the
creation operator of an excitation where the three ingoing channels
are implied and form an outgoing excitation in channel 0| out.
The corresponding quantum mode is the analogous Hawking mode. The
outgoing modes are related to the incoming ones via the scattering
matrix S(ω):

    [   ĉ_0;   ĉ_1; ĉ_2^†;       ]
    =
    [ S_00 S_01 S_02; S_10 S_11 S_12; S_20 S_21 S_22;      ] [   b̂_0;   b̂_1; b̂_2^†;        ],

where all the ω dependencies have been omitted for legibility. 
The modes b̂_2 and ĉ_2 are particular in the sense that they have a negative norm and should be quantized inverting the usual role of the creation and annihilation operators <cit.> in order to satisfy the standard Bose commutation relations. The mode ĉ_2 is analogous to what is called the partner is the context of Hawking radiation. We denote the mode associated with ĉ_1 the companion; Lorentz invariance prevents such a mode to exist in a gravitational black hole, but it is unavoidable in an analog system.

The fact that the outgoing operators fulfil the canonical commutation relations implies that the scattering matrix S(ω) obeys the skew-unitarity condition

    S^†η S =η =S η S^† ,
            η= diag(1,1,-1)  .

For ω>Ω the mode with subscript 2 (the partner) disappears because the corresponding ingoing and outgoing channels do (cf. Fig. <ref>) and the S-matrix becomes 2× 2 and unitary. In this case the vacuum of the outgoing modes (the ĉ's) is identical to the vacuum of the incoming ones (the b̂'s) and the analogous Hawking effect disappears.
The value of the corresponding threshold energy is

    Ω=q^* V_d-ω_ B,d(q^*)   
        q^*ξ_d=(-2+m_d^2/2+
    m_d/2√(8+m_d^2))^1/2.


It has been shown in Ref. <cit.> that the three-mode system describing the analogous black hole horizon can be modeled by an optical setup simply composed of a parametric amplifier and a beam splitter, as depicted in Fig. <ref>.

Entanglement is localized in the two-mode squeezed state f_1|f_2. It is then dispatched by mean of a beam splitter which performs a non local transformation from the effective modes (f_0,f_1,f_2)
onto other effective modes (e_0,e_1,e_2) which connect to the physical outgoing modes (c_0,c_1,c_2) by a simple local linear unitary Bogoliubov 
transformation (LLUBO), see details in Appendix <ref>. We stress that the configuration depicted in Fig. <ref> captures the essence of the analogue black hole configuration and is generic: A similar model has been used in Ref. <cit.> for describing an optical system containing a pair white-black hole analog and then studied on general grounds in <cit.>.
Albeit the fact that we allow the down-conversion process to violate Lorentz invariance (since we 
enforce the Bogoliubov dispersion relation) the two-mode squeezed state built on f̂_1 and f̂_2 is the closest possible analogue to an ideal black hole. The additional beam splitter is inherent to any analogue model. The transmission coefficient of this beam splitter plays the role of the grey-body factor which is relevant in the gravitational context, see discussion in <cit.>.



 §.§ Density matrix and measure of bipartite entanglement


For describing analog black holes physics one should consider observables operating in the Fock space of the outgoing modes, as would 
an observer located outside the horizon of a gravitational black hole. However, the physical implementation of the analog black hole configurations
presented in Appendix <ref>
is realized in the vacuum  of the incoming modes.
This mismatch is the origin of the quantum evaporation process, as first presented by Hawking <cit.>.
An ingoing vacuum mode of frequency ω (which we denote as |0_ω⟩^ in) relates to an outgoing vacuum mode |0_ω⟩^ out 
through

    |0_ω⟩^ in = 1/|S_22|   
    e^(X_02 ĉ_0^†  + X_12 ĉ_1^†)  ĉ_2^†   |0_ω⟩^ out,

where the coefficients X_i2 (i=0 and 1) are expressed in terms of entries of the S-matrix (<ref>):
X_i2(ω)=S_i2(ω)/S_22(ω),
and

    |0_ω⟩^ out=|0_ω⟩_0^ out⊗
    |0_ω⟩_1^ out⊗ |0_ω⟩_2^ out,

where |0_ω⟩_j^ out is the vacuum
of operator ĉ_j(ω), with j∈{0,1,2}. 
 Eq. (<ref>) can be rewritten in terms of the number states of quasi-particles of type c_j 

    |n_ω⟩_j=1/√(n!)(ĉ_j^†(ω))^n
     |0_ω⟩_j^ out,

as

    |0⟩^ in = 1/|S_22|∑_μ,ν=0^∞√(μ+ νμ) X_ 02^μ X_ 12^ν
    |μ⟩_ 0|ν⟩_ 1|μ+ν⟩_ 2.

It will be often the case below that, as in expression (<ref>), we drop 
the explicit ω-dependence for legibility. 

Relationships (<ref>) and  (<ref>) ensure that the system is in a pure three modes Gaussian state which is fully characterized by its 6× 6 covariance matrix σ(ω) with entries

    σ_ℓ m≡1/2 ⟨ξ̂_ℓ ξ̂_m +
      ξ̂_m  ξ̂_ℓ⟩ - ⟨ξ̂_ℓ⟩ ⟨ξ̂_m ⟩,

where operator ξ̂_ℓ (or ξ̂_m) is one of the 6 components of the column vector 
ξ̂ = √(2) 
(q̂_0,p̂_0,q̂_1,p̂_1,q̂_2,p̂_2)^ T, where 

    q̂_j(ω)=1/√(2)(ĉ_j+ĉ_j^†)    p̂_j(ω)=i/√(2)(ĉ_j^†-ĉ_j).

In expression (<ref>) and in all the following, the
averages ⟨⋯⟩ are performed over the density
matrix of the system.  This density matrix is simply
ρ=|0 ⟩^ in  ^ in⟨ 0| in the ideal case just described.
We also consider more realistic situations where some incoherent
excitations are present and the system is not in a pure state. A
simple manner to account for this situation would be to assume that
the system is in a thermal state. This is however impossible because
the analog configurations depicted in Appendix <ref>
are thermodynamically unstable. A way to circumvent this problem has
been proposed in Refs. <cit.>. It consists in postulating
that the system was initially in thermal equilibrium at temperature T
with a constant density and velocity (n_u and V_u, respectively)
and that the flow has been adiabatically modified by slowly ramping the
appropriate external potential, eventually reaching the black 
hole configuration
of interest. This situation, although idealized, it less schematic
than the zero excitation regime. It emulates the experimental
situation of Refs. <cit.> if the system is
considered to have been initially in equilibrium in the frame attached to the
flowing condensate. In this case one has (for j∈{0,1,2})

    n̅_j(ω)≡⟨b̂^†_j(ω) b̂_j(ω) ⟩
    = n_ th[ ω_ B,α(  q_j| in(ω)  ) ],

where n_ th(ω) is the thermal Bose occupation
distribution at temperature T and energy ħω, whereas
ω_ B,α(q_j| in) is the Bogoliubov
dispersion relation (<ref>), with α = u if j=0 and
α=d if j=1 or 2. The functions q_j| in(ω)
appearing in expression (<ref>) are pictorially defined in
Fig. <ref>. For instance q_2| in(ω) is the function
that, to a given angular frequency ω∈[0,Ω], associates a
wave-vector along the 2| in dispersion branch[
  q_2| in(ω)∈[q^*,q_0] with q^*=q_2| in(Ω)
  and q_0=q_2| in(0), see Fig. <ref>.]. We loosely refer
to the cases where n̅_j (j∈{0,1,2}) is equal to the
right-hand side (r.h.s.) of (<ref>) as “finite
temperature” situations. In the simplest configuration, denoted as
“zero temperature” the system is in the pure state
|0⟩^ in and the n̅_j's are all equal to zero.  We
will not specify the values of the n̅_j's in the following, so
that the formulas we give are generally valid, even in situations
where the occupation numbers should not be given by formulas of the
type of Eq. (<ref>). However, for illustrative purposes,
all the figures of the paper are plotted in specific temperature
cases.


Several observables have been proposed to theoretically evaluate the
bipartite entanglement in the context of analogue gravity, such as the
Cauchy-Schwarz criterion
<cit.>,
the generalized Peres-Horodecki parameter
<cit.>, the logarithmic negativity
<cit.>, the
entanglement entropy <cit.>, the
entanglement of formation <cit.>, the Gaussian contangle
<cit.>.  In the present work, for reasons exposed in
Appendix <ref>, we chose to evaluate the bipartite
entanglement between modes i and j by means of a quantity
Λ^(i|j)(ω) ∈ ]-∞,1] which is a monotonous
measure of entanglement that has been denoted as the “PPT measure” in
Ref. <cit.>.  States i and j are separable when
Λ^(i|j)<0, which is always the case when (i|j)=(0|1): the
companion and the Hawking modes are not entangled. The two other
couples of modes, (0|2) and (1|2), are always entangled at T=0
for all ω. Their entanglement decreases with increasing
temperature by an amount specified by the PPT measure which reads
explicitly (i=0 or 1) <cit.>:

    Λ^(i|2)(ω)=     
        - ⟨ĉ_i^† ĉ_i ⟩
        - ⟨ĉ_2^† ĉ_2 ⟩
        +
        √((⟨ĉ_i^† ĉ_i ⟩
        - ⟨ĉ_2^† ĉ_2 ⟩)^2 +
        4 |⟨ĉ_i  ĉ_2 ⟩|^2
        )  .

As clear from the above expression, a key ingredient for characterizing entanglement is the determination of average values of different combination of two creation or annihilation operators of the outgoing 
modes. At the experimental level, this determination could reveal difficult for quantities such as ⟨ĉ_0  ĉ_2 ⟩ for instance. Steinhauer proposed a possible way to extract this information from the knowledge of the density-density correlation function <cit.>. As stressed in Refs. 
<cit.> this method needs to be used with more care than initially thought, but is indeed a possible 
manner to obtain the information. At the theoretical level, it is a straightforward matter to compute the expectation values of products of the ingoing creation and annihilation operator (quantities such as (<ref>) for instance). From there, expression (<ref>) makes it possible to compute the equivalent expressions for the 
outgoing operators, which are the quantities of interest. The relevant formulae are given in Appendix <ref>, Eqs. (<ref>).




§ BIPARTITE NONLOCALITY


Most Bell-like inequalities proposed in the context of continuous variables hinge on a discretization process <cit.>. Indeed, since the set of outcomes for a given observable is typically unbounded in continuous variable systems, it seems a priori difficult to derive an upper bound of the expectation value of a Bell-type observable. It is not impossible, though: relying on the  Fine–Abramsky–Brandenburger theorem <cit.>, 
the authors of Refs. <cit.> were able to derive continuous Bell inequalities for continuous and unbounded observables.
However, to our knowledge, no practical (theoretical or experimental) use of such fully continuous approaches has been successfully implement so far.
A method to circumvent this difficulty is based on discretization and consists of using a dichotomic binning of the outcome results. In this case, by using observables which rely on continuous measurements but can only take a finite number of outcomes, one can construct a Bell inequality similar to those derived for discrete variables. A well-known example of such observables are the so-called pseudo-spin operators: they live in an infinite-dimensional Hilbert space  but the outcome of their measurement is either -1 or +1.

Due to its practicability, we chose to follow this discretization approach to derive Bell-like inequalities in the context of BEC
analogue black holes. In this work we use the GKMR (Gour, Khanna, Mann and Revzen) pseudo-spins introduced in Ref. <cit.>:
To an outgoing mode j (j∈{0,1,2}) of energy ω is associated a Hermitian vectorial operator Π̂^(j)(ω) with Cartesian coordinates


    Π̂^(j)_x(ω) = ∫_0^+∞ dq  (
    |q⟩_j _j⟨ q|-|-q⟩_j _j⟨ -q|
    ), 
       Π̂^(j)_y(ω) = i ∫_0^+∞ dq (
    |q⟩_j _j⟨ -q|-|-q⟩_j _j⟨ q|
    ), 
       Π̂^(j)_z(ω) = ∫_-∞^+∞ dq  
    |q⟩_j _j⟨ -q|,


where |q⟩_j is the eigenstate  associated to the eigenvalue q of the position operator q̂_j(ω) (<ref>).
The operators (<ref>) anti-commute with each others and all square to unity. They verify the expected spin commutation relations, such as

    [Π̂_x^(j),Π̂_y^(j)]=2   i Π̂_z^(j),

and similar relations upon circular permutations of the indices x, y and z. 
In Appendix <ref> we recall the properties of the pseudo-spin operator and of its eigenvalues which are useful in the following.  

These pseudo-spin operators have been studied in contexts similar to ours in Refs. <cit.>.
Compared with other pseudo-spin operators such as for instance those introduced by Banaszek and Wodkiewicz
<cit.>, the GKMR spins (<ref>) have the advantage of having simple Wigner transforms, which makes the computation of their expectation values over Gaussian states relatively easy, as detailed in Appendix <ref>. 
We stress here that different choices of spin representation lead to different values of averages of the Bell operators (<ref>) <cit.>. From this remark naturally results that, generally speaking, the observables (<ref>), (<ref>) and (<ref>)  we use below are witnesses of nonlocality: Their violation of Bell-type inequalities is a sufficient but not necessary test of nonlocal behavior.
As a final remark, we note here an unforeseen benefit of the use of the pseudo-spin operators (<ref>): it will be shown in  Sec. <ref> that the structure 
of the zero temperature ground state of the analog system (the vacuum |
0⟩^ in)
is most easily analyzed in terms of a combination of eigenstates of operator Π̂_x.

Equipped with the pseudo-spin operators (<ref>) we can define a CHSH Bell operator <cit.> measuring the correlations between the emitted quasi-particles of type i and j:

    ℬ̂^(i|j)(ω) =     (a+ a')·Π̂^(i)⊗b·Π̂^(j)
    
       +    
    (a-a')·Π̂^(i)⊗b'·Π̂^(j) ,

where a, a', b and b' are unit vectors.
Given a unit vector n it is easily checked that (n·Π̂^(j))^2=1, meaning that the Hermitian operator n·Π̂^(j) has eigenvalues ± 1.
It then follows from direct inspection that, if one believes in local realism, one should expect that a measure of the operator ℬ̂^(i|j) yields a result ± 2 <cit.>.
The tenant of local realism is thus violated when ⟨ℬ̂^(i|j)⟩>2, whereas Cirel’son bound <cit.> imposes
⟨ℬ̂^(i|j)⟩≤ 2√(2).  Since the modes 0 and 1 are not entangled the quantity ⟨ℬ̂^(0|1)⟩ is always lower than 2 and its computation is of no interest to us. For attempting to violate as much as possible Bell inequality one should consider the modes 0 and 2 (or 1 and 2) and look for an arrangement of the four measurement directions a, a', 
b and b' which maximizes ⟨ℬ̂^(i|2)⟩ (i=0 or 1). This procedure is explained in Appendix <ref> and makes it possible to analytically compute the quantity 

    B^(i|2)(ω)≡ _a,a',b, b'⟨ℬ̂^(i|2) (ω) ⟩.

The corresponding explicit expression is given in Eq. (<ref>). 



The values of B^(0|2) and B^(1|2) are plotted as functions of ω for different temperatures in Fig. <ref> for a waterfall configuration with a downstream Mach number [In this theoretical model configuration, fixing the value of m_d determines all the other dimensionless parameters: m_d=m_u^-2=n_u/n_d=V_d/V_u, see Ref. <cit.>.] 
m_d≡ V_d/c_d=2.9, the same as in the Technion 2019 experiment <cit.>. We also plot for comparison the values of the 
corresponding PPT measures  Λ^(0|2)  and Λ^(1|2), as defined by (<ref>). The figure illustrates that, as well known, 
entanglement is necessary but not sufficient for violating Bell inequality. Also, the amount by which the Bell inequality is 
violated is not proportional to the amount of entanglement. This is clearly seen, for instance, by comparing the values of 
B^(0|2) and Λ^(0|2) at T=0: the maximum violation of Bell inequality is not achieved for the maximal entanglement. 
One can also notice that the violation of Bell inequality is much less resilient to temperature than is the entanglement. These 
features can be most easily understood in the framework of the optical model represented in Fig. <ref>. They originate from the
the dilution of entanglement between the three modes caused by the beam splitter, as discussed in Appendix <ref>.

Fig. <ref> indicates that, in the waterfall configuration we consider (with m_d=2.9),
the entanglement is lower and the violation of Bell inequality less significant for the correlations among 
modes 1 and 2 than for the correlations among modes 0 and 2. However, this is not always the case. This point is illustrated in Fig. <ref>  which displays, for all waterfall configurations, the zero 
temperature values of max_ω B^(i|2) and max_ωΛ^(i|2) (for 
i=0 and 1) as functions of the upper Mach number m_u (the maximisation 
is performed at fixed m_u, for ω∈[0,Ω]). All the possible waterfall 
configurations are considered since m_u spans the whole interval [0,1]. In this figure we aim 
at evaluating the largest amount of entanglement and nonlocality reached in each configuration. This is the reason why
we plot the maximum values taken by the quantities 
Λ^(i|2) and B^(i|2) over the energy interval [0,Ω] [where Ω is defined 
by Relation (<ref>)] since this is only for energies in this interval that spontaneous emission of 
quasi-particles occurs.

For instance, the situation depicted in Fig. <ref> (m_d=2.9) corresponds in Fig. <ref> to the point m_u=0.587 since for waterfall configurations m_u and m_d are related by 
(<ref>). And indeed, Fig. <ref> shows that for this value of m_u the maximum over ω of B^(0|2) is 2.25, and the one of B^(1|2) is 2, as observed in Fig. <ref>.

Fig. <ref> shows that for values of m_u larger than 0.6, the entanglement is mainly concentrated between modes 2 and 0, 
i.e., between the Hawking quantum and the partner. This is indicated by the fact that both the PPT measure and the CHSH parameter significantly point to nonseparability and nonlocality between these two modes. 
For m_u≤ 0.2 instead, the figure shows that entanglement is concentrated between the partner and the companion (modes 2 and 1, respectively). 

A plot similar to the one of Fig. <ref>, but where the
quantities are evaluated at finite temperature, enables to evaluate
the resilience of entanglement and non-separability to an increase of
temperature.  This check is performed in Fig. <ref> which
shows that, whereas at T=0.2  gn_u the PPT measure is not
dramatically affected, the CHSH parameters B^(0|2) and B^(1|2)
no longer show evidences of violation of Bell inequality, except in
the (1|2) sector for waterfall configurations with
m_u≲ 0.15 and, in a lesser extent, in the (0|2) sector for
m_u≳ 0.85.

Figs. <ref> and <ref> illustrate a specific feature of analog models: there is no recipe enabling to qualify one of the positive norm modes (in our specific case, mode 0 or 1) as 
unessential. For instance, it is incorrect to study the system discarding a priori the companion (mode 1) from one's analysis. Fig. <ref> shows that, at zero temperature, this is allowed in some regions of parameters, but incorrect in others. Fig. 
<ref> even shows that only the companion-partner correlations display 
violation of Bell inequality above a certain temperature. In this instance it is the Hawking mode which is unessential. 
It is therefore important to give a proper account of all the modes involved in our analog system.

Another way to consider the same problem is to use the equivalent model depicted in Fig. <ref>: it is well known that a 
nondegenerate optical parametric amplifier generates an EPR state when the squeezing parameter tends to infinity (see, e.g., Ref. <cit.>). As discussed in Appendix <ref> this is the case for the two mode squeezed state (<ref>) when ω→ 0.
The zero energy transmission coefficient Γ_0 of the effective beam splitter 
(pictorially defined in  Fig. <ref>) tends to 0 or 1 when the upper Mach number m_u tends to 0 or 1, respectively 
[see e.g., Eq. (<ref>) which holds for the waterfall 
configuration]. In these 
two limits the system forms an EPR state between two of the three outgoing modes (and the third one can be omitted). This EPR state involves, either the Hawking quantum 
and the partner (when m_u → 1) or the companion and the partner (when m_u → 0).
This is the reason why the entanglement and nonlocality bounds are reached in this two limits in Fig. <ref>, and similarly in Figs. <ref> and <ref> of Appendix <ref> for the delta peak and flat profile configurations, respectively.

We conclude this section by noticing that, from a quantum information perspective,  
the fact that reduced bipartite states of a tripartite system are entangled 
is of no particular significance per se. It is however important for future experimental studies of analog systems to determine for which configurations, and to quantify to which extent, the three-mode
acoustic Hawking emission is bipartite
entangled and nonlocal. Besides, we will see in the next section that this resilience of entanglement to partial tracing acquires a particular significance when examining the exact nature of the long wavelength components of the three-mode state |0⟩^ in (the ground state of the system). 




§ TRIPARTITE NONLOCALITY
 

Equipped with the same pseudo-spin operators as the ones defined in Sec. <ref> one can define a three-mode Bell operator of a type similar to the two-mode one (<ref>). This operator measures the correlations between the outgoing quasi-particles of type i, j and k. It is defined as  <cit.>:

    𝒮̂^(i|j|k)(ω)  =  
     12 ℬ̂^(i|j)⊗c'·Π̂^(k)
    +
     12 ℬ̂'^(i|j)⊗c·Π̂^(k),

where ℬ̂'^(j|k)(ω) is the same as ℬ̂^(j|k)(ω) defined in (<ref>) with the primes reversed, and c and c'  are normalized vectors
(as well as  a, a', b, and b' involved in the definition of  ℬ̂^(j|k)).
Expanding expression (<ref>) shows that 𝒮̂^(i|j|k) is invariant upon a permutation of its indices, provided the names of the unit vectors (a,b,c) and (a',b',c')
undergo the same permutation. In the following we arbitrarily chose the order (i,j,k)=(0,1,2).


 Similarly to what occurs for the two-mode operator, the principle of local realism, if correct, should predict that
a measure of the operator 𝒮̂^(0|1|2) yields results ± 2.  The system violates this principle
when the average of the operator (<ref>) is larger than two. This is often referred to as violation of 
Svetlichny inequality. It is important to note that the observable 𝒮̂^(0|1|2) is
specially designed in order to be sensitive to genuine tripartite nonlocality (see discussions in Refs. <cit.> and references therein): a tripartite system can involve nonlocal correlations between any of its bipartitions and still not violate the Svetlichny inequality ⟨𝒮̂^(0|1|2)⟩ <2. A simple example of a system which displays two mode nonlocality but does not pass the Svetlichny test is presented in Appendix <ref> [see Eq. (<ref>) and the discussion below].

For attempting to reach a maximum violation of Svetlichny (i.e., three-mode Bell) inequality it is necessary to choose a particular arrangement of the vectors a, 
a', b, b', c and  c' which maximizes the expectation value
 ⟨𝒮̂^(0|1|2)⟩. 
Do do so, we resort to a genetic algorithm which is presented in Appendix <ref> and numerically determines the quantity

    S^(0|1|2)(ω) ≡ _a,a',b, b', c, c'⟨𝒮̂^(0|1|2) (ω) ⟩.

It is shown in Appendix <ref> 
that the tripartite parameter
S^(0|1|2)(ω) is bounded from above by 2√(2) 
[this is Eq. (<ref>)]
and that this bound is reached at ω=0 and T=0 [see Eq. (<ref>)].
The behavior of S^(0|1|2)(ω) at zero temperature is displayed in Fig. <ref> for different realizations of the waterfall configuration. Although the bound 2√(2) is reached in the long wave length limit for all the configurations, even a weak temperature is able to destroy the 
signal of nonlocality as illustrated in Fig. <ref>.


This is connected to the loss of purity of the finite temperature system and
can be understood analytically, again in the long wave length limit, as discussed in Appendix <ref>. This sensitivity to a small finite temperature probably precludes the experimental observation of tripartite nonlocality by means of the observable (<ref>). However, although the zero temperature behavior is certainly difficult to observe, it is rich of fundamental insight on the nature of the state of the system, as we now discuss.

As argued in Appendix <ref>, the fact that, at zero temperature and ω=0 the system reaches the tripartite upper bound (S^(0|1|2)=2√(2), see Fig. <ref>) mathematically stems from the fact that, under these conditions the system is in a pure state and displays perfect correlations in the following expectation values [see Eqs. (<ref>), (<ref>) and (<ref>)]: 

    ⟨Π̂_y^(0)⊗Π̂_y^(1)⊗Π̂_z^(2)⟩    = 
    ⟨Π̂_y^(0)⊗Π̂_z^(1)⊗Π̂_y^(2)⟩
        = - ⟨Π̂_z^(0)⊗Π̂_y^(1)⊗Π̂_y^(2)⟩
        =
    ⟨Π̂_z^(0)⊗Π̂_z^(1)⊗Π̂_z^(2)⟩ = 1.

All the other expectation values of products of three components of the Π̂ operator are zero. Since all the operators in (<ref>) have only  ± 1 as eigenvalues, this means that
each term  must reach its extremal value
and the vacuum mode of zero frequency 
|0_ω=0⟩^ in must therefore be an eigenstate of the operators
appearing in (<ref>): 


    Π̂_y^(0)⊗Π̂_y^(1)⊗Π̂_z^(2)
    |0_ω=0⟩^ in
    = + |0_ω=0⟩^ in , 
       Π̂_y^(0)⊗Π̂_z^(1)⊗Π̂_y^(2) 
    |0_ω=0⟩^ in
    = + |0_ω=0⟩^ in , 
       Π̂_z^(0)⊗Π̂_y^(1)⊗Π̂_y^(2)
    |0_ω=0⟩^ in
    = - |0_ω=0⟩^ in , 
       Π̂_z^(0)⊗Π̂_z^(1)⊗Π̂_z^(2)
    |0_ω=0⟩^ in
    = + |0_ω=0⟩^ in .


Relations (<ref>)
define a state which exhibits the GHZ paradox, contradicting local hidden variable theories by means of a single measurement; cf. e.g., the discussion in <cit.>. In the same line, the above noted fact that our system saturates the upper bound  2√(2) of the Svetlichny parameter is a characteristic also shared by GHZ states. 
Another resemblance lies in the fact that the GHZ states 
possess maximal tripartite
entanglement, quantified by the residual tangle <cit.> while similarly, in the domain of temperature and energy where Eqs. (<ref>) hold, the residual contangle <cit.> which measures genuine tripartite entanglement in our continuous variable system diverges <cit.>.
Despite these clear similarities, there is a significant difference between the state
|0_ω=0⟩^ in and the GHZ states commonly considered in quantum information theory: they have different behaviors upon partial tracing.
As well known, taking the partial trace over one of the three modes of a GHZ state yields an unentangled mixed state.
At variance, partial tracing
the state |0_ω=0⟩^ in over modes 0 or 1 leads to an entangled two-mode mixed state, 
as discussed in Sec. <ref> (see, e.g., the values Λ^(0|2)(ω=0)
and Λ^(1|2)(ω=0) in Fig. <ref>). As we now argue, the explanation of this conundrum lies in the fact that, whereas the 
GHZ states are usually build on qubits, the state we consider is an 
infinite sum of degenerate GHZ states of a continuous variable system. 

To tackle this issue it is convenient to expand the state  
|0_ω=0⟩^ in over the eigenstates of the operators Π̂_x^(j) (j=0, 1 or 2). As discussed in Appendix <ref> these eigenstates 
can be 
written as |x^±_n⟩_j: they are
labeled by their eigenvalue (± 1) plus another integer index (n in the above expression) associated to the infinite degeneracy of both eigenvalues. 
Indeed, at variance with what occurs for a regular spin operator, the projection of the pseudo-spin (<ref>)
over a given axis (here Π̂_x) has infinitely degenerate eigenvalues, i.e.,  
there exists an infinite number of mutually orthogonal eigenstates with the same eigenvalue (+1 or -1):

    ∀ n  ∈ℕ,   Π̂_x^(j) |x^±_n⟩_j
        = ± |x^±_n⟩_j,
  
whereas

    ∀ (n,m)  ∈ℕ^2,     _j⟨ x^±_n|x^±_m⟩_j
    =δ_n,m.

This can be shown by directly constructing the eigenstates of Π̂_x from the number states, see Eqs. (<ref>) and (<ref>).
The expansion of |0_ω=0⟩^ in  over the complete basis formed by these states is of the form:

    |0_ω=0⟩^ in
    =∑_σ_0,σ_1,σ_2
    l,m,n
    C^σ_0,σ_1,σ_2_l,m,n
    |x^σ_0_l,x^σ_1_m,x^σ_2_n⟩.

In the summation appearing in the above expression, σ_0, σ_1 and σ_2=± whereas (l,m,n)∈ℕ^3 and we dropped the indices j=0, 1 or 2 of the kets for legibility. 

It follows from relations
(<ref>) and
(<ref>) that:

    C_l,m,n^-σ_0,-σ_1,-σ2
        =-σ_0σ_1      C_l,m,n^-σ_0,-σ_1,-σ_2
    
        =-σ_0σ_2      C_l,m,n^-σ_0,-σ_1,-σ_2
    
        =σ_1σ_2      C_l,m,n^-σ_0,-σ_1,-σ_2.

This imposes that the only nonzero coefficients in expansion (<ref>) are those for which σ_1=σ_2=-σ_0. For such coefficients Eqs. (<ref>) and (<ref>) impose that

    C^+–_l,m,n=C^-++_l,m,n≡ C_l,m,n.

Expansion (<ref>) thus simplifies to

    |0_ω=0⟩^ in =∑_l,m,n
    C_l,m,n(     
    |x^+_l,x^-_m,x^-_n⟩ 
    +
       
    |x^-_l,x^+_m,x^+_n⟩),

which shows that the vacuum |0_ω=0⟩^ in
of the b̂_j(ω=0) operators (the ingoing ground state) is an 
infinite sum of degenerate GHZ states. It is this property which enables the reduced state obtained after partial tracing over one 
mode to remain entangled despite the GHZ nature of the system. Indeed, tracing over mode 
0 for instance leads to a reduced density matrix

    _(0)( 
    |0_ω=0⟩^ in  ^ in⟨ 0_ω=0|
        )    =
    ∑_m,n,μ,ν C^(0)_m,n,μ,ν   ( 
    |x^-_m x^-_n⟩ ⟨ x^-_μ x^-_ν | 
         +
    |x^+_m x^+_n⟩ ⟨ x^+_μ x^+_ν | ).

In this expression the kets and bras concern modes 1 and 2 (mode 0 has been traced out) and

    C^(0)_m,n,μ,ν=∑_l=0^∞ C^*_l,μ,ν C_l,m,n.

If the eigenvalues +1 and -1 of operator Π̂_x were non-degenerate, this would impose m=μ and n=ν in expression (<ref>) and the corresponding reduced state would be clearly separable. Nothing similar occurs in our situation, and indeed the reduced states are typically entangled, as shown  in Sec. <ref>.

Since we now understand the exact GHZ nature of the zero temperature and ω=0 state of the system, it is of interest to quantify to what extent this feature persists at finite 
temperature and finite energy. To this aim, we use the genetic algorithm presented in Appendix <ref> to
compute the optimum of the Mermin parameter <cit.>

    M^(0|1|2)(ω) ≡ _a,a',b, b', c, c'|⟨ℳ̂^(0|1|2) (ω) ⟩|,

where

    ℳ̂^(0|1|2)(ω)=     
    -a·Π̂^(0)⊗b·Π̂^(1)⊗c·Π̂^(2)
       +a·Π̂^(0)⊗b'·Π̂^(1)⊗c'·Π̂^(2)
       +a'·Π̂^(0)⊗b·Π̂^(1)⊗c'·Π̂^(2)
       +a'·Π̂^(0)⊗b'·Π̂^(1)⊗c·Π̂^(2).

We note here that 𝒮̂^(0|1|2)=12 ℳ̂^(0|1|2) + 12 ℳ̂'^(0|1|2), where ℳ̂'^(0|1|2) is the same as
ℳ̂^(0|1|2) with the prime reversed <cit.>.

The largest possible value of the Mermin parameter (<ref>) is 4. This upper bound is reached for a state verifying the relations (<ref>), such as a GHZ state or the low wavelength 
component of the ground state of our system as just 
discussed[From (<ref>) it is clear that at zero temperature and ω=0 the maximum (<ref>) is obtained when a= b'= c'= e_z and a'= b= c= e_y
and reaches the upper bound 4.].  Indeed, at T=0, M^(0|1|2)(0)=4 for all black hole configurations. This is illustrated in Fig. <ref>

in which M^(0|1|2) is plotted as a function of ω for different values of m_u, at zero and finite temperature (T=0 and T=0.1  gn_u). 
If the departure of M^(0|1|2) from its upper bound is taken as an indication of how much the system differs from a GHZ state, this suggests that the GHZ character of the state is restricted to the low energy and low temperature sector.
We stress however that such a criterion is only indicative: The  Svetlichny and Mermin parameters provide useful bounds, but not proper measures. This is illustrated by their poor effectiveness for evaluating genuine tripartite entanglement: 
It is known <cit.> that if


    M^(0|1|2)    >2√(2),
      
    S^(0|1|2)    >2,


the system exhibits genuine three-mode entanglement. Figs. <ref> and <ref> show that this type of entanglement is certainly reached at T=0 for low energy. However, the criteria (<ref>) are here too restrictive, since the computation of the residual contangle done in Ref. <cit.> demonstrates that at T=0 the system is genuinely tripartite entangled for all energies and all configurations. This general conclusion certainly cannot be reached by 
inspecting in Figs. <ref> and <ref> in which domain of energy and for which value of m_u the criteria (<ref>) are met.
However, these criteria are valid and can be useful, for instance at finite temperature. In this case the state of the system is not pure, the evaluation of the residual contangle appears to be very difficult, and the criterion (<ref>) is the only way we know how to demonstrate genuine tripartite entanglement, which, as can be inferred from the tendency displayed in Fig. <ref>, is reached at low T, low ω and small values of m_u.

Finally, the Mermin parameter (<ref>) is also an interesting witness of nonlocality.
A local hidden variable theory predicts that it  should
verify the Mermin-Klyshko inequality M^(0|1|2)≤ 2.
It can be seen from Fig. <ref> that this inequality is violated at zero temperature for all waterfall configurations. For T>0 instead,
M^(0|1|2)(0)=0.
The Mermin parameter has the same behaviour as the Svetlichny parameter, for the same reason: At finite temperature the state of the system is no longer pure, and in this case all the expectation values of products of three components of the pseudo-spin tend to zero in the long wavelength limit, see discussion in Appendix <ref>.
However, contrarily to what occurs for the Svetlichny parameter, there are accessible finite temperature situations where the Mermin parameter is larger that the nonlocality threshold M^(0|1|2)=2 (compare Figs. <ref> and <ref>). In that respect, the Mermin parameter may even reveal more useful that the  
CHSH parameter (<ref>). For instance, at T=0.1 g n_u, for m_u=0.3 the largest value of M^(0|1|2) is 2.19 (as seen from Fig. <ref>), higher than the largest values reached by both CHSH parameters B^(0|2) and B^(0|1) in the same situation (2 and 2.017, respectively).
 
 


§ CONCLUSION


In this work we conducted a systematic theoretical  study of violation of bipartite and tripartite Bell inequalities in an analogue of black hole realized in the flow 
of a quasi-one dimensional Bose-Einstein condensate. There is a reasonable hope to witness bipartite entanglement in such systems, and also possibly bipartite nonlocality whereas the observation of 
signatures of genuine tripartite nonlocality would presumably be more difficult. The reason is that violation of Svetlichny inequality does not resist much to an 
increase in temperature. In that respect it is worth underlying that, from the three black hole configurations we have considered, the waterfall configuration is the one for which the signature of nonlocality is the less sensitive to an increased temperature: compare for instance the finite temperature signal of Fig. <ref> with the one of Figs. <ref> and <ref>.

An important ingredient of our theoretical treatment is the assumption of Gaussianity which, in our case, results from the Bogoliubov approach.  This is certainly a sensible first order hypothesis, but recent studies <cit.> showed that nonlinear effects might significantly affect quantum emission processes in the context of analogue physics. It would then be of great interest to evaluate in detail the amplitude and the relevance of nonlinear back reaction effects on acoustic Hawking radiation in the black hole analogues we consider.

Within the Gaussianity assumption, the main theoretical tool of our study is the covariance matrix. As indicated in Sec. <ref>, indirect techniques have been used for measuring some of the entries of this matrix in BEC physics. Other systems, such as exciton-polaritons in microcavities, or more generally setups involving nonlinear light, could enable other types of measurements, with possibly more direct access to the phases of the averaged quadratures. 
Our approach can be adapted with little modification to such systems. For the quantities we are interested in, a key element is the long wavelength behavior of the dispersion relation (and of the corresponding S matrix). For instance the GHZ nature of the long wavelength modes is generic in any system 
with a sonic-like low energy dispersion relation. The observables studied in the text might however be modified in a system with a different long wave-length behavior, for instance in the presence of massive modes, as is the case in the polariton context, but also in coherently coupled two-component BECs.

A fundamental outcome of our work 
stems from a flaw in the analogy. Our analogue system breaks Lorentz invariance, but
paradoxically this apparent shortcoming (characteristic of 
analogue physics <cit.>) turns into an advantage: It opens the possibility of tripartite entanglement and nonlocality, thanks to
the existence not only of a Hawking and a partner mode, but also of a third mode we denote as the companion
This peculiar nonlocal tripartite configuration,
together with
the continuous nature of the degrees of freedom, 
induces a surprising property of the system: the
long wavelength quantum modes consist in a superposition
of  degenerate GHZ
states which,  at variance with GHZ states build on
qubits, remains entangled after partial tracing.
We expect this feature to be generic, thus suggesting that condensed matter analogues may indeed open new
prospects of robust information protocols and provide efficient platforms to study
the flow of multipartite quantum information.



  We thank M. J. Jacquet, J. Martin and M. Walschaers for inspiring
  discussions.  We also gladly acknowledge C. E. Lopetegui for
  relevant remarks concerning continuous-variable systems. We are
  particularly grateful to A. Buchleitner for thorough exchanges on
  many of the aspects addressed in this work.  This project has
  received support by the DFG funded Research Training Group “Dynamics
  of Controlled Atomic and Molecular Systems” (RTG 2717).







§ DIFFERENT BLACK HOLE CONFIGURATIONS


The background classical field Φ of Eq. (<ref>) is solution of a Gross-Pitaekskii equation which is the stationary and classical version of Eq. (<ref>): 

    -ħ^2/2m∂_x^2Φ
        +[U(x) + g   n-μ]Φ =0 ,

where n=|Φ|^2. An analogue black hole horizon is realized if the flow is subsonic far upstream (for x→ -∞ in our convention) and takes the form of a downstream supersonic plane wave (for x>0). The corresponding classical field behaves as

    Φ(x)={[    √(n_u)exp( i k_u x) ϕ_u(x)                                                       x≤ 0,; √(n_d)exp( i k_d x)exp( iβ_d)                                                       x≥ 0. ].

In this expression lim_x→-∞ϕ_u(x)=exp( iβ_u), where β_u is a constant, as is β_d. In all the text the indices u and d refer to the upstream and downstream flow. For instance
n_u and n_d are the upstream and downstream asymptotic densities. Also
k_α = m V_α/ħ (α=u or d), where V_u (V_d) is
the asymptotic upstream (downstream) flow velocity. The asymptotic velocities of sound
c_u and c_d are defined by m c_α^2=g_α n_α, where[We consider a possible
position-dependent nonlinear coefficient for being able to treat the flat profile
configuration of Ref. <cit.>.]
g_u,d=lim_x→-∞,+∞g(x). The
healing lengths and Mach numbers are defined as ξ_α=ħ/(m c_α) and
m_α=V_α/c_α, respectively.
The form of the flow pattern is specified by the value of the above parameters and 
by the function ϕ_u(x). In the remaining of this subsection we briefly describe the 3 configurations studied in the present work and refer to <cit.> for a detailed presentation 
of these configurations
with the precise values of the corresponding parameters.

The first configuration we consider is an idealized one, introduced in
Ref. <cit.> which we denote as “flat profile”. It
consists in a constant uniform plane wave flow with ϕ_u(x)=1,
β_d=0, n_u=n_d≡ n_0 and k_u=k_d≡ k_0. Such a
configuration can realize an analog black hole in the presence of a
piece-wise nonlinear parameter g(x)=g_uΘ(-x)+g_dΘ(x),
supplemented by an also piece-wise potential U(x), as explained in
Refs. <cit.>. In the two other configurations
we consider g is a constant.

The second configuration is denoted as “waterfall”: The potential
U(x) is a step function and the upstream profile is half a dark
soliton. It is close to the experimental realization of the Technion
group <cit.> and catches important aspects
of the density-correlation pattern reported in <cit.>,
with nonetheless some caveats, see <cit.>.

The third and last configuration is denoted as “delta peak”: the
potential U(x) is a repulsive Dirac distribution located at x=0
and the upstream profile is a part of a dark soliton, but not exactly
one half of it, as is the case in the waterfall configuration.

Both configurations are depicted in Fig. <ref>. In this figure the region x>0 is shaded in order to remind that it corresponds to the interior of the analog black hole. It is however important to recall that the precise location of the horizon separating the interior and the exterior of an analog black hole is ill-defined, see, e.g., the discussion in Sec. II.A of Ref. <cit.>.

Note that, while the upstream and downstream Mach numbers can be fixed independently in the flat profile configuration (with the only constrain that m_u<1<m_d), these two quantities are not independent in the waterfall and delta peak configurations: 


    m_d=m_u^-2     , 
    m_d/m_u= ( 
     -1+√(1+8 m_u^-2)/ 2)^3/2





§ COVARIANCE MATRICES


The 6× 6 covariante matrix σ defined in Eq. (<ref>) can be
written in terms of 2× 2 submatrices σ_i and
ε_ij:

    σ=[     σ_0    ε_01    ε_02; ε_01^ T     σ_1    ε_12; ε_02^ T ε_12^ T     σ_2 ],

where

    σ_i =
    [            ⟨  2q̂_i^2 ⟩ ⟨q̂_i p̂_i + p̂_i q̂_i⟩; ⟨q̂_i p̂_i + p̂_i q̂_i⟩             ⟨ 2p̂_i^2 ⟩;                         ]
    ,

and

    ε_ij =
    2  [ ⟨q̂_i q̂_j⟩ ⟨q̂_i p̂_j⟩; ⟨p̂_i q̂_j⟩ ⟨p̂_i p̂_j⟩;             ].

In the case we consider, the 
definitions (<ref>) and the
specific form of the transformation (<ref>)
result in the following expressions:

    σ_i =(1 + 2 ⟨ĉ_i^†ĉ_i ⟩) 1_2,


    ε_01 =2
    [  Re⟨ĉ_0  ĉ_1^†⟩ -Im⟨ĉ_0  ĉ_1^†⟩;  Im⟨ĉ_0  ĉ_1^†⟩  Re⟨ĉ_0  ĉ_1^†⟩ ],

and (for i=0 or 1)

    ε_i2  = 2
    [  Re⟨ĉ_i  ĉ_2 ⟩  Im⟨ĉ_i  ĉ_2 ⟩;  Im⟨ĉ_i  ĉ_2 ⟩ -Re⟨ĉ_i  ĉ_2 ⟩ ].

In the limit of zero temperature the system is in a three mode pure
Gaussian state. Its covariance matrix (<ref>) can accordingly be
brought by LLUBOs (local linear unitary Bogoliubov transformations) to
a “standard form” in which matrices σ_i are proportional to
the identity [as they already are, cf. (<ref>)] and matrices ε_ij are
diagonal <cit.>.  After this operation
the matrices ε_ij take the following form <cit.>

    ε_01
    = 2  |⟨ĉ_0  ĉ_1^†⟩|  1_2  , 
      ε_i2 = 2   |⟨ĉ_i  ĉ_2 ⟩|  σ_z,

where i∈{0,1} and σ_z is the third Pauli matrix.

The situation at finite temperature is less simple. The system is in a
mixed state with no special symmetry, and the covariance matrix
(<ref>) cannot be put in a standard form where the matrices
ε_ij are all diagonal <cit.>. However, the situation simplifies again if
one is interested in bipartite entanglement only (say, between modes i
and j). In this case one should trace out the third mode (let's denote it by k)
which simply amounts to remove from the total covariance matrix (<ref>) the
two rows and two columns where index k appears. The remaining
4 × 4 covariance matrix associated with the reduced two-mode
state reads

    σ^(i|j)=
    [     σ_i    ε_ij; ε_ij^ T     σ_j ],

where the 2× 2 blocks are the same as the ones in (<ref>).
This reduced covariance matrix
can always be brought by LLUBOs to its standard form
<cit.>, in which the matrix ε_ij takes again
the form (<ref>) with here an average taken over the finite
temperature state, as explained in the main text, Sec. <ref>.

It has been argued in <cit.> that an efficient measure of
bipartite entanglement was given by the “PPT measure” Λ^(i|j)≡ 1-ν^-_(i|j)
where ν^-_(i|j)   is the lowest simplectic eigenvalue of the partial
transpose of σ^(i|j) [The partial transposition
corresponds to a mirror reflection in phase space which inverts the p_j coordinate, leaving q_i, p_i and q_j unchanged <cit.>. The resulting covariance matrix can be brought by means of a symplectic transform to a diagonal form <cit.>. The corresponding diagonal elements are the simplectic eigenvalues. They are twice degenerate, and in our 4× 4 case there are thus two such eigenvalues: ν^-_(i|j) and ν^+_(i|j), with ν^±_(i|j)∈ℝ^+ and ν^-_(i|j)≤ν^+_(ij).]. The largest entanglement corresponds to Λ^(i|j)=1 while separability is reached 
when Λ^(i|j)<0. This separability condition can be shown to be equivalent to the Peres-Horodecki criterion 
<cit.>. Contrarily to other observables, 
such as the Cauchy-Schwarz criterion  or the generalized Peres-Horodecki parameter which have been often used in the domain, the PPT measure has the advantage of being an entanglement monotone. 
Other observables have been used in the context of analog gravity, which are monotonous measures of entanglement, such as the 
entanglement entropy, the entanglement of formation or the logarithmic negativity, but they all have some drawbacks: the state 
our our system is mixed in a finite temperature situation, which discards the entanglement entropy as a possible measure. The 
entanglement entropy generalizes for mixed states to the entanglement of formation <cit.>, but this quantity is 
not easily determined in non symmetric two-mode Gaussian states such as the ones we consider[In the situations we consider the reduced state of modes i and j is nonsymmetric,
since in general the mixedness a_i and a_j (defined in (<ref>)) are not equal.]. 
The logarithmic negativity shares with the  entanglement of formation the drawback of possibly violating monogamy inequalities <cit.>. This is a relatively mild drawback in the context of evaluating bipartite entanglement, but becomes prohibitive in the tripartite context. The Gaussian contangle was introduced in Ref. <cit.> as a quantity which has none of the previous deficiencies. It has been studied in the gravitational context <cit.> and also
in analog gravity <cit.>, but it is not of very practical use as its determination requires numerical minimisation of a complex expression. The PPT measure we use in the present study was introduced in Ref. <cit.> as a measure which mimics many aspects of the Gaussian contangle but is much simpler to evaluate. Its evaluation is 
as simple as the one of the above quoted quantities albeit it 
shares none of their drawbacks. The explicit expression of Λ^(i|2)(ω) characterizing the coupling between mode i=0 or 1 and mode j=2 is given in Eq. (<ref>) of the main text.

As clear from the above expressions (<ref>), (<ref>) and (<ref>), the theoretical evaluation of the components of the covariance matrix relies on the computation of averages of two creation or annihilation operators of the outgoing modes. The relevant expressions are determined from Eq. (<ref>) and read (i=0 or 1):

    ⟨ĉ^†_i ĉ_i⟩ =    
        |S_i0|^2 n̅_0 + |S_i1|^2 n̅_1 + |S_i2|^2 (1+n̅_2),         
    ⟨ĉ_2^†ĉ_2⟩=    |S_20|^2 n̅_0 
    + |S_21|^2 n̅_1 + |S_22|^2 (1+n̅_2)-1,
    ⟨ĉ_0 ĉ_1^†⟩=    S_00S_10^* n̅_0 + S_01S_11^* n̅_1 + S_02S_12^* (1+n̅_2),
    ⟨ĉ_i ĉ_2⟩ =    S_i0S_20^* n̅_0 
    + S_i1S_21^* n̅_1 + S_i2S_22^* (1+n̅_2),

where the n̅_j's (j=0, 1 or 2) are the occupation numbers of the incoming modes, see Eq. (<ref>). 
The coefficients of the S matrix appearing in the above formulae can be determined numerically as explained in Ref. <cit.>. 
Expressions (<ref>) are useful for computing the PPT measure of entanglement (<ref>) and the CHSH
(<ref>) and Svetlichny (<ref>) parameters.
As an illustration we now indicate how to compute the PPT measure Λ^(i|2) at zero temperature. Using here the definition (<ref>) as a short hand notation one gets from (<ref>):

    Λ^(i|2)(ω) =    √([|S_i2|^2+ sinh^2(r_2)]^2 + 4 |S_i2|^2)
        -|S_i2|^2-sinh^2(r_2).

It was shown in Ref. <cit.> that the ratio |S_i2|^2/sinh^2(r_2) tends to a constant when ω→ 0. Let's denote Γ_i the value of this constant (i=0 or 1) [Γ_0 and Γ_1 are the low energy limits of the transmission and reflection coefficients (γ_0 and γ_1, respectively) of the beam splitter involved in the effective optical model depicted in Fig. <ref>, see Appendix <ref>. In the notations of Appendix <ref>, Γ_0=lim_ω→ 0cos^2θ.]. A simple expansion of (<ref>) shows that

    lim_ω→ 0Λ^(i|2) =
        2 Γ_i/1+Γ_i.

In the waterfall configuration <cit.>

    Γ_0=4 m_u/(1+m_u)^2,

and from relation (<ref>) it follows that Γ_0+Γ_1=1. Hence, for the waterfall configuration


    lim_ω→ 0Λ^(0|2) =
        8 m_u/1+6m_u+m_u^2, 
       lim_ω→ 0Λ^(1|2) =
        (1- m_u)^2/1+m_u^2.


The maximum value of Λ^(0|2) is always reached at ω=0, thus the T=0
numerically determined value max_ωΛ^(0|2) plotted in Fig. 
<ref> is identical to (<ref>). The maximum value of the PPT 
measure of entanglement between modes 1 and 2 is reached at ω=0 only for m_u≲ 0.18: the 
numerically determined value max_ωΛ^(1|2) plotted in Fig. 
<ref> is thus identical to (<ref>) in this range of values of m_u.



§ AN ANALOGUE OPTICAL SYSTEM


As discussed in Ref. <cit.>, the entanglement in the system can be localized by a transformation involving effective modes f̂_0,f̂_1,f̂_2 schematically represented in Fig. <ref>. These modes are related to the physical outgoing modes by:

 

    f̂_0     = -sinθ ê_0 + cosθ ê_1, 
    f̂_1     = cosθ ê_0 + sinθ ê_1, 
    f̂_2   =ê_2,


 where

    ê_0 = S_02^*/|S_02| ĉ_0,   ê_1 = S_12^*/|S_12| ĉ_1,   ê_2 = S_22/|S_22| ĉ_2,

with

    cosθ=|S_02|/sinh r_2,   sinθ=|S_12|/sinh r_2,

and

    r_2(ω)= arsinh√(|S_22|^2-1).

It has been shown in Ref. <cit.> that an incoming vacuum mode of frequency ω relates to that of the effective f modes by

    |0_ω⟩^ in = exp[r_2 ( f̂_1^† f̂_2^† - f̂_1  f̂_2 )]  |0_ω⟩^f,

indicating that the T=0 state of the system is, as far as the f_1 and f_2 modes are concerned, a two-mode squeezed vacuum with squeezing parameter r_2(ω). The f_1 and f_2 modes are mixed by
a beamsplitter (see Fig. <ref>) of transmission and reflection coefficients cos^2θ and sin^2θ, respectively. It has been argued in Ref. <cit.> that 
the long wavelength limit of the transmission coefficient,
Γ_0=lim_ω→ 0cos^2θ, plays the role of the grey-body factor of the analog black hole.


From definitions (<ref>) and (<ref>) and expressions (<ref>) it is a simple matter
to evaluate the following averages

    ⟨f̂_0^†f̂_0 ⟩ = ∑_i=0^1
    |S_12S_0i-S_02S_1i|^2  n̅_i/sinh^2 r_2,
    
       ⟨f̂_1^†f̂_1 ⟩ =
    cosh^2 r_2  n̅_01 + sinh^2 r_2 (1+n̅_2), 
       ⟨f̂_2^†f̂_2 ⟩ =
    sinh^2 r_2  n̅_01 + cosh^2 r_2 (1+n̅_2) -1,
    
       ⟨f̂_1 f̂_2 ⟩ = 
    cosh r_2 sinh r_2 (n̅_01+n̅_2 +1).

The occupation numbers n̅_j (j=0, 1 or 2)  in the above expressions are defined in (<ref>) and use has been made of the shorthand notation

    n̅_01≡|S_20|^2/sinh^2 r_2 n̅_0 +
    |S_21|^2/sinh^2 r_2 n̅_1.

It is quite informative to quantify the entanglement between the two
effective modes f_1 and f_2 and the amount by which the
corresponding squeezed state violates Bell inequality. At T=0 all
the relevant quantities can be expressed in terms of the squeezing
parameter (<ref>) involved in the transformation
(<ref>). For instance, the analogues for the the modes
f̂_1 and f̂_2, of the PPT measure of entanglement and of
the CHSH parameter defined for the modes ĉ_i and ĉ_2
by Eqs. (<ref>) and (<ref>), respectively, read

    Λ^(f_1|f_2)=1-exp(-2 r_2),

and <cit.>    

    B^(f_1|f_2)=
        2√(1+4/π^2arctan^2[sinh(2 r_2)]).


The finite temperature values of these quantities can be obtained by
replacing the ĉ operators by the f̂'s in expressions
(<ref>) and (<ref>) and using formulae
(<ref>).  They are represented in Fig. <ref> as
functions of ω for a specific black hole configuration
(waterfall with m_d=2.9). As anticipated, it appears that the
violation of Bell inequality is here much more resilient to
temperature than for the true Hawking-partner pair or the
companion-partner pair (compare with Fig. <ref>). This is
explained by the fact that the beam splitter distributes and, so to
say, dilutes the entanglement. This is clear at T=0: in this case
the f_0 mode is empty [cf. (<ref>)] and the system is, as
far as the f modes are concerned, in a pure two-mode squeezed vacuum
state. At ω=0 for instance, the system is maximally entangled
(the squeezing parameter r_2→∞) and both
Λ^(f_1|f_2) and B^(f_1|f_2) reach their upper bounds (1
and 2√(2) respectively).  At variance, if working with the true
outgoing modes described by the ĉ (or equivalently the
ê) operators, for studying two-mode entanglement it is
necessary to trace the (occupied) third one. The resulting density
matrix is mixed, and in this case, entanglement and violation of Bell
inequality are two different things, as clear from Fig. <ref>.


Another way to tackle this issue consists in expressing the PPT measure of entanglement (<ref>) between the partner (mode 2) and mode i (i=0 or 1) in terms
of the parameters of the equivalent optical system. Denoting as 
γ_0(ω)=cos^2θ and 
γ_1(ω)=sin^2θ the transmission and reflection coefficient of the beam splitter makes it possible to write, at zero temperature, the PPT measure  (<ref>) under the form

    Λ^(i|2)(ω)    = (1+γ_i) sinh r_2×
       {√(cosh^2 r_2 -
    (1-γ_i/1+γ_i)^2)
    -sinh r_2 }.

Algebraic manipulations then show that, because γ_i≤ 1, Λ^(i|2) is always lower that Λ^(f_1|f_2). 
This means that the entanglement between modes i and  2 is always lower than the entanglement between the modes issued from the parametric down converted represented in Fig. <ref>.
The equality is reached when γ_i=1, and, since γ_0+γ_1=1, in this case the other channel (let's denote it as =1-i) is not entangled with mode 2: Λ^(|2)=0. So, indeed, the effect of the beam splitter is, so to say, to dispatch the entanglement of the effective squeezed modes f_1 and f_2 between modes c_0, c_1 and c_2.

After discussing two modes entanglement it is also interesting to briefly consider tripartite entanglement in the effective optical system depicted in Fig. <ref>. This system has been designed – by means of a procedure called entanglement localization <cit.> – in such a way that 
it concentrates entanglement between the effective squeezed modes f_1 and f_2.
Therefore, it should be expected that, in the f basis, the quantity ⟨𝒮̂^(f_0|f_1|f_2)⟩, which is sensible to a genuine tripartite nonlocality, will never reach values above two. Indeed, in this basis, at zero temperature one has

    ⟨Π̂_r^(f_0)⊗Π̂_s^(f_1)⊗Π̂_t^(f_2)⟩ = ⟨Π̂_r^(f_0)⟩⟨Π̂_s^(f_1)⊗Π̂_t^(f_2)⟩.

Since

    ⟨Π̂_z^(f_0)⟩ = 1,   ⟨Π̂_x^(f_0)⟩ =
    ⟨Π̂_y^(f_0)⟩ = 0,

it is legitimate to take, for evaluating expression (<ref>) the vector a_+ = ± e_z (taking a_- = ±e_z instead doesn't change the result (<ref>) below). This leads to

    max_θ⟨𝒮̂^(f_0|f_1|f_2)⟩ = 
     |⟨   b·Π̂^(f_1)⊗c' ·Π̂^(f_2)⟩
    
     +    ⟨b' ·Π̂^(f_1)⊗c·Π̂^(f_2)⟩| .

Then, noticing that at zero temperature 
⟨Π̂^(f_1)_z⊗Π̂^(f_2)_z⟩ =1
at all energies and that the vector b, b', c and c' are not constrained with respect to each other, one can take these vectors to be + e_z, thus obtaining

    S^(f_0|f_1|f_2)(ω) T=0 = 2.


Of course the value of S^(f_0|f_1|f_2)
decreases
when the temperature increases. This quantity is thus always lower than 2, indicating,  as expected, that there is 
no genuine tripartite nonlocality between the effective f modes.

From  Eqs. (<ref>), (<ref>) and also from the facts that 
⟨Π̂_s^(f_1)⊗Π̂_t^(f_2)⟩=0
if s≠ t and that, at ω=0 and T=0,

    ⟨Π̂_x^(f_1)⊗Π̂_x^(f_2)⟩=
    -⟨Π̂_y^(f_1)⊗Π̂_y^(f_2)⟩=1,

it is easily found that, at ω=0, the optimized Mermin parameter (<ref>) of the f modes is

    M^(f_0|f_1|f_2)(0) T=0 = 2.

This shows that, at variance with the c modes, the effective f modes do not violate the Mermin-Klyshko inequality and certainly do not exhibit the GHZ paradox. Results (<ref>) and (<ref>) were expected: since the f modes do not exhibit tripartite entanglement they should violate none of the inequalities (<ref>).

It follows from (<ref>) that when m_u=0 or 1, at ω=0, cosθ=0 or 1, respectively. 
In this case Eqs. (<ref>), (<ref>) 
and (<ref>) indicate that
the f modes are connected to the c modes by LLUBOs:
the tripartite character of the true system thus disappears in these two limiting cases, even when T=0. However our results indicate that these 
two limits are singular, since the GHZ character of the system is observed (at ω=0 and T=0) for all m_u∈]0,1[.



§ EIGENSTATES OF THE PSEUDO-SPIN OPERATORS


In this appendix we list the properties of the eigenstates of the pseudo-spin operators (<ref>) which are useful in the main text. We will consider a given one of the three modes (j=0, 1 or 2), always the same, and will omit the associated label (j) in order to lighten the notations. In a similar way, we do not write the ω-dependence which is implicit in all this appendix.

From the definition (<ref>) it follows that the matrix elements of Π̂_z between two number states are

    ⟨ n | Π̂_z |m ⟩ =
    [(-)^n+(-)^m] ∫_0^∞ dq ψ_n(q)ψ_m(q),

where ψ_n(q)=⟨ n|q⟩
=
(2^n n! √(π) )^-1/2exp(-q^2/2) H_n(q) is a normalized Hermite function (H_n a Hermite polynomial). The prefactor in the r.h.s. of 
(<ref>) imposes that n and m have the same parity. Thus ψ_n and ψ_m are both even or odd functions of q, which makes it possible to extend by symmetry the integration range in the r.h.s. of (<ref>) to the whole real axis, leading to

    ⟨ n | Π̂_z | m ⟩=(-)^nδ_n,m,

and thus

    Π̂_z = ∑_n=0^∞(
        |2n⟩⟨ 2n| - |2n+1⟩⟨ 2n+1|
        ):

the eigenstates of Π̂_z are the number states, they have eigenvalue ± 1 depending on their parity. In the following we will denote them as

    |z^+_n⟩=|2n⟩,      |z^-_n⟩=|2n+1⟩.

Contrarily to usual spins 1/2, the two eigenvalues (here +1 and -1) are infinitely degenerate. This property is shared by any projection of the pseudo-spin operator; it is in particular true for the operator Π̂_x. We denote the eigenstates of Π̂_x associated to the eigenvalue ± 1  by
|x^±_n⟩, they
can be constructed 
by rotating the
eigenstates |z^±_n⟩ of Π̂_z by angles ±π/2 around the y axis:

    |x^±_n⟩ = R̂_y(±π/2) 
       |z^±_n⟩ =
       1/√(2)(1∓Π̂_y)
       |z^±_n⟩,

where R̂_y(θ)=exp(-i/2θ Π̂_y) is the operator of rotation of angle θ around the y axis. The fact that
|x^±_n⟩ is an eigenstate of Π̂_x with eigenvalue ± 1 is easily checked by direct application of Π̂_x to the left of expression (<ref>) and use of the relations

    Π̂_r Π̂_s = i  ε_rst Π̂_t,

where ε_rst is the totally antisymmetric Levi-Civita symbol and (r,s,t)=x, y or z. By the same token it is also easily proven that

    Π̂_y |x^±_n⟩ =
        ∓ i |x^∓_n⟩,     Π̂_z |x^±_n⟩ =
        |x^∓_n⟩.






§ COMPUTING EXPECTATION VALUES OF OBSERVABLES


The use of the technique of Wigner transform is particularly well suited for
determining the different contributions of the pseudo-spin operator to the CHSH measures (<ref>) and (<ref>). The reason is twofold: first, we consider a Gaussian state and for computing the required averages we thus just need to evaluate Gaussian integrals weighted by the Wigner transforms of the spin operators and secondly, the Wigner transforms of the pseudo-spins (<ref>) are simple enough that the relevant integrals can be evaluated analytically.

We will consider n=3 (and also n=2) modes Gaussian states and the computation of the average of an operator Â is performed in an abstract dimensional phase space according to 

    ⟨Â⟩ = ∫ d^nq d^np 
        W_ρ̂(q,p)  
        W_Â(q,p) .

In this expression W_ρ̂ is the Wigner transform of the density matrix. In the Gaussian case we consider it can be computed from the knowledge of the covariance matrix <cit.>:

    W_ρ̂(q,p) = 1/π^n √( detσ)exp{-12 ξ^ T σ^-1ξ},

where σ is the total (or reduced, as appropriate) covariance
matrix defined in Eq. (<ref>). In this expression, n=3
and
ξ = √(2)  (q_0,p_0,q_1,p_1,q_2,p_2)^ T in the three-mode case. In the reduced two-mode case,
n=2 and one should remove from the expression of ξ
the entries corresponding to the subscript of the traced mode.  The
Wigner transform (<ref>) of a Gaussian state is non-negative,
and it was originally considered impossible to violate Bell’s
inequality under such conditions <cit.>. This was latter
proven incorrect <cit.>. In particular,
Revzen et al. <cit.> proved that observables can be
associated to the violation of Bell inequality over a Gaussian state
provided their Wigner transform takes values different from the
eigenvalues of the associated quantal operator. We will see that the
pseudo-spins (<ref>) we consider belong to this class of
observables, denoted as “improper” in Ref. <cit.>.

The other term involved in the integration (<ref>) is the Wigner
transform W_Â of the operator Â. It is a
function defined in phase space by means of an integration over
“representation space”

    W_Â(q,p) =
    ∫ d^nzexp{i p·z}⟨q-12 z | Â |
    q+12 z⟩.
  
In this expression q (as well as z and p) is a vector in an abstract n=3 dimensional space (with basis e_0,e_1,e_2). 
The kets involved in (<ref>) are of the type |Q⟩=|
Q_0⟩_0 ⊗ |Q_1⟩_1 ⊗ |Q_2⟩_2 where Q
=Q_0 e_0 + Q_1 e_1 + Q_2 e_2 and
|Q⟩_j is the eigenstate of operator q̂_j associated to the eigenvalue Q (j=0, 1 or 2). In the case of a reduced two modes Gaussian state 
n=2 and the vector Q is two-dimensional: the component associated to the traced mode disappears, as it also does in  |Q⟩.

For evaluating expectation values such as those appearing in
Eqs. (<ref>) and (<ref>) we need to compute the Wigner
transforms of the components of the pseudo-spin operators. These are
to be evaluated in a 2 dimensional phase space, since operator
Π̂^(j) concerns a single mode (mode
j). The result is independent of j and reads:

    W_Π̂_x(q,p)=      (q), 
    
    W_Π̂_y(q,p)=     i δ(q)
    ∫_-∞^∞ dx   (x) exp{-2 i p x}
    
    =   δ(q)  P(1/q), 
    
    W_Π̂_z(q,p)=    πδ(q)δ(p) ,

where P denotes the principal value. It is shown in the main
text that the eigenvalues of the projections along a given axis of the
pseudo-spin operator Π̂ are ± 1. The above
expressions of the Wigner transforms thus demonstrate that, contrarily
to Π̂_x, operators Π̂_y and Π̂_z are
“improper” in the sense of Ref. <cit.>: they may be
involved in violation of Bell inequality even for a state with a
non-negative Wigner transform, such as the Gaussian state we consider.

For evaluating the expectation values appearing in (<ref>) we
need to compute integrals such as

    ⟨Π̂^(i)_z ⊗Π̂^(2)_z⟩
    =∫    d^2q d^2p 
    W_Π̂_z(q_i,p_i)×
        W_Π̂_z(q_2,p_2)
    W_ρ̂(q,p)  ,

where q=(q_i,q_2) and p=(p_i,p_2) (i=0 or 1). 
The explicit calculation yields

    ⟨Π̂^(i)_z ⊗Π̂^(2)_z⟩
    =1/a_i a_2 -  4|⟨ĉ_iĉ_2⟩|^2 ,

where 

    a_j=2 ⟨ĉ^†_jĉ_j⟩ + 1

is known as the local mixedness of mode j (j=0,1 or 2). The other terms involved in the determination of the CHSH parameter (<ref>) can be evaluated by the same technique. The computation is similar to the one figuring in the Appendix of Ref. <cit.>. One finds that 
all expectation values for which the index z appears a single time cancel at all temperatures. The values of the other non-zero averages are  (i=0 or 1):

    ⟨Π̂^(i)_x ⊗Π̂^(2)_x⟩ =    2/πarctan2  ⟨ĉ_iĉ_2⟩/√(a_i a_2 - 
        4( ⟨ĉ_iĉ_2⟩)^2), 
    ⟨Π̂^(i)_y ⊗Π̂^(2)_y⟩ =    -1/A_i2⟨Π̂^(i)_x ⊗Π̂^(2)_x⟩, 
    ⟨Π̂^(i)_x ⊗Π̂^(2)_y⟩ =    2/π a_2arsinh2  ⟨ĉ_iĉ_2⟩/√(A_i2),
    ⟨Π̂^(i)_y ⊗Π̂^(2)_x⟩ =    2/π a_iarsinh2  ⟨ĉ_iĉ_2⟩/√(A_i2),

where A_i2 is defined in Eq. (<ref>) below, and we recall that the expression of the quantity ⟨ĉ_iĉ_2⟩ is given in (<ref>).

For studying the Svetlichny observable it is necessary to evaluate averages involving the Cartesian coordinates of three pseudo-spins, of the type

    T_rst(ω) ≡⟨Π̂^(0)_r ⊗Π̂^(1)_s
        ⊗Π̂^(2)_t⟩,

where r, s and t∈{x,y,z}. These quantities are zero if z is not one of the indices or appears exactly twice. 
In order to write down tractable expressions in the other cases, it is convenient to introduce new compact notations:

    A_01=    a_0a_1- 4 |⟨ĉ_0ĉ_1^†⟩|^2,
    
    A_i2i<2=    a_ia_2- 4 
    |⟨ĉ_iĉ_2⟩|^2,


    Z_0=-2 a_0 ⟨ĉ_1ĉ_2⟩^*
    + 4 ⟨ĉ_0ĉ_1^†⟩⟨ĉ_0ĉ_2⟩^*
    ,
       Z_1=-2 a_1 ⟨ĉ_0ĉ_2⟩^*
    + 4 ⟨ĉ_0ĉ_1^†⟩^* 
    ⟨ĉ_1ĉ_2⟩^*
    ,
       Z_2=-2 a_2 ⟨ĉ_0ĉ_1^†⟩^*
    + 4 ⟨ĉ_1ĉ_2⟩⟨ĉ_0ĉ_2⟩^*
    ,

and

    δ=     a_0 a_1 a_2 +16  {⟨ĉ_0ĉ^†_1⟩⟨ĉ_1 ĉ_2⟩⟨ĉ_0 ĉ_2⟩^*
        }
           -4 a_0 |⟨ĉ_1 ĉ_2⟩|^2
        -4 a_1 |⟨ĉ_0 ĉ_2⟩|^2
        - 4 a_2 |⟨ĉ_0ĉ^†_1⟩|^2.

δ is the square of the determinant of the 6× 6 covariance matrix (<ref>). At T=0 the system is in a pure state and δ=1 for all values of ω <cit.>, whereas δ= O(1/ω) at finite temperature, as can be shown on the basis of the low energy expansion of the matrix elements of the S-matrix given in Ref. <cit.>.

The explicit theoretical evaluation of the quantities (<ref>), (<ref>), (<ref>) and (<ref>) is easily done from Eqs. (<ref>).
A long computation shows that the averages (<ref>) which are non zero can be expressed 
in terms of these quantities:

    T_zxx=-2/π a_0arctan  Z_0/√(A_01A_02-(  Z_0)^2)
    ,
        T_xzx=-2/π a_1arctan  Z_1/√(A_01A_12-(  Z_1)^2)
    ,
        T_xxz=-2/π a_2arctan  Z_2/√(A_02A_12-(  Z_2)^2)
    ,


    T_zyy=-a_0/δ T_zxx,   T_yzy=-a_1/δ T_xzx,
        T_yyz=+a_2/δ T_xxz,


    A_02 T_zxy=A_01 T_zyx=   arsinh( Z_0/√(a_0δ)),
    
     A_01 T_yzx=A_12 T_xzy=   arsinh( Z_1/√(a_1δ)),
    
     -A_12 T_xyz=A_02 T_yxz=   arsinh( Z_2/√(a_2δ)),

and

    T_zzz=1/δ.

The low energy behavior of the quantities (<ref>) is dictated by the one of the local mixednesses a_0, a_1 and a_2 which diverge as 1/ω for all temperature. As a result

    lim_ω→ 0 T_zxx=lim_ω→ 0 T_xzx
        =lim_ω→ 0 T_xxz∀ T=0.

On the other hand, the behavior of the quantities (<ref>) depends on the behavior of the (square root δ of the) 
determinant of the covariance matrix which is unity at T=0 <cit.>. In this case 

    T_zzz(ω)T = 0=1,

and

    lim_ω→ 0 T_zyy=lim_ω→ 0 T_yzy
        =lim_ω→ 0 T_yyzT = 0= 1.

The reason for this behavior is that, at T=0, δ=1 whereas
the arguments of all the arctan terms in (<ref>) diverge as ω^-1/2, as shown by detailed inspection based on
Eqs. (<ref>), (<ref>), (<ref>) and the asymptotic expression of the coefficients of the S matrix given in <cit.>.

Alternatively, at finite temperature δ diverges at low energy (as 1/ω) and the limits (<ref>) all cancel

    lim_ω→ 0 T_zyy=lim_ω→ 0 T_yzy
        =lim_ω→ 0 T_yyzT  0= 0,

as also does lim_ω→ 0 T_zzzT  0=0.



§ ANALYTIC MAXIMIZATION OF ⟨ℬ̂^(I|2)⟩
 

In this appendix we present the maximisation of the expectation value of the CHSH operator (<ref>). We only state the results useful for the main text but do not detail the procedure because it is well known, see e.g., Refs. <cit.>.
In particular, it can be shown that the CHSH parameter B^(i|2)(ω) defined in  (<ref>) reads

    B^(i|2)=2√(λ_1+λ_2)  ,

where λ_1 and λ_2 are the two largest eigenvalues of 
matrix ^ T  T  T where T(ω) is the 3× 3 matrix with entries 

    T_rs
    =⟨Π̂^(i)_r ⊗Π̂^(2)_s⟩ ,

with (r,s)∈{x,y,z}^2. The results presented in Appendix <ref> show that in all the cases we consider the matrix
T(ω) is block diagonal:

    T=[ T_xx T_xy    0; T_yx T_yy    0;    0    0 T_zz;      ],

and thus

    ^ T  T  T=[      A      C      0;      C      B      0;      0      0 T^2_zz;        ],

with

    A    =  T_xx^2
    +  T_yx^2,  
    B = T_yy^2
    +  T_xy^2,
     
    C   =  T_xx T_xy+ T_yy T_yx.

The largest eigenvalues of matrix ^ T  T  T are
T_zz^2 and

    1/2( A+B + √((A-B)^2+4 C^2)).

They can be computed from Expressions (<ref>), (<ref>) and (<ref>). Then Eq. (<ref>) determines the value of the CHSH parameter (<ref>). 


One may also chose a different strategy for attempting 
to maximise the expectation value of the operator ℬ^(i|2). Instead of choosing to work in the basis of the c-
modes, one may perform a LLUBO for attempting to simplify the form of the covariance matrix.
As already stated in Appendix <ref>, in the bipartite case, the 4× 4 covariance
matrix (<ref>) associated with the reduced two-mode state (i|2) can be brought by LLUBOs to a standard form where the matrix ε_i2 are diagonal <cit.>. Working in this basis does not alter the entanglement properties of the system (they remain unaffected compared to that of the c-modes), but makes the computations easier and may improve the signal of nonlocality. In this basis the result (<ref>) is not affected and expressions
(<ref>) modify to

    T_xx
    =     2/πarctan2  | ⟨ĉ_iĉ_2⟩ |/√(A_i2)= -A_i2 T_yy, 
     T_xy =      T_yx = 
    0,

where A_i2 is defined in Eq. (<ref>). The matrix T is thus diagonal and expression (<ref>) is equal to T_xx^2. Equation (<ref>) then reads


    B^(i|2) =      2 √( T_xx^2+ T_zz^2)
    
    =     
    2  √(4/π^2arctan^2(2|⟨ĉ_iĉ_2⟩|/√(A_i2))
    + 1/A_i2^2)  .


We note here that the difference between the result  (<ref>) evaluated in the c-mode basis and expression (<ref>) is always small. The reason is that, in the c-mode basis, the off diagonal  entries of the upper left blocks of matrix (<ref>) and (<ref>) are always small compared to the diagonal ones, because at all temperature 
| ⟨ĉ_iĉ_2⟩|≪ 
| ⟨ĉ_iĉ_2⟩|.
 However, our numerical checks always demonstrate a small increase of the Bell parameter (<ref>) compared to the one evaluated using (<ref>) in the c-mode basis. We thus present our numerical results in Figs.  <ref>, <ref>, <ref>, 
<ref>, <ref> and <ref> using formula (<ref>). To summarize, the result (<ref>) should be considered as an optimized value of the witness of nonlocality B^(i|2).
We note here that we do not have a general proof of the better efficiency of the method which consists in using the basis in which the covariance matrix is in its standard form, but we believe that a general mathematical result of this type would be quite useful.





§ NUMERICAL MAXIMIZATION OF ⟨𝒮̂^(0|1|2)⟩
 

As explained in Sec. <ref>, the maximal value of the three-mode Bell operator ⟨𝒮̂^(0|1|2)⟩ can be found by optimizing the orientation of the unit vectors 
a, 
a', b, b', c and  c'. 
However, solving this optimization problem analytically proves challenging due to the need to maximize a function depending on 12 real parameters. Indeed, the orientation of each of the six previous normalized vectors in the 3-dimensional physical space corresponds to two degrees of freedom, leading in total to twelve parameters. 

Consequently, we resort to a numerical method to evaluate the maximal violation of Bell inequalities and determine the corresponding optimal orientations for the vectors a, a', b, b', c, and c'. More explicitly, we use a genetic algorithm which has proved very efficient for optimizing 
a function over a large parameter space <cit.>. This algorithm is based on natural selection: the code starts with a 
random set of solutions (in our case each of them consists of 12 
parameters) which form all together what we call a population. We then compute, for each set of vectors 
(a, a', b, b', c, c'), the expectation value 
⟨𝒮̂^(0|1|2)⟩ by means of the technique exposed in Appendix <ref>. For instance the contribution of the term
 ⟨a·Π̂^(0)⊗b·Π̂^(1)⊗c·Π̂^(2)⟩ to 
⟨ℬ̂^(0|1|2)⟩
 can be evaluated from the knowledge of the terms T_rst defined in Eq. (<ref>):

    ⟨a·Π̂^(0)⊗b·Π̂^(1)⊗c·Π̂^(2)⟩=∑_r,s,t a_r b_s c_t  T_rst,

where the sum runs over the indices (r,s,t)∈{x,y,z}^3. At variance with the bipartite case, at finite temperature it not not possible to find a LLUBO enabling to cast the 
6× 6 covariance matrix (<ref>) under a standard form where all the 
ε_ij matrices are diagonal (see the discussion in Appendix <ref>). The expectation values are thus computed in the natural basis of the c-modes where the values of the T_rst coefficients are given by  Eqs. (<ref>), (<ref>), (<ref>) 
and (<ref>).

At each step of the algorithm the code computes the expectation value 
⟨𝒮̂^(0|1|2)⟩ for a set of vectors (a, a', b, b', c, c') and ranks the members of the population by computing a fitness scaling function, a kind of selection rule: only the members of the population with the lowest fitness value will be retained – they are called the parents – and used to generate new sets of solutions, called the children. Then, at the next step, the selection rules are applied to the children, some of them become parents in turn and engender a new generation. The algorithm stops when all children look like their parents, or, in other words, when

    |v_ℓ+1 - v_ℓ| < δ,

where v_ℓ = (a, 
a', b, b',
c,c') is the set of solutions at step ℓ of the algorithm, and δ is the chosen convergence precision fixed prior the beginning of the selection process.

In our case, the fitness scaling function is simply the opposite of the average ⟨𝒮̂^(0|1|2)⟩ of the three-mode Bell operator. Trying to obtain the lowest fitness score is thus equivalent to maximize the Bell operator. For a given set of vectors v_ℓ at step ℓ, the next generation is computed as follows: v_ℓ+1 = v_ℓ + w_ℓ, where w_ℓ is a random weight which controls the mutations between the parents v_ℓ and the children v_ℓ+1, and which tends to decrease when the code starts to converge. For a detailed presentation of the algorithm we refer to Ref. <cit.>.
Note finally that the procedure just presented is also used for determining the optimized Mermin parameter (<ref>).



§ TRIPARTITE CIREL'SON BOUND AND ANALYTIC MAXIMIZATION OF ⟨𝒮̂^(0|1|2)(Ω=0)⟩
 

In this section we first study the upper bound of the quantity S^(0|1|2)(ω) defined in (<ref>) and then study the possibility of maximisation of the average ⟨𝒮̂^(0|1|2)⟩ of the operator (<ref>) by an appropriate choice of the measurement directions a, a', b, b', c and c'.

Denoting as Â=a·Π̂^(0), 
Â'=a'
·Π̂^(0),
B̂=b·Π̂^(1), etc. makes it possible to write the square of the tripartite Svetlichny operator (<ref>) as

    4 (𝒮̂^(0|1|2))^2=    8 + 
        {Â,Â'}⊗{B̂,B̂'}⊗{Ĉ,Ĉ'}
        - 2 [Â,Â']⊗[B̂,B̂']⊗1^(2)
        - 2 1^(0)⊗ [B̂,B̂']⊗ [Ĉ,Ĉ']
        - 2 [Â,Â'].⊗1^(1)⊗ [Ĉ,Ĉ'],

where [·,·] and {·,·} denote the commutator ans the anticommutator, respectively.
From the SU(2) algebra of the pseudo-spins it is easily proven that


    [Â,Â']    =2 i (a×a')·Π̂^(0), 
    {Â,Â'}    =2 (a·a') 1^(0),


with similar formulae for the quantities B̂, B̂' and
Ĉ, Ĉ'.
Since the operators Â, B̂ and Ĉ operate in different Hilbert spaces, one can assume without loss of generality that all the vector products of type (<ref>) appearing in (<ref>) are colinear with e_z.
In this case all the unit vectors a, a', b, b', c and c' lie in the xy plane.
Denoting as θ_a the angle between a' and a, θ_b the angle between b' and b, and θ_c the angle between c' and c (all these angles being in [0,π]), one gets

    (𝒮̂^(0|1|2))^2=   2+
    2cosθ_acosθ_bcosθ_c 
        +2sinθ_asinθ_b Π̂_z^(0)⊗Π̂_z^(1)⊗1^(2)
        +2sinθ_bsinθ_c 1^(0)⊗Π̂_z^(1)⊗Π̂_z^(2)
        +2sinθ_asinθ_c Π̂_z^(0)⊗1^(1)⊗Π̂_z^(2).

It is a simple matter to check that when θ_a, θ_b and θ_c run through [0,2π]^3 the eigenvalues of the operator appearing in the right hand side of the above formula are all comprised in [0,8]. The extremal values 0 and 8 are reached for (θ_a,θ_b,θ_c)=(0,0,π) and (π/2,π/2,π/2), respectively. It then follows that for any choice of the measurement directions
a, a', b, b', c and c'

    ⟨𝒮̂^(0|1|2)⟩^2 ≤⟨ (𝒮̂^(0|1|2))^2 ⟩≤ 8,

and thus the tripartite entanglement parameter (<ref>)
verifies

    S^(0|1|2)(ω)≤ 2√(2).

This is the tripartite equivalent of Cirel'son bound.

In the remaining of this Appendix we consider a related but somehow different problem: In order to violate as much as possible the Svetlichny inequality 
S^(0|1|2)<2 we aim at choosing 
measurement directions a, a', b, b', c and c' which 
maximize the expectation value of 𝒮̂^(0|1|2). This can be done numerically as explained in Appendix <ref>. We show here that this maximization can also be performed analytically in a particular instance.
Setting a·a' = cos 2θ with θ∈ [ 0, π/2 ], we define two new unit vectors a_± by

    2cosθ a_+ = a + a',  
    2sinθ a_- = a - a',

where by definition a_+·a_- = 0.
In terms of the new vectors (<ref>) the average of the three-mode operator (<ref>) reads

    ⟨𝒮̂^(0|1|2)⟩ =
    (W+Z) cosθ + (X-Y) sinθ,

where

    W     =   ⟨a_+·Π̂^(0)⊗b·Π̂^(1)⊗c'·Π̂^(2)⟩, 
    
    X     =   ⟨a_-·Π̂^(0)⊗b' ·Π̂^(1)⊗c' ·Π̂^(2)⟩, 
    
    Y     =   ⟨a_-·Π̂^(0)⊗b·Π̂^(1)⊗c·Π̂^(2)⟩,
    
    Z     = ⟨a_+·Π̂^(0)⊗b' ·Π̂^(1)⊗c·Π̂^(2)⟩.

It is convenient to introduce temporarily notations

    W+Z= Acosδ,   X-Y= Asinδ,

which make it possible to cast (<ref>) under a simple form

    ⟨𝒮̂^(0|1|2)⟩ =     Acos (θ - δ).

The maximum value of this expression is A and thus

    _θ⟨𝒮̂^(0|1|2)⟩ =
     √((W + Z )^2+(X-Y )^2).


The next step of the maximization procedure is easily performed at zero temperature and zero energy (ω=0).
The reason is that, as illustrated in a typical case by (<ref>), the explicit expressions of the W, X, Y, Z coefficients in 
(<ref>) and
(<ref>) involve combinations of terms of the type T_rst [as defined in Eq. (<ref>)] 
that take on particularly simple values at T=0 and ω=0.
Firstly, 
all the T_rst with at least one component along x vanish. Therefore, the contribution of the components along x of the different vectors involved in (WXYZ) cancels. It is thus enough to perform the maximization only considering vectors a_±, b, b', c and c' 
lying in the y,z plane. 
A second simplification stems from the fact that, in the y,z plane, all the T_rst's with an odd number of y cancel at T=0 and ω=0. In the waterfall and delta peak configurations[The situation is slightly different in the case of a flat profile configuration, where T_yyz = -1 and T_zzz =  T_zyy =  T_yzy = 1, with all other coefficients also vanishing. Once this modification is accounted for, the maximization procedure yields the same final result.] the only non-zero coefficients are (see the discussion in Appendix <ref>)

    T_zzz =  T_yyz =  T_yzy = -  T_zyy = 1.

In this case, the particular choice

    a_+ = b = c' =e_z,   a_- = - b' = - c =e_y,

plugged in Eqs. (<ref>) leads to

    W+Z= T_zzz+ T_zyy=2,
        X-Y=- T_yyz+ T_yzy=2.

and expression (<ref>) then shows that the upper bound (<ref>) is reached. It is therefore not possible that another arrangement of vectors
 a, a', b, b', c and c'
reaches a higher value and thus

    S^(0|1|2)(ω=0) T = 0=  2√(2).

The situation is completely different at finite temperature. In this case, as explained in Appendix <ref> all the averages of the type
(<ref>) cancel when ω→ 0. It then follows that the quantities (<ref>) behave in the same way and thus

    S^(0|1|2)(ω=0) T  0=  0.




§ RESULTS FOR OTHER TYPES OF ANALOG BLACK HOLES


The results we have obtained have been presented in the main text for a waterfall configuration. The reason is that this is the only one which has been realized experimentally so far <cit.>. 
For completeness – and also for emphasizing the typicality of the results presented in the main text –  we present in this appendix some equivalent results for the delta peak and flat profile configurations defined in Appendix <ref>. 



It is in particular of interest to plot the equivalents of
Fig. <ref> for these alternative configurations. This is
done in Figs. <ref> and <ref>. In the flat profile
configuration the value of m_d is a free parameter. Thus, for
comparing the results of the flat profile configuration with
Fig. <ref> we impose m_d=1/m_u^2, as is the case for the
waterfall configuration, cf. Eq. (<ref>). Such a procedure
is not required (nor possible) for the delta peak configuration where
fixing m_u unambiguously determines m_d, see (<ref>).
Fig. <ref> shows that the delta peak configuration has a
special feature: the EPR state formed by the companion and the partner
(see Sec. <ref>) is only reached at extremely low values of
the upstream mach number m_u, whereas bipartite nonseparability
between the Hawking quantum and the partner (modes 0 and 2,
respectively) is significant in a wide range of values of m_u
(roughly speaking, for m_u ≳ 0.2). Despite this peculiarity,
it is fair to say that Figs. <ref> and <ref> both
display the same general trend as Fig. <ref>, supporting the
idea that the behavior discussed in the main text is generic. We do
not produce finite temperature figures equivalent to
Fig. <ref> in order not to overload the paper and because,
as expected, the delta peak and the flat profile configurations behave
similarly to the waterfall configuration when the temperature is
increased.  There is however a quantitative change which is worth
noticing: in these two alternative configurations the violation of
bipartite Bell inequality does not persist as much as for the
waterfall configuration when T is increased.


 
It is also quite interesting to study genuine tripartite nonlocality in the delta peak and flat profile configurations by computing the zero temperature value of the Svetlichny parameter S^(0|1|2)(ω) and by then presenting
the equivalents of Fig. <ref>.


This is done in Figs. <ref> and <ref>. Here
again, the phenomenology is the same as the one discussed in the main
text: There is a clear signal of nonlocality at T=0, in a domain of
energy typically more extended than for the waterfall configuration,
but as is the case for the waterfall configuration this signal does
not persist at small finite temperature (the corresponding plots are
not shown in order not to overload the paper).




Finally, it is also important to address the GHZ character of the long
wave length modes in the delta peak and flat profile
configurations. As argued in Sec. <ref> of the main text,
both configurations display the GHZ paradox at ω=0 and verify
M^(0|1|2)(0)=4 at zero temperature.  The behavior of the optimized
Mermin parameter M^(0|1|2)(ω) (<ref>) is represented at
finite and zero temperature in Figs. <ref> and
<ref> which correspond to the delta peak and flat
profile configuration, respectively. At zero temperature the signal
(<ref>) of genuine tripartite entanglement is more pronounced
for the delta peak and flat profile configurations that for the
waterfall. This was also the case for the signal (<ref>) of
genuine tripartite nonlocality.  However, this signal, although more
noticeable at T=0, is less resilient to an increase of temperature
than for the waterfall configuration. The deleterious effect of
temperature is particularly pronounced for the flat profile
configuration, see Fig. <ref>. However, although it is
clear that each configuration bares its specificities, it is also
clear than the general trend is the same: the departure from the GHZ
signal M^(0|1|2)(0)=4 increases at high energy. At finite
temperature, the GHZ behavior is lost. The Mermin parameter
nonetheless provides a signal of nonlocality if larger than 2. For
that matter, the signal is more pronounced at T≠ 0 for the
waterfall configuration than for the two others.





