Vacuum State of the Dirac Field in de Sitter Space and Entanglement Entropy

We compute the entanglement entropy of a free massive Dirac field between two causally disconnected open charts in de Sitter space. We first derive the Bunch-Davies vacuum mode functions of the Dirac field. We find there exists no supercurvature mode for the Dirac field. We then give the Bogoliubov transformation between the Bunch-Davies vacuum and the open chart vacua that makes the reduced density matrix diagonal. We find that the Dirac field becomes more entangled than a scalar field as $m^2/H^2$ becomes small, and the difference is maximal in the massless limit.


Introduction
Since when Einstein-Podolsky-Rosen (EPR) pointed out in 1935, quantum entanglement has fascinated many physicists because of its counterintuitive nature that a local measurement on a particle may affect the outcome of a local measurement on a distant particle instantaneously [1]. After Aspect et al. convincingly tested the quantum nature of entanglement by measuring correlations of linear polarizations of pairs of photons [2,3], much attention has been paid to this genuin quantum property in various research areas including quantum information theory, quantum communication, quantum cryptography, quantum teleportation and quantum computation.
Turning our eyes on cosmology, in de Sitter space where the universe expands exponentially, any two of mutually separated region eventually becomes causally disconnected. This is most conveniently described by spanning open universe coordinates on two open charts in de Sitter space. The positive frequency mode functions of a free massive scalar field for the Euclidean vacuum (the Bunch-Davies vacuum) that have support on both regions were derived in [4]. Using them, several studies have been made on the quantum entanglement, particularly, entanglement entropy [5], on negativity [6], which is a measure of entanglement for any mixed states involved in subsystems, and on quantum discord [7].
Quantum entanglement between two causally disconnected regions in de Sitter space was first studied by Maldacena and Pimentel [5]. They showed that the entanglement entropy, which is a measure of quantum entanglement, of a free massive scalar field between two disconnected open charts is non-vanishing. Motivated by this result, there have been several attempts to test the idea of multiverse by studying long range correlations of various states that quantum entanglement naturally gives rise to [8,9,10].
In order to gain some insight into relativistic quantum information, qunatum entanglement between causally disconnected regions in flat space was investigated by Fuentes-Schuller and Mann [11] by making use of the Rindler coordinates. They studied the entanglement between two causally disconnected modes of a free scalar field as viewed by two relatively accelerated observers in Rindler space. It was found that the entanglement is degraded from the perspective of accelerated observers in flat space and, in particular, the entanglement disappears for an infinitely accelerated observer. Interestingly, however, Alsing et al. showed that unlike bosonic fields fermionic fields always remain entangled even in the limit of infinite acceleration [12]. Then Datta showed that quantum discord, which is a measure of all quantum correlations including quantum entanglement, never disappears in this limit [13].
The two open charts of de Sitter space are analogous to the Rindler wedges in flat space in the sense that an observer in the region described by one of the charts has no access to the field modes in the other causally disconnected region. Therefore it is of interest to see how the spacetime curvature will affect these results obtained in flat space.
In this paper, we first study the Bunch-Davies vacuum of the Dirac field in open charts by extending the previous work on scalar fields [4]. Then using thus obtained spinor mode functions, we compute the entanglement entropy.
The paper is organized as follows. In section 2, we consider the Dirac equation and its conserved currents in the open chart. In section 3, we derive the mode functions in each of

The Dirac equation in open chart
The open de Sitter space in the case of a free massive scalar field is studied in detail in [4]. In Figure 1, the de Sitter space and its Penrose diagram are depicted. We extend the study [4] to the case of a Dirac field in this section. We do not include dynamical gravity as in [4].

