Fermionic vacua and entanglement in hyperbolic de Sitter spacetime

We point out a non-trivial feature of the vacuum structure of free massive or massless Dirac fields in the hyperbolic de Sitter spacetime. Here we have two causally disconnected regions, say $R$ and $L$ separated by another region, $C$, appropriate to discuss the long range correlations of two superhorizon separated observers. There are local modes of the Dirac field having supports individually either in $R$ or $L$ as well as global modes found via analytically continuing the $R$ modes to $L$ and vice versa. However, we show that unlike the case of a real scalar field, the analytic continuation does not preserve the orthogonality of the resulting global modes. Accordingly, we need to orthonormalise these modes by making their suitable linear combinations prior to the field quantisation, in order to preserve the canonical anti-commutation relations. The most general form of this orthonormalisation is achieved via introducing a spacetime independent continuous parameter, $\theta_{\rm RL}$, resulting in a one-parameter family of global vacua. We emphasise that unlike the case of the so called de Sitter $\alpha$-vacua, introducing such parametrisation is mandatory in our present scenario, in order to preserve the natural canonical structure. Using these vacua, we next investigate both entanglement and R\'enyi entropies by tracing over the states belonging to either $R$ or $L$ region and demonstrate their variations with respect to the parameter $\theta_{\rm RL}$.


