Bipartite and tripartite entanglement in a Bose-Einstein acoustic black hole
Mathieu Isoard, Nadia Milazzo, Nicolas Pavloff, Olivier Giraud
BBipartite and tripartite entanglement in a Bose-Einstein acoustic black hole
Mathieu Isoard,
1, 2
Nadia Milazzo,
1, 3
Nicolas Pavloff, and Olivier Giraud Universit´e Paris-Saclay, CNRS, LPTMS, 91405, Orsay, France Physikalisches Institut, Albert-Ludwigs-Universit¨at Freiburg, D-79104 Freiburg, Germany. Institut f¨ur theoretische Physik, Universit¨at T¨ubingen, 72076 T¨ubingen, Germany (Dated: February 12, 2021)We investigate quantum entanglement in an analogue black hole realized in the flow of a Bose-Einstein condensate. We argue that the system is described by a three-mode Gaussian state andconstruct the associated covariance matrix at zero and finite temperature. We study the associatedbipartite and tripartite entanglement measures and determine the best experimental configurationfor their observation.
I. INTRODUCTION
Analogue gravity aims at providing platforms makingit possible to conduct laboratory studies of phenomenaat the interface between general relativity and quantumphysics, such as Hawking radiation [1] and black hole su-perradiance [2], for which in the gravitational context di-rect observation is not possible and no complete theoryexists. It has also been suggested that analogue mod-els can bring some insight on the information loss para-dox [3, 4]. The concept has now broadened so as to in-clude experimental tests of physical effects of relevance incosmological scenarii, such as dynamical Casimir effect,Kibble-Zurek mechanism, Zakharov oscillations, Hubblefriction... see, e.g., [5] and references therein.However, in order to reach scientifically meaningful re-sults based on the study of an analogue model, it is im-portant to precisely characterize the experimental systemsupporting the analysis and to correctly circumscribe thephenomenon under scrutiny. The present work aims atfollowing this line of research in the case of an analogueof event horizon realized in a Bose-Einstein condensed(BEC) ultracold atomic vapor. The use of a BEC as ananalogue model has been first suggested by Garay et al. [6], followed by many others. This motivated Steinhauerand his group to develop and then ameliorate an exper-imental setup making it possible to realise an acoustichorizon in a quasi one-dimensional BEC [7–10]. Con-comitantly, the theoretical study of this system by meansof a Bogoliubov decomposition has been first suggested in[11], then gradually refined [12–14] until a point where adetailed comparison with experiments has been possible[15]. In this line, a crucial question is the quantum na-ture of the Hawking radiation observed in [9, 10]: is thephenomenon simply triggered by noise or does it corre-spond to spontaneous quantum emission as in Hawking’soriginal scenario ? A natural test of the latter consists indemonstrating entanglement of the Hawking pair.Entanglement detection and characterization has at-tracted a great deal of effort in the past two decades, asit has been identified as a key resource for quantum in-formation processing [16]. A quantum state is entangledif it is not separable, i.e., if it cannot be written as aconvex sum of product states [17]. One of the simplest necessary separability criteria is given by the positivityof the partial transpose (PPT), first proposed for discretevariables [18, 19] and extended to the continuous case in[20]. A wealth of entanglement measures were discussedin the literature, both for discrete and continuous vari-ables. For bipartite pure states, quantitative measures ofentanglement include the entanglement entropy (whichcan be shown to be unique if some additional natural re-quirements are imposed), or the concurrence [21]. In themixed state case, it is possible to construct ’good’ en-tanglement measures in many different ways, which areinequivalent in the sense that they lead to different order-ings of entangled states [16]. A possible way is to extendmeasures for pure states via a convex roof construction:entanglement of a mixed state is then defined by a mini-mization over all its possible pure state decompositions.For instance, entanglement entropy generalizes for mixedstates to the entanglement of formation [22].A striking difference between classical correlations andquantum entanglement is that the latter is monogamous.This means that a particle which is maximally entan-gled with a partner cannot be entangled with a thirdparty, or in other words that any amount of entangle-ment shared with a particle limits the entanglement thatcan be shared with another particle. In the case ofthree qubits, this limitation to bipartite entanglementwas expressed in [23] through an inequality that must besatisfied by an entanglement measure called the concur-rence, or more precisely by its square, the tangle. Thismonogamy inequality was later generalized to an arbi-trary number of qubits [24], and to continuous variables,as we now discuss.In the case of continuous variables, to which the sit-uation of black hole analogues pertains, Gaussian statesare the most natural objects with which one is led todeal. From a qualitative point of view, entanglementcan be detected by the PPT criterion, which is a neces-sary and sufficient separability condition for 1 × N -modeGaussian states [20] (the three-mode case, which is rel-evant to our situation, was investigated in [25]). Quan-titatively, entanglement can be measured by the loga-rithmic negativity, which measures by which amount thePPT criterion is violated [26]. For continuous variablesit is generally highly difficult to make use of the con- a r X i v : . [ c ond - m a t . qu a n t - g a s ] F e b vex roof construction, both analytically and numerically,as the optimization has to take place over all pure statedecompositions. To circumvent this issue, Gaussian en-tanglement of formation was defined in [27], restrictingthe convex roof construction to Gaussian pure state de-compositions. This quantity provides an upper boundfor the entanglement of formation, and is more amenableto calculations. In [28] it was shown that Gaussian en-tanglement of formation and entanglement measured bynegativity are inequivalent measures. In [29], it was pro-posed to construct a specific measure of entanglement,the contangle (continuous tangle), defined as the convexroof extension of the square of the logarithmic negativity,so that the monogamy inequality expressed by this mea-sure also holds for Gaussian states. It was generalized to N -mode Gaussian states in [30]. The amount by whichboth sides of the monogamy inequality differ provides anestimate for multipartite entanglement.Several theoretical works have addressed the issue ofentanglement in analogue gravity systems. Most of them[31–37] discuss qualitative measures such as the Peres-Horodecki or Cauchy-Schwarz criterion. In the presentwork we follow Refs. [38–40] and focus on quantitativemeasures. To conduct this computation, the specifici-ties of BEC physics are important to take into account.In particular, dispersive effects complexify the standardHawking quantum/partner picture by introducing newpropagation channels. Accordingly the system is de-scribed by a three-mode Gaussian state, and we buildbipartite and tripartite entanglement measures from thecovariance matrix of bipartitions of the system.The paper is organized as follows. In Section II wereview the formalism describing acoustic black holes. InSection III we review the basics of Bogoliubov transfor-mations and apply it to our situation. The descriptionof Gaussian states appearing in the scattering processsesinvolved in black hole analogues is discussed in SectionIV. Section V is dedicated to the investigation of two-mode and three-mode entanglement in Gaussian states,and more specifically to the Gaussian states we are con-sidering here. The case of a finite-temperature setting isexamined in Section VI, and concluding remarks are pre-sented in Section VII. Some technical points are given inthe appendixes. In Appendix A we recall some propertiesof Bogoliubov transformations. In Appendix B we givesome useful explicit expressions of the covariance matrix.In Appendix C we recall the low frequency behavior ofthe coefficients describing the scattering of linear wavesby the acoustic horizon. Appendix D details the con-struction making it possible to localize intrication in oursystem. In Appendix E we establish a formula makingit possible to compute the Gaussian contangle at finitetemperature. x U ( x ) = − U Θ( x ) n ( x ) = | Φ( x ) | n u → ← n d V u V d FIG. 1. Waterfall configuration. The flow is directed from leftto right. The downstream classical field Φ( x >
0) is exactlya plane wave of density n d and velocity V d . Φ( x <
0) is thehalf profile of a dark soliton which is asymptotically a planewave of density n u and velocity V u , see Eqs. (3). The farupstream flow is subsonic, with a velocity 0 < V u < c u and thedownstream flow is supersonic with a velocity V d > c d > c α = ( gn α /m ) / is the speed of sound in region α = u (upstream) or d (downstream). The region x > II. ANALOGUE BLACK HOLE IN BECS
We consider a stationary flow of a one-dimensionalBEC which is upstream subsonic and downstream super-sonic. This ”transsonic” configuration mimicks a blackhole since acoustic excitations generated in the down-stream supersonic region are expected to be dragged bythe flow, and not detected in the upstream region.
A. The background flow
The quantum operator ˆΨ is decomposed into a classicalfield Φ (describing the stationary flow of the condensate)supplemented by small quantum fluctuations (describedby an operator ˆ ψ ) according toˆΨ( x, t ) = exp( − iµt/ (cid:126) ) (cid:16) Φ( x ) + ˆ ψ ( x, t ) (cid:17) , (1)where µ is the chemical potential . The function Φ issolution of a classical Gross-Pitaevskii equation, with theaddition of an external potential U ( x ) used to implementthe transsonic flow: µ Φ( x ) = − (cid:126) m ∂ x Φ + (cid:0) g | Φ | + U ( x ) (cid:1) Φ , (2)where g > The validity of this decomposition is discussed, e.g., in [36] operator ˆ ψ describes the quantum fluctuations on top ofthe background Φ.Several configurations realizing an analogue black holehave be proposed in the past [12, 41–43]. The approachwe use in this work is valid in a general setting, butfor the sake of illustration we will present numerical re-sults for the so-called ”waterfall configuration” [14] whichhas been experimentally realized in [9, 10] and whichhas been shown to lead to a significant violation of theCauchy-Schwarz criterion in [36]. In this configuration U ( x ) = − U Θ( x ), where U > n d and velocity V d > x >
0) and half a dark soliton in theupstream region ( x <
0) with asymptotic density n u andvelocity V u >
0, meaning thatΦ( x >
0) = √ n d exp ( imV d x/ (cid:126) ) exp( i β d ) , Φ( x → −∞ ) = √ n u exp ( imV u x/ (cid:126) ) exp( i β u ) , (3)where β u and β d are constant phase factors. This settingis illustrated in Fig. 1 (see details in [14]). B. Elementary excitations
Since the far upstream and downstream backgroundflows are uniform, the elementary excitations which forma basis set for the quantum operator ˆ ψ are plane waves inthese two regions, with dispersion relations of Bogoliubovtype (see, e.g., [44]):( ω − q V α ) = ω B ,α ( q ) , α = u or d , (4)where V u and V d are the upstream and downstream ve-locities, and ω B ,α is the Bogoliubov dispersion relation ω B ,α ( q ) = c α q (cid:112) ξ α q / , (5) c α = ( gn α /m ) / being the speed of sound and ξ α = (cid:126) / ( mc α ) the ”healing length”, in the far upstream regionif α = u and in the downstream region if α = d . Theleft-hand side of Eq. (4) includes a Doppler shift causedby the velocity V α of the background.It will be useful in the following to define the quantities m α = V α c α , α = u or d , (6)known as the upstream ( α = u ) and downstream ( α = d )Mach numbers. It was shown in [14] that in the water-fall configuration which we use below to exemplify ourresults, the following relations hold V d V u = n u n d = 1 m u = m d = (cid:18) ξ d ξ u (cid:19) = (cid:18) c u c d (cid:19) , (7)which means that this configuration is uniquely charac-terized once the value of m u , say, is fixed. The flow being upstream subsonic ( V u < c u , i.e., m u <
1) and downstream supersonic ( V d > c d , i.e., m d > | in and 0 | out. In the downstream supersonic regionthere are four branches: 1 | in, 1 | out, 2 | in and 2 | out, thelast two branches being limited to ω ∈ [0 , Ω], where Ω isthe frequency at which these two branches coalesce, andwhose value is given byΩ = q ∗ V d − ω B ,d ( q ∗ ) with q ∗ ξ d = (cid:18) − m d m d (cid:113) m d (cid:19) . (8)For future convenience (see Sec. VI) we define functions q | in ( ω ), q | in ( ω ) and q | in ( ω ) as the reciprocal of theBogoliubov dispersion relation (4) along some of thesebranches; q | in ( ω ) and q | in ( ω ) are defined for ω > q | in ( ω ) only for ω ∈ [0 , Ω]. A number of previous works[14, 15, 31, 36, 43] followed the convention introduced in[13], in which indices u , d
1, and d ω lower than thethreshold Ω, to a specific scattering process of elemen-tary excitations onto the analogue event horizon. Forinstance, a wave issued from the interior region along thechannel identified as 1 | in in Fig. 2 is transmitted to theexterior along the 0 | out channel and reflected back alongthe 1 | out and 2 | out channels. The corresponding (com-plex) transmission and reflection amplitudes are denotedas S ( ω ), S ( ω ) and S ( ω ), respectively. They areobtained by imposing matching conditions at the hori-zon, as explained in Ref. [14]. The quantum boson op-erator corresponding to this whole process is denoted asˆ b ( ω ). Similarly, a wave incident along the 0 | in channel istransmitted towards the interior of the black hole alongchannels 1 | out (amplitude S ( ω )) and 2 | out (amplitude S ( ω )) and reflected along 0 | out (amplitude S ( ω )); thecorresponding quantum mode is associated with opera-tor ˆ b ( ω ) . A third mode describes the scattering ofa wave issued from the channel 2 | in onto the outgoingchannels 0 | out, 1 | out and 2 | out. The channels labelled2 | in and 2 | out are particular, in the sense that they havea negative norm (i.e., a negative energy in the rest frameof the fluid). We shall not discuss this subtelty furtherhere, and refer the reader to [45–48]. As a result, themode initiated by the incoming channel 2 | in should be In the terminology we use, it is important to make a distinctionbetween the quantum modes and the propagation channels: amode corresponds to a whole process typically involving one orseveral incoming channels and one or several outgoing channels. q [a.u.] ← | in 1 | out → | out → ← | in( q ∗ , Ω) − q q q [a.u.] ω [ a . u . ] | in →← | out FIG. 2. Graphical representation (in arbitrary units) of the positive frequency part of the dispersion relation (4) in thefar upstream (left plot) and downstream (right plot) regions. The downstream region (grey background) is the interior ofthe analogue black hole while the upstream region (white background) is the exterior. In the upstream region, any given ω (represented by a horizontal dashed line) corresponds to two channels of propagation denoted as 0 | in and 0 | out. In thedownstream region it corresponds to four or two solutions, depending if ω is larger or smaller than Ω. The arrows indicate thedirection of propagation of the corresponding waves, and the channels are labelled 1 or 2, with an additional ”in” (or ”out”)indicating if the wave propagates towards (or away from) the horizon. quantized using an operator ˆ b † ( ω ), i.e. inverting the roleof the creation and annihilation operators used for thetwo other modes. Only in this way do the propagatingmodes behave as bosons satisfying the usual commuta-tion relations (cid:104) ˆ b i ( ω ) , ˆ b † j ( ω (cid:48) ) (cid:105) = δ i,j δ ( ω − ω (cid:48) ) (cid:104) ˆ b i ( ω ) , ˆ b j ( ω (cid:48) ) (cid:105) = (cid:104) ˆ b † i ( ω ) , ˆ b † j ( ω (cid:48) ) (cid:105) = 0 , (9)for i and j ∈ { , , } . Another consequence is thatthe 3 × S ( ω ) whose elements are the S ij ( ω ) obeys a skew-unitarity relation [13] : S † ηS = η = SηS † , η = diag(1 , , − . (10)For ω > Ω the situation is drastically different: the chan-nels 2 | in and 2 | out disappear (cf. Fig. 2), as well as theoperator ˆ b ( ω ), and the S -matrix becomes 2 × b modes as ”incoming” since they cor-respond to scattering processes initiated by a single waveincident along one of the three ”in” channels directed to-wards the horizon: 0 | in, 1 | in and 2 | in. One could equiv-alently choose to work with ”outgoing modes” [11] de-scribing processes each resulting in the emission of a sin-gle wave along one of the three ”out” channels 0 | out,1 | out and 2 | out. We denote the corresponding quantumoperators as ˆ c ( ω ), ˆ c ( ω ) and ˆ c ( ω ). They relate to theincoming operators via [13] ˆ c ˆ c ˆ c † = S S S S S S S S S ˆ b ˆ b ˆ b † , (11)where for legibility we omit the ω dependence of all theterms. The definition (11), together with the property(10), ensures that the ˆ c operators obey the same commu-tation relations (9) as the ˆ b operators and thus describebosonic quasiparticles. In the setting we consider, the analogue of the Hawk-ing radiation spectrum is the number of excitations emit-ted per unit time and per unit frequency into the sub-sonic region ( x < c † ( ω )ˆ c ( ω ) over the state vector. From relation (11) onesees that this current is non zero when the state vector isthe vacuum | (cid:105) b of incoming modes: b (cid:104) | ˆ c † ( ω )ˆ c ( ω ) | (cid:105) b = | S ( ω ) | ; this is the analogous Hawking effect. Themode associated with operator ˆ c is thus denoted theHawking outgoing mode. The other outgoing modes, as-sociated with operators ˆ c and ˆ c are denoted the com-panion and the partner, respectively.As can be seen from expression (11), the outgoing oper-ators ˆ c and ˆ c † are expressed as a combination of the ingo-ing annihilation and creation operators ˆ b and ˆ b † . There-fore, it is possible to associate a Bogoliubov transforma-tion with our analogue system. This is the aim of thenext section. III. BOGOLIUBOV TRANSFORMATIONS
Bogoliubov transformations are linear transformationsof creation and annihilation operators that preserve thecanonical commutation rules [49]. In the context of quan-tum field theory in curved spacetime, these transforma-tions are at the heart of the Hawking process; indeed,since they mix annihilation and creation operators theycan give rise to spontaneous emission of particles fromvacuum [50–55]. This mixing of operators also occursfor analogue black holes, as is clear from Eq. (11). Thisway of viewing the emergence of the analogue Hawkingradiation through a Bogoliubov transformation makes adirect connection with the gravitational case: as shownby Hawking in 1974 [50, 51], one of the parameters in-volved in the Bogoliubov transformation, the so-called β -coefficient, is directly related to the number of particlescreated by black holes. In our case, we can derive such aparameter and compare its properties with Hawking’s β -coefficient; in particular, through this approach, we willbe able to question the thermality of the analogue Hawk-ing radiation (see section IV C). Furthermore, identifyingthe Bogoliubov transformation will be an important stepto understand and study the entanglement properties ofthe analogue Hawking radiation (see section V).The section is divided into two parts. First, we con-sider an arbitrary (but unitary) Bogoliubov transforma-tion and derive its properties. Then, we apply these re-sults to the particular case of analogue black holes inBECs starting from expression (11). A. General setting
We start by briefly recalling some well-known factsconcerning unitary Bogoliubov transformations [45, 49,56]. Some useful intermediate results are given in Ap-pendix A.Let us consider N boson operators ˆ b , . . . , ˆ b N satisfyingthe usual commutation relations [ˆ b i , ˆ b † j ] = δ i,j . Definingthe column vector b = (ˆ b , . . . , ˆ b N , ˆ b † , . . . , ˆ b † N ) T , (12)the Bose commutation relations can be rewritten as[ b i , b j ] = (cid:101) J ij , with (cid:101) J = (cid:18) N − N (cid:19) , (13)where N is the N × N identity matrix. A (unitary) Bo-goliubov transformation is a linear transformation map-ping the operators ˆ b i onto new operators ˆ c i definedthrough c i = N (cid:88) j =1 T ij b j , or equivalently c = T b . (14)For unitary Bogoliubov tranformations c has the form c = (ˆ c , . . . , ˆ c N , ˆ c † , . . . , ˆ c † N ) T , (15)i.e., that c i + N = c † i . In this case the matrix T admitsthe block decomposition T = (cid:18) α ∗ − β ∗ − β α (cid:19) , (16)where α and β are N × N matrices. Operators ˆ c i and ˆ b i can then be related by a unitary operator T such thatˆ c i = T † ˆ b i T, (17)whose explicit construction from matrix T is detailed inAppendix A.In general, the transformation T in (14) mixes cre-ation and annihilation operators, so that the vacua | (cid:105) b and | (cid:105) c , defined byˆ b i | (cid:105) b = 0 , and ˆ c i | (cid:105) c = 0 , i ∈ { , . . . , N } , (18) differ. These vacua are related via the identity | (cid:105) b = T | (cid:105) c , (19)as is clear from the fact that ˆ b i T | (cid:105) c = T ˆ c i | (cid:105) c = 0.Defining the N × N matrix X = − β ∗ α − and usingthe decomposition (A7), it is possible to write Eq. (19)under the explicit form | (cid:105) b = 1(det α ) e (cid:80) i,j X ij ˆ c † i ˆ c † j | (cid:105) c . (20)A simple example of a Bogoliubov transformation is theone leading to two-mode squeezed states [57, 58]. For areal squeezing parameter r , a two-mode squeezed state isobtained by applying the two-mode squeezing operator T = e r ( b † b † − b b ) (21)to the vacuum state. The corresponding Bogoliubovtransformation is of the form (16) with N = 2 and α, β given by α = (cid:18) cosh r
00 cosh r (cid:19) , β = − (cid:18) r sinh r (cid:19) . (22) B. Bogoliubov transformation in a transsonic BEC
In the case described in Sec. II of a transsonic flowrealised in a BEC, b and c correspond to sets of ingo-ing and outgoing modes. The associated column vectors b = (ˆ b , ˆ b , ˆ b , ˆ b † , ˆ b † , ˆ b † ) T and c = (ˆ c , ˆ c , ˆ c , ˆ c † , ˆ c † , ˆ c † ) T are related by Eq. (11). One can express this relationequivalently as c = T b , with T a Bogoliubov transfor-mation of the form (16) with α = S ∗ S ∗ S ∗ S ∗
00 0 S , β = − S ∗ S ∗ S S , (23)where for legibility we do not write the ω -dependence ofthe scattering amplitudes. This yields X = 1 S S S S S . (24)From relation (10) one can show that det α = | S | , andthus (20) takes the simple form | (cid:105) b = 1 | S | e ( X ˆ c † + X ˆ c † ) ˆ c † | (cid:105) c . (25)A word of caution is in order here. The case we con-sider in the present section is different from the discus-sion of the previous section III A because the modes arehere continuously distributed along the energy axis, seeSec. II B. A natural way to set up a framework encom-passing the two situations consists in discretizing the en-ergies with a small mesh ∆ ω and to define coarse-grainedoperators ˆ B i,p = 1 √ ∆ ω (cid:90) ω p +1 ω p d ω ˆ b i ( ω ) , (26)and ˆ C i,p = 1 √ ∆ ω (cid:90) ω p +1 ω p d ω ˆ c i ( ω ) , (27)where i ∈ { , , } , p ∈ N and ω p = p ∆ ω . It is easy tocheck that these operators obey the standard Bose com-mutation rules, such as [ ˆ B i,p , ˆ B † j,q ] = δ i,j δ p,q for instance.If ∆ ω is small compared to the typical scale of variationof the elements of the S -matrix, then the ˆ C i,p and theˆ B j,p are related by a relation analogous to (11): ˆ C ,p ˆ C ,p ˆ C † ,p = S S S S S S S S S ˆ B ,p ˆ B ,p ˆ B † ,p , (28)where the S ij should be evaluated at ω p . Thus the rela-tion (25) should be replaced by | (cid:105) b = 1 (cid:81) ∞ p =0 | S ( ω p ) | × e (cid:80) ∞ p =0 ( X ( ω p ) ˆ C † ,p + X ( ω p ) ˆ C † ,p ) ˆ C † ,p | (cid:105) c . (29)This remark being made, in the following we favor legibil-ity over formal rigor: We will continue to write relationsof the type (25), instead of the more rigourous but cum-bersome Eq. (29), keeping in mind that the correctionof ”naive” expressions – such as Eqs. (30), (31) or (36)below – is straightforward.From (25), if we define the Fock state basis of quasi-particles of type c by | n (cid:105) i = 1 √ n ! (ˆ c † i ) n | (cid:105) c , (30)where i is the mode number, then the explicit expansionof the vacuum | (cid:105) b reads | (cid:105) b = 1 | S | ∞ (cid:88) n,n (cid:48) =0 (cid:115)(cid:18) n + n (cid:48) n (cid:19) X n X n (cid:48) | n (cid:105) | n (cid:48) (cid:105) | n + n (cid:48) (cid:105) . (31)It is convenient for future use in sections IV D and V C tointroduce a new set operators e = (ˆ e , ˆ e , ˆ e , ˆ e † , ˆ e † , ˆ e † ) T .By writing S ij ( ω ) = v ij ( ω ) e i ϕ ij ( ω ) , v ij ≥ , ≤ i, j ≤ , (32)we define the operators ˆ e , ˆ e and ˆ e asˆ e = e − i ϕ ˆ c , ˆ e = e − i ϕ ˆ c , ˆ e = e i ϕ ˆ c (33) (note the + sign in front of ϕ ). This defines a local uni-tary Bogoliubov transformation, as it does not mix anni-hilation and creation operators. In particular | (cid:105) e = | (cid:105) c .Using the notations of section III A, this transformationcan be cast in the form e = R c , (34)where R = diag (cid:0) e − i ϕ , e − i ϕ , e i ϕ , e i ϕ , e i ϕ , e − i ϕ (cid:1) . (35)Then, using expression (24) and this new set of creationand annihilation operators e , Eq. (25) becomes | (cid:105) b = 1 v e v − ( v ˆ e † + v ˆ e † ) ˆ e † | (cid:105) e . (36) IV. THREE-MODE GAUSSIAN STATESA. Gaussian states
In order to set up notations we start by reviewing theformalism for Gaussian states (see [59] for a review).Gaussian states are states whose Wigner function is aGaussian. A Gaussian state ρ can be entirely describedby its first and second moments. We define the covari-ance matrix σ of ρ as the real symmetric positive-definitematrix σ ij ≡ (cid:104) ˆ ξ i ˆ ξ j + ˆ ξ j ˆ ξ i (cid:105) − (cid:104) ˆ ξ i (cid:105) (cid:104) ˆ ξ j (cid:105) , (37)where ˆ ξ i are components of the vector ξ = √ q , ˆ p , . . . , ˆ q N , ˆ p N ) T of quadratures relative to mode i , defined so that [ˆ q i , ˆ p j ] = i δ i,j . In the definition (37)and in all the following the averages (cid:104)· · · (cid:105) are taken overthe density matrix rho characterizing the state of the sys-tem, which, in the simpler case, is the projector onto thevacuum state | (cid:105) b . We shall discuss in section VI how togeneralize to a finite-temperature configuration.The commutation relations between the ˆ ξ i can be ex-pressed as [ ˆ ξ i , ˆ ξ j ] = 2 i J ij , ∀ i, j ∈ { , . . . N } with J = N ⊕ J i , J i = (cid:18) − (cid:19) . (38)Entanglement properties of a quantum state are un-changed by local unitary (LU) operations, so that meanvalues of position and momentum operators can be set to0. A Gaussian state is then entirely specified by its co-variance matrix, which can be rewritten in terms of 2 × σ = σ ε · · · ε N ε T . . . . . . ...... . . . . . . ε N − N ε T N · · · ε T N − N σ N , (39)with ε ij = 2 (cid:18) (cid:104) ˆ q i ˆ q j (cid:105) (cid:104) ˆ q i ˆ p j (cid:105)(cid:104) ˆ p i ˆ q j (cid:105) (cid:104) ˆ p i ˆ p j (cid:105) (cid:19) (40)and σ i = (cid:18) (cid:10) q i (cid:11) (cid:104){ ˆ q i , ˆ p i }(cid:105)(cid:104){ ˆ q i , ˆ p i }(cid:105) (cid:10) p i (cid:11) (cid:19) , (41) { ., . } denoting the anticommutator. A covariance matrix σ satisfies the inequality σ + i J ≥ , (42)which is a consequence of the canonical commutation re-lations and positivity of the density matrix [60, 61]. Inparticular, σ is a positive matrix. B. Transformations of Gaussian states
Partial trace of a Gaussian state is particularly simple.The covariance matrix of the reduced state is simply ob-tained by discarding the lines and columns correspondingto the modes over which the partial trace is done (see,e.g., [62]). For instance, the two-mode state obtainedfrom (39) by tracing out modes 3 to N has covariancematrix σ ij = (cid:18) σ i ε ij ε T ij σ j (cid:19) , (43)where the 2 × ρ ( i ) of mode i is a single-mode Gaussian state entirely specified by thecovariance matrix σ i .Let us now turn to the transformation of the covari-ance matrix under a Bogoliubov transformation. Forthe vector of creation and annihilation operators b , wedenote by ξ b the corresponding vector of position andmomentum operators ξ b = √ q , ˆ p , . . . , ˆ q N , ˆ p N ) T withˆ q j = (ˆ b j + ˆ b † j ) / √ p j = i (ˆ b † j − ˆ b j ) / √
2. We thus have ξ b = U b , with U = √ . . . √ . . . − i √ . . . i √ . . . √ . . . √ . . . − i √ . . . i √ . . . . . . . . . (44)a 2 N × N unitary matrix. Similarly, ξ c = U c . TheBogoliubov transformation c = T b then entails that ξ c = S T ξ b with S T = U T U † . It can be proved thatthe matrix S T ∈ Sp(2 N, R ) is real and symplectic [61].Since this transformation is linear, we get from Eq. (37)that a Gaussian state in mode b is a Gaussian state inmode c with covariance matrix σ (cid:48) = S T σ S T T . (45) As guaranteed by Williamson theorem [63], it isalways possible to find a symplectic transform thatbrings σ to a canonical diagonal matrix σ (cid:48) =diag( ν , ..., ν N , ν , ..., ν N ), which is unique up to the or-dering of the ν j . The ν j are called the symplectic eigen-values of σ . They can be directly obtained from the eigen-values of the matrix J σ , which are given by ± iν j [64]. Interms of the ν j the uncertainty relation (42) immediatelyyields ν j ≥ j = 1 , ..., N . (46) C. Thermal states
The symplectic eigenvalues have an appealing physi-cal interpretation. Indeed, they can be related with themean particle number of a thermal state. Recall that ageneric (single-mode) thermal state is a state whose den-sity matrix in the Fock space spanned by vectors | n (cid:105) isof the form ρ th ( a ) = 2 a + 1 ∞ (cid:88) n =0 (cid:18) a − a + 1 (cid:19) n | n (cid:105) (cid:104) n | , (47)with a some parameter. Denoting as ˆ n the correspondingnumber operator, since ¯ n ≡ (cid:104) ˆ n (cid:105) = tr( ρ th ˆ n ) = ( a − a is simply related with the mean parti-cle number as a = 2¯ n + 1. The state ρ th ( a ) is in fact aGaussian state with 2 × σ th = a .This means that a single-mode covariance matrix in di-agonal form describes a thermal state with mean particlenumber ¯ n = ( a − / ν = a .Another way of representing a thermal state is to set¯ n = sinh r , which yields ρ th ( a ) = 1cosh r ∞ (cid:88) n =0 (tanh r ) n | n (cid:105)(cid:104) n | , a = cosh(2 r ) . (48)The vacuum state is a thermal state ρ th (1) with meanoccupation numbers ¯ n j = 0. The purity of ρ th ( a ) canbe readily calculated from (47); it reads tr[ ρ th ( a )] =1 / cosh(2 r ) = 1 /a . The quantity a being the inverse ofthe purity of a single-mode reduced density matrix, it isreferred to as the local mixedness [65].More generally [66], an N -mode Gaussian state witharbitrary covariance matrix σ can be brought to a prod-uct of thermal states ρ ν = N ⊗ j =1 ρ th ( ν j ) . (49)Indeed, if S is the symplectic transformation that di-agonalizes σ as diag( ν , ν , ..., ν N , ν N ) = S σ S T , then itcan be realized on the Gaussian state by a unitary evolu-tion generated by a quadratic Hamiltonian (see e.g. [61]).Therefore Williamson’s theorem ensures that any Gaus-sian state can be decomposed into a thermal state whosemean occupation number in mode j is obtained from thesymplectic eigenvalue ν j as ¯ n j = ( ν j − / ν j ≥ ρ ν ) = 1 / √ det σ . D. Vacuum as a three-mode Gaussian state
The Bogoliubov transformation associated with thescattering process (11) leads to a three-mode Gaussianpure state, given by (31). The covariance matrix of thevacuum | (cid:105) b is the identity matrix . Applying (45)we thus get that the covariance matrix of state (31) is S T S T T . Using the explicit expression of T derived from(23), and the explicit expression (44) for U , we obtainthe 6 × × σ i and ε ij , i ∈ { , , } , simply read σ i = (cid:16) (cid:104) ˆ c † i ˆ c i (cid:105) (cid:17) (50)and ε i = (cid:18) (cid:104) ˆ c i ˆ c (cid:105) (cid:104) ˆ c i ˆ c (cid:105) (cid:104) ˆ c i ˆ c (cid:105) − (cid:104) ˆ c i ˆ c (cid:105) (cid:19) ,ε = (cid:18) (cid:104) ˆ c ˆ c † (cid:105) (cid:104) ˆ c ˆ c † (cid:105)− (cid:104) ˆ c ˆ c † (cid:105) (cid:104) ˆ c ˆ c † (cid:105) (cid:19) . (51)In the two-mode and three-mode cases, it is known[67, 68] that all pure Gaussian states can be broughtby LLUBOs (local linear unitary Bogoliubov transfor-mations) to a standard form where matrices σ i are pro-portional to the identity and matrices ε ij are diagonal.In order to get such a standard form, we use the set ofoperators ˆ e j related with the ˆ c j by e = R c , where R has been defined in (35). Using the results of IV B andapplying Eq. (45), the covariance matrix of state (36) inmode e is σ = S R S T S T T S T R . (52)Note that S R is a rotation operator. Indeed, one caneasily find that S R = diag { R ( ϕ ) , R ( ϕ ) , R ( − ϕ ) } ,where R ( φ ) = (cid:18) cos φ sin φ − sin φ cos φ (cid:19) . (53)Then, one proves that the covariance matrix σ defined by(52) is in the standard form, with σ i given by (50) and ε ij = 2 |(cid:104) ˆ c i ˆ c † j (cid:105)| = 2 v i v j , i, j = 0 , ε i = 2 |(cid:104) ˆ c i ˆ c (cid:105)| σ z = 2 v i v σ z , i = 0 , , (54)where σ z is the third Pauli matrix. Following the no-tation introduced in section IV C for thermal states, wedefine real parameters r i ≥ a i ≥ n i = (cid:104) ˆ c † i ˆ c i (cid:105) = sinh ( r i ) = a i − . (55) . . . . . . . ¯ h ω/ ( g n u ) . . . . . . . . a i Ω ω c a a a FIG. 3. Local mixedness a i ( ω ) [see Eqs. (56)] for each mode 0,1 and 2 as functions of the dimensionless quantity (cid:126) ω/ ( g n u ),for a waterfall configuration with m u = 0 .
59. The frequency ω c indicates the turning point above which a becomes lowerthan a . The upper-bound frequency Ω corresponds to thevanishing of the mode 2 [see Eq. (8)]. To be completely accurate, we recall that the operatorsˆ c i = ˆ c i ( ω ) all depend on the energy (cid:126) ω of the elemen-tary excitations. Therefore, the above defined quantities a i also depend on ω . They can be written explicitly asfunctions of the scattering matrix coefficients: a ( ω ) = 1 + 2 | S ( ω ) | a ( ω ) = 1 + 2 | S ( ω ) | a ( ω ) = − | S ( ω ) | (56)(see Appendix B). From the solution of the scatteringproblem in the waterfall configuration, we calculate thescattering amplitudes S ij ( ω ) following [14]. This makesit possible to compute the three local mixednesses a , a and a as functions of the frequency. In particular wehave a ( ω ) + a ( ω ) = a ( ω ) + 1 , (57)which stems from relations (56) and (10). Figure 3shows the associated curves. Here, these coefficients arecomputed for a waterfall configuration with downstreamMach number m d = 2 .
9, which is the one for which theexperiment of [10] has been realized. In our case, thiscorresponds to an upstream Mach number m u = 0 . ω c the lowest of the three parametersis a ; above this frequency, the minimum value becomes a . The value of this frequency is determined numeri-cally and is equal to ω c = 0 . g n u / (cid:126) for m u = 0 .
59. Weobserve that the ratio ω c / Ω (where Ω is the frequency(8) at which mode 2 vanishes and also depends on m u )decreases when m u decreases. The local mixednesses a , a and a go to 1 when ω → Ω, which means that thepopulations of all modes vanish.Using Eqs. (56) one may rewite expressions (50) and(54) in terms of the a i as σ i = a i , i = 0 , , ε ij = √ a i − (cid:112) a j − , i, j = 0 , i (cid:54) = j ) ε i = √ a i − √ a + 1 σ z , i = 0 , . (58)The 6 × c -modestate, but the covariance matrix associated with state(36); since c and e only differ by phases, the entangle-ment properties are the same. When considering entan-glement in Sec. V we will therefore use the standard form(58). In the case of a pure three-mode Gaussian state,the three local symplectic invariants a i fully determinethe entanglement content of any given bipartition [68].As we shall see in Section V, the blocks of the covariancematrix σ in the form of expressions (58) will allow us tocompute the amount of bipartite and tripartite entangle-ment.As mentioned in Section IV B, σ i is the covariance ma-trix of the reduced state ρ ( i ) of mode i . Given its diagonalform, one gets from Section IV C that ρ ( i ) is a thermalstate with local mixedness a i . It can also be consideredas a reduced state of a two-mode squeezed state withsqueezing parameter r i [58, 69, 70]. In this respect, thestudy of the reduced state ρ (0) is of particular interestin the context of analogue gravity, since the number ofemitted quanta in the 0 mode gives access to the Hawk-ing radiation spectrum. In the context of general rela-tivity, this spectrum is exactly Planckian [50, 51], with atemperature which is called the ”Hawking temperature”, T H . For the analogue model we consider, dispersive effectsignificantly affect this result. Indeed, if one defines aneffective temperature T eff such that:¯ n = 1exp( (cid:126) ω/T eff ) − , (59)one finds from (48), (55) and (56) that T eff is frequency-dependent: T eff ( ω ) = (cid:126) ω r ( ω )] = (cid:126) ω ln[1 + | S ( ω ) | − ] . (60)Figure 4 represents the corresponding frequency depen-dence of the effective temperature for a waterfall con-figuration with m u = 0 .
59. We note here that similarresults have been obtained numerically in [12]. In thelong wavelength limit ω → V. ENTANGLEMENT IN THREE-MODEGAUSSIAN STATESA. Bipartite entanglement
The criterion usually used to detect entanglement inbipartite systems is the Peres-Horodecki (or PPT) crite- . . . . . . . ¯ h ω/ ( g n u ) . . . . T e ff / ( g n u ) T H FIG. 4. Effective temperature T eff defined in Eq. (60) as afunction of the frequency ω (blue curve). The dashed red lineis the Hawking temperature T H given by (C5). rion [18, 19]. It is a necessary and sufficient separabilitycondition for bipartite 1 × ( N − i | jk or, after tracing out mode k , bipartitions i | j . This criterion states that a state ρ is separable if and only if its partial transpose ρ PT withrespect to the first mode (mode i in the above notation)is positive. Partial transposition of an N -mode Gaussianstate is equivalent to mirror reflection in phase space forthe Wigner function [20]. The covariance matrix of ρ PT is given by σ PT = Λ σ Λ , with Λ = σ z ⊕ N − . (61)According to Eq. (46), the positivity of ρ PT is equivalentto positivity of the symplectic eigenvalues of σ PT , whichwe denote by ν PT j . Thus, the necessary and sufficientseparability criterion becomes ν PT j ≥ , j = 1 , . . . , N. (62)In our case N = 3. Let us investigate bipartite entan-glement of two-mode states obtained by tracing out thethird one. As discussed in section IV B, the covariancematrix associated with the two-mode state i, j obtainedby tracing out mode k is σ ij given by (43). Its symplecticeigenvalues ν ± are given by2 ν ± = ∆ ij ± (cid:113) ∆ ij − σ ij (63)with ∆ ij = det σ i + det σ j + 2 det ε ij [71]. The symplecticeigenvalues ν PT ± of σ PT ij associated with the partial trans-pose are given by2 ( ν PT ± ) = ∆ PT ij ± (cid:113) (∆ PT ij ) − σ ij , (64)with ∆ PT ij = det σ i + det σ j − ε ij . For a two-modestate, the PPT criterion is in fact equivalent to condition ν PT − ≥ ν PT + ≥ ν PT − inour case. Note that, again, since the local mixednessesappearing in Eq. (58) depend on the frequency ω , thelowest symplectic eigenvalue ν PT − also depends on ω .By using the fact that a i ≥ i = 0 , , ν PT − ≥ |
1, independently of the frequency. There-fore, the reduced state of modes 0 | ν PT − ofthe reduced covariance matrices σ and σ are lowerthan 1, which implies that the reduced state of modes0 | | ω , see for instance Fig. 8(a),where the blue curve represents 1 − ν PT − ( ω ) computed forthe reduced state of modes 0 | (cid:40) |(cid:104) ˆ c i ˆ c j (cid:105)| > (cid:104) ˆ c † i ˆ c i (cid:105) (cid:104) ˆ c † j ˆ c j (cid:105) , if i ∈ { , } , j = 2 , (cid:104) ˆ c i ˆ c † j (cid:105)| > (cid:104) ˆ c † i ˆ c i (cid:105) (cid:104) ˆ c † j ˆ c j (cid:105) , if i (cid:54) = j ∈ { , } . (65)Using Eqs. (54) and (58), one finds (cid:104) ˆ c † i ˆ c i (cid:105) (cid:104) ˆ c † j ˆ c j (cid:105) =sinh r i sinh r j , |(cid:104) ˆ c i ˆ c (cid:105)| = sinh r i cosh r and |(cid:104) ˆ c ˆ c † (cid:105)| = sinh r sinh r . Therefore, when consider-ing the bipartition 0 |
1, one concludes immediately thatthe second inequality of (65) is never true; one has in-stead the equality |(cid:104) ˆ c ˆ c † (cid:105)| = (cid:104) ˆ c † ˆ c (cid:105) (cid:104) ˆ c † ˆ c (cid:105) for all fre-quencies ω . Therefore, the reduced state 0 | | |
2, since tanh( r ) < r > | |
2: these statesare always entangled.However, the “Cauchy-Schwarz criterion” does not giveany clue about the amount of entanglement shared byeach bipartition. Indeed, as we will discuss in sectionVI, contrary to the PPT criterion, in an experimentalsetup for which the temperature of the system cannot beexactly equal to zero, a stronger violation of the Cauchy-Schwarz inequality does not necessary imply a greateramount of entanglement.
B. Tripartite entanglement
1. Monogamy inequality
Monogamy is a fundamental property of entanglementcorrelations. It can be described by monogamy inequal-ities, which in the case of a tripartite system with sub- systems labelled by ( i, j, k ) takes the form E ( i | jk ) − E ( i | j ) − E ( i | k ) ≥ E ( A | B ) is a proper measure of bipartite entan-glement between subsystems A and B (nonnegative onseparable states and monotonic under (G)LOCC). Thisinequality expresses the fact that the total amount ofentanglement that can be shared between i and j andbetween i and k is upper bounded by the amount ofentanglement between i and jk taken as a whole. Theleft-hand side of inequality (66) provides a quantifier ofgenuine tripartite entanglement.Not all entanglement measures satisfy a monogamy in-equality. However, it is possible to find and constructproper measures of entanglement which satisfy these re-lations, both in the qubit case and in the continuous-variable case. In the case of qubits, the monogamy in-equality holds for entanglement measured by the squareof the concurrence. For Gaussian states a measure sat-isfying (66) was constructed in [29]; it is called the con-tangle E τ and it corresponds to the squared logarithmicnegativity. For an arbitrary pure state ρ = | ψ (cid:105)(cid:104) ψ | withcovariance matrix σ p it is defined as E τ ( σ p ) = (ln (cid:107) ρ PT (cid:107) ) , (67)where (cid:107) ˆ O (cid:107) = tr (cid:112) ˆ O † ˆ O is the trace norm.The state considered in our case is a pure three-modeGaussian state; thus, any bipartition i | jk is a pure state,for which the term E ( i | jk ) in (66) can be computed easily(see Section V B 2 below). On the contrary, the two otherterms of (66) correspond to reduced two-mode states,which are mixed. The squared logarithmic negativitycan be extended to mixed states by taking the infimumover all convex decompositions of ρ in terms of pure states {| ψ i (cid:105)} . In order to get a quantity more amenable to com-putations, the Gaussian contangle G τ was defined by re-stricting this convex-roof construction to decompositionsover pure Gaussian states only. The Gaussian contanglecan be expressed as G τ ( σ ) = inf σ p ≤ σ E τ ( σ p ) , (68)where the notation σ p ≤ σ means that the matrix σ − σ p is positive semidefinite. It is an upper bound to thetrue contangle E τ obtained from unrestricted pure-statedecompositions, but for pure states both coincide.For three qubits the residual tangle (or three-way tan-gle) E ( i | jk ) − E ( i | j ) − E ( i | k ) provides a measure of tripartiteentanglement. It has an explicit expression [23], which issymmetric in the three qubits. The corresponding quan-tity in the continuous case is no longer symmetric in thethree modes. One can however define a permutation-invariant quantity by minimizing it over all permutationsof the modes [29]. This measure of tripartite entangle-ment shared among Gaussian modes was called residualcontangle [68]. Its explicit expression reads G res τ = G ( i | j | k ) τ = min i,j,k (cid:16) G ( i | jk ) τ − G ( i | j ) τ − G ( i | k ) τ (cid:17) . (69)1
2. Pure-state contangle
Let us consider first a bipartition i | jk . For a pure state,the Gaussian contangle G ( i | jk ) τ coincides with the truecontangle E ( i | jk ) τ . In general, for a multimode Gaussianstate | ψ (cid:105) with covariance matrix σ p ( p stands for pure)and generic bipartition i . . . i N − | i N ( N = 3 in our case),the squared logarithmic negativity can be written as [29] E i ...i N − | i N τ ( σ p ) = (cid:88) j : ν PT j < ln ν PT j , (70)where ν PT j are the symplectic eigenvalues associated withthe partial transpose state ρ PT .It is actually possible to write Eq. (70) in terms ofthe local mixedness a i N associated with mode i N . In-deed, for any covariance matrix σ associated with apure multimode Gaussian state and generic bipartition i . . . i N − | i N , there exists a local symplectic transfor-mation S such that [72] S σ S T = N − ⊕ σ sq , (71)where σ sq is the covariance matrix of a two-modesqueezed state and reads a i N (cid:113) a i N − a i N − (cid:113) a i N − (cid:113) a i N − a i N − (cid:113) a i N − a i N . (72) In the case of a tripartite system ( N = 3), a direct proofof Eqs. (71)–(72), as well as explicit expressions of thesymplectic matrix S for each bipartition , 12 | | |
2, can be found in Appendix D.The symplectic eigenvalues of σ PT (corresponding totaking the partial transpose with respect to mode i N ) arethen readily obtained from (64), using the form (72); thesymplectic eigenvalue 1 has degeneracy 2 ( N − e ± r iN , with twofolddegeneracy. They can be related to the local mixedness a i N through the relations a i N = cosh(2 r i N ). Equation(70) then gives E i ...i N − | i N τ ( σ p ) = arcsinh (cid:16)(cid:113) a i N − (cid:17) = 4 r i N , (73)which only depends on the local mixedness of mode i N and has a simple expression in terms of r i N . We willmake calculations more explicit for our system in thenext section.
3. Residual contangle
Equation (73) provides an explicit expression for thefirst term G ( i | jk ) τ in (69). For a pure three-mode Gaus-sian state, an explicit expression of G ( i | j ) τ and G ( i | k ) τ can also be obtained. Indeed, in this specific case, any re-duced two-mode state saturates the uncertainty relation(42) and belongs to a class of states called Gaussian leastentangled mixed states (GLEMS). For GLEMS, one has[68] G ( i | j ) τ = arcsinh [ (cid:113) m GLEMS ( a i , a j , a k ) − , (74)where m GLEMS can be explicitly calculated as a functionof the three local mixednesses, as shown in Appendix E.In our case (see Eqs. (E24)–(E27)), we obtain G (0 | τ = 0 ,G ( j | τ = arcsinh (cid:18)
21 + a k (cid:113) ( a + 1) ( a j − (cid:19) = arcsinh (cid:32) |(cid:104) ˆ c j ˆ c (cid:105)| (cid:104) ˆ c † k ˆ c k (cid:105) (cid:33) , (75)with j = 0 , k = 1 or j = 1 , k = 0.Let us now introduce the quantity G res( i ) τ = G ( i | jk ) τ − G ( i | j ) τ − G ( i | k ) τ , (76)such that the residual contangle is given by G res τ = min i ∈{ , , } (cid:104) G res( i ) τ (cid:105) . (77)Using (73) and (75), Eq. (76) yields G res(0) τ = arcsinh (cid:18)(cid:113) a − (cid:19) − arcsinh (cid:18)
21 + a (cid:112) ( a + 1) ( a − (cid:19) , (78) G res(1) τ = arcsinh (cid:18)(cid:113) a − (cid:19) − arcsinh (cid:18)
21 + a (cid:112) ( a + 1) ( a − (cid:19) , (79)and G res(2) τ = arcsinh (cid:18)(cid:113) a − (cid:19) − arcsinh (cid:18)
21 + a (cid:112) ( a + 1) ( a − (cid:19) − arcsinh (cid:18)
21 + a (cid:112) ( a + 1) ( a − (cid:19) . (80)The residual contangle only depends on the three localmixednesses a , a and a (this is no longer true at fi-nite temperature, see Sec. VI B and Appendix B). Theminimum over all possible permutations of i, j and k inEq. (77) can be obtained by choosing as reference mode i the one with smallest local mixedness [68].We can then compute the residual Gaussian contanglefor our three-mode Gaussian state using the a i given in2 . . . . . . . ¯ h ω/ ( g n u ) . . . . G r e s τ Ω ω c . . . ¯ h ω/ ( g n u ) − G r e s ( ) τ − G r e s ( ) τ × − ω c Ω FIG. 5. Residual tangles G res(0) τ (blue), G res(1) τ (green), G res(2) τ (red). The frequency ω c is the value above which a becomeslower than a (see Fig. 3), and coincides with the point abovewhich G res(0) τ < G res(1) τ . Inset is the difference G res(0) τ − G res(1) τ (cyan). The upper bound frequency Ω corresponds to thevanishing of the mode 2. (56). The results for m u = 0 .
59 are shown in Fig. 5.Tripartite entanglement naturally emerges from quantumfluctuations around a sonic horizon and diverges whenthe energy goes to zero. This divergence always comesfrom the first term in Eqs. (78), (79) and (80). Indeed,this term diverges as ln ω (see discussion in Sec. V B 4).On the other hand, from the results found in Ref. [14],it may be proven that G ( j | τ given by expressions (75)for j = 0 , m u < j = 0 , k = 1 or j = 1 , k = 0, G ( j | τ = ω → arcsinh (cid:18) | F F j || F k | (cid:19) , (81)where the explicit expressions of the constant coeffi-cients | F i | , i ∈ { , , } are given in Appendix C. Itmeans in particular that the entanglement of bipartitions j | , j ∈ { , } remains finite at zero energy, while the tri-partite entanglement becomes infinite. Then, for higherfrequencies, the residual tangle decreases rapidly to zeroand vanishes at the upper-bound frequency Ω.Moreover, we show in the inset of Fig. 5 (cyan curve)that while at low frequency G res(1) τ < G res(0) τ the situationis reversed for ω > ω c , i.e., when a < a (the differenceis anyway quite small). We note that this result maybe different for Mach numbers different from the value m u = 0 .
59 we consider here. In particular, based on theestimate (82) below, one can show that, when m u < . G res(0) τ becomes the contribution whichminimizes (77). These results have been previously presented in [73].
FIG. 6. Measure of tripartite entanglement G res τ as a functionof the (dimensionless) energy (cid:126) ω/ ( g n u ) and of the upstreamMach number m u ∈ [0 . , .
95] defined in Eq. (6). The pinkcurve corresponds to the upper bound frequency Ω (8). Fora fixed value m u mode 2 only exists for a frequency ω lowerthan Ω( m u ). Beyond this value the tripartite system { , , } no longer exists; this is the reason why the correspondingarea is left blank. The horizontal red line corresponds tothe value m u (cid:39) .
59, corresponding to m d = 2 . (cid:82) G res τ d ω over frequencies ω ∈ [0 , Ω] for each value of m u .The red dot pinpoints the numerical estimate of the integralfor the specific value m u (cid:39) .
59, while the blue dot locatesthe maximum of the black curve reached for m u = 0 .
14. Thelight blue horizontal line on the left graph corresponds to thisvalue.
4. Experimental perspectives
The waterfall model we use has proven to provide afairly good description of the experimental setting [15].In this section we use the relevance of our model to assesswhat is the best choice of parameters for an experimentalmeasure of tripartite entanglement.Figure 6 displays the amount of genuine tripartite en-tanglement G res τ expected in our one-dimensional ana-logue black hole for a given frequency (horizontal axis)and a specific upstream Mach number m u (vertical axis).For most cases, as proved by this two-dimensional graph,the entanglement is indeed shared among the Hawking,the partner and the companion quanta. Therefore, thecompanion plays an important role in the distribution ofentanglement within the emitted quanta.From Fig. 6 one sees that the amount of tripartite en-tanglement is maximal for m u = 0 .
14, in the sense thatthe integral of G res τ ( ω ) over all frequencies is maximal forthis value of upstream Mach number. This specific valueof m u is indicated by the light blue horizontal cut on thegraph. It has been determined by numerical integrationof the residual contangle (77). One can also obtain ananalytic estimate of this value of m u , as we now explain.3From the low-frequency behavior (C1) of the componentsof the S -matrix involved [through Eq. (56)] in (78), (79)and (80), one obtains the following expression for thelow-frequency residual contangle: G res τ ( ω ) (cid:39) ω → min i ∈{ , , } (cid:20) ln (cid:18) | F i | (cid:126) ω/g n u (cid:19)(cid:21) . (82)The value of m u for which G res τ ( ω ) in (82) is the largestis thus simply the value for which the mimimum of | F | , | F | and | F | reaches a maximum. From the analyticexpressions (C2), (C3) and (C4) of these coefficients oneobtains m u = 0 .
17. Although this value has been deter-mined using a different criterion than the numerical es-timate m u = 0 .
14 plotted in Fig. 6 (the former is basedon the low ω behavior and the latter on the integratedsignal) the fact that both values are quite close confirmstheir relevance. C. Entanglement localization
The tripartite entanglement of our system can be con-centrated in a two-mode state by applying a local lin-ear Bogoliubov transformation [74, 75]; this is calledentanglement localization. This transformation can beobtained by means of the symplectic transformation S given by (71). To this mapping (71) between σ and σ (cid:48) = ⊕ σ sq one can associate the Bogoliubov transfor-mation T = U † S U [see Eq. (45)]. The modes e definedin (33) (which coincide with the modes c up to a phase)are mapped through this Bogoliubov transformation tonew modes f . The Bogoliubov transformation from e to f = ( ˆ f , ˆ f , ˆ f , ˆ f † , ˆ f † , ˆ f † ) T is denoted T e → f , and thus wehave f = T e → f e . This transformation is such that thetripartite entanglement e | e | e gets completely localizedin a two-mode squeezed state.Let us consider in turn the different cases. If we con-sider bipartitions ij | k = 12 | | e ,as derived explicitly in Appendix D [see in particularEq. (D22)], the new operators ˆ f i and ˆ f correspond to amixing of annihilation and creation operators ˆ e i , ˆ e andˆ e † i , ˆ e † . In the case of bipartition ij | k = 01 | e ,entanglement can also be localized but without mixingannihilation and creation operators. The correspond-ing Bogoliubov transformation is given by Eq. (D25)and corresponds to a change of basis from { ˆ e , ˆ e , ˆ e } to { ˆ f , ˆ f , ˆ f } given byˆ f = − sin θ ˆ e + cos θ ˆ e , (83a)ˆ f = cos θ ˆ e + sin θ ˆ e , (83b)ˆ f = ˆ e , (83c)where (see Appendix D 4)cos θ = sinh r sinh r and sin θ = sinh r sinh r . (84) Parametric DownConversion e e f f e f = Pump Beamsplitter
Partner
Analogue black hole
CompanionHawking
FIG. 7. Schematic representation of an optical process equiv-alent to the Hawking emission in the transsonic BEC sys-tem we consider. Entanglement is localized in the two-modesqueezed state f | f . The mode f being empty is representedby a dashed line. The transformation leading to entanglement localizationis thus particularly simple in the case of bipartition 01 | | (cid:105) b = T | (cid:105) f , T = e r ( ˆ f † ˆ f † − ˆ f ˆ f ) . (85)The operator T is a two-mode squeezing operator of theform (21) between ˆ f and ˆ f , with squeezing parame-ter r defined in (55). Note that the modes ˆ e and ˆ e that are combined in (83a) and (83b) are those of pos-itive norm (that is, the ones corresponding to positive-frequency branches in Fig. 2); this leads to a squeezedstate between the only mode of negative norm (mode 2)and a combination of the modes of positive norm (modes0 and 1), exactly as occurs in the gravitational case [76].To summarize, the tripartite entanglement in our sys-tem can be unitarily localized by linearly combiningmodes ˆ e and ˆ e as in Eqs. (83a)–(83b) to obtain modeˆ f , which forms a two-mode squeezed state with ˆ f . Be-sides, using the definition (84) and noticing that (cid:104) ˆ e † ˆ e (cid:105) = |(cid:104) ˆ c † ˆ c (cid:105)| one obtains (cid:104) ˆ f † ˆ f (cid:105) = sin θ (cid:104) ˆ c † ˆ c (cid:105) + cos θ (cid:104) ˆ c † ˆ c (cid:105)− θ cos θ |(cid:104) ˆ c † ˆ c (cid:105)| = (sin θ | S | − cos θ | S | ) = 0 . (86)This means that mode f is not occupied. This comesas no surprise since the corresponding local mixednessis equal to 1 in the transformed covariance matrix givenby (71), which entails from (55) that the mean particlenumber is equal to 0. One can thus schematically de-scribe the Bogoliubov transformation (83) operating inour analogue black hole by means of the equivalent opti-cal setup represented in Fig. 7: non-degenerate paramet-ric down-conversion in a nonlinear crystal creates a two-mode squeezed state. One of the modes is the partnerˆ f = ˆ e . The other one, ˆ f , is directed to a beam-splitter4that generates the two other outgoing channels ˆ e andˆ e , respectively the companion and the Hawking mode.We note that various theoretical proposals and experi-mental works have addressed the issue of generating andmeasuring tripartite entangled states for continuous vari-ables, based on different setups of nonlinear optical para-metric oscillators [75, 77–79]. The analogue black hole weconsider here is another such setup. It is quite peculiar inthe sense that genuine tripatite entanglement is realizedalthough two of the outgoing modes (0 and 1) are notentangled. VI. FINITE TEMPERATURE
We previously considered the zero-temperature case,where all the averages (cid:104)· · · (cid:105) in Sec. IV A are taken overthe vacuum state | (cid:105) b . In the present section we study afinite temperature system. A. Finite-temperature states
Because of the existence of negative-energy modes, thetranssonic flow we consider is energetically unstable andcannot support a thermal state. However, one can definea finite temperature configuration [12, 13, 36] as follows:One considers a BEC initially flowing at a constant ve-locity V u (with uniform density n u ), and at thermal equi-librium at temperature T BEC in the frame moving alongwith the fluid. Then the potential U ( x ) of Eq. (2) isslowly ramped up until the system reaches the configura-tion described in Sec. II A and Fig. 1. At the end of thisadiabatic branching process one can define an occupa-tion number ¯ n i ( ω, T BEC ) for each of the incoming modesˆ b i . As explained, e.g., in [36], for a fixed frequency ω these occupation numbers are given by (cid:104) ˆ b † i ( ω )ˆ b i ( ω ) (cid:105) = ¯ n i ( ω, T BEC ) = n th (cid:2) ω B ,α (cid:0) q i | in ( ω ) (cid:1)(cid:3) , (87)where n th ( (cid:36) ) = (exp( (cid:36)/T BEC ) − − is the thermal Boseoccupation distribution. In this expression ω B ,α ( q ) is theBogoliubov dispersion relation (5), with α = u if i = 0and α = d if i = 1 or 2, and the functions q i | in ( ω ) aredefined below Eq. (8).The regime in which the separation (1) between a clas-sical field and quantum fluctuations is valid and wherethe Bogoliubov treatment of the fluctuations applies hasbeen denoted as the “weakly interacting quasiconden-sate regime” in [80]. It is valid up to a temperature T BEC (cid:39) g n u [15], where g is the coefficient of the nonlin-earity in the Gross-Pitaevskii equation (2). For typicalexperimental parameters g n u (cid:39) T BEC , the vacuum state | (cid:105) b isreplaced by a product of thermal states of b modes givenby ρ a b = ⊗ i =0 ρ th ( a bi ) , a bi = 1 + 2¯ n i ( ω, T BEC ) (88)where ¯ n i is given by Eq. (87). We recall that we usethe term ”thermal” to designate that state in a loosesense, since, as explained in the beginning of this section,the occupation numbers (87) do not correspond to anequilibrium distribution in the transsonic configurationwe consider. The covariance matrix associated with thisstate is given by σ th = diag( a b , a b , a b , a b , a b , a b ). Afterthe Bogoliubov transformation c = T b , the covariancematrix becomes σ (cid:48) th = S T σ th S T T (see Eq. (45)). The2 × σ i and ε ij in the block decomposition (39)of σ (cid:48) th are given by expressions (50) and (51) where theaverages (cid:104) . . . (cid:105) should be replaced by (cid:104) . . . (cid:105) th = Tr [ ρ a b . . . ] . (89)In particular, σ i = a i, th , (90)where a i, th = 1 + 2 (cid:104) ˆ c † i ˆ c i (cid:105) th , i ∈ { , , } , (91)is the corresponding local mixedness [compare to (50)and to the first of Eqs. (58)].We conclude this short section by noting that the opti-cal analogue proposed in Fig. 7 remains relevant at finitetemperature. The difference with the zero-temperaturecase is just the occupation number of the f -modes: theynow acquire an incoherent contribution. In particular theoccupation (cid:104) ˆ f +0 ˆ f (cid:105) is no longer zero as in Eq. (86). Thissuggests a possible experimental study of the effects oftemperature on tripartite entanglement: one could real-ize the optical setup of Fig. 7, send a non-coherent beamalong the mode f , and evaluate the associated effect onentanglement in the system. B. Detection of entanglement
Contrary to the zero-temperature case, σ (cid:48) th is associ-ated with a mixed state with no special symmetry, andit cannot be put in a standard form where the matrices ε ij are all diagonal [68]. In this section we thus restrictour study to bipartite entanglement. In this case, the4 × ij can always be brought by LLUBOs to itsstandard form [67]. One easily proves that matrices ε ij have the same form as those in the zero-temperature case,namely ε ij = 2 |(cid:104) ˆ c i ˆ c † j (cid:105) th | , i, j = 0 , , i (cid:54) = jε i = 2 |(cid:104) ˆ c i ˆ c (cid:105) th | σ z , i = 0 , . (92)5As a consequence, the lowest symplectic eigenvalue as-sociated with the partial-transposed reduced two-modestate ij takes the same form as in the zero-temperaturecase. Eq. (64) still holds, and in particular2 ( ν PT − ) = ∆ PT ij − (cid:113) (∆ PT ij ) − σ ij , (93)with heredet σ = (cid:16) a , th a , th − |(cid:104) ˆ c ˆ c † (cid:105) th | (cid:17) , ∆ PT = a , th + a , th − |(cid:104) ˆ c ˆ c † (cid:105) th | , det σ i = (cid:0) a i, th a , th − |(cid:104) ˆ c i ˆ c (cid:105) th | (cid:1) , i = 0 , PT i = a i, th + a , th + 8 |(cid:104) ˆ c i ˆ c (cid:105) th | , i = 0 , . (94)Note that the above expressions only involve moduli ofmean values, so that we could equivalently use opera-tors ˆ e i instead of ˆ c i since the transformation defined byEqs. (34) and (35) is diagonal. The explicit form of thequantities appearing in Eqs. (94) is given in AppendixB. At variance with the zero-temperature case they donot depend only on the local mixednesses. They shouldbe experimentally accessible through the measurement ofthe structure form factor and of real space density corre-lations [81], meaning that the PPT criterion can be usedto experimentally detect entanglement.The PPT criterion asserts that the bipartite state isentangled iff 1 − ν PT − > . (95)In the following we denote this quantity as the ”PPTmeasure”. It is of particular interest to focus on the bi-partition 0 | ω of the elementary excitations and for different tempera-tures ranging from 0 to 1.5 g n u (blue curves). In eachplot the dashed blue curves display the same quantity atzero-temperature for comparison.It is instructive to compare the conclusions drawn fromthe study of the PPT measure with those obtained us-ing the criterion of violation of the Cauchy-Schwarz in-equality (65). According to this criterion, the analogueHawking–partner pair 0 | CS ≡ |(cid:104) ˆ c ˆ c (cid:105) th | − ( a , th −
1) ( a , th − > . (96)In the following we denote ∆ CS as the ”Cauchy-Schwarzparameter”. It is represented by the red curves in Figs. 8which confirm the results obtained with the PPT cri-terion: the blue and red curves are positive in the sameregion and cross zero exactly at the same frequency. Thismeans, as expected, that both criteria lead to the samequalitative result for entanglement detection. However,as we shall see in Sec. VI C, they lead to different quan-titative estimation of the amount of entanglement. The analogue Hawking pair is entangled in the rangeof frequencies for which inequalities (95) and (96) hold.This corresponds to the grey shaded regions bounded bytwo vertical black dot-dashed lines in Figs. 8. The rangeof parameters over which entanglement can be observeddecreases when the temperature of the Bose gas increases.As shown in Fig. 8(f) thermal effects destroy the entan-glement between the Hawking pair for T BEC (cid:38) . g n u .Therefore the temperature of the experimental systemshould not exceed this limiting value to be able to ob-serve entanglement.Figure 9(a) represents the PPT measure 1 − ν PT − ofthe Hawking pair 0 | T BEC = 0 . g n u fordifferent configurations parameterized by the upstreamMach number m u . As already seen in Fig. 8, which cor-responds to the specific case m u = 0 .
59, a finite temper-ature reduces the range of frequencies for which entan-glement occurs. One observes in this new plot that theentanglement of the Hawking pair persists for a largerfraction of the available frequency domain when the pa-rameter m u is closer to unity. This is in agreement withthe results obtained in [33]; it was noticed that not onlythe temperature T BEC destroys the entanglement of theanalogue Hawking pair, but also that a strong “coupling”of mode 1 with the other modes can affect their entan-glement. This coupling is measured through the squaredmodulus of the scattering matrix coefficients | S ( ω ) | and | S ( ω ) | . One finds numerically (and analyticallyin the low- ω sector [14]) that these two quantities de-crease when m u increases. This exactly corresponds tothe results presented in Fig 9(a): when m u increases,the coupling between 0-1 and 1-2 decreases, and indeedleads to a stronger violation of PPT criterion for a largerfraction of frequencies. However, it is important to notethat the previous discussion is only valid at low enoughtemperatures. This is illustrated in Fig. 9(b): for a tem-perature as large as T = 1 . g n u the region where thepair is entangled greatly diminishes and entanglementonly survives at moderate values of m u (at variance withthe conclusion of the above discussion). Likewise, at thistemperature, even in the region where entanglement ispresent, the PPT measure is significantly lower than inthe equivalent regions in Fig. 9(a).It is also interesting to study the entanglement of theHawking pair, not as a function of the absolute temper-ature, but as a function of the Hawking temperature T H whose m u -dependent expression is reminded in Eq. (C5).The Hawking temperature, which has been defined in sec-tion IV D, is a measure of the intensity of the one-bodyHawking signal: there is no obvious reason why entangle-ment between modes should disappear when T BEC > T H .This is indeed what is observed in Fig. 10: entangle-ment persists in sizeable regions even when T BEC = 10 T H .We can conclude from the above discussion that whereasentanglement persists for temperatures noticeably largerthan T H , it is significantly reduced when T BEC becomeslarger than the chemical potential g n u .6 . . . . . . . . . . . . . (a) T BEC = 0 G ( | ) τ , ∆ C S , − ν P T − . . . . . . . . . . . . . (b) T BEC = 0.3 g n u . . . . . . . . . . . . . (c) T BEC = 0.5 g n u G ( | ) τ , ∆ C S , − ν P T − . . . . . . − . . . . (d) T BEC = 1.0 g n u . . . . . . ¯ h ω/ ( g n u ) − . − . . . (e) T BEC = 1.5 g n u G ( | ) τ , ∆ C S , − ν P T − . . . . . . ¯ h ω/ ( g n u ) − . − . . . (f ) T BEC = 1.8 g n u FIG. 8. Evolution of the PPT measure 1 − ν PT − (blue), of the Cauchy-Schwarz parameter ∆ CS (red) and of the Gaussian contangle G (0 | τ (green) for the bipartite system 0 | (cid:126) ω/ ( g n u )and for different temperatures of the system, denoted by T BEC , ranging from 0 (a) to 1 . g n u (f). All the plots are obtainedfor an upstream Mach number m u = 0 .
59. The dashed blue curves in graphs (b)-(f) correspond to the zero-temperature valueof 1 − ν PT − . The gray areas indicate the range of frequencies for which the bipartite system is entangled [see text and Eqs. (95)and (96)]. The purple dots locate the upper-bound frequency Ω at which mode 2 vanishes. (a)(b) FIG. 9. PPT measure 1 − ν PT − of the Hawking pair 0 | m u and of thefrequency ω , for temperatures (a) T BEC = 0 . g n u and (b) T BEC = 1 . g n u . The pink curve corresponds to the upper-bound frequency Ω (8). For a fixed value m u , i.e., along ahorizontal cut on the graph, mode 2 only exists for a frequency ω lower than Ω( m u ) (see Fig. 2). The dashed black curvecorresponds to 1 − ν PT − = 0 and thus delimits the region wherethe analogue Hawking pair is entangled. C. Measurement of entanglement
The violation of the Cauchy-Schwarz inequality is oftenused to study the entanglement between the elementaryexcitations in the context of analogue gravity [31–33, 35–37]. However, while this criterion tells us whether thebipartite system is entangled or not, the Cauchy-Schwarzparameter ∆ CS it is not a good measure of the amountof entanglement at finite temperature.To clarify this point, we compute the amount of en-tanglement at finite temperature for the bipartition 0 | G (0 | τ defined inEq. (68). This computation is slightly more difficult herethan in the zero-temperature case, where we obtained FIG. 10. Same as Fig. 9 for T BEC = 10 T H . Eq. (75). In the presence of temperature the reducedtwo-mode state 0 | G (0 | τ = arcsinh (cid:26)(cid:113) min θ [ m ( θ )] − (cid:27) , (97)where m ( θ ) is explicitly given by Eq. (E23). In Figs. 8we represent by a green solid line the value of G (0 | τ inthe range of frequencies for which the system is entangled(the minimum over the angle θ in Eq. (97) is obtained nu-merically). The results for G (0 | τ confirm that the PPTand Cauchy-Schwarz criteria correctly determine the re-gion where entanglement exists.As expected, entanglement decreases as the tempera-ture increases. In the zero-temperature case [Fig. 8(a)],both 1 − ν PT − and ∆ CS vary in the same way as G (0 | τ .The situation at finite temperature is different: whilethe quantities 1 − ν PT − and G (0 | τ appear to behave in thesame way, being increasing and decreasing in the sameregions and having a maximum at the same value of ω ,this is not the case for ∆ CS whose maximum is shiftedwith respect to the two others, see Figs. 8(b)-(f).In order to illustrate this phenomenon, in Fig. 11 werepresent ∆ CS plotted as a function of G (0 | τ for sev-eral temperatures. These are parametric curves obtainedfrom expressions (96) and (97), ω playing the role of theparameter. Except at T BEC = 0, ∆ CS is not a monotonousfunction of G (0 | τ , as demonstrated by the closed loopswith regions of negative slope observed for each finitetemperature. Another way to note the same point is toremark that the maximal violation of Cauchy-Schwarzinequality (∆ CS maximal) is not reached when G (0 | τ ismaximal. This confirms that the parameter ∆ CS is notan entanglement monotone.8 . . . . . . . G (0 | τ . . . . . ∆ C S T BEC = 0 T BEC = 0 . T BEC = 0 . FIG. 11. Evolution of ∆ CS given by Eq. (96) as a func-tion of the measure of bipartite entanglement G (0 | τ givenby expressions (97) and (E23), for the same set of temper-atures as in Fig. 8, ranging from T BEC = 0 (blue curve) to T BEC = 1 . g n u (red curve), with m u = 0 .
59. When possible,the corresponding temperature for each curve is indicated onthe graph (we dropped the factor g n u for readability). Notethat for T BEC >
0, the curves describe a loop.
In Fig. 12 we underline the difference between the be-haviors of the Cauchy-Schwarz parameter and the PPTmeasure by plotting 1 − ν PT − as a function of G (0 | τ . The . . . . . . . G (0 | τ . . . . . . − ν P T − T BEC = 0 T BEC = 0 . T BEC = 0 . T BEC = 0 . T BEC = 1
FIG. 12. Evolution of 1 − ν PT − given by Eq. (93) as a func-tion of the measure of bipartite entanglement G (0 | τ givenby expressions (97) and (E23), for the same set of temper-atures as in Fig. 8, ranging from T BEC = 0 (blue curve) to T BEC = 1 . g n u (red curve), with m u = 0 .
59. When possible,the corresponding temperature for each curve is indicated onthe graph (we dropped the factor g n u for readability). Foreach temperature the common maximal value of 1 − ν PT − and G (0 | τ is marked with a point. difference with Fig. 11 is striking. For each temperature,1 − ν PT − is a monotonous increasing function of G (0 | τ . Itis not easily seen in the figure, but for finite T BEC the relation between the two quantities is not one to one: foreach G (0 | τ there are two (close) values of 1 − ν PT − whichcoalesce at the common maximum of the two quantities,marked with a point on Fig. 12. This confirms withoutambiguity that the PPT measure is still an entanglementmonotone at finite temperatures. We also note that allthe curves in Fig. 12 almost superimpose, meaning thatrelation between the two quantities 1 − ν PT − and G (0 | τ isvery weakly dependent on temperature, which makes thePPT measure an even better candidate for quantifyingentanglement.To conclude this section, we would like to insist on thepositive aspects of using PPT measure in future experi-mental studies of analogue black hole configurations: (i)the calculations of Sec. VI B show that the computationof the lowest symplectic eigenvalue requires essentiallythe knowledge of the the same quantities (94) as ∆ CS [compare Eqs. (93) and (96)]. It means that the valueof ν PT − is experimentally accessible and can be measuredfrom the density correlations along the acoustic blackhole in the same manner as discussed in Section VI B;(ii) contrary to Cauchy-Schwarz criterion, it is possible toassociate to the PPT criterion a ’good’ measure of entan-glement, whatever the temperature of the system is; (iii) ν PT − has a much simpler expression than G (0 | τ in termsof the local mixednesses. As a last remark, we empha-size that the use of the symplectic spectrum to constructmeasures of entanglement has already been proposed in[38–40] for several analogue models. We confirm in thispaper its relevance for quantifying entanglement in thetranssonic flow of a BEC. VII. CONCLUSION
In the present work, we have investigated entanglementproperties of Gaussian modes emitted from an analogueblack hole realized in the flow of a Bose-Einstein conden-sate. These Gaussian modes result from the scattering ofelementary excitations onto the black hole horizon, whichinduces emission of an analogue Hawking radiation.Gaussian states are entirely characterized by their firstand second moments. Thus, their entanglement proper-ties can be expressed in terms of their covariance ma-trix. Using results developed in the field of continuous-variable entanglement, we have characterized bipartiteentanglement in our system via the Gaussian contan-gle; for tripartite entanglement we used the Gaussianresidual entanglement defined from monogamy inequal-ities. We quantified the bipartite entanglement at zeroas well as at finite temperature. We quantified tripar-tite entanglement at zero temperature, and identifiedthe best configuration for its experimental observation: G res is larger for waterfall configurations with moder-ate upstream Mach number ( m u (cid:39) .
15) and at smallfrequencies. We illustrated the entanglement propertiesof the system by means of an equivalent simple optical9setup. We also showed that, quite counter-intuitively,while there is no bipartite entanglement between two ofthe outgoing modes (tracing out the third one), there isnevertheless genuine tripartite entanglement between thethree modes.Experimentally, various quantum optics setups can re-alize tripartite entangled states. The quantities we con-sidered here should be experimentally accessible from themeasure of density correlations [9, 10, 81]. From the pointof view of quantum information, our work underlines theusefulness of results on continuous-variable entanglementrestricted to Gaussian states, as those naturally emergefrom our setting. The advantage of working with Gaus-sian states is of course manifold: these states can be de-scribed in terms of their covariance matrix only; in thecase of mixed states, the convex roof construction re-stricted to decompositions onto Gaussian pure states al-lows to derive analytic expressions for the entanglementmeasure, and thus to characterize the entanglement con-tent of the various partitions of the state.Extensions of the present work include the investiga-tion of zero-norm modes, which were shown in [15] toplay an important role in the correct quantum descrip-tion of Bogoliubov excitations . The study of the influ-ence of thermal effects on the amount of tripartite en-tanglement also constitutes a natural continuation of thepresent study.During the completion of this work, we became awareof the preprint [82] which studies tripartite entanglementin an analogue system thanks to the residual contangle,as done in the present paper. The model in [82] corre-sponds to a ”subluminal” dispersion relation, whereas inthe BEC case we consider, the dispersion is rather ”su-perluminal”.
ACKNOWLEDGMENTS
We thank M. Jacquet for inspiring discussions on en-tanglement in analogue gravity. We also thank A. As-pect, T. Bienaim´e, Q. Glorieux and C. Westbrook forfruitful exchanges.
Appendix A: Bogoliubov transformations
In this Appendix we detail some intermediate stepsuseful for establishing the results presented in Sec. III A.The form (12) and (15) of vectors b and c implies thatthe 2 N × N matrix T defining the unitary Bogoliubovtransformation (14) has a block structure given by (16).In order that the ˆ c i defined in Eq. (14) satisfy bosoniccommutation relations, the matrix T must verify T (cid:101) J T T = (cid:101) J , (A1)where (cid:101) J is defined in Eq. (13). Eq. (A1) means that T belongs to the symplectic group Sp(2 N, C ). As a conse-quence, one has T − = − (cid:101) J T T (cid:101) J = (cid:18) α T β † β T α † (cid:19) . (A2)Condition (A1) can be reexpressed in terms of the two N × N matrices α and β as α α † − β β † = N , α β T − β α T = 0 α † α − β T β ∗ = N , α T β ∗ − β † α = 0 . (A3)The matrix T being symplectic, it can be written as T = e (cid:101) J Q (A4)with Q a 2 N × N symmetric matrix.The unitary operator T relating operators ˆ c i and ˆ b i according to ˆ c i = T † ˆ b i T, (A5)is defined as T = e b T Q b , (A6)as can be shown by using the Baker-Campbell-Hausdorffformula [45]. Note that using Eqs. (14) and (A4) one has b T Q b = c T ( T − ) T Q T − c = c T Q c .It is possible to show [49, 56, 83] that T can be uniquely decomposed into the product T = (det α ) − / exp N (cid:88) i,j =1 X ij ˆ c † i ˆ c † j exp N (cid:88) i,j =1 Y ij ˆ c † i ˆ c j exp N (cid:88) i,j =1 Z ij ˆ c i ˆ c j , (A7)where X , Y , Z are N × N matrices defined by X = − β ∗ α − , e − Y T = α, Z = α − β. (A8)The interest of the decomposition (A7) lies in the fact that all annihilation operators have been put to the right.Therefore, when applied to the vacuum | (cid:105) c , T only acts through matrix X . This directly yields Eq. (20).0 Appendix B: Explicit expression of the covariance matrix
In this Appendix we give a useful formula for the covariance matrix, present explicit expressions necessary forevaluating its finite-temperature form, and discuss their zero-temperature limit.The decomposition (32) makes it possible to write the covariance matrix σ = S T S T T of Sec. IV D under the form σ = v v v cos( φ ) − v v sin( φ ) 2 v v cos( φ ) 2 v v sin( φ )0 1 + 2 v v v sin( φ ) 2 v v cos( φ ) 2 v v sin( φ ) − v v cos( φ )2 v v cos( φ ) 2 v v sin( φ ) 1 + 2 v v v cos( φ ) 2 v v sin( φ ) − v v sin( φ ) 2 v v cos( φ ) 0 1 + 2 v v v sin( φ ) − v v cos( φ )2 v v cos( φ ) 2 v v sin( φ ) 2 v v cos( φ ) 2 v v sin( φ ) − v v v sin( φ ) − v v cos( φ ) 2 v v sin( φ ) − v v cos( φ ) 0 − v , (B1)where v ij and ϕ ij are defined in Eq. (32) and φ ij = ϕ i − ϕ j .Also, for explicitly computing the finite temperature entanglement properties studied in Sec. VI B [see Eqs. (94)and (96)] one uses the formulae: (cid:104) ˆ c ˆ c † (cid:105) th = S S ∗ (1 + ¯ n ) + S S ∗ (1 + ¯ n ) + S S ∗ ¯ n , (cid:104) ˆ c † i ˆ c i (cid:105) th = | S i | ¯ n + | S i | ¯ n + | S i | (1 + ¯ n ) , i = 0 , , (cid:104) ˆ c i ˆ c (cid:105) th = S i S ∗ (1 + ¯ n ) + S i S ∗ (1 + ¯ n ) + S i S ∗ ¯ n , i = 0 , , (cid:104) ˆ c † ˆ c (cid:105) th = | S | (1 + ¯ n ) + | S | (1 + ¯ n ) + | S | ¯ n , (B2)where the quantities ¯ n ¯ n and ¯ n are defined in Eq. (87), and, as in Eq. (B1), we do not write the explicit ω dependences for legibility.At zero temperature the above equations reduce to (cid:104) ˆ c ˆ c † (cid:105) = S S ∗ + S S ∗ = S ∗ S , (cid:104) ˆ c † i ˆ c i (cid:105) = | S i | , i = 0 , , (cid:104) ˆ c i ˆ c (cid:105) = S i S ∗ + S i S ∗ = S i S ∗ , i = 0 , , (cid:104) ˆ c † ˆ c (cid:105) = | S | + | S | = − | S | , (B3)where use has been made of property (10). UsingEqs. (B3) and expressions (56), one may show that thefinite-temperature components (90) and (92) of the co-variance matrix reduce at T BEC = 0 to the form (58) asthey should. However, at finite temperature, Eq. (B2)holds instead of (B3), implying that, contrarily to thezero-temperature case, the covariance matrix, its sym-plectic eigenvalues, and thus the entanglement propertiesof the system, do not depend only on the local mixed-nesses.
Appendix C: Long wavelength limit of the scatteringamplitudes
In the long wavelength limit, the S i coefficients of the S -matrix (11) behave as S i ( ω ) = F i (cid:114) gn u (cid:126) ω + O ( ω / ) , i ∈ { , , } , (C1)where the F i are dimensionless constant coefficients. Forthe waterfall configuration we consider here, analytic ex- pressions of their moduli have been determined in [14]: | F | = 2 m u (1 − m u ) (1 + m u ) (1 + m u ) (1 + m u + m u ) , (C2) | F | = 12 (1 − m u ) (1 + m u ) (1 + m u ) (1 + m u + m u ) , (C3)and | F | = 12 (1 − m u ) (1 + m u + m u ) , (C4)where m u is the upstream Mach number. From the low-frequency behavior of the scattering coefficients it hasalso been possible in [14] to evaluate the analogue Hawk-ing temperature of the waterfall configuration. The resultreads T H gn u = 12 (1 − m u ) (2 + m u )(1 + 2 m u ) . (C5) Appendix D: Entanglement localization in atripartite system
In this Appendix we present the specifics of the processof entanglement localization discussed in Sec. V C. Let σ be a covariance matrix associated with a pure three-modeGaussian state. We want to determine the explicit formof the symplectic matrix S which transforms σ accordingto S σ S T = ⊕ σ sq , (D1)1where σ sq is the covariance matrix of a two-modesqueezed state [see Eq. (72)].
1. General form of the symplectic matrix
Consider a bipartition ij | k . The covariance matrix as-sociated with the subsystem k is denoted as σ k and theone associated with subsystem ij reads σ ij = (cid:18) σ i ε ij ε T ij σ j (cid:19) . (D2)The whole covariance matrix associated with the tripar-tite system is then σ = σ i ε ij ε ik ε T ij σ j ε jk ε T ik ε T jk σ k . (D3)Consider the case where the covariance matrix is in itsstandard form (58), i.e. σ k = a k and either σ ij = a i c a i − cc a j − c a j , c = √ a i − (cid:112) a j + 1 (D4)for bipartitions ij | k = 02 | | ε i = cσ z )or σ ij = a i c a i cc a j c a j , c = √ a i − (cid:112) a j − | ε = c ). The dif-ference in the sign in front of c between (D4) and (D5)is actually of great importance and leads to two differenttypes of symplectic transformations in Eq. (D1). Notethat we also impose a i < a j in (D4); in fact, the order ofthe local mixednesses does not matter in (D5), as shallbe clear at the end of this section.The symplectic eigenvalues σ k = a k are ν k = a k .Using Williamson theorem, we can bring σ ij to a diagonalmatrix( σ ij ) (cid:48) = S ij σ ij ( S ij ) T = diag { ν i , ν i , ν j , ν j } , (D6)where we ordered the symplectic eigenvalues such that ν i < ν j . Easy calculations lead to S ij = a b a ηbηb − a b − a , (D7)with a = − (cid:115) a j ν j − a i ν i ν j − ν i , b = (cid:115) a i ν j − a j ν i ν j − ν i . (D8) and η = − ij | k = 02 | |
0, and η = 1 for bipartition 01 |
2. The coefficients a and b satisfythe identity a + ηb = a i + ηa j ν i + ην j = 1 . (D9)The last equality is valid only if S ij is a symplectic ma-trix.Expressions (D7) and (D8) are valid for any covariancematrix σ ij of the form (D4)-(D5). In our case, we canfurther simplify these expressions using the purity con-straint of the three-mode Gaussian state under consider-ation. Indeed, one can easily prove that for any reducedtwo-mode states ij of a pure three-mode Gaussian state,∆ ij = det σ ij + 1 = det σ k + 1 [68]. Therefore, consid-ering the reduced state jk , Eq. (63) immediately gives ν i = 1 and ν j = (cid:112) det σ ij = a k , which imply from thelast equality of (D9) that a i + ηa j = ν i + ην j = 1 + ηa k .This expression is true for the case (D4) iff a i < a j , be-cause η = −
1. For (D5), the order is not importantbecause η = 1. Thus, (D8) simplifies to a = − (cid:115) ( a j − η ) ( a k + η ) a k − ,b = (cid:115) ( a i −
1) ( a k + η ) a k − . (D10)
2. Standard form
The symplectic matrix defined by S = S ij ⊕ S k , (D11)with S k = a k , transforms the covariance matrix (D3)to σ (cid:48) = S σ S T = (cid:18) σ (cid:48) ij KK T σ (cid:48) k (cid:19) , (D12)with K some matrix and σ (cid:48) ij = a k
00 0 0 a k , σ (cid:48) k = (cid:18) a k a k (cid:19) . (D13)Following [72], we first notice that − ( J σ ) = , (D14)where we recall that J is given by Eq. (38). The previ-ous expression is not difficult to prove: the Williamsontheorem ensures the existence of a symplectic matrix O mapping the covariance matrix σ to the identity (for apure state all the symplectic eigenvalues are equal toone); thus, one obtains − ( J σ ) = − J O O T J O O T = .Note that one also has − ( J σ (cid:48) ) = . Then, inserting2expression (D12) in Eq. (D14) and using the fact that σ (cid:48) ij and σ (cid:48) k are diagonal, it is easy to prove the followingconditions for the matrix K: ( σ (cid:48) ij ) − J ij K J k K T = ( σ (cid:48) k ) − J k K T J ij K = − σ (cid:48) ij K + J ij K J k σ (cid:48) k = 0 , (D15)where J k is defined in Eq. (38) and J ij = J i ⊕ J j . Thelast condition in (D15) implies that a given coefficientK mn (cid:54) = 0 iff ( σ (cid:48) ij ) mm = ( σ (cid:48) k ) nn , i.e., if and only if thesymplectic eigenvalue on the row m of the subsystem ij matches with the one on the column n of the subsystem k . Therefore, one can rewrite expression (D12) in theform σ (cid:48) = (cid:18) (cid:101) σ (cid:19) , with (cid:101) σ = (cid:18) a k (cid:101) K (cid:101) K T a k (cid:19) , (D16)where we introduced a new 2 × (cid:101) K. Then, notic-ing that − ( J ij (cid:101) σ ) = , one obtains (cid:40) (cid:101) K J k (cid:101) K T = (1 − a k ) J k , J k (cid:101) K J k = (cid:101) K . (D17)The above conditions lead to (cid:101) K = (cid:18) a √ λ − a −√ λ − a − a (cid:19) , (D18)where λ = (cid:112) a k −
1. Then, given that σ is in its standardform [meaning that ε ik and ε jk are diagonal matrices, seeEq. (D3)] and remembering that S ij is given by expres-sion (D7), one can see easily that (cid:101) K must be diagonal;therefore, a = λ . As a result, one finds that (cid:101) σ = (cid:18) a k (cid:112) a k − σ z (cid:112) a k − σ z a k (cid:19) , (D19)which exactly corresponds to the covariance matrix of asqueezed state with squeezing parameter r k >
0, withcosh(2 r k ) = a k (see Eq. (72)). This statement endsthe proof: the symplectic matrix S given by expression(D11), with S ij explicitly written in equations (D7) and(D8) and S k = , indeed lead to the transformation(D1).
3. Bipartitions | and | For bipartitions ij | k = 02 | | e , thesymplectic transformation (D11) involves the matrix S ij (with j = 2) given by Eq. (D7) with η = −
1. Using theidentity (D9), we introduce a parameter γ , such that a =cosh γ and b = sinh γ , where a and b are the coefficientsof the S ij matrix. In this case, one finds S i = − cosh γ γ − cosh γ − sinh γ − sinh γ γ
00 sinh γ γ , (D20) with cosh γ = cosh r / cosh r k , sinh γ = sinh r i / cosh r k ,computed from expressions (D10).To this symplectic transformation S σ S T at the level ofthe covariance matrix corresponds a Bogoliubov trans-formation f = T e → f e , with T e → f = U † S U, (D21)(see Eqs. (44) and (45)). Using the explicit expression ofthe symplectic matrix (D20), one finds ˆ f i = − cosh γ ˆ e i + sinh γ ˆ e † j ˆ f j = − sinh γ ˆ e † i + cosh γ ˆ e j ˆ f k = ˆ e k . (D22)Therefore, entanglement can be localized in the subsys-tem f | f k with k = 0 or 1 only through a Bogoliubovtransformation which mixes annihilation and creation op-erators.
4. Bipartition | For the bipartition ij | k = 01 | S is givenby Eq. (D7) with η = 1. One finds S = − sin θ θ − sin θ θ cos θ θ
00 cos θ θ , (D23)with cos θ = sinh r / sinh r , sin θ = sinh r / sinh r , us-ing again the identity (D9) and expressions (D10). Theassociated Bogoliubov transformation T e → f = U † ( S ⊕ S ) U (D24)leads to the new set of operators ˆ f = − sin θ ˆ e + cos θ ˆ e ˆ f = cos θ ˆ e + sin θ ˆ e ˆ f = ˆ e , (D25)where, as in the previous subsection, f and f are newcombinations of modes e and e , and f = e . Here,there is no mixing of annihilation and creation operatorsand the matrix S ⊕ S is unitary. Appendix E: Computation of the finite-temperatureGaussian contangle
In this Appendix we explain how to obtain expression(E23) used in Eq. (97) for evaluating the Gaussian con-tangle at finite temperature. We could not find a deriva-tion of this formula in the literature, and since the explicitform given in [28] has some missprints, we find it usefulto give the whole proof, following the same path as in[28]. For a general (mixed or pure) two-mode Gaussian3state, a measure of bipartite entanglement is given bythe Gaussian contangle G τ ( σ ) defined in Eq. (68). It hasbeen proven in [28] that finding the infinimum over pureGaussian states amounts to minimize m ( x , x , x ) = 1 + x det Γ , (E1)with det Γ = x − x − x , where x , x , and x mustbelong to the following cones x = a + b − (cid:115) ( x − c + ) + (cid:18) x − a − b (cid:19) ,x = a + b d + (cid:115)(cid:16) x + c − d (cid:17) + (cid:18) x + a − b d (cid:19) , (E2)where a , b , c + and c − are the coefficients of the covariancematrix σ written in the standard form and associatedwith a given two-mode Gaussian state: σ = a c + a c − c + b c − b . (E3)In Eqs. (E2), d = a b − c − . The minimum of expression(E1) is located at the intersection of both cones (E2)[28]. Therefore, in the following, we aim at finding thisintersection, which corresponds to an ellipse. To findthe equation of this ellipse, we first make a change ofcoordinates (Lorentz boost): x (cid:48) = γ ( x − v x ) ,x (cid:48) = γ ( x − v x ) ,x (cid:48) = x , (E4)with v = a − ba + b d + 1 d − < d > ,γ = ( a + b ) ( d − (cid:112) ( a d − b ) ( b d − a ) . (E5)We find after simplifications: (cid:40) ( x (cid:48) − α ) − ( x (cid:48) − β ) − ( x (cid:48) − γ ) = 0 , ( x (cid:48) − α ) − ( x (cid:48) − β ) − ( x (cid:48) − γ ) = 0 , (E6)with α = γ ( a − b )2 (cid:18) a + ba − b − v (cid:19) , β = c + ,α = γ ( a − b )2 d (cid:18) a + ba − b + v (cid:19) , β = − c − d , (E7)and γ = γ = − γ ( a − b ) d − . (E8) Note that α and α simplify to α = 2 a b d − a − b (cid:112) ( a d − b ) ( b d − a ) ,α = d (cid:0) a + b (cid:1) − a b d (cid:112) ( a d − b ) ( b d − a ) . (E9)Let us now make another change of variables x (cid:48)(cid:48) = x (cid:48) − L + ,x (cid:48)(cid:48) = x (cid:48) − H + ,x (cid:48)(cid:48) = x (cid:48) − γ = x (cid:48) − γ , (E10)with L + = α + α a b ( d − d (cid:112) ( a d − b ) ( b d − a ) ,H + = β + β c + d − c − d . (E11)This leads to (cid:40) ( x (cid:48)(cid:48) − L − ) − ( x (cid:48)(cid:48) − H − ) − x (cid:48)(cid:48) = 0 , ( x (cid:48)(cid:48) + L − ) − ( x (cid:48)(cid:48) + H − ) − x (cid:48)(cid:48) = 0 , (E12)with L − = α − α (cid:112) ( a d − b ) ( b d − a )2 d ,H − = β − β c + d + c − d . (E13)Looking at Eqs. (E12), one sees that the changes of co-ordinates (E4) and (E10) make it possible to eliminateone variable ( x (cid:48)(cid:48) ) and to symmetrise the equations. Notethat both cone tops belong to the plane x (cid:48)(cid:48) = 0.The intersection of the cones (E12) is now simple tofind. By combining equations (E12), one can first elimi-nate x (cid:48)(cid:48) to find the relation between x (cid:48)(cid:48) and x (cid:48)(cid:48) : x (cid:48)(cid:48) = H − L − x (cid:48)(cid:48) . (E14)Inserting the previous relation in one of the equations(E12) yields (cid:18) − H − L − (cid:19) x (cid:48)(cid:48) + x (cid:48)(cid:48) = L − − H − , (E15)which exactly corresponds to the equation of an ellipse.Let us define the angle θ such that x (cid:48)(cid:48) = H − cos θ, x (cid:48)(cid:48) = L − cos θ,x (cid:48)(cid:48) = (cid:113) L − − H − sin θ. (E16)At this stage, we have everything to express Eq. (E1)only in terms of the parameter θ and coefficients of thecovariance matrix. Since the Lorentz boost preserves the4relations between both cones, one can find the minimumof the function m in the basis ( x (cid:48) , x (cid:48) , x (cid:48) ), that is to say m = 1 + x (cid:48) x (cid:48) − x (cid:48) − x (cid:48) . (E17)Using Eqs. (E10), (E11), (E13) and (E16), one finds x (cid:48) = H − cos θ + L + ,x (cid:48) = L − cos θ + H + = 12 d [ c + d − c − + (cid:112) ( a d − b ) ( b d − a ) cos θ (cid:105) ,x (cid:48) = (cid:113) L − − H − sin θ + γ . (E18)This gives x (cid:48) − x (cid:48) − x (cid:48) = α α − β β − γ − L − H + − H − L + ) cos θ − γ (cid:113) L − − H − sin θ. (E19)After some simplifications, the first right-hand side termof Eq. (E19) reads α α − β β − γ = a + b + 2 c − c + d . (E20) Expanding the coefficient of − cos θ in the second right-hand side term of Eq. (E19) leads to2 ( L − H + − H − L + ) = α β − β α = (cid:8) a b c − + ( a + b ) c + c − + c − (cid:2) a (1 − b ) + b (cid:3) − a b c + ( a + b − (cid:9) × (cid:104) d (cid:112) ( a d − b ) ( b d − a ) (cid:105) − . (E21)The coefficient of − sin θ in the last right-hand side termof Eq. (E19) reads2 γ (cid:113) L − − H − = − a − b d × (cid:115) − ( c + d + c − ) ( a d − b ) ( b d − a ) . (E22)The last step consists of inserting expressions (E20),(E21) and (E22) in Eq. (E19); then, Eq. (E19) in ex-pression (E17). This leads to the final result m ( θ ) = 1 + 12 d (cid:104)(cid:112) ( a d − b ) ( b d − a ) cos θ + c + d − c − (cid:105) × (cid:40) ( a + b + 2 c − c + ) − cos θ a b c − + ( a + b ) c + c − + c − (cid:2) a (1 − b ) + b (cid:3) − a b c + ( a + b − (cid:112) ( a d − b ) ( b d − a )+ ( a − b ) sin θ (cid:115) − ( c + d + c − ) ( a d − b ) ( b d − a ) (cid:41) − . (E23)The explicit expressions of a , b , c + , c − and d = a b − c − used for evaluating G (0 | τ in Eq. (97) are a = a , th , b = a , th , c + = 2 |(cid:104) ˆ c ˆ c (cid:105) th | , c − = − c + and d = a , th a , th − |(cid:104) ˆ c ˆ c (cid:105) th | .Let us now find the explicit expressions of the Gaussiancontangles G ( j | τ given by (74), i.e., to compute G ( j | τ = arcsinh (cid:26)(cid:113) min θ [ m ( θ )] − (cid:27) , (E24)where m ( θ ) is given by Eq. (E23) and j = 0 ,
1. First, bynoticing that any reduced two mode state of a pure threemode Gaussian state belongs to the class of GLEMS [68],expression (E23) simplifies to [28] m GLEMS ( θ ) = 1 + ( A cos θ + B ) d [( g −
1) cos θ + g + 1] , (E25) where g = √ det σ , with σ given by (E3), A = c + d + c − and B = c + d − c − . Using our notations and theexplicit expression of the covariance matrix written inthe standard form (58), for a given bipartition j |
2, onehas a = a , b = a j , c + = − c − = (cid:112) a j − √ a + 1, d = g = a k ; we recall that j = 0 or 1 and that the remaining(third) mode (1 or 0) is denoted as k . One proves inthis case that the minimum over θ in expression (E25) isreached when [28]cos θ (cid:63) = − a k . (E26)Inserting this expression in Eq. (E25) leads to m GLEMS ( θ (cid:63) ) = (cid:18) − a j + a k a k (cid:19) . (E27)5Using this result in Eq. (E24) yields immediately expres- sions (75), remembering that a j + a k = a + 1. [1] W. G. Unruh, Experimental black-hole evaporation?,Phys. Rev. Lett. , 1351 (1981).[2] T. Torres, S. Patrick, A. Coutant, M. Richartz, E. W.Tedford, and S. Weinfurtner, Rotational superradiantscattering in a vortex flow, Nat. Phys. , 833 (2017).[3] P. Chen and G. Mourou, Accelerating plasma mirrors toinvestigate the black hole information loss paradox, Phys.Rev. Lett. , 045001 (2017).[4] S. Liberati, G. Tricella, and A. Trombettoni, The infor-mation loss problem: An analogue gravity perspective,Entropy , 940 (2019).[5] M. J. Jacquet, S. Weinfurtner, and F. K¨onig, The nextgeneration of analogue gravity experiments, Phil. Trans.R. Soc. A. , 20190239 (2020).[6] L. J. Garay, J. R. Anglin, J. I. Cirac, and P. Zoller, Sonicanalog of gravitational black holes in Bose-Einstein con-densates, Phys. Rev. Lett. , 4643 (2000).[7] O. Lahav, A. Itah, A. Blumkin, C. Gordon, S. Rinott,A. Zayats, and J. Steinhauer, Realization of a sonic blackhole analog in a Bose-Einstein condensate, Phys. Rev.Lett. , 240401 (2010).[8] J. Steinhauer, Observation of self-amplifying Hawking ra-diation in an analogue black-hole laser, Nat. Phys. ,864 (2014).[9] J. Steinhauer, Observation of quantum Hawking radia-tion and its entanglement in an analogue black hole, Nat.Phys. , 959 (2016).[10] J. R. M. de Nova, K. Golubkov, V. I. Kolobov, andJ. Steinhauer, Observation of thermal Hawking radiationand its temperature in an analogue black hole, Nature(London) , 688 (2019).[11] U. Leonhardt, T. Kiss, and P. ¨Ohberg, Theory of elemen-tary excitations in unstable Bose-Einstein condensatesand the instability of sonic horizons, Phys. Rev. A ,033602 (2003).[12] J. Macher and R. Parentani, Black-hole radiation inBose-Einstein condensates, Phys. Rev. A , 043601(2009).[13] A. Recati, N. Pavloff, and I. Carusotto, Bogoliubov the-ory of acoustic Hawking radiation in Bose-Einstein con-densates, Phys. Rev. A , 043603 (2009).[14] P.-E. Larr´e, A. Recati, I. Carusotto, and N. Pavloff,Quantum fluctuations around black hole horizons inBose-Einstein condensates, Phys. Rev. A , 013621(2012).[15] M. Isoard and N. Pavloff, Departing from thermality ofanalogue Hawking radiation in a Bose-Einstein conden-sate, Phys. Rev. Lett. , 060401 (2020).[16] R. Horodecki, P. Horodecki, M. Horodecki, andK. Horodecki, Quantum entanglement, Rev. Mod. Phys. , 865 (2009).[17] R. F. Werner, Quantum states with Einstein-Podolsky-Rosen correlations admitting a hidden-variable model,Phys. Rev. A , 4277 (1989).[18] A. Peres, Separability criterion for density matrices,Phys. Rev. Lett. , 1413 (1996).[19] M. Horodecki, P. Horodecki, and R. Horodecki, Separa-bility of mixed states: necessary and sufficient conditions, Phys. Lett. A , 1 (1996).[20] R. Simon, Peres-Horodecki separability criterion for con-tinuous variable systems, Phys. Rev. Lett. , 2726(2000).[21] S. Hill and W. K. Wootters, Entanglement of a pair ofquantum bits, Phys. Rev. Lett. , 5022 (1997).[22] C. H. Bennett, D. P. DiVincenzo, J. A. Smolin, and W. K.Wootters, Mixed-state entanglement and quantum errorcorrection, Phys. Rev. A , 3824 (1996).[23] V. Coffman, J. Kundu, and W. K. Wootters, Distributedentanglement, Phys. Rev. A , 052306 (2000).[24] T. J. Osborne and F. Verstraete, General monogamyinequality for bipartite qubit entanglement, Phys. Rev.Lett. , 220503 (2006).[25] G. Giedke, B. Kraus, M. Lewenstein, and J. I. Cirac, Sep-arability properties of three-mode gaussian states, Phys.Rev. A , 052303 (2001).[26] G. Vidal and R. F. Werner, Computable measure of en-tanglement, Phys. Rev. A , 032314 (2002).[27] M. M. Wolf, G. Giedke, O. Kr¨uger, R. F. Werner, andJ. I. Cirac, Gaussian entanglement of formation, Phys.Rev. A , 052320 (2004).[28] G. Adesso and F. Illuminati, Gaussian measures of entan-glement versus negativities: Ordering of two-mode Gaus-sian states, Phys. Rev. A , 032334 (2005).[29] G. Adesso and F. Illuminati, Continuous variable tan-gle, monogamy inequality, and entanglement sharing inGaussian states of continuous variable systems, New J.Phys. , 15 (2006).[30] T. Hiroshima, G. Adesso, and F. Illuminati, Monogamyinequality for distributed gaussian entanglement, Phys.Rev. Lett. , 050503 (2007).[31] J. R. M. de Nova, F. Sols, and I. Zapata, Violation ofCauchy-Schwarz inequalities by spontaneous Hawking ra-diation in resonant boson structures, Phys. Rev. A ,043808 (2014).[32] X. Busch, I. Carusotto, and R. Parentani, Spectrum andentanglement of phonons in quantum fluids of light, Phys.Rev. A , 043819 (2014).[33] X. Busch and R. Parentani, Quantum entanglement inanalogue Hawking radiation: When is the final state non-separable?, Phys. Rev. D , 105024 (2014).[34] S. Finazzi and I. Carusotto, Entangled phonons in atomicBose-Einstein condensates, Phys. Rev. A , 033607(2014).[35] D. Boiron, A. Fabbri, P.-E. Larr´e, N. Pavloff, C. I. West-brook, and P. Zi´n, Quantum signature of analog Hawk-ing radiation in momentum space, Phys. Rev. Lett. ,025301 (2015).[36] A. Fabbri and N. Pavloff, Momentum correlations as sig-nature of sonic Hawking radiation in Bose-Einstein con-densates, SciPost Phys. , 019 (2018).[37] A. Coutant and S. Weinfurtner, Low-frequency analogueHawking radiation: The Bogoliubov-de Gennes model,Phys. Rev. D , 025006 (2018).[38] B. Horstmann, R. Sch¨utzhold, B. Reznik, S. Fagnocchi,and J. I. Cirac, Hawking radiation on an ion ring in thequantum regime, New J. Phys. , 045008 (2011).[39] D. E. Bruschi, N. Friis, I. Fuentes, and S. Weinfurtner, On the robustness of entanglement in analogue gravitysystems, New J. Phys. , 113016 (2013).[40] M. J. Jacquet and F. Koenig, The influence of spacetimecurvature on quantum emission in optical analogues togravity, SciPost Phys. Core , 5 (2020).[41] R. Balbinot, A. Fabbri, S. Fagnocchi, A. Recati, andI. Carusotto, Nonlocal density correlations as a signatureof Hawking radiation from acoustic black holes, Phys.Rev. A , 021603 (2008).[42] I. Carusotto, S. Fagnocchi, A. Recati, R. Balbinot, andA. Fabbri, Numerical observation of Hawking radiationfrom acoustic black holes in atomic Bose–Einstein con-densates, New J. Phys. , 103001 (2008).[43] I. Zapata, M. Albert, R. Parentani, and F. Sols, ResonantHawking radiation in Bose–Einstein condensates, New J.Phys. , 063048 (2011).[44] L. P. Pitaevskii and S. Stringari, Bose-Einstein Conden-sation and Superfluidity , International Series of Mono-graphs on Physics (Oxford University Press, Oxford,United Kingdom, 2016).[45] J.-P. Blaizot and G. Ripka,
Quantum Theory of FiniteSystems (MIT Press, Cambridge, Mass, 1986).[46] A. L. Fetter, Theory of a dilute low-temperature trappedBose condensate, in
Bose-Einstein Condensation inAtomic Gases, Proceedings of the International School“Enrico Fermi”, Course CXL , edited by M. Inguscio,S. Stringari, and C. E. Wieman (IOS Press, Amsterdam,1999) pp. 201–263.[47] C. Barcel´o, L. J. Garay, and G. Jannes, Two faces ofquantum sound, Phys. Rev. D , 044042 (2010).[48] S. J. Robertson, The theory of Hawking radiation in lab-oratory analogues, J. Phys. B: At. Mol. Opt. Phys. ,163001 (2012).[49] R. Balian and E. Br´ezin, Nonunitary Bogoliubov trans-formations and extension of Wick’s theorem, Nuovo Ci-mento B , 37 (1969).[50] S. W. Hawking, Black hole explosions?, Nature (London) , 30 (1974).[51] S. W. Hawking, Particle creation by black holes, Com-mun. Math. Phys. , 199 (1975).[52] B. S. DeWitt, Quantum field theory in curved spacetime,Phys. Rep. , 295 (1975).[53] P. C. W. Davies and S. A. Fulling, Radiation from movingmirrors and from black holes, Proc. R. Soc. Lond. A ,237 (1977).[54] R. M. Wald, On particle creation by black holes, Com-mun. Math. Phys. , 9 (1975).[55] L. Parker, Probability distribution of particles created bya black hole, Phys. Rev. D , 1519 (1975).[56] K. Takayanagi, Utilizing group property of Bogoliubovtransformation, Nucl. Phys. A , 17 (2008).[57] B. L. Schumaker and C. M. Caves, New formalism fortwo-photon quantum optics. II. Mathematical foundationand compact notation, Phys. Rev. A , 3093 (1985).[58] G. S. Agarwal, Quantum optics (Cambridge UniversityPress, 2012).[59] C. Weedbrook, S. Pirandola, R. Garc´ıa-Patr´on, N. J.Cerf, T. C. Ralph, J. H. Shapiro, and S. Lloyd, Gaussianquantum information, Rev. Mod. Phys. , 621 (2012).[60] R. Simon, E. C. G. Sudarshan, and N. Mukunda,Gaussian-Wigner distributions in quantum mechanicsand optics, Phys. Rev. A , 3868 (1987).[61] R. Simon, N. Mukunda, and B. Dutta, Quantum-noisematrix for multimode systems: U( n ) invariance, squeez- ing, and normal forms, Phys. Rev. A , 1567 (1994).[62] G. K. Giedke, Quantum information and continuous vari-able systems , Ph.D. thesis, Universit¨at Linz Bibliothek,4040 Linz-Auhof (2001).[63] J. Williamson, On the algebraic problem concerning thenormal forms of linear dynamical systems, Amer. J.Math. , 141 (1936).[64] A. Serafini, Multimode uncertainty relations and separa-bility of continuous variable states, Phys. Rev. Lett. ,110402 (2006).[65] G. Adesso, A. Serafini, and F. Illuminati, Quantificationand scaling of multipartite entanglement in continuousvariable systems, Phys. Rev. Lett. , 220504 (2004).[66] G. Adesso, A. Serafini, and F. Illuminati, Extremal en-tanglement and mixedness in continuous variable sys-tems, Phys. Rev. A , 022318 (2004).[67] L.-M. Duan, G. Giedke, J. I. Cirac, and P. Zoller, Insep-arability criterion for continuous variable systems, Phys.Rev. Lett. , 2722 (2000).[68] G. Adesso, A. Serafini, and F. Illuminati, Multipartite en-tanglement in three-mode gaussian states of continuous-variable systems: Quantification, sharing structure, anddecoherence, Phys. Rev. A , 032345 (2006).[69] K. Mølmer, A. Perrin, V. Krachmalnicoff, V. Leung,D. Boiron, A. Aspect, and C. I. Westbrook, HanburyBrown and Twiss correlations in atoms scattered fromcolliding condensates, Phys. Rev. A , 033601 (2008).[70] M. Perrier, Z. Amodjee, P. Dussarrat, A. Dareau, A. As-pect, M. Cheneau, D. Boiron, and C. I. Westbrook, Ther-mal counting statistics in an atomic two-mode squeezedvacuum state, SciPost Phys. , 2 (2019).[71] A. Serafini, F. Illuminati, and S. D. Siena, Symplectic in-variants, entropic measures and correlations of Gaussianstates, J. Phys. B: At. Mol. Opt. Phys. , L21 (2004).[72] A. Botero and B. Reznik, Modewise entanglement ofgaussian states, Phys. Rev. A , 052311 (2003).[73] M. Isoard, Theoretical study of quantum correlations andnonlinear fluctuations in quantum gases , Ph.D. thesis,Universit´e Paris-Saclay (2020).[74] A. Serafini, G. Adesso, and F. Illuminati, Unitarily lo-calizable entanglement of Gaussian states, Phys. Rev. A , 032349 (2005).[75] D. Daems, F. Bernard, N. J. Cerf, and M. I. Kolobov,Tripartite entanglement in parametric down-conversionwith spatially structured pump, J. Opt. Soc. Am. B ,447 (2010).[76] N. D. Birrell and P. C. W. Davies, Quantum Fields inCurved Space (Cambridge University Press, 1982).[77] T. Aoki, N. Takei, H. Yonezawa, K. Wakui, T. Hiraoka,A. Furusawa, and P. van Loock, Experimental creationof a fully inseparable tripartite continuous-variable state,Phys. Rev. Lett. , 080404 (2003).[78] A. S. Villar, M. Martinelli, C. Fabre, and P. Nussen-zveig, Direct production of tripartite pump-signal-idlerentanglement in the above-threshold optical parametricoscillator, Phys. Rev. Lett. , 140504 (2006).[79] A. S. Coelho, F. A. S. Barbosa, K. N. Cassemiro, A. S.Villar, M. Martinelli, and P. Nussenzveig, Three-colorentanglement, Science , 823 (2009).[80] P. Deuar, A. G. Sykes, D. M. Gangardt, M. J. Davis,P. D. Drummond, and K. V. Kheruntsyan, Nonlocal paircorrelations in the one-dimensional Bose gas at finitetemperature, Phys. Rev. A , 043619 (2009).[81] J. Steinhauer, Measuring the entanglement of analogue Hawking radiation by the density-density correlationfunction, Phys. Rev. D , 024043 (2015).[82] Y. Nambu and Y. Osawa, Tripartite entanglementof Hawking radiation in dispersive model (2021), arXiv:2101.11764 [gr-qc].[83] X. Ma and W. Rhodes, Multimode squeeze operators andsqueezed states, Phys. Rev. A41