Some measures for fermionic entanglement in the cosmological de Sitter spacetime

We investigate two measures of quantum correlations and entanglement, namely the violation of the Bell-Mermin-Klyshko inequalities and the quantum discord, for Dirac fermions in the cosmological de Sitter background of dimension four. The BMK violation is focussed on the vacuum for two and four mode squeezed states and it is shown to increase with the number of modes. For the quantum discord, we investigate a maximally entangled in-vacuum state. Qualitative similarities as well as differences of our results with that of different coordinatisations of de Sitter in the context of scalar and fermionic field theory is discussed.

The relativistic sector, where particle pair creation can take place is always interesting in this context, for such pairs turn out to be entangled. The most well studied case of quantum entanglement in this context corresponds to the Rindler left-right wedges (or the maximally extended near horizon geometry of a nonextremal black hole) [21,22,23,24,25,26,27,28,29,30,31] (also references therein). Due to pair creation, the entanglement or quantum correlation between two maximally entangled Bell pair gets degraded in the Rindler frame, e.g. [26]. See e.g. [32] and references therein for a discussion on entanglement in the context of Schwinger pair creation. We further refer our reader to [33] and references therein for a discussion on entanglement in the soft sector of quantum electrodynamics.
Another well motivated sector of the relativistic entanglement is the cosmological scenario, important chiefly because such study might provide us insight about the initial state as well as the geometry of the early inflationary universe. Such studies might involve investigations like the vacuum entanglement entropy, various other measures of quantum entanglement between initially entangled states for both bosonic and fermionic fields as well as the quantum decoherence of cosmological perturbations and their possible observational consequences, e.g. [34,35,36,37,38,39,40,41,42,43,44,45,46,47,48,49,50,51] and references therein.
In this work we wish to investigate two measures -namely the violation of the Bell-Mermin-Klyshko inequality [2,3,10,11,12] and the quantum discord [16,17] for Dirac fermions in the cosmological de Sitter background.
The Bell inequality [2] (see also [20] and references therein) is a measure of non-locality for a two-partite quantum system. Later such inequality was extended to multipartite systems [3,10,11,12], altogether regarded as the Bell-Mermin-Klyshko (BMK) inequalities. In the nonlocal regime of quantum mechanics BMK inequalities may be violated, thereby clearly distinguishing quantum effects from that of any local classical hidden variables. As the number of partite is increased in a system, the upper bound of the BMK violation also increases, e.g. [13,14]. Given two subsystems, on the other hand, quantum discord is a suitable measure of all correlations including entanglement between them [16,17]. Accordingly, even if there is no entanglement for a mixed state, the quantum discord can be non-vanishing. The key ingredient of the computation of discord is the quantum mutual information between the subsystems. One also needs to optimise over all possible measurements performed on one of the subsystems.
The basic computational tools we shall use in this paper can be seen in [41,42] and references therein. Our study on the BMK violation will be focused on the vacuum, whereas for the quantum discord we shall work on some maximally entangled initial Bell state, which will give us insight about correlations between non-vacuum states. In [41], the quantum discord corresponding to a maximally entangled state for two scalar fields was investigated in the hyperbolic de Sitter background. In [42], the infinite BMK violation was demonstrated for a massless scalar field in a cosmological background which is de Sitter and radiation dominated respectively in the past and future. We shall compute these two measures for massive Dirac fermions in the cosmological de Sitter background in order to see how much similar or dissimilar the result is with the already existing results.
The rest of the paper is organised as follows. In the next section we construct the relevant two and four mode squeezed states using the Bogoliubov relations discussed in Appendix A. We compute the BMK violation for the two and four mode squeezed states in Section 3. Computation of the discord can be seen in Section 4. Finally we discuss the results and conclude in Section 5.
We work in spacetime dimension four and set c = 1 = .