Introduction
The de Sitter (dS) spacetime is a solution of the Einstein equations with a positive cosmological constant and without any other matter field. It is a maximally symmetric spacetime with a constant positive curvature, having various physical interests. First, owing to the recent observed phase of the accelerated expansion of our universe, there is a strong possibility that our universe is endowed with a small but positive cosmological constant, Λ. Second and perhaps more importantly, the high degree of homogeneity and isotropy of the current universe at large scales indicate that our early universe had to go through an inflationary phase. For modelling such inflationary scenario, dS is a very well motivated candidate [1]. Moreover being maximally symmetric, it is comparatively easier to do explicit computations in dS.
The dS spacetime, being maximally symmetric and positive curvature, could also be derived by compactifying the Minkowski spacetime of one higher dimension as follows [2] (also references therein). If we are interested in (3 + 1)-dimensions, we consider a (4 + 1)-dimensional Minkowski spacetime with coordinates (T, X 1 , X 2 , X 2 , X 4 ) and choose a hypersurface where H 2 = Λ/3. After making appropriate coordinatisation on it, the resulting (3+1)-dimensional spacetime turns out to be dS. The five translations out of the 15 generators of the Poincarré group of the Minkowski spacetime are sacrificed during this compactification, making the resulting dS spacetime maximally symmetric, with SO(4, 1) being the isometry group. Accordingly, the dS can be covered by a single chart, having topology R 1 × S 3 with S 3 being the spatial slicing. The global covering of the dS contains both contracting and expanding parts. In cosmological applications however, one is instead interested in its expanding segment, leading to the flat spatial slicing, covering only the half of the global dS. Now, due to the accelerated expansion of the dS, not even the entire of the cosmological dS could be causally connected. This leads to the static coordinatisation, valid only up to the radius H −1 , The radius H −1 is known as the cosmological event horizon, sets a natural causal boundary of the universe an observer located at R = 0 can see. The cosmological event horizon is a non-degenerate Killing horizon having thermal properties analogous to that of a black hole, though it is observer dependent owing to the maximal symmetry of dS [3,4,5]. Also just like a black hole horizon, the t−r part of the near cosmological horizon metric can be cast into a Rindler form. One can then use the standard replica trick method [6] to investigate the entanglement of states of a quantum field theory residing inside and outside the horizon and it would turn out to be proportional to the horizon area. It is also worth mentioning that one can have stationary black hole solutions embedded in a de Sitter universe, known as the de Sitter black holes. Such black hole spacetimes are naturally associated with two temperatures for particle creation corresponding to the two Killing horizons, leading to complex vacuum structures, e.g. [3,7].
The entanglement entropy mentioned above is not long range at the leading order. This can be due to the pile up of the accessible microstates of the field in an infinitesimal width across the Killing horizon [8]. Is there a natural set up to investigate the long range quantum correlation of two observers in the dS space, not only causally separated but so by a large distance, say of the order of the superhorizon size? This question was first addressed in [9] using the coordinatisation of [10,11], known as the hyperbolic dS spacetime describing two casually disconnected expanding bubbles, can be realised as follows. Imagine a spacetime region, causally connected initially is expanding in the global dS. After 'sufficiently' long time, certainly not all part of it will remain causally connected. Such description can be achieved via appropriate analytic continuation of the Euclidian dS into the Lorentzian sector to show that the entire dS splits into three causally disconnected regions given by Eq. (2), Eq. (3). The relationship between their coordinatisation with that of the Euclidian dS is given by Eq. (1). Fig. 1 represents the Penrose-Carter diagramme of this geometry. The chief difference of this with the static one is -in this case we do not have a Rindler like wedge to separate two causally disconnected regions R and L, but now an entire superhorizon region C sits in between them. Since the spatial foliations in R and L are 3-hyperboloids, they are also known as the hyperbolic or open dS spacetime.
In this work we wish to investigate the entanglement properties of the fermionic vacua in the hyperbolic dS. Let us first briefly review the case of a real scalar field [9]. We first define orthonormal local basis mode functions having supports either in R or in L, with definite positive or negative frequency behaviour in the asymptotic past. We then make a field expansion using them and define the local vacuum as a direct product state between the vacua of R and L. However, if there is any correlation between the local states, clearly there must be some mode functions having supports in both R and L. These modes are obtained by relaxing the earlier condition of rigid support and then by analytically continuing them from one region into the other along a complex path dictated by Eq. (1) [10,11]. These global modes turn out to be orthogonal as well and after normalising, we make an alternative field expansion using them and define a suitable global vacuum. The comparison of the two mode expansions yields a Bogoliubov relation, yielding in turn an expansion of the global vacuum in terms the local states. It turns out from this expansion that the states belonging to R and L are entangled. One then finds the matrix representation of the reduced density operator ρ R found by tracing out the states belonging to either R or L and the entanglement entropy density is given by −Tr(ρ R ln ρ R ). The full entanglement entropy is then found by integrating over the spatial volume and regularising the volume integration by introducing a suitable cut off. The entanglement entropy thus found will not be proportional to any area. This procedure will be more explicit in the due course of this paper.
In fact a lot of effort has been given to explore such long range quantum correlations in the dS spacetime so far, in various coordinatisations, e.g. [12]- [28]. They not only involve the computations of the entanglement or the Rényi entropy (e.g. [29] and references therein), but also studies of other measures like Bell's inequality violation, entanglement negativity and discord, in the Bunch-Davies or more general α-vacua (e.g. [30,31,32] and references therein). We also refer our reader to [33] for a study on the effect of the background fields on the entanglement entropy in dS. See e.g. [34,35] for discussions on holographic entanglement entropy in dS. In particular, possible observable effects of quantum entanglement in various primordial perturbations has been suggested in [36,37,38,39,40].
Most of the references cited above focus on the scalar field theory and discussions on other spin fields seems really sparse. In [41], the entanglement properties of a Dirac fermion was discussed and certain qualitative differences with a scalar field were pointed out. See [38,42,43] for aspects of fermionic entanglement in the cosmological dS. See also [44,45] for discussions on fermionic entanglement in the Rindler spacetime.
We wish to point out here a non-trivial feature associated with the vacuum structure of the Dirac fermion in the hyperbolic dS. After reviewing the local and global mode solutions of the free Dirac equation in Section 2, we show in Section 3 that unlike the case of a scalar field [10], the analytic continuation procedure does not preserve the orthogonality of the global mode functions, Eq. (29). However, in order to preserve the canonical anti-commutation structure we must work with orthonormal modes only, e.g. [46] and hence we need to orthogonalise them by making suitable linear combinations. We argue in Section 3 that a priori such orthogonalisation must be sufficiently general and hence we introduce a one parameter (θ RL ) family of orthonormal modes in Eq. (32). Using these modes along with the local ones, we quantise the fermion field and find out the Bogoliubov relations compatible with the canonical anti-commutation structure in Section 4.1. Since the global vacuum Eq. (48), depends upon the parametrisation, so will the expectation value of any observable quantity. We compute the entanglement and the Rényi entropies respectively in Section 4.2 and Section 4.3 and depict their variation with respect to θ RL , e.g. in Fig. 2. Finally we conclude with some possible future directions in Section 5.
We would like to emphasise that such one parameter freedom is analogous but not the same as the de Sitter α vacua, e.g. [31], for in our case introducing the parametrisation was completely mandatory to preserve the canonical structure and this seems a unique feature of the Dirac fermions in the hyperbolic covering of dS only (c.f., discussions below Eq. (32)). To the best of our knowledge this has not been reported earlier.
We shall work with the mostly positive signature of the metric in (3 + 1)-dimensions and will set c = G = = κ B = 1 throughout.

The Dirac modes in hyperbolic dS
In this section we shall review the basic geometry of the hyperbolic dS, the local or the Bunch-Davies Dirac modes and their analytic continuation to form the global modes. The detail of these can be seen in [10,11,41,47,48].

