Entangled de Sitter from Stringy Axionic Bell pair I: An analysis using Bunch Davies vacuum

In this work, we study the quantum entanglement and compute entanglement entropy in de Sitter space for a bipartite quantum field theory driven by axion originating from ${\bf Type~ IIB}$ string compactification on a Calabi Yau three fold (${\bf CY^3}$) and in presence of ${\bf NS5}$ brane. For this compuation, we consider a spherical surface ${\bf S}^2$, which divide the spatial slice of de Sitter (${\bf dS_4}$) into exterior and interior sub regions. We also consider the initial choice of vaccum to be Bunch Davies state. First we derive the solution of the wave function of axion in a hyperbolic open chart by constructing a suitable basis for Bunch Davies vacuum state using Bogoliubov transformation. We then, derive the expression for density matrix by tracing over the exterior region. This allows us to compute entanglement entropy and R$\acute{e}$nyi entropy in $3+1$ dimension. Further we quantify the UV finite contribution of entanglement entropy which contain the physics of long range quantum correlations of our expanding universe. Finally, our analysis compliments the necessary condition for the violation of Bell's inequality in primordial cosmology due to the non vanishing entanglement entropy for axionic Bell pair.


Introduction
The information encoded in entanglement entropy is very powerful to distinguish quantum states with long range correlation as appearing in the context of condensed matter physics [1,3,58]. It plays a significant role in the field of quantum information [4][5][6][7][8] and cosmoloogy [9][10][11][12][13][14][15][16][17][18][19][20][21][22][23][24][25]. It is also useful in quantum field theory to identify the specific nature of the long range quantum correlations that we generally get using the standard vacuum states like the Chernikov-Tagirov, Bunch-Davies, Hartle-Hawking or the most general α vacuum [26][27][28]. Till date quantum entanglement has been one of the most mysterious and at the same time fascinating features of foundational aspects of quantum mechanics. This is because, it is believed that performing a local measurement may instantaneously affect the outcome of the measurement beyond the physical light cone under consideration. This is interpreted as the violation of causality and commonly known as the Einstein-Podolsky-Rosen (EPR) paradox [29]. However, it is important to note that neither quantum information gets transferred in such a physical measurement, nor it changes the causality constraints. Nevertheless, quantum entanglement may have applications to explain the origin of various physical phenomena [30][31][32]. The prime example is the well known Schwinger effect in de Sitter [33,34] space describing pair production of particles in presence of uniform electric field [19]. In this specific context, the pair production occurs spontaneously with a certain finite coordinate separation in field space. Such physical states of the particle pairs are required to be correlated quantum mechanically.
In the present work we consider an axion field theory obtained from compactification of Type IIB string theory on a Calabi-Yau three fold (CY 3 ) but in presence of NS5 brane. So that the model consist of effective potential which breaks the shift symmetry associated with the axion field in a controled manner. This model has been analyzed as a candidate for driving inflation [35][36][37] as well as a quintessence model for the late time acceleration of our universe [38]. Furthermore, in a recent analysis it was shown that this model captures an interesting phenomena of violation of Bell's inequality in primordial cosmology [16][17][18] by formation of EPR pair. In the present analysis we will strengthen the formation of EPR pair by demonstrating that the entanglement entropy of axionic EPR pair is non zero, which confirms the existence of superhorizon long range correlation in primordial cosmology. This connection also will be helpful in future to provide an algorithm to contrast various models of inflation [40][41][42][43]. Additionally, we we compute Rényi entropy in 3 + 1 D de Sitter space and its connection with quantum entanglement in presence of axionic model.
In fig. (1), we have schematically shown the logic involved in computing the entanglement entropy in de Sitter space. The plan of the rest of the paper is as follows. In Section 2 we briefly review the strategy for computing the entanglement entropy in de Sitter space. In Section 3 we follow the above strategy for the afore mentioned stringy axionic model. In Section 4 we present our conclusion and future prospects. Finally, in Appendix we give the detailed derivation of entanglement entropy and Rényi entropy in the present context.

Computational strategy: Brief review
Here we briefly review the method to derive entanglement entropy in 3+1 dimensional de Sitter space following ref. [15] and references therein. We consider a closed surface Σ in a hypersurface where time is frozen. Consequently, the space-like hypersurface is divided into an interior and exterior region which are identified as RI (which is L in open chart) and RII (which is R in open chart) region in the representative schematic diagram shown in fig. (2). If we represent the total Hilbert space of the system under consideration by H, then in the present context one can use the following approximate bipartite decomposition [44], given by 1 , H = H INT ⊗ H EXT . This implies that H can be written as a direct product of two Hilbert spaces H INT and H EXT which is consistent with the requirement of local QFT. In this context, the building blocks of H INT and H EXT are the localized modes in the RI and RII region respectively. This allows us to construct the density matrix for the internal RI region by tracing over all degrees of freedom in the external RII region and is given by: (2.1) Here the vacuum state |BD is the Bunch-Davies vacuum. Using the well known Von Neumann entropy formula, the entanglement entropy in de Sitter space can be expressed in terms of the density matrix as: which quantifies how much a given quantum state is quantum mechanically entangled. In this paper we compute this expression explicitly and establish its relationship with Bell's inequality violation in cosmology [16][17][18].
Here we consider Σ to be a closed surface S 2 which has radius R S 2 with the restriction, R S 2 >> R dS = H −1 , where H is the Hubble parameter and R dS is the de Sitter radius. This surface S 2 separates the spatial slice into the interior (L) and exterior (R) region which is shown in fig. (2). The reduced density matrix, which is a key ingredient for computing entanglement entropy, is obtained by tracing over the exterior (R) region. Also it is important to note that the total entanglement entropy can be expressed as a sum of UV divergent and UV finite contribution as: When the S 2 is taken at the boundary of 3 + 1 D de Sitter space then SO(1, 3) symmetry 2 plays significant role to define the entanglement entropy for a given principal quantum number. We note that entanglement entropy in de Sitter space is completely different from that of the dynamical gravity [45][46][47][48][49][50][51][52][53][54], which we are not considering. In 3 + 1 dimensions, the UV-divergent contribution is appearing due to the local quantum field theoretic effects and it takes the following simplified form [15,19,30,31]: where UV is the short distance lattice UV cut-off of the bipartite local quantum field theory under consideration, A ENT represents the proper area of the entangling region and c i ∀i = 1, 2, 3, 4 are the numerical coefficients which carries different physical significance in the present computation. For an example, in Eqn (2.4), the first term represent area dependent contribution to the UV divergent part of the entropy, which is originated from entanglement between particles in the vicinity of the entangling surface under consideration. Further taking flat space limit, where the Hubble parameter H → 0, the UV-divergent part of the entropy of the bipartite local quantum field theoretic system for scalar field (axion) can be expressed as 3 :