Fermionic squeezed states
Based upon the discussion of Appendix A, we shall construct below the two-and four-mode fermionic squeezed states, to be useful for our purpose. Corresponding to the field quantisations Eq. (53), Eq. (58), we define the 'in' and 'out' vacua as, The Bogoliubov relations of Eq. (59) show that the 'in' vacuum can be expressed as a squeezed state over all the 'out' states, where |0 ± k,s, out respectively represent a particle and an antiparticle vacuum.
We shall work here with a specific value of the spatial momentum, k, e.g. [41]. We also note from Eq. (59) that the helicities do not mix in the Bogoliubov transformations. Thus due to the various anti-commutation relations, the squeezed state expansion corresponding to different s values in Eq. (1) will just factor out. This permits us to go for another simplification -to restrict ourselves to a specific s value as well. In other words, we shall work with a subspace of |0 in corresponding to specific k and s. Thus instead of Eq. (1), we work with (after normalisation) where we have suppressed the index s, since we are restricting ourselves to any single value of it. We shall further comment on the more general helicity summed state at the end of Section 3.2. The above is called a two-mode squeezed state.
The notion of the two-mode squeezed state can easily be extended if we include more than one fermionic fields, say Ψ 1 (x), Ψ 2 (x), . . . , each quantised in a way described in Appendix A and further mix these particle species via some Bogoliubov transformations (see [42] for discussions on scalar field theory, also [52,53]). Let us consider two fermionic fields, Ψ 1 and Ψ 2 with their 'in' vacuum |0 in 1 and |0 in 2 respectively, where a in,1 , b in,1 and a in,2 , b in,2 are the annihilation operators corresponding to Ψ 1 and Ψ 2 respectively. The combined 'in' vacuum for these two field system is then |0 in 1 ⊗ |0 in 2 . Let us suppose that these two fields are correlated via a simple mixing transformation as, in ( k, s) = Γ k a in,2 ( k, s) where Γ k and ∆ k are Bogoliubov coefficients satisfying, |Γ k | 2 + |∆ k | 2 = 1. It is easy to check that the operators defined above satisfy the canonical anti-commutation relations. We assume that such squeezing between different field species is weak, i.e., ∆ k Γ k 1 Let us denote the vacuum state corresponding to the new operators in Eq. (4) by |0 , From Eq. (4), |0 can then be expanded as, We focus as earlier on a specific k and s value. Suppressing the index s, making the expansion of the exponential only up to the second order owing to the weak squuezing, and after normalising, Eq. (6) takes the form where A k and B k depend upon Γ k and ∆ k and |A k | 2 + |B k | 2 = 1.
The 'in' states appearing on the right hand side of Eq. (7) can further be expanded according to Eq. (2). If we take the rest mass of both the fields to be the same, the Bogoliubov coefficients (α k , β k ) corresponding to these two fields are then same as well, cf. Appendix A. This yields, known as the four mode squeezed state. Note that the above construction is similar but not exactly the same as the de Sitter α-vacua (e.g. [55]). This is because the latter mixes the modes of a single quantum field whereas here we have mixed two different fields themselves. It is also clear that the construction of such squeezed states goes beyond two fields. For example, with three fermionic fields one can construct a six mode squeezed state.
Having constructed the necessary states, we are now ready to go into computing the BMK violation and the quantum discord.

The BMK inequalities
The construction of the Bell and the Bell-Mermin-Klyshko (BMK) operators for fermions is similar to that of the scalar field theory [42] (also references therein). Let us consider two sets of non-commuting observables defined respectively over the Hilbert spaces H X and H Y : {X, X ∈ H X } and {Y, Y ∈ H Y }. We assume that these are spin-1/2 operators along specific directions, such as X = n i σ i , X = n i σ i , where σ i 's are the Pauli matrices and n i , n i are unit vectors on the three dimensional Euclidean space. The eigenvalues of these operators are ±1. The Bell operator B ∈ H X ⊗ H Y is defined as In theories with local classical hidden variables we have the so called Bell's inequality, B 2 ≤ 1 and | B | ≤ 1 [3]. However, this is violated in quantum mechanics. Indeed, from Eq. (9) we have (suppressing the tensor product sign), where I is the identity operator. Using the commutation relations for the Pauli matrices, one gets B 2 ≤ 2 or | B | ≤ √ 2, obtaining a violation of Bell's inequality, where the upper bound is known as the maximum violation [4].
The above construction can be extended to multipartite systems, known as the Mermin-Klyshko or collectively the Bell-Mermin-Klyshko inequalities. The relevant operator is recursively defined as, where we have defined All the operators correspond to spin-1/2 systems as earlier.
In classical hidden variable theories we have O n = ±O n , yielding the Mermin-Klyshko inequalities whereas in quantum mechanics one has [12,13,14], Thus the BMK inequality will be violated for multipartite states, n ≥ 2. We shall mention another relevant form of the BMK inequality at the end of the next subsection.