The separation of variables
The dS with hyperbolic spatial slicing is obtained by analytic continuation of the Euclidean sphere S 4 to the Lorentzian sector [10,11]. Such analytic continuation gives rise to three distinct and causally disconnected regions of the global dS spacetime, say R, L and C. The relationship between the Lorentzian and Euclidean coordinates (τ, ρ) of these three regions is given by, The metrics in the three regions read, where dΩ 2 is the metric on S 2 and r and t are scaled by H −1 . We shall briefly review now the solutions of the free Dirac equation in the hyperbolic dS spacetime, Eq. (2), first obtained in [47,48]. Note that there exists no normalisable solution in the region C, Eq. (3) [41]. The Dirac equation in a generally covariant form reads where the spin covariant derivative is defined as where a, b stand for the local Lorentz indices. The Ricci rotation coefficients ω's are given by, e.g. [49] : ω µ a b = −e b ν ∂ µ e a ν − Γ λ µν e a λ and Σ ab = γ a , γ b /4. We also have in accordance with the mostly positive metric signature, so that, [γ µ , γ ν ] + = 2g µν 1 4×4 The Dirac equation expanded in the background of Eq. (2) reads where we have suppressed the region indices, f = sinh t and the 'tilde' stands for the spatial section of the metric i.e., / . Taking now the ansatz, where each of the above entires are two component spinors, we obtain from Eq. (7), Combining the equations above, we obtain the single second order equation for φ t , The equation for φ b is analogous, except one has a minus sign instead of plus in front of the second term on the left hand side of the above equation.

The local, Bunch-Davies mode functions in R and L regions
Let us first focus on region R of Eq. (2). Even though the metric in the region L is formally similar, note from Eq. (1) that t L = −t R − iπ. Since the Dirac operator is linear in derivative, this will give rise to some alterations in the solutions for the mode functions in L. Also, the modes we are going to find below will be assumed to have supports in the respective regions only and hence we shall call them local, as opposed to the global ones to be derived in Section 2.3 and Section 3. We make the following ansatz for the variable separation, where z = cosh t and the two component spinor χ pjm (r, θ, φ) is purely spatial, defined on the hyperboloid of Eq. (2), satisfying [47,48], where p is a real, positive and continuous eigenvalue and (j, m) stand for the eigenvalues of the spin-1/2 weighted spherical harmonics. The symbol '(−)' in χ basically stands for the minus sign appearing on the right hand side of the above equation. We shall see below that there exists another spinor orthogonal to χ (−) , to be denoted by χ (+) for which the equation analogous to Eq. (12) has a plus sign on the right hand side.
The two spinors χ (±) form a complete orthonormal set of eigenfunctions of the operator / ∇, with the normalisation being fixed as where the 'dagger' denotes the hermitian conjugation, s = ± and h is the determinant of the metric of the 3-hyperboloid in Eq. (2). For our present purpose we do not need to use anything more explicit on χ (±) . Plugging Eq. (11) into Eq. (10) one now obtains The general solution of the above equation reads, One can find analogous solutions for g(z) appearing in Eq.
It is easy to see by using the asymptotic form of the hypergeometric functions [51], that u p and v p realise positive frequency solutions in the distant past t → 0 (z → 1). Thus their complex conjugations become negative frequency solutions.
We shall now choose our basis modes by ascribing to them definite positive and negative frequency behaviour existing in the asymptotic past. Also, once we choose the upper component of Eq. (11), the lower component is automatically dictated by Eq. (9). Putting these in together, the two possible mode functions in R, respectively positive and negative frequency are given by As we mentioned earlier, we also need to consider the solutions obtained by changing the spatial part of the ansatz from χ . This corresponds to taking p → −p in Eq. (12). We have Having defined these modes, we need to further make sure that they are orthogonal with respect to the Dirac inner product, defined as where the indices (a) and (b) collectively characterises the nature of the spinor, including the various eigenvalues. The above equation is analogous to Eq. (13) with the difference of having the time dependent function coming from a t = const. section of Eq. (2). Also, since we are dealing with local R modes now, the spacelike hypersurface above will be located in that region only. Since the inner product is independent of time e.g. [50], the integral of Eq. (18) could be evaluated with any suitable choice of time. In particular, by choosing t → 0 it is easy to see for the four modes of Eq. (16), Eq. (17) that a) the inner product of a mode with itself equals δ(p − p ) δ jj δ mm , whereas (using also Eq. (13)) b) with any two different modes vanish. This establishes the orthonormality of our chosen local mode functions.
Having thus obtained the basis mode functions in the R, we now move on to L. The coupled first order equations in this region becomes, Note the alterations of plus and minus signs compared to Eq. (9), as we indicated at the beginning of Section 2.2. This implies that the changes : u p → v p and v p → −u p would give us the correct modes in this region. We thus respectively find the analogues of Eq. (16), Eq. (17) in the region L, and The orthonormality of these modes can be checked as earlier. Also, since these modes are local, the orthogonality between the R and L-modes follows trivially. This completes our discussion on the local mode functions. Since by construction Eq. (16), Eq. (17), Eq. (20) and Eq. (21) have definite positive and negative frequency behaviour in the distant past, they can be regarded as the local Bunch-Davies modes for the Dirac field.