7)
2 In arbitrary D dimensions entanglement entropy in de Sitter space is invariant under SO(1, D) isometry group. 3 Here we in the flat space limit H → 0 we get: lim H→0 ln ( UV H) = 2 ln ( UV ) + ∆, (2.5) where ∆ characterizes the following finite contribution in the flat space limit, which is given by: In this context the coefficients c 2 and c 3 contribute in the flat limiting result of the entropy of the system under consideration. Also it is important to note that, the coefficients c 2 and c 3 are short distance lattice UV cut-off independent and therefore universal in nature. On the other hand, the last term in Eqn (2.4), explicitly incorporates the effect of bulk extrinsic and intrinsic curvature contribution in de Sitter background. In our computation of entanglement entropy in de Sitter space we only restrict our discussion within the UV-finite contributions which contain the physics of long range correlations of the quantum mechanical state. Most important fact is that, the total entanglement entropy in de Sitter space is invariant under SO(1, 4) isometry group and consequently both UV-divergent and UVfinite contributions remain invariant under the same group. However, for the computation we only use the connected subgroup of SO (1,4), which is SO (1,3) in this context. For more details see ref. [15]. Additionally, in the present context in the subhorizon scale, where kc S η >> −1 4 , plays crucial role as the long range entanglement is valid within this scale. In the superhorizon scale, where kc S η << −1 or more precisely at late time scale i.e. at the limiting case η → 0, long range contributions in the de Sitter entanglement entropy freezes. However, in this limit some additional contribution appears which only contributes to short distance scale and the associated entanglement can be expressed in a local form in (quasi) de Sitter space. Consequently, at late time scale, η → 0 of the UV-finite contributions in the entanglement entropy can be written as [15,19]: Here the UV finite part gets contributions from two different sources. First part being proportional to the area of entangling surface of S 2 and the proportionality factor is dependent on a overall coefficient and Hubble parameter H. Second part is proportional to the logarithm of the product of area of S 2 and Hubble parameter H and the proportionality factor is identified to be the interesting part of the entanglement entropy which we compute explicitly from the prescribed setup.
Additionally it is important to note the overall characteristic coefficient as appearing in the first term in UV finite contribution is not completely independent. Once we know the coefficient of the logarithmic term, one can determine the other coefficient in terms of this contribution. Further taking the flat space limit, H → 0, the UV-finite contributions to the entanglement entropy can be expressed as 5 : Here k is the momentum scale of the corresponding Fourier modes in De Sitter space, cS is the speed of sound and η is the conformal time. 5 Here we use the flat space limit H → 0 we get: where Θ characterizes the following contribution in the flat space limit, which is given by: which is in general divergent but can be regularized by putting the previously introduced lattice UV cut-off UV as: where H UV is the harmonic number of order UV .
Further, we introduce a new parameter A C representing the entangling area in comoving coordinates, which is defined as, A C = A ENT H 2 η 2 . In terms of newly introduced parameter the UV-finite part of the entanglement entropy can be written as: where the co-efficient of the logarithmic contribution, c 6 , characterizes the long range correlations. In our computation our prime objective is to compute the expression for c 6 in de Sitter space. On general ground, c 6 can be expressed as a sum of two possible contributions of the extrinsic curvature of the entangling surface (Σ ≡ S 2 ) in comoving coordinate system as [55]: where g 1 and g 2 are two undetermined co-efficients and K AB is the extrinsic curvature of entangling surface. Here we note that, the co-efficient c 6 is identified as the interesting part of the UV − finite contribution to the entanglement entropy, as represented by S intr , which we compute in this paper.
In the next section,we compute the entanglement entropy for axionic Bell pairs using the results of this section.  Figure 2. Schematic diagram for the geometrical construction of the bipartite system. First we consider S 2 with radius R >> R dS = H −1 , which at late global time τ global can mapped to half of a S 3 with boundary S 2 at the equator. For simplicity we can also use hyperbolic slices where the interior of the S 3 mimics the role of the "left" (L) slices. Corresponding Penrose diagrams are drawn for completeness [15].

Entanglement for axionic Bell pair
Let us first discuss about the underlying symmetries of the present setup in de Sitter space. For this let us consider, a surface S 2 defined by the relation, where R S 2 is the radius of the spherical surface S 2 and x i ∀i = 1, 2, 3 are the flat coordinates. To consider a spherical surface sufficiently large compared to the size of the horizon one need to consider here R S 2 >> η, where η is the conformal time defined as, η = dt/a(t). Here a(t) is the scale factor which is appearing in the spatially flat k = 0 FLRW metric in 3 + 1 dimensions, written in the conformally flat form as: (3.1) We note that for R S 2 >> η one can suppress the conformal time coordinate and we achieve this by fixing the value of R S 2 in the η → 0 limit. This surface in this context is invariant by an SO(1, 3) subgroup, which is the SO(1, 4) quasi de Sitter isometry group in left region. Technically this is commonly known as left invariant in which we perform our rest of the computation presented in this paper. Also for our purpose we use de Sitter global coordinate system in which the equal time slices are three sphere S 3 . In this context the entangling surface is two sphere S 2 , which is the equator of three sphere S 3 . Further using quasi de Sitter isometry group as η = 0 one can consider a map which transform S 2 to the equator of three sphere S 3 on the boundary of 3 + 1 D quasi de Sitter space. However, this gives rise to divergent contributions, which can be regularized by transforming back to the S 2 to a surface at global time τ global .
In the present situation we also use hyperbolic slices which are geometrically characterized by the equation, 5 j=1 Y 2 j = H −2 , where for Euclidean signature we can write down the following parametric equations: cos τ E cos ρ E for j = 5.

(3.2)
Heren j ∀j = 1, 2, 3 are the j-th component of the unit vector defined in R 3 .
In this hyperbolic slices the metric with Euclidean signature can be expressed as: where dΩ 2 2 is defined in the two sphere S 2 . Further we perform the following analytic continuation in the fifth coordinate: and redefine the rest of the coordinates as, X k = Y k ∀k = 1, 2, 3, 4. Consequently, in the Lorentzian signature one can write the following geometrical equation: On the other hand, in the Lorentzian signature one can consider three different region, which are characterized as: R : (3.7) L : Further, in hyperbolic slices the metric with Lorentzian signature can be expressed for the three regions as: In fig. (2) schematic diagram for the geometrical construction and underlying symmetries of the bipartite system.

Wave function of axion in an open chart
In this section our prime objective is to compute the expression for entanglement entropy in de Sitter space driven by an axion field. This axion field originates from the R-R sector of Type IIB string theory compactified on a Calabi-Yau three fold (CY 3 ) in presence of NS 5 brane. For details of the compactification and the resulting effective action, in 3 + 1 dimensions, for the axion field, see refs. [35][36][37][38][39]. Using this effective action we discuss that such stringy axions can be treated as necessary ingredient to study quantum entanglement in the context of primordial cosmology 6 . Let us start with the aforementioned effective action for axion field: where φ is the axion field and the corresponding effective potential can be expressed as: 13) 6 From the view of quantum information theory violation of Bell's inequality is more appropriately connected to non-locality. Violation of Bells inequality (CHSH) in bipartite quantum systems necessarily requires quantum entanglement, but the reverse statement is not always true [56,57]. In a bipartite quantum system described by a specific class of state, for instance say for the Werner states with non vanishing entanglement entropy can be easily connected to the violation in Bell's inequality [58]. However, without complete understanding of the quantum state being considered in such analysis, it is not possible in general to determine exactly the extent of Bell's inequality violation.
where µ 3 is the mass scale, f a is the axion decay constant and we have defined a new parameter b In this case scale Λ G is given by, Λ G = m SU SY L 3 / √ α g s e −cS inst , where S inst is the instanton action, c ∼ O(1) is a constant factor, m SU SY is the supersymmetry breaking scale, α is the Regge slope parameter, g s is the string coupling constant and L 6 is the volume factor in string units. Using this action we compute the wave function of axion in open chart. In fig. (3) we have shown the schematic behaviour of the axion potential. For the computation of the wave function here we consider the following two possibilities:

Case I:
In the first approximation of the computation we restrict ourself up to the linear term of the effective potential: (3.14) Here the linear term of the effective potential can be treated as a source term in the equation of motion where the axion has no mass.

Case II:
On the other hand, in the small field limit φ << f a , where analytic solution of the equation of motion is possible, one can approximate the non-perturbative contribution as, cos φ fa ≈ . Consequently in small field limit φ << f a one can approximate the total effective potential for axion as: Here we define the effective mass of the axion as: where the axion decay consatnt in general follow a time dependent profile and for our purpose we use the following choice, f a = 100 − 80 1+ ln η ηc 2 H, which was used in refs. [16][17][18] to demonstrate Bell's inequality violation in primordial cosmology. It is important to mention here that, the one point function from axionic fluctuation in primordial cosmology is exactly proportional to the parameter m axion /f a and such contributions are apeparing in Case II, which we further use for the computation of entanglement entropy to prove the existance of EPR Bell pairs in cosmology.
Further using Eqn (3.12) one can write down the following equation of motion for the axion for the above mentioned two cases:

17)
For Case II : where 2 is the 3 + 1 D D'Alembertian operator which can be expressed in open chart as [59]: Also the exact time dependence in the scale factor a(t) for an open chart of the de Sitter background is given by: Additionally, the Laplacian operatorL 2 H 3 in hyperbolic slice H 3 can be written as [59]: which satisfies the following eigenvalue equation: Here Y plm (r, θ, Φ) are the eigenfunctions ofL 2 H 3 which can be written using method of separation of variable here one can decompose Y plm (r, θ, Φ) into a radial and angular part as: where Y lm (θ, Φ) is the spherical harmonic function and P The full solution of the equations of motions for the Case I and Case II can be expressed in the following forms: where U Λ (t, r, θ, Φ) forms a complete basis of mode function labeled by index Λ which is fixed by the three quantum numbers p, l and m and an index σ = ± for R and L hyperboloid in the present context.
It is important to note that equation of motions for the Case I and Case II are inhomogeneous differential equations with constant source and for this the total solutions of these equations can be expressed as: where the time dependent part of the wave function χ p,σ (t) forms a complete set of positive frequency function which can be written as a sum of complementary part (χ (c) p,σ (t)) and particular integral part (χ (p) p,σ (t)): Here our prime objective is to explicitly compute both of the contributions and use them further in the computation of entanglement entropy. Further, one can write down the following equation of motion for the complementary part of the time dependent contribution in the total solution as: where D t is the differential operator in the present context, which can be expressed for the above mentioned two cases as: for Case II. (3.28) The final solution for the complementary time dependent part for Case I and Case II are given by: Case I : for L.

(3.29)
Case II : for L, (3.30) where in Case II we introduce a parameter ν, which is defined for axion as: • In the result for Case II, in the limit m axion << H the parameter ν is given by, ν = 3/2. For this value of ν, the solutions for Case I and Case II coincide. For this reason one can treat Case I as a special situation of Case II.
• Similarly in the limit m axion = √ 2H we have ν = 1/2. This is an important class of solution where axion is conformally coupled. Using this class of solution one can get back the flat space limiting result of the entanglement entropy as the metric for de Sitter space is conformally flat.
• In the case where m axion < √ 2H parameter ν is lying within the range, 1/2 < ν < 3/2 and this is also an important class of solution which is commonly known as the solution in the low mass region.
• Further in the limit m axion = H parameter ν is given by, ν = 5/2 and this is a restrictive class of solution where axion mass is exactly equal to the Hubble friction.
• For m axion >> H, we have, ν = i m 2 axion H 2 − 9 4 ≈ im axion /H. This is an important class of solution in the present context corresponding to high mass. Note that, for √ 2H < m axion < 3H/2 the parameter ν is lying within the window, 0 < ν < 1/2.
• Additionally it is important to note that, the final solution for the complementary time dependent part is symmetric under the exchange of the sign of the quantum number p i.e. χ • In this context the overall normalization constant of the time dependent complementary part of the solution is fixed by the following Klien-Gordon inner product [59], χ p,σ (t) = N pσ δ σσ , where N pσ is the overall normalization constant, which is given by: for Case II. (3.32) Further, one can write down the following equation of motion for the particular integral part of the time dependent contribution in the total solution as: where D t is the differential operator as given in Eqn (3.28) and J = µ 3 . For Case II depending on the time dependence of the axion effective mass as well as the time dependent axion decay constant, source function can be time dependent or may be constant. But for general consideration we assume that the source function is time dependent. The final solution for the time dependent "particular integral" part is given by: where G σ (t, t ) is the Green's function for axion field, given by: pn,σ (t) represent the solution for the complementary part with slight modification due to the replacement of p by p n ∀ n = 0, · · · , ∞ and is represented by: Case I : for L.

(3.36)
Case II : for L. (3.37) Further, using the results obtained for the total solution of the EOM given in Eqn (3.24) we expand the field in terms of creation and annihilation operators in Bunch-Davies vacuum, which is given by: where the Bunch-Davies vacuum state is defined as: a σplm |BD = 0 ∀σ = (+1, −1); p = 0, · · · , ∞; l = 0, · · · , p − 1, m = −l, · · · , +l. (3.39) In this context, the Bunch-Davies mode function U σplm (r, t, θ, φ) can be expressed as: After substituting Eqn (3.40) in Eqn (3.38) we get the following expression for the wave function: Below we use the following sets of shorthand notations to denote various associated Legendre polynomials as appearing for the solutions of the wave functions: Complementary solution : for Case II.

(3.42)
Particular solution : for Case II. where the index q = R, L for the right and left handed hyperbolic region in the open chart. Thus, the total time dependent part of the solution can be simplified to the following form: where we use the following redefined symbol: Here p n ∀n = 0, · · · , ∞, σ = ±1 and q = R, L. Also the normalization constant N p for the complementary part is defined as, N p = 2 sinh πp N pσ , where N pσ is defined in Eqn. (3.32). For the particular part if we replace p by p n then the normalization constant N pn can be defined as, N pn = 2 sinh πp n N pnσ . In Eqn (3.44), expansion coefficient functions for complementary solution (α σ q ,β σ q ) for q = R, L are defined as: (3.47) Exapnsion coefficients for particular solution (ᾱ σ q ,β σ q ) can be obtained by replacing p by p n in Eqn. (3.46) aand Eqn. (3.47). Before going to the further details, let us first analyze Eqn (3.44): • The full solution is valid for p = p n ∀n = 0, · · · , ∞, which is the necessary condition to solve the inhomogeneous differential equation using Green's function technique.
• If the source is absent, then the solution for the complementary part is perfectly consistent with ref. [15].
Further using Eqn (3.44) we write the following matrix equation in a compact form as: where for the complementary part of the solution we define the following matrices: Similarly for the partucluar solution we also define following matrices: where σ = ±1, q = R, L and I, J = 1, 2, 3, 4.
On the other hand the redefined normalization constant for the particular part of the solution N p,(n) can be expressed as, N p,(n) = 2 sinh πp n N pnσ p 2 − p 2 n . Further using Eqn (3.48) the Bunch-Davies mode function can be written as: where a I = (a σ , a † σ ) represents a set of creation and annihilation operator. On the other hand we define: where a (c) σ,n ) are the set of creation and annihilation operators which act on the complementary and particular part respectively. Thus, the operator contribution for the total solution is: where by inverting Eqn (3.52) we have expressed: a (c) (3.54) We define the following inverse matrices: where σ = ±1, q = R, L and I, J = 1, 2, 3, 4. The entries of these inverse matrices are given by:  Similarlyγ jσ,n andδ * jσ,n can be obtained by replacing p by p n in Eqn. (3.56) and Eqn. (3.57). Further we introduce set of rules which are useful to differentiate the operations of these operators in the complementary and particular solution. These set of rules are: • Rule I: The operator corresponding to the complementary solution is insensitive to the particular • Rule II: The operator corresponding to the particular solution is insensitive to the complementary solution i.e. a (p) For simplification we quantify the components of a I = (a σ , a † σ ) for σ = ±1 as: Further the Bunch-Davies vacuum state can be expressed in terms of the direct product of R and L vacua by using Bogoliubov transformation as: Here we assume that the Hilbert space corresponding to the Bunch-Davies vacuum state H BD can be decomposed into two separable portions as, H BD := H R ⊗ H L , where H R and H L are the Hilbert space corresponding to R and L vacuum state. In the present context operatorÔ is given by: where the coefficients m ij andm ij,n will be determined later. In Eqn (3.60) we will keep only linear terms in the exponential for simplicity. Note that R and L vacuum states can be expressed as: with (c) and (p) representing the complementary and particular part respectively. Further in principle one can express the particular part as: Also the annihilation operator satisfy: as well as the following commutation relations: Further one can write the annihilation of Bunch-Davies vacuum in terms of the annihilations of R and L vacua as: where neglecting contribution from the higher powers of creation operators, A (q) s ∀s = 1, 2, 3, 4 are defined as: This implies that the following condition holds good: Thus, the complementary part and particular part vanish independently as the solutions are independent of each other. Consequently, we get the following constraints: Further using Eqn (3.74) and Eqn (3.75), the mass matrices corresponding to the complementary part and particular part can be expressed as: Next substituting the right hand side of the above mentioned equations explicitly the entries of the mass matrices can be expressed for i, j = R, L as: Entries of the matrixm ij,n can sismilarly written by replacing p by p n in Eqn. (3.77). Before further discussion here we point out few important features: • For the Case I we observe that for the complementary and particular part of the solution But for Case II we find that (3.79) which is non vanishing for 0 < ν < 3/2 and ν > 3/2. A special case appear for ν = 3/2, where the result vanishes and one can get back the result of Case I. Another important observation is that simultaneously it is not possible to fix p = 0, ν = 3/2 for complementary solution and p n = 0, ν = 3/2 for particular solution, as both of them give divergent contribution in the diagonal component of the mass matrix.
• For the Case I we see that for the complementary and particular part of the solution But for Case II we find that (3.81) Here also, Case II coincides with Case I for ν = 3/2. Additionally, the non vanishing off diagonal components for both the cases indicate the signature of quantum entanglement, which finally gives rise to a non vanishing entanglement entropy. In what follows, we will explore this possibility in detail.
• Creation operators of b oscillators can be redefined by absorbing the overall phase contribution e iθ .
• The eigenvalues for the complementary part and particular part of the solution are given by: Their explicit expressions are: Case I :