Computing the violation
Let us first compute the BMK violation for the two-mode squeezed state defined in Section 2. Following [9,42], we introduce a pseudospin operator S, where n ≡ (sin θ cos φ, sin θ sin φ, cos θ) is a spatial unit vector, S ± are the ladder operators and n.S has eigenvalues ±1. We have the operations over the orthonormal states |0 , |1 (corresponding to a fix value of the helicity s), Since Eq. (14) is defined on an Euclidean plane, we can use its rotational invariance to set φ = 0 in Eq. (14), so that the expectation value of the Bell operator, Eq. (11), in the two-mode squeezed state (Eq. (2)) becomes where S, in Eq. (11) and F (θ 1 , θ 2 ) is given by The other F 's can be found by replacing θ 1 , θ 2 suitably above. The violation of the BMK inequality can be examined for various values of the angles as well as the Bogoliubov coefficients. For example for θ 1 = 0, θ 1 = π/2 and θ 2 = −θ 2 , we obtain which is plotted in Fig. 1. In particular, the maximum violation B 2 ∼ √ 2, is achieved for the maximum value of the Bogoliubov coefficient, |β k | ∼ 0.707 (cf., Appendix A). Note also from Eq. (59) that the variation of the Bogoliubov coefficient correspond to the variation of the rest mass only. Thus for a given particle species, the value of the Bell violation does not change with respect to the momentum k.
As Eq. (13) indicates, the upper bound of the Bell violation can be increased by going to four or higher mode squeezed states. We now wish to demonstrate such BMK violation for a four-mode squeezed state discussed in Section 2. Using Eq. (11), the relevant BMK operator can be written as where O i = n i · S and O i = n i · S, for i = 1, 2, 3, 4 and n i 's are unit spacelike vectors on the Euclidean 3-plane as earlier. We shall compute the expectation value of the above operator in the state Eq. (8). Denoting the first operator appearing within the square bracket on the right hand side of the above equation by E 4 , we find We obtain after some algebra, where we have denoted for the sake of brevity, The expectation values corresponding to the other operators in Eq. (19) can be found by suitably permuting their angular arguments in Eq. (20), Eq. (21) and Eq. (22). Fig. 2 depicts the expectation value of Eq. (19) with respect to the angular variable θ. 1 We have taken A k ∼ 0.975 and B k ∼ 0.224 in the expression of the four-mode squeezed state, Eq. (8). Comparison with Fig. 1 shows that the BMK violation has increased, in agreement with Eq. (13).
Let us now recall that for a general many particle quantum state which can be grouped into entangled and non-entangled parts, the BMK violation can also be written in a more illuminating form [15], where N represents the total number of partite states and where K 1 denotes the number of single partite states which are not entangled with other (N − 1)-partite states. L = M l=2 K l , with K l being the number of groups consisting of l entangled partite states and M stands for the largest number of such states in a group. It follows that N = M l=1 lK l . Eq. (23) indicates that the upper bound of the BMK violation increases with the number of modes in the squeezed states discussed in Section 2 as follows. For a two-mode squeezed state, we have N = 2, L = 1, K 1 = 0, so that Z = 1. For a four mode squeezed state on the other hand, we have N = 4, L = 1, K 1 = 0 and hence Z = 3. Following [42] as of the scalar field, it is then easy to argue that if we include more than one k value into the state, for the two-mode squeezed state the BMK violation does not increase any further, but for the four mode squeezed state it increases without bound, eventually leading to an infinite violation of the BMK inequality. We shall not go into any detail of this. Finally, we note that we have worked with states with a single value of the helicity s, Eq. (2), Eq. (8). It is easy to see that the above results will remain unchanged even if we work with a squeezed state where the s values are summed over, as follows. As we have argued below Eq. (1), any squeezed state expansion that sums over s, will factor out between two normalised sub-sectors corresponding to those s values. On the other hand, since the pseudospin operator Eq. (14) acts on states only with a specific s value, it is clear that the bra and ket for the squeezed state expansion corresponding to the other s value will just give unity, while computing the expectation values like Eq. (16).
Finally, we note the qualitative similarity between the BMK violations for a scalar field [42] with that of the Dirac fermions. The BMK violation discussed above basically probes the vacuum state of the cosmological de Sitter spacetime. We also wish to discuss in this work some correlation properties (both classical and quantum) associated with the maximally entangled states. The quantum discord is one such suitable measure, which we investigate below using the two-mode squeezed state, Eq. (2).