Analytic continuation and the global modes
Eq. (1) indicates that the charts covering the R and L regions would be diffeomorphic to each other via some complex coordinate transformations. This allows us, as opposed to the local modes, to write down modes belonging to one region into another [41], reviewed below. Modes thus found which have supports in both R and L (recalling there exist no normalisable mode in the region C for Dirac fermions, as opposed to a scalar field [41]) are called the global mode functions. Since the coordinate transformation between these two regions are complex, writing an R mode in terms of the L coordinates and vice versa would essentially involve an analytic continuation of the hypergeometric functions Eq. (15) and their complex conjugations on the complex plane. The spatial functions χ (±) however, do not undergo any formal changes under this procedure. The argument of a complex coordinate z will always be taken to be its principle value (−π < arg(z) < π) below. For u p (z R ) in Eq. (15), we extend the functions analytically from R (z R ≥ 1) to C (−1 ≤ z R ≤ 1) by taking z R − 1 = e −πi (1 − z R ). Then we change the variable by z R = −z L in the C region. Then we continue analytically from C to L (z L > 1) by performing 1 − z L = e πi (z L − 1), to obtain In order to cast the above expression into the form of our initial mode functions, we recall now the identity, e.g. [51] F Substituting this into Eq. (22) with z ≡ (1 + z L )/2 and also using where we have written Similarly, The opposite procedure, i.e. the analytic continuation L → R yields exactly formally similar results. We list below the resulting (unnormalised) global mode functions originating from the R → L analytic continuation, Likewise we have for the L → R analytic continuation, Note that any specific mode looses its characteristic of being purely positive or negative frequency under the analytic continuation. Since each of the above mode functions have supports in both the regions, when normalised, they are supposed to be the global versions of the local modes appearing in Eq. (16), Eq. (17), Eq. (20) and Eq. (21). However, we shall see below that the modes of Eq. (27) and Eq. (28) do not form an orthogonal set under the global inner product, in contrast to the scalar field theory [10]. Hence we cannot simply treat the normalised versions of the above modes as our global basis modes.