86)
λ ±,n =m RR,n ±m RL,n = e iθ √ 2 e −pnπ (cos πν ± i sinh p n π) √ cosh 2πp n + cos 2πν .    However, this is not a physical basis as it is extremely complicated to trace over the R and L contributions when the Bunch-Davies vacuum state is represented using Eqn (3.60) and Eqn (3.61).
In fig. (4(a)) and fig. (4(b)), we have shown the behaviour of the eigenvalue (λ ± ) vs the momentum p for fixed values of the mass parameter |ν| in the small and large axion mass limits. Here we cover the mass parameter range 0 < |ν| < 3/2 for both the cases. In fig. (4(a)) for the small mass limiting range for large values of the momentum p the two branches of solution for |ν| = 0 and |ν| = 1/2, 3/2 coincides with each other and in the small values of the momentum p the two branches of solution can be separately visualized. Exactly opposite situation appears when we consider the large mass limit in fig. (4(b)). Further in fig. (4(c)) and fig. (4(d)), we have explicitly shown the behaviour of the eigenvalue (λ ± ) vs mass parameter |ν| for the fixed values of the the momentum p in the small and large mass limits. For fixed value of p (p = 1 and p = 2) in the large mass limit, we initially get decaying behaviour and then after |ν| = 3/2 they saturate. On the other hand, for same values of p, in the large mass limit, we get oscillating behaviour. Combined effects for large and small mass limits are plotted in fig. (4(e)), where we have explicitly shown that for ν 2 < 0 (for p = 1 and p = 2) we get distingushable behaviour of the mass eigen values. Once ν 2 → 0 magnitude of the eigevalues increse and finally when ν 2 > 0, both of the plots show apeariodic oscillations.
To find a suitable basis where we can trace over all contributions from R and L region we need to perform another Bogoliubov transformation incorporating new sets of operators, which are given by: with the following constraints: Using the above set of operators one can write the Bunch-Davies vacuum state in terms of new Bogoliubov transformed basis represented by R and L vacuum state as: where we have introduced a new operatorQ, defined in the basis as: where γ p and Γ p,n are defined as the coefficients of the complementary and particular part of the operator, which play significant role to construct the density matrix as well as the entanglement entropy. In this context the exponential of the operator eQ is approximated as: The overall normalization factor N p is defined as: Here due to the second Bogoliubov transformation the direct product of the R and L vacuum state is connected with the direct product of the new R and L vacuum state as: Before going to the further details let us first mention the following useful commutation relations of the creation and annihilation operators of the R and L vacua as given by: In this context, the operations of creation and annihilation operators defined on the Bunch-Davies vacuum state are appended bellow: Further, one can express the annihilation operators after second Bogoliubov transformation in terms of the annihilation operator obtained after first Bogoliubov transformation as: Here the matrices G I J and G (n) I J are defined as: where the entries of the matrices are given by: Further using Eqn (3.88) and Eqn (3.89), in Eqn (3.99) and Eqn (3.100), we get the following sets of homogeneous equations for the two cases:

107)
For Case − II : Similarly for Case II for the diagonal and off diagonal terms we get: 3. In the Case II, if we consider γ p and Γ p,n are purely imaginary i.e. γ * p = −γ p , Γ * p,n = −Γ p,n , and if we set, v * =v, u * =ū, V * n =V n , U * n =Ū n . Consequently four sets of equation reduces to two sets of homogeneous equations for the two cases the normalization |u| 2 − |v| 2 = 1 and |U n | 2 − |V n | 2 = 1 are explicitly imposed for the complementary and particular part of the solution.
Finally, the non trivial solutions obtained from these system of equations for the two cases can be expressed as:   Further substituting the entries of the mass matrices one can further simplify the result as: In fig. (5(a)) and fig. (5(b)), we have explicitly shown the behaviour of |γ p | vs the momentum p for the fixed values of the mass parameter |ν| in the small and large mass limits. Here we cover the mass parameter range 0 < |ν| < 3/2 for both of the cases. In fig. (5(a)) and fig. (5(b)), for the small mass limiting range for large values of the momentum p the two branches of solution for |ν| = 0 and |ν| = 1/2, 3/2 coincide with each other and in the small values of the momentum p the two branches of solution can be separately visualized. Further in fig. (5(c)) and fig. (5(d)), we have shown the behaviour of |γ p | vs mass parameter |ν| for the fixed values of the the momentum p in the small and large mass limits. For fixed value of p at p = 1 and p = 2 in the large mass limiting situation, we initially get decaying behaviour and then after |ν| = 3/2 both of the results obtained for p = 1 and p = 2 saturates. On the other hand, for same values of p in the large mass limit, we get oscillating behaviour. Combined effects for large and small mass limits are plotted in fig. (5(e)), where we have explicitly shown that for ν 2 < 0 for p = 1 and p = 2 we get distinguishable behaviour of the mass eigen values. Once ν 2 → 0 magnitude of the eigevalues increase and finally when ν 2 > 0 both of the plots show apeariodic behaviour.