The metric in the open chart is
where the indices (µ , ν) run from 1 to 4 and H −1 is the curvature radius of the de Sitter space. The second line is the metric in the tetrad system and the indices (A, B) and (a, b) run from 1 to 4 and 1 to 3, respectively. The relation between curved and flat metrics is given by g µν = η AB e µ A e ν B where the tetrad field satisfies e A = e µ A dx µ . We defined e a = sinh tẽ a .
We choose 4 × 4 gamma matrices in accordance with our metric signature (−, +, +, +) such that where I andγ a are the 2 × 2 unit matrix and gamma matrices respectively. The gamma matrices satisfy γ A , γ B = 2η AB .
The Dirac equation for a four dimensional massive spinor Ψ is then expressed as Here commutators are Σ AB = 1 4 γ A , γ B and the spin connection is given by If we use Eq. (2.1), the above equation becomes where a prime denotes the derivative with respect to the time t, f (t) = sinh t, m is the mass of Ψ and/ ∇ =γ a∇ a . If we define 1 where Ω is the three-dimensional angle, then the Dirac equation becomes Combining Eqs. (2.6) with (2.7), the equation for φ + is given by Once we obtain a solution for φ + , we can get a solution for φ − by using Eq. (2.7).
The Dirac equation gives its conserved Noether current expressed as Here, the Dirac adjoint is defined byΨ ≡ Ψ † γ 0 . Then the charge density (µ = 0) is given by where we used (γ 0 ) 2 = −1. We use the charge density for mode functions to have orthonormality relations below.

The mode function in each R or L region
In this section, we consider the mode functions in the region R or L. Since these two regions are completely symmetric, the argument in this section can be applied to both R and L regions although we don't specify the region R or L below.

Positive frequency mode
We can separate variables according to [14]: where the three-dimensional spinors χ Once we obtain a solution for φ p (t), we can get a solution for φ − (t, Ω) = ϕ p (t)χ The positive frequency solution Ψ + is found to be where we defined z ≡ cosh t and the subscript ↑ indicates spin-up, where p > 0. Note that φ p and ϕ p realize positive frequency in the distant past t → 0 (z → 1).
This solution is normalized such that where the Dirac inner product for the mode function is given by at z → 1.

Negative frequency mode
The negative frequency mode function of Eq. (2.8) is given by separating variables as in [14] φ where the three-dimensional spinors χ (+) p m (Ω) have two components and satisfỹ / ∇χ where p is positive and continuous and the spinors χ Thus, the negative frequency mode function is obtained by replacing Eqs. (3.6) and (3.7) by p → −p. We find where p > 0 and * denotes complex conjugation. Here φ * p and −ϕ * p realize negative frequency at t → 0 (z → 1). This solution is normalized in the same way as Eq. (3.8).
Note that the negative frequency mode is not simply the complex conjugate of the positive frequency mode but a minus sign is necessary for the lower 2 components.

Spin-down solutions for the positive and negative mode functions
If we take negative p and interchange φ + and φ − in Eq. (2.5), we find the Dirac equation (2.4) does not change. Such solutions correspond to spin-down solutions.
The spin-down solution corresponding to the positive frequency mode Ψ + ↑ in Eq. (3.5) is given by .
(3.15) Note that the above solution is not only exchanging φ and ϕ in Eq. (3.5) but also χ p m . Similary, the spin-down solution for the negative frequency mode Ψ − ↑ in Eq. (3.10) is .
The general solution is expressed as a linear combination of the spin-up and down positive and negative frequency mode functions obtained in this section.

Analytic Continuation
In this section, we pick up the positive frequency mode functions corresponding to the Euclidean vacuum (the Bunch-Davies vacuum). To do this, we require that they are analytic when they are continued to the lower hemisphere of the Euclidean de Sitter space.

The relation between the Lorentzian and the Euclidean coordinates
The open chart is obtained by analytic continuation of the Euclidean sphere S 4 . The relation between the Lorentzian coordinate of three parts and the Euclidean coordinate is as follows.

The analytic continuation of the positive frequency mode function
We first ignore the normalization and focus on the time dependent part of Eq. (3.5) in region where the argument z R means that Ψ has support only in region R and vanish in region L. In order to construct positive frequency mode functions corresponding to the Euclidean vacuum (the Bunch-Davies vacuum), we need to extend Ψ analytically from region R to region L. The analytic continuation from the region R to L with the branch cut [−1, 1] goes through from Imz < 0 to Imz > 0.
To do it, we first extend Ψ Then we change the variable by z R = −z L in the C region. Since the hypergeometric function becomes singular at z L = 1, we use some formulae in Appendix A. Then we extend Ψ where we have defined . (4.6) Note that Eq. (4.5) has support only in region L and vanishes in region R due to the causally disconnected nature of region R and L. Note also that for the case of a massless Dirac field (m = 0), A = 0 and B = e −πp . In the limit of p = 0, A = sinh m H π and B = cosh m H π. From Eq. (4.5), we find that the solution in the region L can then be given by 2 where the argument z L means Ψ has support only in region L and vanishes in region R.
Thus if we extend Ψ analytically from region L to region R, the above solution becomes .
(4.8) 2 One may think that the solution could be another candidate solution in the region L.
However, the solution does not satisfy the orthonormality relation.
By using the solutions given in Eqs. (4.4), (4.5), (4.7) and (4.8), we find that the above solution satisfies the orthonormality relation in the far past z → 1.