Constructing the global orthonormal modes
The normalisation integration for the global modes looks formally similar to Eq. (18), with the difference that we need to perform the integration now on a suitably chosen global spacelike hypersurface existing both in R and L. As in the scalar field theory [10], we choose the global hypersurface to be a union of three different pieces in Fig. 1 -one t R = const. hypersurface in R and one t L = const. hypersurface in L, smoothly connected via another spacelike hypersurface existing in C. Since there exists no normalisable mode function in C [41], the global normalisation integral then reduces to two local pieces in R and L.
The global inner product then can be explicitly computed using the expressions of Eq. (27), Eq. (28).
As an example, let us consider the pair {Ψ +L (−) , Ψ −R (−) }. The global inner product for them reduces to, where the suffix 'G' stands for global. The inner products on the right hand side are of course, local. We can express Ψ −R (−) in terms of z L (for the first) and Ψ +L (−) in terms of z R (for the second) via Eq. (27) and Eq. (28). Then the local inner products can be evaluated explicitly using Eq. (13) and Eq. (18). We find explicitly (after suppressing the various δ-functions arising from Eq. (13) for the sake of brevity), where λ 2 is given by Eq. (25). It can further be explicitly checked in the similar manner that the eight global fermionic modes of Eq. (27), Eq. (28) can be grouped into four pairs such that the members of any pair are not orthogonal (although inter-pair orthogonality is satisfied) with respect to the global inner product, given by Note that there is no pair belonging to either only R or only L. Also λ 2 = 0 irrespective of whether the field is massive or massless. In any case, it is evident that we cannot use the modes appearing in Eq. (27), Eq. (28) as our global basis modes, for they would lead to a non-preservation of the canonical anti-commutation relations when used as basis of field expansion, e.g. [46] and references therein. Hence we need to make suitable linear combinations between the members of each pair of Eq. (29), to orthogonalise them. This can be achieved via the standard Gram-Schmidt procedure in which one vector is held fixed while the other is suitably 'rotated' to make it orthogonal to the first. Thus from the first pair of Eq. (29) for example, we may construct a new pair, say Ψ , Ψ −R (−) , orthogonal under the global inner product, where Clearly, such orthogonalisation procedure is never unique, for we could also have constructed a different orthogonal pair, say Ψ, Ψ +L (−) , with Note also from Eq. (27), Eq. (28) that when a mode of one region is continued analytically into another, it becomes a superposition of both positive and negative frequency modes. Thus in order to accommodate sufficient generality into our orthogonalisation scheme, we must 'rotate' both the solutions simultaneously in an equal footing. Also, since such 'rotations' could be done in an infinite number of ways, we may characterise them by introducing a continuous parametrisation. Thus, after explicitly evaluating the inner products in the expressions of Ψ , Ψ appearing above, we introduce two new mode functions where ∆θ 1 and ∆θ 2 are one-parameter coefficients given by and The angle θ RL does not depend upon any spacetime coordinate. It is easy to check that the two mode functions defined in Eq. (30) are indeed orthogonal under the global inner product. Moreover, since θ RL is a continuous parameter, we have infinite number of choices of such orthogonal pairs. In particular, for θ RL = 0, and π/2 respectively, we recover our earlier pairs, Ψ, Ψ +L In the same spirit, by choosing appropriate linear combinations in three other cases in Eq. (29) one can generate orthogonal pairs of global mode functions. The full set of eight orthonormal global modes is given by, where we have written for the normalisations, Using Eq. (29), the orthonormality of these modes under the global inner product can explicitly be checked at once. This orthogonalisation of the global modes was not performed in [41].
We wish to emphasise a couple of things before proceeding any further. First, note once again that Eq. (29) dictates that we must make linear combinations between the members of a given pair (such as in Eq. (30)) in order to form an orthogonal pair, for the orthogonality between any two modes picked up from two different pairs is satisfied. Thus we could have taken four different θ RL 's to obtain the eight global orthonormal mode functions. However, introducing such abundance of parameters would certainly weaken the predictability of the theory and hence we have kept the parametrisation minimal.
Second, it is clear from Eq. (27) and Eq. (28) that the mixing of positive and negative frequency solutions in Eq. (32) depends upon θ RL . Accordingly, when we use Eq. (32) to make a field quantisation, the vacuum corresponding to them, or the global vacuum, will depend upon the parametrisation θ RL . Thus we may expect physical observable quantities to depend upon the parametrisation. Since it does not seem a priori that there is any reason to prefer any particular value of θ RL , the actual value or any bound on it can only be determined via the comparison of the theoretical prediction with the observational data.
Finally, note that such one-parameter family of vacua invariant under the dS isometries is already known, i.e. the so called α-vacua, e.g. [30,31,32] and references therein 1 . The α-vacua are constructed out of suitable linear superpositions of the basis modes corresponding to the Bunch-Davies vacuum via a continuous parameter independent of the spacetime coordinates, analogous to what we have done here. However, it should be emphasised that the modes related to the Bunch-Davies vacuum indeed furnish a perfectly valid description of the theory and the α-vacua, even if it turns out to be phenomenologically more preferable during future cosmological observations, are not necessary on the grounds of any fundamental consistency in defining the field theory. In our case however, we have seen that introducing θ RL was absolutely necessary for obtaining the most general orthogonal global modes, which we may treat (after normalising) as our global basis mode functions.
Our analysis seems to thus indicate that such one-parameter freedom for the fermionic field in the hyperbolic dS spacetime, Eq. (2), Eq. (3) is an intrinsic feature. Note that such feature is not present for fermions in dS spacetime with flat spatial slicing e.g. [31], nor it is present for a real scalar field in the hyperbolic slicing [10]. Since this one parameter freedom essentially originates from the analytic continuation process described in Section 2.3, such freedom does not seem to be present for a complex scalar field as well, for its basis mode functions will be identical to that of a real scalar field. We also note that for a massive vector field the basis modes are identical to that of the scalar [53] and hence this particular feature will be absent for it as well. Thus the natural emergence of such intrinsic one parameter freedom seems to be a novel characteristic of fermions only, in the hyperbolic dS spacetime. This is one of the chief results of this paper and to the best of our knowledge, it has not been reported earlier.
Having defined the appropriate orthonormal global mode functions, in the next section we shall make expansion of the field operator both in terms of the local and global modes and will obtain the relevant Bogoliubov transformation coefficients. Since the hyperbolic dS is a natural set up to discuss the nonlocal and long range correlations of two superhorizon separated observers, we are specifically interested in computing the entanglement entropy and the Rényi entropy for the fermionic vacua. Such investigation will provide us the effect of the parameter θ RL into the non-locality, along with the rest mass of the field.

Computation of the entanglement and the Rényi entropies
The chief goals of this section are a) to discuss the canonical field quantisation in terms of both local and global modes derived in Section 2.2 and Section 3 and to find out the relevant Bogoliubov coefficients and b) using these coefficients we shall express the global vacuum as a squeezed state over the local vacuum. From this we shall compute the reduced density operator by partially tracing the global density operator over either R or L region. This will lead to the numerical computations of the entanglement and the Rényi entropies.