Quantum discord
Let us first define the quantum discord [16,17]. In classical information theory we have the mutual information between two random variables X and Y , defined as where are the Shannon entropies with respective probabilities P (X) and P (Y ), and is the joint Shannon entropy with the joint probability P (X, Y ) for both the variables X and Y . One can relate the joint probability to the conditional probability as where P (Y |X) is the probability of Y if X is given with probability P (X). Thus the joint entropy H(X, Y ) can be rewritten as Hence we rewrite Eq. (25) as where is regarded as the conditional entropy. It represents the average over X of the Shannon entropy of Y , given X. The above construction is purely classical. For a quantum system on the other hand, the Shannon entropy gets replaced by the von Neumann entropy, S(ρ) = −Trρ log ρ, where ρ is the density operator. Also, P (X, Y ), P (X) and P (Y ) get replaced respectively by the whole system density operator ρ X,Y , the reduced density operator of the subsystem X (ρ X = Tr Y ρ X,Y ) and the reduced density operator of the subsystem Y (ρ Y = Tr X ρ X,Y ).
The notion of the conditional probability P (Y |X) in quantum mechanics requires projective measurements, achieved via a complete set of projection operators, Π i = |ψ i ψ i |, for all i. The density operator of Y after we have measured X is given by, where p i = Tr X,Y ( i ρ X,Y i ). This allows us to define a quantum analogue of the conditional entropy, where S(ρ Y |i ) is the von Neumann entropy for the density operator defined in Eq. (29). The 'min' appearing in the above expression corresponds to measurements which disturb the system the least. This minimises the influence of the projectors. Putting these all in together, we write down the quantum analogues of Eq. (25), Eq. (28), respectively as However, it is clear that the classically equivalent expressions Eq. (25), Eq. (28) of the mutual information need not necessarily be equal in the quantum scenario, for they involve different measurement procedure and quantum measurement on one subsystem can affect the other. Accordingly, the quantum discord is defined as [16,17], We wish to compute the above quantity in the cosmological de Sitter background we are interested in, for a maximally entangled 'in' state, Using the expression for the two-mode squeezed state Eq. (2), the above state can be expanded in terms of the 'out' states. We shall consider two cases below. In the first case the states denoted by the momentum l will be held intact. In another case we shall consider their squeezed state expansion as well. This will help us to probe the correlations between the 'in-out' and the 'out-out' sectors. Accordingly as the first case, using Eq. (2) for the modes k, we rewrite Eq. (33) as The reduced density matrix ρ XY (X ≡ l, Y ≡ k), after tracing out over − k and suppressing the level 'out' for the sake of brevity is given by, We also have, Using Eq. (35) and Eq. (36), we now compute the von Neumann entropies, In order to compute the conditional entropy S(Y |X), Eq. (30), which needs a minimisation over the projective measurements, we take the usual projection operators [16,41,24], wherex 2 1 +x 2 2 +x 2 3 = 1, representing spatial unit vectors. Note that the above operates only on the l sector of Eq. (35). Since the relevant Hilbert space is 2-dimensional, we have taken two projectors which are orthogonal to each other and follows the identity 2 ± = ± . Using now Tr(AB) = Tr(BA), we get from Eq. (29) after some algebra, In the usual parametrisation,x 1 = sin θ cos φ,x 2 = sin θ sin φ,x 3 = cos θ, the conditional entropy is given by Eq. (30) and is found to be independent of the azimuthal angle φ, where the suffix 'min' stands for the minimisation with respect to θ. The quantum discord, Eq. (32), is given by, For different values of the Bogoliubov coefficients, D θ appearing above minimises at θ = π/2, as shown in Fig. 3, Fig. 4. For the maximum value |β k | ≈ 0.707 in the de Sitter, we have D π/2 ≈ 0.377. Note also that D θ is never vanishing, not even for |β k | → 0, in which limit it equals log 2.   59), the variation of the Bogoliubov coefficient correspond to the variation of the rest mass only, and hence for a given species of field, |β k | or the quantum discord is fixed with respect to its spatial momentum k.
As we mentioned earlier, we shall now examine the second scenario where both l and k sectors in Eq. (33) undergo the squeezed state expansion, so that We define the reduced density operator ρ XY (X ≡ l, Y ≡ k) by tracing out over the (− l, − k) sector, where we have suppressed the level 'out' as earlier for the sake of brevity. We find after some algebra, the following von Neumann entropies, Since the l sector of Eq. (43) is two dimensional by the virtue of the Pauli exclusion principle, the conditional entropy, Eq. (30), for this system can still be computed as earlier using the projectors of Eq. (38). 2 We find where p ± and C ± are given by and the suffix 'min' in Eq. (45) refers to minimisation with respect to the polar angle θ. Note that alike Eq. (40), the above is independent of the azimuthal angle φ. Even though the θ dependence of Eq. (46) looks different compared to Eq. (40), the discord, D θ = S(ρ X ) − S(ρ X,Y ) + S(Y |X), still minimises at θ = π/2, as depicted in Fig. 5. For the Bogoliubov coefficients close to their maximal values, |β l | = |β k | ≈ 0.7, we have D π/2 ≈ 0.146, which is less compared to our previous case, Fig. 3. This shows the degradation of correlations in the 'out-out' states compared to the 'in-out' states. We also have plotted the variation of the discord with respect to the two Bogoliubov coefficients in Fig. 6. As earlier, the discord is never vanishing.  The above analysis was done by tracing out the (− k, − l) subsectors in Eq. (42). One can also trace out the other subsectors and analogously compute the quantum discord between them. We expect qualitatively similar results to hold. We also note the qualitative similarity of Fig. 3, Fig. 5 with the earlier studies on scalar and fermion fields in the Rindler spacetime [27,24,25] and as well as on the scalar field in the hyperbolic de Sitter spacetime [41].
Finally, we wish to mention very briefly the logarithmic negativity [7,8] for this system. The log-negativity is a measure of pure quantum correlations, useful in particular, for mixed ensembles. In order to compute this, one needs to take the partial transpose of the matrix representation of Eq. (35), whose eigenvalues are 1/2, 1/2, ±|α k | 2 /2. The log-negativity is given by, Since we do not have any extreme squeezing limit (|β k | → 1, |α k | → 0) for our case, cf. Appendix A, the above quantity is always non-vanishing. Similar conclusion holds for the mixed density operator of Eq. (43) as well. This is qualitatively different for the scalar field theory in the context of transition from the de Sitter to radiation dominated era [41], where such extreme squeezing may happen for which the log-negativity vanishes, even though the quantum discord remains non-vanishing. For our case both log-negativity and the quantum discord are always non-vanishing.

Discussions
In this work we have computed the BMK violation and quantum discord for massive Dirac fermions in the cosmological de Sitter background, respectively in Section 3 and Section 4. For the BMK violation, we have focussed on the vacuum corresponding to the two-and four-mode squeezed states whereas for the latter we have used a maximally entangled Bell state as our 'in' state. Our motivation was to see how the results differ subject to different coordinatisation of the de Sitter as well as the spin of the field. We thus first note the qualitative similarities of the variations with respect to parametrisation angles (Fig. 1, Fig. 2, Fig. 3, Fig. 5) with a scalar field theory discussed respectively in the context of transition from de Sitter to radiation dominated era [42] and the hyperbolic de Sitter spacetime [41]. However, the chief qualitative difference for the cosmological de Sitter from them or the non-inertial frame (e.g. [21]), or even the static de Sitter spacetime (e.g. [50]) is that the Bogoliubov coefficients for the 'in' and 'out' vacua in this case is independent of the spatial momentum or the total energy, Eq. (59). Similar feature is seen for a scalar field theory in the global or the cosmological de Sitter spacetime, e.g. [58,59]. Thus as we have indicated in e.g. Fig. 6, the variation of the Bogoliubov coefficients correspond to the variation of the rest mass only. In other words, unlike the standard R-L entanglement, for a given particle species, any of our results in the cosmological de Sitter spacetime is fixed. We also note that there exists no extreme squeezing limit (|β k | → 1) for the states constructed in Section 2 and we always have |β k | 0.707, as discussed at the end of Appendix A. Due to this reason, as discussed at the end of the preceding section, the logarithmic negativity is never vanishing. This is qualitatively different from the scenario reported in [41] for a scalar field theory where the log-negativity can indeed vanish, indicating the complete decay of quantum entanglement due to particle creation but the discord, being a measure of all correlations, survive.
Investigation of these results and also the decoherence properties in the presence of background primordial electric and magnetic fields seems interesting, due to the Schwinger pair creation mechanism. We hope to return to this issue in a future work.

Acknowledgement
The research of SB is partially supported by the ISIRD grant 9-289/2017/IITRPR/704. Majority of HG's work was done when he was a Masters student at Indian Institute of Technology Ropar, India, and he also acknowledges the same grant for partial support. HG's current research is supported under the CSIR SPMF scheme of CSIR, India.

A The Dirac mode functions and Bogoliubov coefficients
In this Appendix we shall briefly review the solutions of the Dirac equation in the cosmological de Sitter background. The Greek indices appearing below will stand for the spacetime whereas the Latin ones will stand for the local Lorentz frame.
The de Sitter metric in the conformally flat form reads, where H = Λ/3 with Λ being the cosmological constant, and the conformal time η varies from, −∞ < η < 0. We shall also require for our purpose the standard cosmological time t, given by The Dirac equation in a curved background reads, where D µ is the spin covariant derivative e.g. [54], where Σ ab = γ a , γ b /4 and the Ricci rotation coefficients ω's are given by The Gamma matrices appearing in Eq. (49) can be expanded in terms of the tetrad as, γ µ = e µ a γ a and they satisfy the anti-commutation relation [γ µ , γ ν ] + = 2g µν I 4×4 Making the choice of the tetrad, e µ a ≡ Hη diag (1, 1, 1, 1) , the Dirac equation in the de Sitter background becomes, Ψ(x) can be quantised in terms of the conformal time as, where the temporal part of the Bunch-Davies mode functions are given by, e.g. [55], where β s = ±1 is the helicity of the mode corresponding respectively to s = 1, 2. φ (s) k and χ (s) k are eigenvectors of the helicity operatork · γ, with eigenvalues β s . They are given by, wherek x ,k y are momentum unit vectors. Also, ν is a constant given by ν (kη) are the Hankel functions of the first and second kind respectively. It is easy to check using their asymptotic properties [57] that u(η → −∞) ∼ e −ikη and v(η → −∞) ∼ e ikη in Eq. (54) and thus respectively behave as positive and negative frequency solutions in the asymptotic past. The orthonormality of the mode functions appearing in Eq. (53) can be checked using Eq. (55) by computing their inner products, say on an η → −∞ hypersurafce. Eq. (54) should be regarded as the fermionic analogue of the Bunch-Davies mode functions [55,56].
The 'out' modes are those which have definite positive and negative frequency characteristics in the asymptotic future, η → 0 − . It is easy to see from the asymptotic expansion of the Hankel function that the conformal time cannot be a good coordinate for this purpose and accordingly we resort to the cosmological time, t. In terms of t, we have from Eq. (54) in the asymptotic future, Accordingly, we choose the temporal part of the 'out' modes as, Thus we have the field quantisation, where the absolute values indicate that the Bogoliubov relations can be defined up to some global phase factors. Note that Eq. (60) is independent of k. This originates from the fact that the frequencies of the 'out' modes, Eq. (57) are solely determined by the rest mass m of the field. Similar thing occurs for scalar field theory in the global and the cosmological de Sitter backgrounds [58,59]. Although we are certain that this result for fermions must have been reported earlier, we were unable to find any reference stating the same.
Recall also that a conformally invariant field like a massless fermion cannot create particles in the conformal vacuum of a conformally flat spacetime such as the de Sitter, e.g. [54]. Thus we have an apparent ambiguity with setting m = 0 in Eq. (60). However, we note that in this case the modes Eq. (57) have no positive and negative frequency characteristics at all. Thus it seems that the m → 0 limit for Eq. (60) is not smooth and one needs to treat the purely massless case separately, retaining explicitly its conformal symmetry by using the conformal vacuum. However, this case will not be relevant for our current purpose. Finally, we note the range of the Bogoliubov coefficient, 0 ≤ |β k | 0.707 where the upper bound correspond to m/H → 0.