The analytic continuation of the negative frequency mode function
We write the time dependent part of Eq. (3.10) as . (4.10) We want to extend Ψ −(R) ↑ analytically from region R to region L. The analytic continuation from the region R to L with the branch cut [−1, 1] goes through from Imz > 0 to Imz < 0.
The procedure is the same as what we did for the positive frequency mode. The solution extended analytically from the region R to L is then give by where A and B are given in Eq. (4.6).
The solution in the region L then can become . (4.12) If we extend Ψ −(L) ↑ from region L to R analytically, the above solution becomes . (4.13) The above solution satisfies the orthonormality relation in the far past z → 1.

The Bunch-Davies vacuum solutions
In the previous section, we found that mode functions obtained in one region are analytically continued to the other region by φ → Aϕ + Bφ * and ϕ → −Aφ + Bϕ * . In this section, we present the Bunch-Davies vacuum solutions of positive and negative frequency mode functions.

Positive frequency mode
The positive frequency mode functions that have support on both region R and L are found The normalization factor N b is computed by (5.5) Note that for a massless Dirac field, N b = 1 + e −2πp .

Negative frequency mode
The negative frequency mode functions that have support on both region R and L are found The normalization factor N b is given by Eq. (5.5).

Supercurvature modes
In this section, we examine if there is a supercurvature mode for the Dirac field. The supercurvature modes exist in the open chart for a massive scalar field [4]. They are known to be able to carry information about the pre-tunneling vacuum state and to be expected that they may explain the dipolar statistical anisotropy [15,16,17]. For a massive vector field, the paper [18] concluded that there is no supercurvature mode in the open chart.
We examine if the mode functions (3.6) and (3.7) for p = ip is normalizable in C region (see Figure 1). The metric in the C region is Note that the volume element is given by √ h = H −3 cos 2 t C cosh 2 r C . By using Eqs. (4.1) and (4.2), we find the relation z R = cosh t R = sin t C , and z R = ±1 correspond to t C = ±π/2. The mode function Ψ †(R) ↑ becomes real in the C region except for the irrelevant overall phase when p is pure imaginary. Since the hypergeometric functions in Eqs. (3.6) and (3.7) converge to unity at z = 1, we can easily read the asymptotic behavior of the mode function Ψ †(R) ↑ at t C = ±π/2 from Eqs. (4.4) and (4.5). If we focus on the behavior around t C = π/2, i.e. z R = 1, the more singular is the upper component proportional to Since the normalization is given by the contribution from the vicinity of t C = π/2 converges if ip is negative. The behavior on the other boundary of the region C can be read from (4.5). For ip < 0, the most singular piece in (4.5) at t C = −π/2 is proportional to Bφ * p (− sin t C ). Thus as long as B does not vanish, the above integral that determines the normalization diverges. Since B never vanishes as easily seen from Eq. (4.6), we conclude that there is no normalizable mode for p 2 < 0. Thus there exists no supercurvature mode.

Entanglement entropy
In this section, we quantize the Dirac field in the open chart in order to discuss quantum entanglement. In order to discuss quantum entanglement from the point of view of, say R region, we need to trace out the degree of freedom of inaccessible L region. Thus we need to change basis to mode functions that have support on either R or L regions. Thus, we consider the Bogoliubov transformation between the Bunch-Davies vacuum and R, L vacua.
We then derive the reduced density matrix in the region R by tracing over the region L and compute the entanglement entropy with the reduced density matrix.