The field quantisation and the Bogoliubov coefficients
Let us first consider the local modes in the R and L region, Eq. (16), Eq. (17), Eq. (20) and Eq. (21) and make the following expansion of the field operator, where s = ±. Note that since the local modes have supports only in the respective regions, the addition of the modes belonging to R and L appearing above has to be understood as a non-local one. The coefficients of the asymptotic negative frequency modes, c R † (s) , d L † (s) have their natural interpretation as the fermionic creation operators for particle and antiparticle in R and L regions respectively. Likewise we may interpret c R (s) , d L (s) as the fermionic annihilation operators in the region R and L respectively. These operators are postulated to satisfy the anti-commutation relations and all other anti-commutators vanish. Eq. (35) can be considered as the quantisation condition corresponding to the local Dirac modes. We may define the local vacua |0 R , |0 L straightforwardly, Likewise, we can also express the Dirac field in terms of the orthonormal global fermionic modes given by Eq. (32). However, the usual form of the field expansion in terms of the (asymptotically) positive and negative frequency modes is not conceivable now, for the global modes are hybrid in nature originating from, as we have seen, the analytic continuation process. This poses the question : how do we identify the creation and annihilation operators? Eq. (25) shows that in the limit m → 0 and p → ∞ , both λ 1 and λ 2 are vanishing, showing in this limit we do not have any analytically continued modes in Eq. (27) and Eq. (28) and accordingly, Eq. (32) reduces to the local mode functions, having well defined positive or negative energy eigenvalues in the asymptotic past. Since in that limiting scenario we do not have any trouble with identifying the creation and annihilation operators, we may make the following expansion of the field Ψ in terms of the global mode functions, where we interpret a 1 , . . . , a 4 and b 1 , . . . , b 4 as the annihilation operators related to the global modes. The global vacuum is defined as We shall see below that equating Eq. (34) with Eq. (37) would lead to Bogoliubov coefficients which will automatically guarantee the desired anti-commutations between a and b operators. Note also that since the Ψ i 's in Eq. (37) depend upon the parametrisation θ RL via Eq. (32), the creation and annihilation operators and hence the vacuum states defined in Eq. (38) will also depend upon θ RL . This will be more explicit from the computations done below. In order to obtain the Bogoliubov coefficients, after equating Eq. (34) and Eq. (37), we successively take eight inner products with both sides with respect to the eight global basis modes defined in Eq. (32). When we take the inner product of a global mode with a local one, we need to use Eq. (27) or Eq. (28) in order to express the former in terms of the local modes. Putting these all in together and using Eq. (13) we obtain for example, for the global mode Ψ 1pjm where we have written α = 2λ 2 ∆θ 1 N 2 b and have suppressed the eigenvalues p, j, m, for the sake of brevity. Similarly, by taking the inner products with the seven other modes Ψ 2 , Ψ 3 . . . , Ψ 8 , we obtain where we have written Using Eq. (35) it can be checked that the global operators satisfy the desired anti-commutation relations, where σ = 1, 2, 3, 4 with all other anti-commutations vanishing. We once again emphasise that had we not properly orthonormalised our global modes, we would not have preserved the canonical anti-commutation structure, essential to preserve the spin-statistics of the field theory. Now, by observing the right hand sides of Eq. (39), Eq. (40), it becomes evident that the set of global operators (a 1 , a 3 , b 2 , b 4 ) and (a 2 , a 4 , b 1 , b 3 ) form two disjoint sectors, for their operator contents on the right hand side are different. This implies that the global vacuum defined in Eq. (38) can also be split into two subspaces, where |0 1 is defined by For the rest of this paper, we will work with the |0 (1) subspace only. As long as we are only concerned with the computation of the entanglement and Rényi entropies, the other subspace, |0 (2) , will produce identical results. From Eq. (39), Eq. (40), it is clear that |0 (1) can be constructed as a squeezed state over the local vacua defined in Eq. (36), where ξ ij 's are four complex numbers, c R (+) |0 (+) (+) L = 0 and N is a normalisation. Also, since the operators c's and d's anti-commute (Eq. (35)), we may further make the following decompositions, . Note that in Eq. (43), we have chosen one particular ordering of the operators and the states. However, a different ordering will always lead to different results owing to the anti-commutation relations. We refer our reader to [45] and references therein on a discussion on possible ambiguities related to such orderings in fermionic field theories.
We may now expand the exponential in Eq. (43) keeping in mind the various anti-commutations and then annihilate |0 (1) by a 2 , a 4 , b 1 , b 3 , to obtain after a lengthy computation, where, The normalisation of the vacuum is given by, Putting these all in together, we may now explicitly write down the global vacuum in terms of the local states, The above expression shows that there will be pair creation in the R and L regions with respect to the global vacuum. Since the states belonging to R and L cannot be factored out, those pairs will be entangled, e.g. [54]. Note also that |0 (1) depends upon θ RL via the coefficients ξ 1 and ξ 2 , Eq. (45), Eq. (46). Now, the vacuum expectation value of any observable say A, with respect to the global vacuum equals Tr (ρ global A), where ρ global = |0 (1) (1) 0| is the density operator corresponding to the global vacuum. This shows as we anticipated earlier, the expectation value will explicitly depend upon the parametrisation θ RL . Certainly this will be true for global excited states as well. However, since the coordinate coverings Eq. (2), Eq. (3) are more apt to describe long range correlations in dS, we shall limit our investigation to such effects only. In particular, starting from ρ global , we shall compute below the reduced density operator and the entanglement and Rényi entropies for |0 (1) .

The entanglement entropy
Let us imagine a local observer residing in R, causally disconnected from L. We start with the full density operator ρ global defined above and trace over the states of the L region inaccessible to that observer, to obtain the reduced density operator in R, We could also have considered an observer instead in L, but that will give identical results. We find a matrix representation of ρ R in terms of the R-state vectors, The entanglement entropy for a mode (with a given p and θ RL value) is given by where λ (i) R 's denote the eigenvalues of Eq. (50). Note that for a bosonic field the density operator, be it global or reduced, is infinite dimensional. Accordingly, in order to trace out the states in L and to find out the eigenvalues of the reduced density matrix, one needs to diagonalise the global density operator first, yielding a reduced density matrix automatically diagonal [9,12]. Similar computations can be seen for fermions in [41]. However, since the the matrix representation Eq. (50) is finite dimensional owing to the Pauli exclusion principle, we have worked without diagonalisation and found out the eigenvalues λ The entanglement entropy per unit volume between R and L is obtained by integrating over all p modes. The final expression of the entanglement entropy is obtained by further integrating the result over the purely spatial section of Eq. (2). Since S(p, m; θ RL ) has no spatial dependence, the integration, being done over a non-compact space, would diverge. Accordingly, one needs to put a cut off at some 'large' radial distance. The resultant regularised volume integral equals 2π, see [9] for details. The regularised entanglement entropy then equals where D(p) is density of states in the momentum space corresponding to the spatial section of Eq. (2) for spin-1/2 fermions, given by D(p) = ( 1 4 + p 2 )/(2π 2 ) [48,55]. The integral in Eq. (52) is convergent and can be calculated numerically.
We have plotted various characteristics of the entanglement entropy in Fig. 2 and Fig. 3 as a function of the parameter ν 2 = 9/4 − m 2 /H 2 subject to different values of θ RL . We may chiefly note the following features.
In Fig. 2, the curves corresponding to θ RL = 0, π/2, are exactly coincident and they show maximal entanglement for all values of ν 2 . The coincidence corresponds to the fact that the coefficients ξ 1,2 (θ RL = 0) = ±ξ 1,2 (θ RL = π/2), in Eq. (45), Eq. (46), and the eigenvalues λ (i) R of the reduced density matrix Eq. (50) are quadratic in ξ 1 and ξ 2 . However, it is easy to see from Eq. (45), Eq. (46) that there is no such general symmetry for θ RL → π/2 − θ RL . Also, while most of the curves in Fig. 2 are monotonic, the curve corresponding to θ RL = π/3 shows extrema. Finally, for any given value of θ RL , the entanglement entropy has its maximum value in the massless limit, ν 2 = 9/4. This might correspond to the fact that in this limit the creation of particle-antiparticle pairs is energetically most favourable. Since such pairs are entangled, Eq. (48), it is perhaps reasonable to expect that the entanglement entropy also gets its maximum value in the massless limit. However, note also that such expectation does not explain the aforementioned extrema in Fig. 2 for θ RL = π/3 in the near massless regime. Thus perhaps we need to investigate other measures of the quantum correlation in order to understand the full implication or at least a more detailed quantification of the parametrisation. The green curve corresponds to θ RL = 0, π/2, the red to θ RL = π/6, the blue to θ RL = π/5, the black to θ RL = π/4 and the pink to θ RL = π/3. We may chiefly note that a) the θ RL = 0, π/2 plots are exactly coincident and they give maximum entanglement entropy for all values of ν 2 b) the pink curve (θ RL = π/3) has extrema, whereas the other plots are monotonic. c) for any given value of θ RL , the entanglement entropy has its maximum value in the massless limit, ν 2 = 9/4. As an aside, we note that one physical quantity of interest is the expectation value of the number operator with respect to the global vacuum, which will give us the number of created particles at a given mode. We shall do it in the massless limit only, in which Eq. (45) and Eq. (46) reduce to ξ 1 = 0, ξ 2 = −e −pπ 1 + e pπ cos 2θ RL e pπ + cos 2θ RL Thus, 1 + e 2pπ e pπ +cos 2θRL For θ RL → 0, π/2, this reduces to the familiar Fermi-Dirac distribution, Away from these values of θ RL , even though the field is massless, the spectrum is not thermal. However this should not be surprising, for two different values of θ RL would correspond to two different global vacuum in Eq. (48), related via some Bogoliubov transformations. Accordingly the spectra of created particles with respect to them would also be different. For α-vacua in the dS spacetime with flat spatial slicing, such non-thermality was also noted earlier in [32].
Nevertheless, it is clear that obtaining a massless thermal spectrum and maximally entangled R, L vacua for θ RL = 0, π/2 (Fig. 2) does not essentially imply that those two values should be preferred. For just like the α-vacua, θ RL could also be related to the initial energy condition associated with the R and L regions. It might also be possible that such one parameter ambiguity is borrowed from the global dS spacetime itself. Accordingly, its actual value should only be determined via some observational effects.