Construction of density matrix
In this subsection our prime objective is construct the density matrix using the Bunch-Davies vacuum state which is expressed in terms of newly obtained sets of annihilation and creation operators in the Bogoliubov transformed frame. Most importantly, we already know that the full Bunch-Davies vaccum state can be expressed as a product of the vacuum state for each oscillator in the Bogoliubov transformed frame. Here each oscillators are labelled by the quantum numbers p, l and m. After tracing over the right part of the Hilbert space we get the following expression for the density matrix for the left part of the Hilbert space as: where the Bunch-Davies vacuum state takes the form: The above two eqns lead to the density matrix for the left part of the Hilbert space as: |Γ p,n | 2r |n, r; p, l, m n, r; p, l, m|, (3.130) where γ p and Γ p,n are defined in the earlier section and have defined the source normalization factor f p as: (3.131) In Eqn. Here we note the following crucial points: 1. The density matrix is diagonal for a given set of the SO(1, 3) quantum numbers p, l, m. This leads to the total density matrix to take the form as: 2. To find out a suitable normalization of the total density matrix, we use the following two results: Consequently using these results we get: This implies that the normalization condition of this total density matrix is fixed by the expression, Trρ L = 1+f p . This is consistent with the ref. [15] where f p = 0. But to maintain always Trρ L = 1 the total density matrix can be redefined by changing the normalization constant as: |γp| 2k |k; p, l, m k; p, l, m| Complementary part |Γp,n| 2r |n, r; p, l, m n, r; p, l, m| This also indicates that in such a situation entanglement is absent among all states which carries non identical SO(1, 3) quantum numbers p, l, m. 8 For an example here one can use the following equivalent ansatz for density matrix: 4. Finally, the total density matrix can be written in terms of entanglement modular Hamiltonian of the axionic Bell pair as, ρ L = e −βH ENT , where at finite temperature T dS of de Sitter space the parameter β is defined as, β = 2π/T dS . This also implies that at finite temperature the modular Hamiltonian can be expressed as 9 : |γ p | 2k |k; p, l, m k; p, l, m| |Γ p,n | 2r |n, r; p, l, m n, r; p, l, m| . (3.142) If we assume that the dynamical Hamiltonian in de Sitter space is represented by entangled Hamiltonian then for a given principal quantum number p the Hamiltonian for axionic Bell pairs can be expressed as: Acting this Hamiltonian on the Bunch-Davies vacuum state we find: where the total energy spectrum of this sytem can be written as: where the energy spectrum corresponding to the complementary and particular part of the wave function are given by: One can also recast the total energy spectrum of this sytem as: where O p is defined as: (3.148) Here in absence of the source term B p,n = 1∀n and consequently ln O p = 0. Note that, 9 If we use the equivalent ansatz for density matrix in presence of axionic source, in that situation the modular Hamiltonian can be expressed as:     for conformally coupled axion (ν = 1/2) and for minimally coupled axion (ν = 3/2) the SO(1, 3) principal quantum number p dependent spectrum can be expressed as: (3.149) In this case, the total energy spectrum is given by: This implies that for conformally coupled axion (ν = 1/2) and for minimally coupled axion (ν = 3/2) the entangled Hamiltonian (H ENT ) and the Hamiltonian for axion (H p ) R×H 3 are equivalent in absence of the linear source term in the effective action.
On the other hand, if we consider any arbitrary mass parameter ν in that case the SO(1, 3) principal quantum number p dependent spectrum can be expressed as: cosh 2πp + cos 2πν ± √ cosh 2πp + cos 2πν + 2 2 , (3.151) In this case, the total energy spectrum is given by: This implies that for with arbitrary parameter ν the entangled Hamiltonian (H ENT ) and the Hamiltonian for axion (H p ) R×H 3 are significantly differ as well in absence of the linear source term in the effective action.
In fig. (6(a)) and fig. (6(b)) we have shown the behaviour of the total energy spectrum for axion with respect to the SO(1, 3) quantum number p for a fixed value of the mass parameter ν without source contribution and with source contribution respectively. From both the plots it is clearly observed that for minimally coupled axion (ν = 3/2) (green) the energy spectrum is linear. Also in this case with source the slope and intercept will change. Also it is important to mention that, for conformally coupled axion (ν = 1/2) the situation is exactly similar and the behaviour overlaps with result obtained for minimally coupled axion (ν = 3/2). Further if we increase the value of the mass parameter to ν = 5 (blue) then it shows small deviation for very small values of p. Next if we go the large mass range, where ν 2 < 0 the energy spectrum shows significant deviation from the linearity for the large values of p. For example, we have considered ν = 3i/2 (orange) and ν = 5i (red) for our analysis. From both the plots it is clear that for ν = 5i (red) deviation from leniearity is more faster than the bevaiour obtained for ν = 3i/2 (orange).
In fig. (6(c)) and fig. (6(d)) we have shown the behaviour of the total energy spectrum for axion with respect mass parameter ν 2 for a fixed value of the SO(1, 3) quantum number p without source contribution and with source contribution respectively. From both the plots it is clearly observed that in presence of the source behaviour of the energy spectrum significantly changes compared to the result obtained for without source.