Canonical quantization
If we take into account the two sets of modes in each R and L region, the Dirac field can be where σ = R or L and we recovered abbreviated subscripts such as p, , m temporally but they are omitted below for simplicity unless there may be any confusion. All other anticommutators vanish. If we focus on a mode p, the Bunch-Davies vacuum is defined by where the + superscript on the ket indicates the spin-down particle and spin-up antiparticle vacua and the − superscript for the spin-up particle and spin-down antiparticle vacua respectively, so that where σ = R and L. Since (a † s ) 2 = (b † s ) 2 = 0, only two states are allowed for particles and antiparticles.

The Bogoliubov transformation
We denote the spin-up (3.5) and down (3.15) positive frequency mode functions in region R by ψ Similarly for the spin-up (4.7) and down (5.4) mode functions in region L, For the negative frequency mode functions, we express the spin-up (3.10) and down (3.16) solutions by and for the spin-up (4.12) and down (5.9) negative frequency mode functions in region L we write In region R, let us denote (c R s , c R † s ) as the annihilation and creation operators for fermions and (d R s , d R † s ) as the annihilation and creation operators for antifermions. The corresponding fermion and antifermion operators in region L are denoted as (c L s , c L † s ) and (d L s , d L † s ). Then the Dirac field can be expanded as These obey the anticommutation relations where σ = R or L. All other anticommutators vanish 3 .
As the Dirac field operator should be the same under this change of basis, we can relate the creation and annihilation operators in the Bunch-Davies vacuum to those in R, L vacua by comparing Eq. (7.1) with (7.9) It follows that We have another set of four relations with their spin-up and down exchanged but it is totally equivalent. Thus, we focus on the relation (7.11) below. Here M is a 4 × 4 matrix where α and β are 2 × 2 matrices respectively defined as Note that A, B and N b are given in Eqs. (4.6) and (5.5). Then, we find where 15) and the components of the matrix M −1 are computed as Since a σ ↓ mixes the spin-down particles and the spin-up antiparticles in region R and L, the Bunch-Davies vacuum for mode p can be regarded as the Bogoliubov transformation of R, L vacua of the form where m ij is an antisymmetric matrix and operators c i and d j satisfy the anticommutation Here the normalization of the Bogoliubov transformation is omitted because we will need another Bogoliubov transformation in Eq. (7.22). Note that we drop spin labels on operators for simplicity because we focus on operators c ↓ and d ↑ as defined by If we apply a σ ↓ in Eq. (7.14) to Eq. (7.17), we have where σ = R and L, and we used the anticommutation relation Eq. (7.2). Using the expressions in Eqs. (7.13) and (7.16), we find where A and B are given in Eq. (4.6). Since A = 0 in the case of masslessness, the density matrix ρ = |0 + BD + BD 0| in terms of Eq. (7.17) becomes diagonal in the |0 + R |0 − L basis. In other cases, however, it is not diagonal and then it is difficult to trace out the L degrees of freedom. Thus, we introduce new operatorsc R andd L and perform a further Bogoliubov with |u| 2 + |v| 2 = 1, so that we obtain the form where a plus sign in front of the square root term is taken to satisfy |γ p | < 1. Note that the second condition of Eq. (7.24) leads to the relations for u * /v * . We find u * /v * = −u/v and we obtain the same γ p in Eq. (7.27). The condition Eq. (7.25) gives the same γ p in Eq. (7.27) and the relation u * /v * = −u/v as well. Note also that in the limit of p = 0, ω = −A/B, ζ = 1/B and ω 2 + ζ 2 = 1 holds, thus we find γ p = 1 independently of m.