The Rényi entropy
Before we end, we wish to briefly discuss the Rényi entropy, a one-parameter generalisation of entanglement entropy, e.g. [29] S q = ln Trρ q 1 − q , q > 0 (54) In the limit q → 0, lim q→0 S q = ln n where n is the dimension of the Hilbert space we are interested in.
Also for q → 1, Eq. (54) reduces to the entanglement entropy as follows. We obtain using L'Hopital's rule, for a normalised density matrix Tr (ρ q ln ρ) Trρ q = −Tr (ρ ln ρ) Finally for q → ∞ we have, where λ max is the largest eigenvalue of the density matrix.
The final Rényi entropy, akin to the expression for the final regularised entanglement entropy, is given by, We have plotted S q (m, θ RL ) in Fig. 4a, Fig. 4b and Fig. 5, as earlier with respect to the parameter ν 2 = 9/4 − m 2 /H 2 , with different values of θ RL . We note chiefly that as a whole, the qualitative nature of the Rényi entropy for different q values we have taken remains the same as that of the entanglement entropy. In particular, a) the values θ RL = 0, π/2 gives maximum Rényi entropy for all values of ν 2 and b) the extrema for θ = π/3 is still present. The green curve corresponds to θ RL = 0, π/2, the red to θ RL = π/6, the blue to θ RL = π/5, the black to θ RL = π/4 and the pink to θ RL = π/3. The green curve corresponds to θ RL = 0, π/2, the red to θ RL = π/6, the blue to θ RL = π/5, the black to θ RL = π/4 and the pink to θ RL = π/3. Since the q value is close to unity, this plot is almost the same to that of the entanglement entropy, Fig. 2, see text for discussions.
This completes our discussion on the computation of the entanglement and the Rényi entropies.

Discussions and outlook
In this paper we have investigated the entanglement and the Rényi entropies between the R and L states of the Dirac field in the hyperbolic dS spacetime, Eq. (2), Eq. (3), as a measure of the long range non-local quantum correlations between these two regions. We have reviewed the derivations of the local or the Bunch-Davies orthonormal mode functions in Section 2.2, Eq. (16), Eq. (17), Eq. (20) and Eq. (21). These modes have only supports in either R or L region. We have further reviewed in Section 2.3, the global mode functions obtained via the analytic continuation of the local modes form one region into another, Eq. (27), Eq. (28). However, we have shown in Section 3 that the global modes thus obtained are not orthogonal and hence we cannot use them to make any field expansion that preserves the canonical anti-commutation structure. Accordingly, we have introduced a continuous, one parameter (θ RL ) family of modes in Eq. (32), orthonormal under the global inner products. We have used these modes to make the expansion of the field operator in Section 4.1 and have found a set of Bogoliubov coefficients, Eq. (40), preserving the canonical anti-commutation structure appropriate for the fermionic field theory. The global vacuum state is expressed with the help of these coefficients as a superposition state of the local states in Eq. (48). It is evident from this expression that the states belonging to R and L are entangled. Accordingly in Section 4.2 and Section 4.3, we have computed the entanglement and the Rényi entropies for the fermionic vacua and have investigated the effect of the parameter θ RL onto them.
The chief findings of these paper could be summarised as follows. First, the natural emergence of the continuous, one parameter family of global modes and vacua. We once again emphasise that a) introducing the parametrisation θ RL was mandatory in order to preserve the canonical structure of the theory in the most general way; b) this certainly makes our parametrisation qualitatively different from the so called dS α-vacua, e.g. [31] and c) it seems such mandatory one parameter freedom is a very unique property of the fermionic field theory in the hyperbolic dS spacetime, when compared to other field theories and other patches of dS (cf., discussions starting below Eq. (32)). Second, we have seen in Section 4.2, Section 4.3 that a) θ RL = 0, π/2 give thermal spectrum for the created massless particles and b) for these two values, plots of both the entropies are not only coincident for all values of field's mass, but also the entropies are maximal compared to other values of the parametrisation, Fig. 2, Fig. 4a, Fig. 4b, Fig. 5. However to the best of our understanding, certainly this does not implies that these two θ RL values should make the preferable choices of the global vacuum state. In fact its actual value should only be determined via comparison of theoretical prediction of some observable with the data.
In other words, it thus seems highly interesting to investigate the effect of θ RL into the other measures of quantum correlations, e.g., the entanglement negativity, the violation of Bell's inequality and the quantum discord etc (see e.g., [15,18,22] and references therein for discussions with scalar field theory) for the fermions. Moreover, we have seen in Fig. 2 -Fig. 5 that parameter values around θ RL = π/3 are qualitatively distinguished from other values. Perhaps these studies would then reveal further interesting properties of this parametrisation (including whether we may at all distinguish the values θ RL = 0, π/2) and will shed some more light on this novel feature of fermions in the hyperbolic dS spacetime. We wish to address these issues in our future works.