Computation of entanglement entropy
In this subsection we derive the expression for entanglement entropy in de Sitter space. In general the entanglement entropy can be written in terms of density matrix using Von Neumann measure as: where the parameter ν is defined earlier. For the Case I and Case II the expression for the entanglement entropy in terms of the complementary and particular part of the obtained solution for a given SO(1, 3) principal quantum number p can be expressed as 10 : Then the quantifying formula for the entanglement entropy in de Sitter space in presence of axion can be expressed as a sum over all possible quantum states which carries SO(1, 3) principal quantum number p. Technically, this sum can be written in terms of an integral as: Consequently, the final expression for the entanglement entropy in 3 + 1 D de Sitter space is given by the following expression: where D 3 (p) characterize the density of quantum states corresponding to the radial functions on H 3 , D 3 (p) = p 2 2π 2 . The volume of the hyperboloid H 3 is denoted by the overall factor V H 3 . This volume is infinite and can be regularize with a large radial cutoff r = L c in the hyperboloid H 3 . Technically, here regularization in volume stands for embedding surface of entanglement from zero to finite conformal time. After performing the regularization the volume of the hyperboloid H 3 for r ≤ L c can be written as 11 where we use the fact that the volume of S 2 is given by, V S 2 = 4π. Here one can interpret each of the term obtained in the regularized volume as: • Here the entangling area A ENT is given by: (3.160) • Consequently the second term is interpreted as logarithm of the entangling area A ENT , which is precisely the following factor: Finally the volume of the hyperboloid H 3 for r ≤ L c can be written in terms of the entangling area A ENT as: For our computation we will not further follow this ansatz of the density matrix. 11 In general for D − 1-sphere, which is denoted as S D−1 volume can be expressed as V S D−1 = 2π D/2 Γ( D In this context the cutoff regulator L c is large and in this limit this is identified with the conformal time as, L c ∼ − ln η. We also define the regularized volume of the hyperboloid H 3 as, V REG Thus the volume of the hyperboloid H 3 for r ≤ L c can be recast in terms of conformal time for large L c limit as: where c 6 is defined as: which represents the long range quantum entanglement in presence of axion. Further substituting the expression for entanglement entropy S(p, ν) computed in presence of axionic Bell pair and integrating over all possible SO(1, 3) principal quantum number, lying within the window 0 < p < ∞, we get the following expression for the co-efficient c 6 : where the integrals I 1 and I 2 can be written in 3 + 1 dimensional space time as: ln |γ p | 2 , (3.167) Here it is important to mention that: • For both Case I and Case II the integral I 2 diverges in 3 + 1 D within 0 < p < ∞. To make it finite we need to regularize this by introducing a cut-off Λ C in x = 2πp as: • For Case I we have two possible solution for γ p . For γ p = e −πp the integral I 1 is always finite for the Case I. But for the other solution of γ p = e πp the integral I 1 is divergent. To make the consistency throughout the analysis we put a cut-off Λ C on the integral I 1 obtained from both of the solutions of γ p . Here we explicitly show that in the final expression for I 1 after taking limit Λ C → ∞ final result is finite for γ p = e −πp and diverges for γ p = e πp .
Here we introduce a change in variable by using x = 2πp and introducing the cut-off Λ C for γ p = e −x/2 we get: where Li n (x)∀n = 2, 3, 4 is the polylogarithm function of n th kind, which is defined as After taking Λ C → ∞ for γ p = e −x/2 we get,lim Λ C →∞ I 1 = 1 90 . Further taking the source free limit in which f p → 0 we get: which is perfectly consistent with the results obtained from ν = 1/2 and ν = 3/2 cases without source in ref. [15]. This also implies that in absence of the source with Λ C → ∞ for this specific branch of solution the final result of the interesting part of the entanglement entropy is not very sensitive to the cut-off of SO(1, 3) principal quantum number p.
On the other hand, for γ p = e x/2 in 3 + 1 D introducing the cut-off Λ C in the rescaled principal quantum number x we get: in which after taking the source free limit in which f p → 0 we get: which is perfectly consistent with the results obtained from ν = 1/2 and ν = 3/2 cases without source for this branch of solution for γ p . In this context, this implies that in absence of the source with Λ C → Large for this specific branch of solution the final result of the interesting part of the entanglement entropy is highly sensitive to the cut-off of SO(1, 3) principal quantum number p.
• For Case II we are dealing with the the situation where due to arbitrary ν dependence the expressions for the solutions of the γ p are very complicated, which are given by, cosh 2πp + cos 2πν ± √ cosh 2πp + cos 2πν + 2 . (3.174) Let us first analyze the representative integral I 1 using both of the solutions obtained for arbitrary ν. Following the previous logical argument here we also put a cut-off Λ C to perform the integral on the rescaled SO(1, 3) quantum number x = 2πp and after performing the integral we will check the behaviour of both of the results. First of all we start with the following integral with " ± " signature, as given by: ln (2G ± (x, ν)) , (3.175) where G ± (x, ν) is defined as: However, it is not possible to solve this integral analytically using the exact mathematical structure. Small mass limiting situations are studied in ν = 1/2 as well in ν = 3/2 case which is appearing in Case I. Here we consider large mass limiting situation which is specifically important to study the physical imprints from Case II. In this large mass limiting situation we devide the total window of the SO(1, 3) principal quantum number p into two sub regions, as given by 0 < p < |ν| and |ν| < p < Λ C and in these consecutive region of interests the two solutions for γ p can be approximated as:  Consequenty, for this large mass limiting situation the regularized specified integral I 1 for the first solution for |γ p | can be written as: (3.179) and for the second solution for |γ p | can be written as: (3.180) Here A(ν), B(ν, Λ C ) abd C(ν, Λ C ) are defined as: (3.183) Further to check the dependency of the cut-off and mass parameter ν in the final results obtained for both of the solution of |γ p | obtained for large mass limiting situation for the range 0 < x < 2π|ν| we take the limit |ν| >> 1, which gives the following results for the first solution for |γ p |: Consequently in absence of source in the large mass limit the interesting part of the entanglement entropy can be written as: For the second solution of |γ p |, we can write: (3.186) and in absence of source in the large mass limit the interesting part of the entanglement entropy can be written as: lim |ν|>>1,fp→0 (3.187) Both of the results obtained for large |ν| limiting situation in absence of the source are consistent with the results obtained ref. [15].
In fig. (7(a)) and fig. (7(b)), we have explicitly shown the behaviour of entanglement entropy in 3 + 1 D de Sitter space in absence (f p = 0) and in presence (f p = 10 −7 ) of axionic source with respect to the mass parameter ν 2 where we consider both small and large mass limiting situation. In both the cases we have normalized the entanglement entropy with the result obtained from ν = 1/2 which is the conformal limiting result. In fig. (7(a)), it is clearly observed that in absence of axionic source in the large mass limiting range where ν 2 < 0, the normalized entanglement entropy S intr /S ν=1/2 saturates to zero asymptotically for large values of ν 2 in the negative axis. In the large mass limiting range actually the interesting part of the entanglement entropy for axion EPR Bell pair can be expressed in terms of its mass for the solution γ p = e −π|m axion /H| as:   which implies less entanglement for large mass limiting situation. Further, if we increase the value of ν 2 and move towards the ν 2 > 0 region then the normalized entanglement entropy S intr /S ν=1/2 will increase to its maximum value S intr /S ν=1/2 max = 1 at ν = 1/2. Further, between 3/2 < ν < 5/2 normalized entanglement entropy S intr /S ν=1/2 shows one oscillation and reaches its maximum value S intr /S ν=1/2 max = 1 again. Similar situation occurs within the interval 5/2 < ν < 7/2 and the same oscillating feature will repeat itself. But all such oscillations are apeariodic in nature and for large positive values of the ν 2 the period of oscillation increases. Similarly in fig. (7(b)), it is clearly depicted that in presence of axionic source in the large mass limiting range where ν 2 < 0, the normalized entanglement entropy S intr /S ν=1/2 saturates to zero asymptotically for large values of ν 2 in the negative axis to a large value, say S intr /S ν=1/2 sat ∼ 0.99975, which is appearing in the representative figure. This is obviously higher than the result obtained from without source.
On the other hand, in presence of axionic source if we increase the value of ν 2 and move towards the ν 2 > 0 region then the normalized entanglement entropy S intr /S ν=1/2 will increase to its maximum value S intr /S ν=1/2 max = 1 at ν = 1/2, which is exactly similar to the result obtained without axionic source. However if we compare both the results then we see that in absence of the axionic source the rate of increase in the normalized value of the entanglement entropy is much higher compared to the result when axionic source is present. This feature clearly indicates that the normalized entanglement entropy for all values of the mass parameter ν 2 is significantly higher in presence of the axionic source compared to the result obtained for without source. Further in between 3/2 < ν < 5/2 normalized entanglement entropy S intr /S ν=1/2 in presence of the axionic source shows one oscillation and reaches its maximum value S intr /S ν=1/2 max = 1 again very fast compared to the result obtained for without source. In presence of source, similar situation occurs in the subsequent intervals and the oscillating feature will repeat itself with increasing apeariodic nature for large positive values of the mass parameter ν 2 . But in presence of the source position of the minima in consecutive oscillations significantly increase compared to the case of without source. Additionally, it is important to note that for ν = 1/2 (red dashed line), ν = 3/2 (green dashed line) and ν = 7/2 (purple dashed line) we get same peak hights (maximum) of the oscillations and we get maximum normalized entanglement entropy at S intr /S ν=1/2 max = 1.
On the other hand, for the solution γ p = e π|m axion /H| the entanglement entropy for axion Bell pair in the large mass limiting range is very large: which implies very large entanglement for large mass limiting situation in a non trivial way. However, in this paper further we have not analyzed this case in detail due to its complicated structure.

Computation of Rényi entropy
In this section we further use the expression for the density matrix to compute the Rényi entropy, which is defined by the following expression 12 : with q > 0. (3.190) From this definition of Rényi entropy one can point out the following characteristics: • Here by taking the limit q → 1 one can produce the entanglement entropy for a given SO(1, 3) principal quantum number p and mass parameter ν i.e. lim q→1 S q (p, ν) = S(p, ν). Here S q (p, ν) characterize the Rényi entropy corresponding to single SO(1, 3) principal quantum number p and mass parameter ν. However, here it is important to note that this statement is only true in absence of any source. In presence of a source one can generally write the following expression: lim q→1 S q (p, ν) = S(p, ν) + ∆S(p, ν, f p ), (3.191) where ∆S(p, ν, f p ) is the additional contribution which is appearing due to axion and vanishes in the limit, lim fp→0 ∆S(p, ν, f p ) = 0.
• Further by taking the limit q → 0 in the definition of Rényi entropy one can produce the Hartley entropy for a given SO(1, 3) principal quantum number p and mass parameter ν i.e. lim q→0 S q (p, ν) = S H (p, ν). This is actually a measure of the dimension of the density matrix and can be written as, S H (p, ν) = ln (dim [H occupied ]) , where H occupied represents the image of the density matrix.
• In the limit q → ∞ Rényi entropy produces the largest eigenvalue of density matrix for a given SO(1, 3) principal quantum number p and mass parameter ν given by: 192) where [ρ L ] max represents the largest eigenvalue of density matrix. Also ∆λ p,ν,fp is the additional contribution arising from the source term and vanishes in the limit, lim fp→0 ∆λ p,ν,fp = 0.
• In absence of the axionic source, the quantifying formula for the Rényi entropy for a specified SO(1, 3) principal quantum number p and mass parameter ν can be expressed as: Further the specific information about the long range quantum correlation in 3+1 dimension is obtained by integrating over the SO(1, 3) quantum number p and a regularized volume integral over hyperboloid, as given by: dp p 2 S q (p, ν), (3.194) which satisfy the constraint, lim q→1 S q,intr = S intr ≡ c 6 . But since our prime objective of this paper is to see the imprints of stringy axion, we have to compute the modified expression for the Rényi entropy for a given SO(1, 3) principal quantum number p and mass parameter ν. In the next subsection we explicitly compute this result, in presence of the source.
The general expression for the Rényi entropy is written in Eqn (3.190), where parameter ν is defined earlier. For the Case I and Case II the expression for the Rényi entropy in terms of the complementary and particular part of the obtained solution for a given SO(1, 3) principal quantum number p can be expressed as: . (3.195) Now to study the properties of the derived result we check the following physical limiting situations as given by: • First of all we take the q → 1 limit for which we get the following simplified expression: (3.196) Further we take the difference between the results obtained from Eqn (3.156) and Eqn (3.196), which we define as: (3.197) which shows that in presence of source, the entanglement entropy and Rényi entropy in the limit q → 1 is not same. Now if we further take the limit f p → 0 then we get, lim fp→0 ∆S(p, ν, f p ) = 0, which is consistent with the result known for free theory.
• Next we take the q → 0 limit for which we get, lim q→0 S q (p, ν) = ln [ ∞ n=1 1] → ∞. This directly implies that the dimension of the image of the density matrix is given by for the present computation as, dim [H occupied ] → ∞, which is perfectly consistent with the fact that the density matrix in the present computation is infinite dimensional.
• Finally, we take the q → ∞ limit for which we get the following simplified expression: lim q→∞ S q (p, ν) ≈ − ln 1 − |γ p | 2 + ln (1 + f p ) . (3.198) This directly implies that the largest eigenvalue of density matrix given by for the present computation as, [ρ L ] max = 1 − |γ p | 2 1 1+fp , which is perfectly consistent with the fact that the density matrix in the present computation is infinite dimensional. Now if we further take the limit f p → 0 then we get, lim fp→0 ∆λ p,ν,fp = 0, which is consistent with the result known for free theory in absence of axion source contribution.
Then the final expression for the interesting part of the Rényi entropy in 3 + 1 D de Sitter space can be expressed as: dp p 2 S q (p, ν). (3.199) which represents the long range entanglement of the quantum state under consideration for the limit q → 1. Further substituting the explicit expression for entanglement entropy S(p, ν) computed in presence of axion and integrating over all possible SO(1, 3) principal quantum number, lying within the window 0 < p < ∞, we get the following expression: where the integrals J 1,q , J 2,q and J 3,q are defined as: . (3.203) Here it is important to note that: • For both Case I and Case II the integral J 2,q diverges in 3 + 1 D. To make it finite we need to regularize this by introducing a change in variable by using x = 2πp and introducing a cut-off Λ C in the rescaled principal quantum number x as given by: (3.204) • For Case I we have two possible solution for γ p . For γ p = e −πp the integrals J 1,q and J 3,q are always finite for the Case I. But for the other solution of γ p = e πp the integrals J 1,q and J 3,q are divergent. To have the consistency throughout the analysis we put a cut-off Λ C on the integrals J 1,q and J 3,q obtained from both the solutions of γ p . Here after introducing the cut-off Λ C in x = 2πp for γ p = e −x/2 we get: (3.205) . (3.206) Taking Λ C → ∞ for γ p = e −x/2 we get: and further taking q → 1, we get, lim Λ C →∞,q→1 J 1,q = 1 90 . Exact analytical form of the integral J 3,q for any q is not computable. But in the limiting case f p → 0 we get, lim Λ C →∞,q→1,fp→0 J 3,q = 0. This implies that in the source free limit where f p → 0 we get: which is consistent with the results obtained from ν = 1/2 and ν = 3/2 cases without source in 3 + 1 D.
On the other hand, for γ p = e x/2 in 3 + 1 D introducing the cut-off Λ C in the rescaled principal quantum number x we get: After taking q → 1 we get: Exact analytical form of the integral J 3,q for any q is not computable. But in the limiting case f p → 0 we get, lim Λ C →Large,q→1,fp→0 J 3,q = 0. This implies that in the source free limit in which f p → 0 we get:, which is perfectly consistent with the results obtained from ν = 1/2 and ν = 3/2 cases without source for this branch of solution for γ p in 3 + 1 D.
• For Case II following the similar procedure we get: where G ± (x, ν) is defined earlier.
Here we consider large mass limiting situation which is specifically important to study the physical imprints from Case II as mentioned earlier.
In this large mass limiting situation we devide the total window of the SO(1, 3) principal quantum number p into two sub regions, as given by 0 < p < |ν| and |ν| < p < Λ C : Consequenty, for this large mass limiting situation the regularized specified integral J 1,q and J 3,q for the first solution for |γ p | in 3 + 1 D can be written as: (3.215) (3.216) and for the second solution for |γ p | can be written as: (3.217) for 0 < x < 2π|ν| for 2π|ν| < x < Λ C . (3.218) Here A(ν) is defined earlier and D(ν, Λ C ) and E(ν, Λ C ) are defined as: Further to check the dependency of the cut-off and mass parameter ν in the final results obtained for both of the solution of |γ p | obtained for large mass limiting situation for the range 0 < x < 2π|ν| we take the limit Λ C → ∞ and q → 1, which gives the following results for the first solution for |γ p | in 3 + 1 D: Further taking |ν| >> 1 limiting case we finally get: On the other hand if we take the source limit f p → 0 then for the integral J 3,q we get, lim Λ C →∞,q→1,|ν|>>1,fp→0 J 3,q = 0. Consequently in absence of source in the large mass limit with q → 1 the interesting part of the Rényi entropy can be written as: For the second solution for |γ p | taking |ν| >> 1 limiting case we get: and lim Λ C →∞,q→1,|ν|>>1,fp→0 J 3,q = 0. Consequently, in absence of source in the large mass limit the interesting part of the Rényi entropy can be written as: Both of the results obtained for large ν limiting situation in absence of the source are consistent with the results obtained in this paper.
In fig. (8(a)) and fig. (8(b)), we have explicitly shown the behaviour of Rényi entropy in 3 + 1 D de Sitter space in absence (f p = 0) and presence (f p = 10 −7 ) of axionic source with respect to the mass parameter ν 2 where we consider both small and large mass limiting situation with cutoff of rescaled principal quantum number fixed at the value Λ C = 10000. In both the cases we have normalized the entanglement entropy with the result obtained from ν = 1/2 which is the conformal limiting result. In fig. (8(a)), it is clearly observed that in absence of axionic source in the large mass limiting range where ν 2 < 0, the normalized Rényi entropy S q,intr /S q,ν=1/2 saturates to zero   Here we set Λ C = 10000. In both the cases we have normalized the entanglement entropy with the result obtained from ν = 1/2 which is the conformal limiting result.   saturates to zero very slowly for the very large mass limiting range with ν 2 < 0. On the other hand, if we increase the value of ν 2 and move towards the ν 2 > 0 region then the normalized Rényi entropy S q,intr /S q,ν=1/2 will increase to its maximum value S q,intr /S q,ν=1/2 max = 1 at ν = 1/2 for all values of q, where q > 0. Further in between 3/2 < ν < 5/2 normalized Rényi entropy S q,intr /S q,ν=1/2 shows one oscillation and reaches its maximum value S intr /S ν=1/2 max = 1 again exactly at ν = 3/2 for all values of q, where q > 0. Here it is important to note that, as the value of q increase the position (height) of the minima of the oscillation decrease and exactly reaches the result with normalized entanglement entropy once we go to the q → 1 limiting situation. Here for q = 0.9 (red) we have represented this feature explicitly. Similar situation occurs within the interval 5/2 < ν < 7/2 and the same oscillating feature will repeat itself for all values of q, where q > 0. But all such oscillations are apeariodic in nature and for large positive values of the mass parameter ν 2 the period of oscillation increases for all values of q, where q > 0. Also in fig. (8(a)), for ν = 1/2 (black dashed line), ν = 3/2 (orange dashed line) and ν = 7/2 (magenta dashed line) we get same maxima of the oscillations at S q,intr /S q,ν=1/2 max = 1. In fig. (8(b)), it is clearly depicted that in presence of axionic source in the large mass limiting range where ν 2 < 0, the normalized Rényi entropy S intr /S ν=1/2 saturates to unity for q = 0.9 (red) or other value for q = 0.5 (green) and q = 0.7 (blue) through a jump at ν = 0 for large values of ν 2 in the negative axis. Except from the jump at ν = 0 in normalized Rényi entropy for q = 0.9 (red) it replicates the limiting situation q → 1 as required to produce normalized entanglement entropy. Such little deviation appears due to the presence of the axionic source term which is explicitly derived earlier.
Most importantly for a given value of q with q > 0 the height of the maxima of oscillation gives the maximum value of the normalized Rényi entropy and its magnitude increases with decreasing value of q which follows the restriction q > 0. On the other hand, in presence of axionic source if we increase the value of ν 2 and move towards the ν 2 > 0 region then the normalized Rényi entropy S q,intr /S q,ν=1/2 will decrease to its minimum value S q,intr /S q,ν=1/2 min = 1 at ν = 1/2, which is the maximum value of normalized Rényi entropy without axionic source. Further in between 3/2 < ν < 5/2 normalized Rényi entropy S q,intr /S q,ν=1/2 in presence of the axionic source shows one oscillation and reaches its minimum value S q,intr /S q,ν=1/2 max = 1 at ν = 3/2. In presence of source contribution minimum value of normalized Rényi entropy is appearing in the subsequent intervals and the oscillating feature will repeat itself with increasing apeariodic nature for large positive values of the mass parameter ν 2 . But in presence of the source position of the minima in consecutive oscillations significantly increase compared to the result obtained without source. Additionally, it is important to note that in fig. (8(b)), for ν = 1/2 (black dashed line), ν = 3/2 (orange dashed line) and ν = 7/2 (magenta dashed line) we get same minimum of the oscillations at S q,intr /S q,ν=1/2 min = 1.
Further in fig. (9(a)) and fig. (9(b)), we have explicitly shown the behaviour of Rényi entropy in 3 + 1 D de Sitter space in absence (f p = 0) and presence (f p = 10 −7 ) of axionic source with respect to the mass parameter ν 2 where we consider both small and large mass limiting situation with different cutoff of rescaled principal quantum number fixed at the small value Λ C = 300. In both the cases we have normalized the entanglement entropy with the result obtained from ν = 1/2 which is the conformal limiting result. In fig. (9(a)), it is clearly observed that in absence of axionic source in the large mass limiting range where ν 2 < 0, the normalized Rényi entropy S q,intr /S q,ν=1/2 saturates to different values depending on the different values of the parameter q as given by, q = 0.1 (black), q = 0.3 (purple), q = 0.5 (green), q = 0.7 (blue), q = 0.9 (red) with large values of ν 2 in the negative axis. On the other hand, if we increase the value of ν 2 and move towards the ν 2 > 0 region then we get the exactly similar features as appearing for fig. (8(a)). In fig. (9(b)), it is clearly depicted that in presence of axionic source in the large mass limiting range where ν 2 < 0, the normalized Rényi entropy S intr /S ν=1/2 saturates at different values for different values of the parameter q for large values of ν 2 in the negative axis. Most importantly for a given value of q with q > 0 the height of the maxima of oscillation increases with increasing value of q with Λ C = 300. On the other hand, in presence of axionic source if we increase the value of ν 2 and move towards the ν 2 > 0 region then the normalized Rényi entropy S q,intr /S q,ν=1/2 will increase to its maximum value S q,intr /S q,ν=1/2 max = 1 at ν = 1/2, which is also the maximum value of normalized Rényi entropy without axionic source with Λ C = 300. Further in between 3/2 < ν < 5/2 normalized Rényi entropy S q,intr /S q,ν=1/2 in presence of the axionic source shows one oscillation and reaches its maximum value S q,intr /S q,ν=1/2 max = 1 at ν = 3/2. Additionally, it is important to note that in fig. (9(b)), for ν = 1/2 (black dashed line), ν = 3/2 (orange dashed line) and ν = 7/2 (magenta dashed line) we get same maximum of the oscillations at S q,intr /S q,ν=1/2 max = 1.
Next, in fig. (10(a)) and fig. (10(b)), we have explicitly shown the behaviour of normalized Rényi entropy in 3 + 1 D de Sitter space in absence (f p = 0) and presence (f p = 10 −7 ) of axionic source with respect to the mass parameter ν 2 where we consider both small and large mass limiting situation alongwith the limiting case q → ∞, which produces the largest eigenvalue of the density matrix. Here it is important to note that, the feature obatined for the q → ∞ case is looks similar to the result obtained for normalized entanglement entropy in de Sitter space in 3 + 1 D. However, the differences are appearing in the height of the minimum of the oscillations obtained from ν 2 > 0 region. It also establishes that the results obtained for q → ∞ and q → 1 are not same. Here also in fig. (10(a)) and fig. (10(b)), for ν = 1/2 (black dashed line), ν = 3/2 (orange dashed line) and ν = 7/2 (magenta dashed line) we get same maxima of the oscillations at S q,intr /S q,ν=1/2 max = 1.