The density matrix
By using Eqs. (7.22) and (7.27), the reduced density matrix in region R is then found to be where the conservation of probability holds, Trρ R = 1. ν 2 S(ν ) Figure 2: Entanglement entropy between two causally disconnected regions for one degree of freedom of a fermion (red) and a scalar field (blue).
As an observer in the region R will observe particles defined by the operatorsc R , the expected number of such particles will be given by where for notational convenience we have defined |1 c 1 d In the case of masslessness, this is expressed by (e 2πp + 1) −1 , which is a thermal state with temperature The entanglement entropy for each mode is calculated to be (7.31) The final entanglement entropy per unit comoving volume between two causally disconnected regions are obtained by integrating over p and a volume integral over the hyperboloid H 3 .
That is, we use the density of states on the hyperboloid [5] S(m) = 2π ∞ 0 dp D(p)S(p, m) .
The density of states for H 3 in the case of the Dirac field is D(p) = ( 1 4 + p 2 )/(2π 2 ) [14,19]. The result is plotted in red line of Figure 2 where ν is defined in Eq. (7.35). We plotted the entanglement entropy for one degree of freedom for comparison with a scalar field, that is, the plot is 1/2 of Eq. (7.32) 4 . Now let us compare the result with the entanglement entropy of a scalar field which is computed by [5] and expressed as where γ is given by cosh 2πp + cos 2πν + √ cosh 2πp + cos 2πν + 2 , (7.34) and a mass parameter is defined by In the case of a massless scalar field (ν = 3/2), we find γ = e −πp , and then the reduced density matrix is found to be thermal ∞ n=0 e −2πpn |n n| , (7.36) with temperature T = H/(2π).
We integrate over p as in Eq. (7.32). The density of states for H 3 in the case of the scalar field is D(p) = p 2 /(2π 2 ) [19]. The result of the scalar field is plotted in blue line in Figure 2.
In Figure 2, we see that the Dirac field gets more entangled than the scalar field as m 2 /H 2 becomes small, and the difference is maximal in the massless limit.

Fermion seems more entangled than scalar in the massless limit
For the massless Dirac field (ν 2 = 9/4), γ p = e −πp which becomes 1 in the limit of p = 0.
Then the entanglement entropy of the Dirac field per each degree of freedom, Eq. (7.31), becomes log 2. 5 Since the density of states of the Dirac field is finite even in the limit of p = 0, 4 As explained below Eq. (7.11) we focus on two degrees of freedom of a fermion since two other degrees of freedom satisfy the same relations. 5 We consider two degrees of freedom of a fermion now as explained below Eq. (7.11) the final entanglement entropy Eq. (7.32) on large scales is finite. For a massless scalar field, on the other hand, the entanglement entropy Eq. (7.33) becomes logarithmically infinite in the limit of p = 0. But the density of states of the scalar field becomes quadratically zero.
Then the entanglement entropy summing over p gives zero in the limit of p = 0. Thus the contribution of the states from large scales to the entanglement entropy becomes large for the massless Dirac field compared to the case of the massless scalar field.
However, one should note that the scalar field entanglement entropy shows a strange behavior as the mass decreases. For m 2 < 2H 2 (ν 2 > 1/4), the entanglement entropy once decreases as the mass decreases. It is known that there exits a supercurvature mode for m 2 < 2H 2 , and its contribution to the long-distance correlation becomes more and more dominant as m 2 /H 2 → 0 [4]. At the moment we have no clue about how the supercurvature mode affects the entanglement. It seems possible that the contribution of the supercurvatue mode, if it could ever be computed, would dominate the entanglement entropy of the scalar field in the small mass limit. This issue needs further studies before we make a firm conclusion.

Summary
We studied the entanglement entropy of a free massive Dirac field between two causally disconnected open charts in de Sitter space. For this purpose, we first derived the Bunch-Davies vacuum mode functions of the Dirac field in the coordinates that respect the open chart. We then gave the Bogoliubov transformation between the Bunch-Davies vacuum and the open chart vacua that makes the reduced density matrix diagonal. We derived the reduced density matrix in one of the open charts (R region) after tracing out the other (L region) and found that the Fermi-Dirac distribution is realized in the limit of masslessness.
We then computed the entanglement entropy of the Dirac field by using the reduced density matrix. We compared the entanglement entropy of one degree of freedom of the Dirac field with that of a scalar field calculated by [5].
We found that the entanglement entropy for the Dirac field gets more entangled than that for a scalar field as m 2 /H 2 becomes small, and the difference is maximal in the massless limit. This is because the contribution of the states from large scales to the entanglement entropy becomes large for the massless Dirac field compared to the case of the massless scalar field. But there is a caveat. In the computation of the entanglement entropy of a scalar field, it is assumed that the supercurvature mode does not contribute. If this assumption were wrong, the entanglement entropy of the scalar field might become larger than that of the Dirac fermion in the small mass limit. In connection with this issue, we also showed that there is no supercurvature mode for the Dirac field.