Conclusion
To summarize, in the present article, we have addressed the following points: • Firstly we have started our discussion with the computational strategy for quantifying the expression for entanglement entropy in de Sitter space for a bipartite quantum field theory of axion with a linear contribution in the effective potential as appearing from string theory using a specific compactification technique. Here we consider two subclass of models for axion, in Case I only linear field contribution is appearing and for Case II additionally a quadratic field combination is appearing in the effective potential for axion in φ < f a regime.
• Further, we have discussed about the basic setup for computing the all possible UV divergent and UV finite contribution of entanglement entropy in 3+1 dimension. We have also discussed the importance of UV finite contribution compared to the UV divergent part to quantify the interesting part of the entanglement entropy.
• Next, we have discussed the geometrical construction and underlying symmetries of the system under consideration in 3 + 1 D de Sitter space. Also we have discussed the string theory origin of axion effective potential and its role in the present context.
• Further we have derived the expression for the wave function of axion which is the solution of the field equation in a hyperbolic slice, commonly known as open chart. We have explicitly shown that the total solution can be written as a sum of a complementary solution and a particular solution, which is a completely new solution in presence of axion source term in the equation of motion.
• Next using the Bunch-Davies initial condition for the choice of the vacuum state we have written the total solution of the wave function in terms of creation and annihilation operators. Using this fact we construct a suitable basis by applying Bogoliubov transformation for the Bunch-Davies vacuum state using which we construct the expression for density matrix by tracing over the exterior region of the prescribed axionic quantum field theory in de Sitter space.
• Further, using the expression for the density matrix and also using the Von Neumann measure of entropy we have derived the new formula for entanglement entropy in 3 + 1 dimensional de Sitter space in presence of axion source and we have also checked the consistency of our result by comparing with the result obtained in ref. [15], which was derived for a free massive scalar field without linear contribution in the effective potential. Here it is important to note that the large mass limiting range we have provided the exact analytical expression for the entanglement entropy. But to analyze the correctness of our derived result and to compare with ref. [15] we have further used numerical techniques to study the behaviour of entanglement entropy with mass parameter ν 2 for Case I and Case II.
• Similarly we have computed the modified expression for the Rényi entropy in presence of axion source in 3 + 1 dimensional de Sitter space using which after taking q → 1 we have found that in presence of axion linear contribution it is not possible to get the exact analytical formula for entanglement entropy using Von Neumann entropy measure. Once we switch off the contribution from the linear term of the axion in the effective potential one can get back the exact formula for entanglement entropy as appearing in ref. [15] by taking q → 1 limit. Our analysis also clearly implies that the definition of Rényi entropy is not universal for any arbitrary structures of effective potential. Only for free theory where no such linear contribution appears, the definition is universal. Here the large mass limiting range we have provided the exact analytical expression for the Rényi entropy. But to analyze the correctness of our result and to compare with ref. [15] we have further used numerical techniques to study the behaviour of entanglement entropy with mass parameter ν 2 for Case I and Case II. Additionally, by using numerical technique and taking q → ∞ we have studied the behaviour of largest eigenvalue of the density matrix with respect to the mass parameter ν 2 for both the cases in presence of axion source.
• In this context it is important to mention here that, the appearance of non vanishing entanglement entropy in de Sitter space in 3 + 1 D directly verifies the correctness of the existence of the one point function in primordial cosmology due to axionic pair. Moreover, it is one of the important facts of cosmological evolution of our universe and temperature fluctuations of the Cosmic Microwave Background (CMB) originated from quantum fluctuations during the initial inflationary era as it acts a theoretical tool using which one can able to break the degeneracy amongst the cosmological predictions of various models.
The future prospects of our work are appended below: • In our analysis we choose Bunch-Davies vacuum state to specify the solution of wave function. Also it is used to determine rest of the quantities related to quantum entanglement. In future one can generalize this study of quantum entanglement in de Sitter space in presence of axion using α vacua [60].
• Using this prescribed methodology one can further derive the expression for two point and any other higher n point correlation function in superhorizon, subhorizon and horizon crossing scale to study the cosmological consequences of quantum fluctuation using Bunch-Davies vacuum as well as α vacua.
• The condition for violation of Bells inequality in the context of primordial cosmology are fixed by the further study of the quantum state in detail. In this connection computation of quantum discord, entanglement negativity and other possible strong measures in presence of axion will is also play important role. We have not addressed all these issues in this paper, which one can address in future for completeness.
Consequently the expression for the entanglement entropy in terms of the complementary and particular part of the obtaned solution can be expressed as:

B Derivation of Rényi entropy
For the Case I and Case II the expressions for the Rényi entropy can be expressed as a sum over four contributions, which can be expressed as: where the four explicit contributions in terms of the functions V m,q (p, ν)∀m = 1, 2, 3 as appearing in the entanglement entropy is given by: |γ p | 2k |k; p, l, m k; p, l, m| |Γ p,n | 2r |n, r; p, l, m n, r; p, l, m| Consequently the expression for the entanglement entropy in terms of the complementary and particular part of the obtaned solution can be expressed as: