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 an axion originating from Type IIB string compactification on a Calabi–Yau three fold (CY3\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathbf{CY^3}$$\end{document}) and in the presence of an NS5\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathbf{NS5}$$\end{document} brane. For this computation, we consider a spherical surface S2\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathbf{S}^2$$\end{document}, which divides the spatial slice of de Sitter (dS4\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathbf{dS_4}$$\end{document}) into exterior and interior sub-regions. We also consider the initial choice of vacuum to be Bunch–Davies state. First we derive the solution of the wave function of the 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 the entanglement entropy and Rényi entropy in 3+1\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$3+1$$\end{document} dimension. Furthermore, we quantify the UV-finite contribution of the entanglement entropy which contain the physics of long range quantum correlations of our expanding universe. Finally, our analysis complements the necessary condition for generating non-vanishing entanglement entropy in primordial cosmology due to the axion.


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][2][3]. It plays a significant role in the field of quantum information [4][5][6][7][8] and cosmology [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 does it change 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 the presence of a 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 the presence of an NS5 brane. So that the model consist of effective potential which breaks the shift symmetry associated with the axion field in a controlled manner. This model has been analysed 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 the interesting phenomenon 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 the 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 the future to provide an algorithm to contrast various models of inflation . Additionally, we compute Rényi entropy in 3 + 1 D de Sitter space and its connection with quantum entanglement in presence of the 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 Sect. 2 we briefly review the strategy for computing the entanglement entropy in de Sitter space. In Sect. 3 we follow the above strategy for the afore mentioned stringy axionic model. In Sect. 4 we present our conclusion and future prospects. Finally, in the appendix we give the detailed derivation of the entanglement entropy and Rényi entropy in the present context.

Computational strategy: brief review
Here we briefly review the method to derive the entanglement entropy in 3 + 1 dimensional de Sitter space following Ref. [15] and the 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 approximate bipartite decomposition [75] given by 1 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 localised 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 ρ = Tr R |BD BD|. (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 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 bound-ary 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] 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 S = S UV-divergent + S UV-finite . (2.3) 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 [76][77][78][79][80][81][82][83][84][85], 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]: 2 In arbitrary D dimensions entanglement entropy in de Sitter space is invariant under SO(1, D) isometry group.
× ln ( UV H ) , (2.4) 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 example, in Eq. (2.4), the first term represent area dependent contribution to the UV divergent part of the entropy, which originates from entanglement between particles in the vicinity of the entangling surface under consideration. Furthermore, 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 3 Here we in the flat space limit H → 0 we get lim H →0 ln ( UV H ) = 2 ln ( UV ) + , (2.5) where characterises the 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 limit 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 Eq. (2.4), explicitly incorporates the effect of bulk extrinsic and intrinsic curvature contribution in de Sitter background. In our computation of the 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. The 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 UV-finite 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 limit 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 an 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.
Furthermore, taking the flat space limit, H → 0, the UV-finite contributions to the entanglement entropy can be expressed as: 5 lim H →0 Furthermore, 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 + finite contributions, (2.13) where the coefficient of the logarithmic contribution, c 6 , characterises 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 [86]: (2.14) where g 1 and g 2 are two undetermined coefficients and K AB is the extrinsic curvature of entangling surface. Here we note that the coefficient 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. 5 Here we use the flat space limit H → 0 we get where characterises the contribution in the flat space limit, which is given by which is in general divergent but can be regularised by putting the previously introduced lattice UV cut-off UV as: where H UV is the harmonic number of order UV .
In the next section,we compute the entanglement entropy for axionic Bell pairs using the results of this section.

Geometrical construction and underlying symmetries
Let us first discuss the underlying symmetries of the present setup in de Sitter space. For this purpose, 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 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 . Furthermore, 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 regularised 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 characterised by the equation, 5 j=1 Y 2 j = H −2 , where for Euclidean signature we can write down the following parametric equations: Heren j ∀ j = 1, 2, 3 are the jth 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 . Furthermore, 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 regions, which are characterised as Furthermore, in hyperbolic slices the metric with Lorentzian signature can be expressed for the three regions as (3.10) (3.11) In Fig. 2 is a schematic diagram for the geometrical construction and underlying symmetries of the bipartite system.

Wave function of the axion in an open chart
In this section our prime objective is to compute the expression for the 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 the presence of an 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]87]. 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 where μ 3 is the mass scale, f a is the axion decay constant and we have defined a new parameter b as, b = 4 G /μ 3 f a . 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 the 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) Footnote 6 continued true [88,89]. 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 easily be connected to the violation in Bell's inequality [2]. 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.
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 . 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 constant in general follows a time dependent profile and for our purpose we use the 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 appearing in Case II, which we further use for the computation of the entanglement entropy to prove the existence of EPR Bell pairs in cosmology.
Furthermore, using Eq. (3.12) one can write down the following equation of motion for the axion for the abovementioned two cases:

17)
For Case II: 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 [90]: 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 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 labelled 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 (3.25) 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. Furthermore, one can write down the following equation of motion for the complementary part of the time dependent contribution in the total solution: 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.29) Case II : where in Case II we introduce a parameter ν, which is defined for axion as (3.31) From the above solutions obtained for Case I and Case II we observe the following characteristic features: • In the above solutions σ can take two values i.e. σ = ±1.
For each values of σ one can write down the final results in two regions -R and L regions of the hyperboloid (H 3 ). • 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 limit result of the entanglement entropy as the metric for de Sitter space is conformally flat.
• In the case where m axion < √ 2H parameter ν lies 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. • Furthermore, 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. • For R and L hyperboloid regions, the solutions have singularities at, ν ≡ ν 1 = − 1 2 − i p and ν ≡ ν 2 = − 1 2 + i p. This implies that to get non-singular solution the magnitude of the axion mass must satisfy the constraint : where N pσ is the overall normalisation constant, which is given by for Case II. (3.32) Furthermore, one can write down the following equation of motion for the particular integral part of the time dependent contribution in the total solution: where D t is the differential operator as given in Eq. (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 might be time dependent or 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 Here, χ (c) p n ,σ (t) represents the solution for the complementary part with a slight modification due to the replacement of p by p n ∀ n = 0, . . . , ∞ and is represented by Case II : Furthermore, using the results obtained for the total solution of the EOM given in Eq. (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 In this context, the Bunch-Davies mode function U σ plm (r, t, θ, φ) can be expressed as After substituting Eq. (3.40) in Eq. (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: (3.42) Particular solution: 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 normalisation constant N p for the complementary part is defined as, N p = 2 sinh π p N pσ , where N pσ is defined in Eq. (3.32). For the particular part if we replace p by p n then the normalisation constant N p n can be defined as, N p n = 2 sinh π p n N p n σ . In Eq. (3.44), expansion coefficient functions for the complementary solution (α σ q , β σ q ) for q = R, L are defined as Expansion coefficients for particular solution (ᾱ σ q ,β σ q ) can be obtained by replacing p by p n in Eqs. (3.46) and (3.47). Before going to the further details, let us first analyse Eq. (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].
Furthermore, using Eq. (3.44) we write the following matrix equation in a compact form: where for the complementary part of the solution we define the following matrices: , Similarly for the particular solution we also define the following matrices: where σ = ±1, q = R, L and I, J = 1, 2, 3, 4.
On the other hand the redefined normalisation constant for the particular part of the solution N p,(n) can be expressed as, N p,(n) = 2 sinh π p n N p n σ p 2 − p 2 n . Furthermore, using Eq. (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) 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: (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 (3.57) Similarlyγ jσ,n andδ * jσ,n can be obtained by replacing p by p n in Eqs. (3.56) and (3.57). Furthermore, we introduce a set of rules which are useful to differentiate the operations of these operators in the complementary and particular solution. These rules are: • Rule I: The operator corresponding to the complementary solution is insensitive to the particular solution i.e. a (c) 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 Furthermore, 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 bŷ where the coefficients m i j andm i j,n will be determined later. In Eq. (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. Furthermore, in principle one can express the particular part as Also the annihilation operator satisfy as well as the following commutation relations: Furthermore, 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: Thus, the complementary part and particular part vanish independently as the solutions are independent of each other. Consequently, we get the following constraints: Furthermore, using Eqs. (3.74) and (3.75), the mass matrices corresponding to the complementary part and particular part can be expressed as (3.76) 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 i j,n can similarly written by replacing p by p n in Eq. (3.77).
Before further discussion here we point out few important features: • For Case I we observe that for the complementary and particular part of the solution But for Case II we find that But for Case II we find that √ 2 e − p n π sinh p n π √ cosh 2π p n + cos 2πν .
Here also, Case II coincides with Case I for ν = 3/2. Additionally, the non-vanishing off-diagonal components for both 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 (3.83) Their explicit expressions are Case I : 84) λ ±,n = ±m RL,n = ± e i(θ+ π 2 ) e − p n π , (3.85) Case II : 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 Eqs. (3.60) and (3.61). In Fig. 4a, b, we have shown the behaviour of the eigenvalue (λ ± ) vs. the momentum p for fixed values . 4 Behaviour of the eigenvalues (λ ± ) in de Sitter space for '+' branch of solution. This figure clearly shows that we cover both large and small axion mass limit situations of the mass parameter |ν| in the small and large axion mass limits. Here we cover the mass parameter range 0 < |ν| < 3/2 for both cases. In Fig. 4a for the small mass limit 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 visualised. Exactly opposite situation appears when we consider the large mass limit in Fig. 4b. Furthermore, in Fig. 4c, d, we have explicitly shown the behaviour of the eigenvalue (λ ± ) vs. mass parameter |ν| for the fixed values of 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. 4e, where we have explicitly shown that for ν 2 < 0 (for p = 1 and p = 2) we get a distinguishable behaviour of the mass eigen values. Once ν 2 → 0 magnitude of the eigenvalues increase 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 (3.88) with the following constraints: (3.91) Using the above set of operators one can write the Bunch-Davies vacuum state in terms of a new Bogoliubov transformed basis represented by the R and L vacuum state as where we have introduced a new operatorQ, defined in the basis aŝ 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 by (3.94) The overall normalisation 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 (3.96) Before going to further details let us first mention the following useful commutation relations of the creation and annihilation operators of the R and L vacua: In this context, the operations of creation and annihilation operators defined on the Bunch-Davies vacuum state are
(3.100) Furthermore, 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 (3.103) Furthermore, using Eqs. (3.88) and (3.89), in Eqs. (3.99) and (3.100), we get the following sets of homogeneous equations for the two cases:   Similarly for Case II for the diagonal and off-diagonal terms we get m RR,n =m LL,n =m * RR,n = √ 2 e − p n π cos πν √ cosh 2π p n + cos 2πν , (3.118) m RL,n =m LR,n = −m * RL,n = e i π 2 √ 2 e − p n π sinh p n π √ cosh 2π p n + cos 2πν .
(3.119) 3. In 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 normalisation |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 For Case−I : For Case−II :  In Fig. 5a, 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. 5a, b, for the small mass limit 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 visualised. Furthermore, in Fig. 5c, d, we have shown the behaviour of |γ p | vs. mass parameter |ν| for the fixed values of 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 limit 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. 5e, 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 eigenvalues 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 vacuum 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: (ρ L ) p,l,m = Tr R |BD BD|, (3.128) where the Bunch-Davies vacuum state takes the form: The above two equations lead to the density matrix for the left part of the Hilbert space as |γ p | 2k |k; p, l, m k; p, l, m| | 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 normalisation factor f p as (3.131) In Eq. (3.130), the states |k; p, l, m and |n, r ; p, l, m are defined in terms of the quantum state in left Hilbert space as 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 (3.133) 2. To find a suitable normalisation of the total density matrix, we use the following two results: Consequently using these results we get This implies that the normalisation condition of this total density matrix is fixed by the expression, Trρ L = 1+ f p . This is consistent with Ref. [15] where f p = 0. But to maintain always Trρ L = 1 the total density matrix can be redefined by changing the normalisation constant as: |γ p | 2k |k; p, l, m k; p, l, m| Complementary part 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. 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 8 For example here one can use the following equivalent ansatz for density matrix: |γ p | 2k |k; p, l, m k; p, l, m| Complementary part (3.139) 9 If we use the equivalent ansatz for density matrix in presence of the axionic source, in that situation the modular Hamiltonian can be expressed as |γ p | 2k |k; p, l, m k; p, l, m| | p,n | 2r |n, r ; p, l, m n, r ; p, l, m| |γ p | 2k |k; p, l, m k; p, l, m| where the total energy spectrum of this system can be written as where the energy spectrum corresponding to the complementary and particular part of the wave function are given by (3.146) One can also recast the total energy spectrum of this system as where O p is defined as (3.148) Here in the absence of the source term B p,n = 1∀n and consequently ln O p = 0. Note that for a conformally coupled axion (ν = 1/2) and for a minimally coupled axion (ν = 3/2) the SO(1, 3) principal quantum number p dependent spectrum can be expressed as In this case, the total energy spectrum is given by (3.150) This implies that for a conformally coupled axion (ν = 1/2) and for a minimally coupled axion (ν = 3/2) the entangled Hamiltonian (H ENT ) and the Hamiltonian for axion (H p ) R×H 3 are equivalent in the 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  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 the absence of the linear source term in the effective action. In Fig. 6a, 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 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 a conformally coupled axion (ν = 1/2) the situation is exactly similar and the behaviour overlaps with the result obtained for minimally coupled axion (ν = 3/2). Furthermore, if we increase the value of the mass parameter to ν = 5 (blue) then it shows a 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 plots it is clear that for ν = 5i (red) deviation from linearity is more faster than the behaviour obtained for ν = 3i/2 (orange). In Fig. 6c, 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 plots it is clearly observed that in the presence of the source the behaviour of the energy spectrum significantly changes compared to the result obtained without source.

Computation of the entanglement entropy
In this subsection we derive the expression for the entanglement entropy in de Sitter space. In general the entanglement entropy can be written in terms of density matrix using Von Neumann measure as Then the quantifying formula for the entanglement entropy in de Sitter space in the presence of the axion can be expressed as a sum over all possible quantum states which carries Footnote 10 continued For our computation we will not further follow this ansatz of the density matrix.

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) characterises 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 regularise with a large radial cut-off r = L c in the hyperboloid H 3 . Technically, here regularisation in volume stands for embedding surface of entanglement from zero to finite conformal time. After performing the regularisation 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 regularised volume as: • Here the entangling area A ENT is given by • Consequently the second term is interpreted as logarithm of the entangling area A ENT , which is precisely the following factor: (3.161) Finally, the volume of the hyperboloid H 3 for r ≤ L c can be written in terms of the entangling area A ENT as (3.162) In this context the cut-off regulator L c is large and in this limit this is identified with the conformal time as, L c ∼ − ln η. We also define the regularised 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 the large L c limit as (3.163) 11 In general for D − 1-sphere, which is denoted as S D−1 volume can be expressed as Finally, summing over all possible principal quantum numbers p of SO(1, 3) we get from Eq. (3.158) the following expression for the entanglement entropy in de Sitter space: S(ν) = c 6 1 4η 2 + ln η , (3.164) where c 6 is defined as , ν), (3.165) which represents the long range quantum entanglement in the presence of axion. Furthermore, substituting the expression for the entanglement entropy S( p, ν) computed in the presence of the axionic Bell pair and integrating over all possible SO(1, 3) principal quantum numbers, lying within the window 0 < p < ∞, we get the following expression for the coefficient c 6 : where the integrals I 1 and I 2 can be written in 3 + 1 dimensional space time as

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 regularise this by introducing a cut-off C in x = 2π p as (3.169) • For Case I we have two possible solution for γ p . For γ p = e −π p the integral I 1 is always finite for 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 the 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 Li n (x) = ∞ k=1 x k k n . After taking C → ∞ for γ p = e −x/2 we get lim C →∞ I 1 = 1 90 . Furthermore, taking the source free limit in which 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 the 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 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 situation where due to arbitrary ν dependence the expressions for the solutions of the γ p are very complicated, which are given by (3.174) Let us first analyse 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 (3.176) However, it is not possible to solve this integral analytically using the exact mathematical structure. Small mass limit situations are studied in the ν = 1/2 and in the ν = 3/2 cases, is appearing in Case I. Here we consider large mass limit situation which is specifically important to study the physical imprints from Case II. In this large mass limit situation we divide the total window of the SO(1, 3) principal quantum number p into two subregions, as given by 0 < p < |ν| and |ν| < p < C and in these consecutive regions of interest the two solutions for γ p can be approximated by |γ p | = e −π |ν| for 0 < p < |ν| e −π p for |ν| < p < C /2π, (3.177) and |γ p | = e π |ν| for 0 < p < |ν| e π p for |ν| < p < C /2π. (3.178) Consequently, for this large mass limit situation the regularised specified integral I 1 for the first solution for |γ p | can be written as (3.179) and for the second solution for |γ p | it can be written as (3.180) Here A(ν), B(ν, C ) and C(ν, C ) are defined as (3.183) Furthermore, to check the dependency of the cut-off and mass parameter ν in the final results obtained for both of the solutions of |γ p | obtained for the large mass limit 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 the absence of the 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 lim |ν| 1 (3.186) and in the absence of the source in the large mass limit the interesting part of the entanglement entropy can be written as Both of the results obtained for the large |ν| limit situation in absence of the source are consistent with the results obtained Ref. [15].
In Fig. 7a, b, we have explicitly shown the behaviour of entanglement entropy in 3 + 1 D de Sitter space in the absence ( f p = 0) and in the presence ( f p = 10 −7 ) of the axionic source with respect to the mass parameter ν 2 where we consider both small and large mass limit situation. In both cases we have normalised the entanglement entropy with the result obtained from ν = 1/2, which is the conformal limit result. In Fig. 7a, it is clearly observed that in the absence of the axionic source in the large mass limit range where ν 2 < 0, the normalised entanglement entropy S intr /S ν=1/2 saturates to zero asymptotically for large values of ν 2 on the negative axis. In the large mass limit 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 the large mass limit situation. Furthermore, if we increase the value of ν 2 and move towards the ν 2 > 0 region then the normalised entanglement entropy S intr /S ν=1/2 will increase to its maximum value S intr /S ν=1/2 max = 1 at ν = 1/2. Furthermore,   Here we set C = 10,000.
In both cases we have normalised the entanglement entropy with the result obtained from ν = 1/2, which is the conformal limit result between 3/2 < ν < 5/2 normalised 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. 7b, it is clearly depicted that in the presence of the axionic source in the large mass limit range where ν 2 < 0, the normalised entanglement entropy S intr /S ν=1/2 saturates to zero asymptotically for large values of ν 2 on 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 the axionic source if we increase the value of ν 2 and move towards the ν 2 > 0 region then the normalised 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 the two results then we see that in the absence of the axionic source the rate of increase in the normalised value of the entanglement entropy is much higher compared to the result when axionic source is present. This feature clearly indicates that the normalised entanglement entropy for all values of the mass parameter ν 2 is significantly higher in the presence of the axionic source compared to the result obtained for without source. Furthermore, in between 3/2 < ν < 5/2 normalised entanglement entropy S intr /S ν=1/2 in the 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 the 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 the 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 (dashed line), ν = 3/2 (green dashed line) and ν = 7/2 (purple dashed line) we get the same peak heights (maximum) of the oscillations and we get maximum normalised 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 limit range is very large: which implies very large entanglement for the large mass limit situation in a non-trivial way. However, in this paper further we have not analysed this case in detail due to its complicated structure.

Computation of the 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 (3.190) From this definition of the 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, ν) characterise 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 the absence of any source. In presence of a source one can generally write the following expression: (3.191) where S( p, ν, f p ) is the additional contribution which is appearing due to axion and vanishes in the limit, lim f p →0 S( p, ν, f p ) = 0. • Furthermore, by taking the limit q → 0 in the definition of the 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 (3.192) where [ρ L ] max represents the largest eigenvalue of density matrix. Also λ p,ν, f p is the additional contribution arising from the source term and vanishes in the limit, lim f p →0 λ p,ν, f p = 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 Furthermore, the specific information as regards the long range quantum correlation in 3+1 dimension is obtained by integrating over the SO(1, 3) quantum number p and a regularised volume integral over hyperboloid, as given by (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 the presence of the source.
The general expression for the Rényi entropy is written in Eq. (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 limit situations: • First of all we take the q → 1 limit for which we get the following simplified expression: (3.196) Furthermore, we take the difference between the results obtained from Eqs. (3.156) and (3.196), which we define as 197) which shows that in the presence of the 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 f p →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, 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 f p →0 λ p,ν, f p = 0, which is consistent with the result known for free theory in the absence of the 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: , ν), (3.199) which represents the long range entanglement of the quantum state under consideration for the limit q → 1.
Furthermore, substituting the explicit expression for the entanglement entropy S( p, ν) computed in the presence of the axion and integrating over all possible SO(1, 3) principal quantum number, lying within the window 0 < p < ∞, we get the following expression: S q,intr = J 1,q + ln 1 + f p J 2,q + J 3,q , (3.200) where the integrals J 1,q , J 2,q and J 3,q are defined as (3.202) (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 regularise 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 Case I. But for the other solution of γ p = e π p the integrals J 1,q and J 3,q are divergent. To have consistency throughout the analysis we put a cut-off C on the integrals J 1,q and J 3,q obtained from the two solutions of γ p . Here after introducing the cut-off C in x = 2π p for γ p = e −x/2 we get C q 4 ln −e − C + 15 3 C q 3 ln −e − C q + 45 2 C q 4 Li 2 e C − 45 2 C q 2 Li 2 e C q − 90 C q 4 Li 3 e C +90q 4 Li 4 e C +90 C qLi 3 e C q − 90Li 4 e C q − π 4 q 4 + π 4 , (3.205) (3.206) 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 limit case f p → 0 we get lim C →∞,q→1, f p →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 (3.209) The exact analytical form of the integral J 3,q for any q is not computable. But in the limit case f p → 0 we get lim C →Large,q→1, f p →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 a similar procedure we get (3.213) where G ± (x, ν) is defined earlier. Here we consider large mass limit situation which is specifically important to study the physical imprints from Case II as mentioned earlier. In this large mass limit situation we divide the total window of the SO(1, 3) principal quantum number p into two sub-regions, as given by 0 < p < |ν| and |ν| < p < C : Consequently, for this large mass limit situation the regularised 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) (3.218) Here A(ν) is defined earlier and D(ν, C ) and E(ν, C ) are defined as   4 Li 4 e C + 2 C qLi 3 e C q − 2Li 4 e C q .
(3.220) Furthermore, to check the dependency of the cut-off and mass parameter ν in the final results obtained for both of the solutions of |γ p | obtained for the large mass limit 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: Furthermore, taking the |ν| 1 limit case we, finally, get lim C →∞,q→1,|ν| 1 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, f p →0  For the second solution for |γ p | taking the |ν| 1 limit case we get lim C →Large,q→1,|ν| 1 Both of the results obtained for large ν limit situation in absence of the source are consistent with the results obtained in this paper.
In Fig. 8a, b, we have explicitly shown the behaviour of Rényi entropy in 3 + 1 D de Sitter space in the absence ( f p = 0) and presence ( f p = 10 −7 ) of the axionic source with respect to the mass parameter ν 2 where we consider both small and large mass limit situation with cut-off of rescaled principal quantum number fixed at the value C = 10000. In both cases we have normalised the entanglement entropy with the result obtained from ν = 1/2, which is the conformal limit result. In Fig. 8a, it is clearly observed that in the absence of axionic source in the large mass limit range where ν 2 < 0, the normalised Rényi entropy S q,intr /S q,ν=1/2 saturates to zero suddenly through a jump at ν = 0 with q = 0.5 (green), q = 0.7 (blue), q = 0.9 (red) for large values of ν 2 on the negative axis. But for q = 0.1 (black) and q = 0.7 (purple) S q,intr /S q,ν=1/2 saturates to zero very slowly for the very large mass limit range with ν 2 < 0. On the other hand, if we increase the value of ν 2 and move towards the ν 2 > 0 region then the normalised 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. Furthermore, in between 3/2 < ν < 5/2 normalised 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  (red) with '+' branch of solution of |γ p | and | p,n |. Here we set C = 300. In both cases we have normalised the entanglement entropy with the result obtained from ν = 1/2, which is the conformal limit result that, as the value of q increases, the position (height) of the minima of the oscillation decrease and exactly reaches the result with normalised entanglement entropy once we go to the q → 1 limit 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. 8a, for ν = 1/2 (black dashed line), ν = 3/2 (orange dashed line) and ν = 7/2 (magenta dashed line) we get the same maxima of the oscillations at S q,intr /S q,ν=1/2 max = 1. In Fig. 8b, it is clearly depicted that in the presence of axionic source in the large mass limit range where ν 2 < 0, the normalised 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 on the negative axis. Except from the jump at ν = 0 in normalised Rényi entropy for q = 0.9 (red) it replicates the limit situation q → 1 as required to produce normalised 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 normalised Rényi entropy and its magnitude increases with decreasing value of q, which follows the restriction q > 0. On the other hand, in the presence of the axionic source if we increase the value of ν 2 and move towards the ν 2 > 0 region then the normalised 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 normalised Rényi entropy without axionic source. Furthermore, in between 3/2 < ν < 5/2 the normalised Rényi entropy S q,intr /S q,ν=1/2 in the 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 the presence of the source contribution the minimum value of the normalised Rényi entropy is appearing in the subsequent intervals and the oscillating feature will repeat itself with increasingly aperiodic nature for large positive values of the mass parameter ν 2 . But in the presence of the source the position of the minima in consecutive oscillations significantly increase compared to the result that is obtained without source. Additionally, it is important to note that in Fig. 8b, for ν = 1/2 (black dashed line), ν = 3/2 (orange dashed line) and ν = 7/2 (magenta dashed line) we get the same minimum of the oscillations at S q,intr /S q,ν=1/2 min = 1.
Furthermore, in Fig. 9a, b, we have explicitly shown the behaviour of the Rényi entropy in 3 + 1 D de Sitter space in absence ( f p = 0) and presence ( f p = 10 −7 ) of the axionic source with respect to the mass parameter ν 2 where we consider both small and large mass limit situation with different cut-off of rescaled principal quantum number fixed at the small value C = 300. In both cases we have normalised the   Here we set C = 10,000. In both cases we have normalised the entanglement entropy with the result obtained from ν = 1/2, which is the conformal limit result entanglement entropy with the result obtained from ν = 1/2, which is the conformal limit result. In Fig. 9a, it is clearly observed that, in the absence of the axionic source in the large mass limit range where ν 2 < 0, the normalised 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 on 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 to those appearing for Fig. 8a. In Fig. 9b, it is clearly depicted that in the presence of the axionic source in the large mass limit range where ν 2 < 0, the normalised Rényi entropy S intr /S ν=1/2 saturates at different values for different values of the parameter q for large values of ν 2 on 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 the axionic source if we increase the value of ν 2 and move towards the ν 2 > 0 region then the normalised 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 normalised Rényi entropy without axionic source with C = 300. Furthermore, in between 3/2 < ν < 5/2 normalised Rényi entropy S q,intr /S q,ν=1/2 in the 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. 9b, for ν = 1/2 (black dashed line), ν = 3/2 (orange dashed line) and ν = 7/2 (magenta dashed line) we get the same maximum of the oscillations at S q,intr /S q,ν=1/2 max = 1. Next, in Fig. 10a, b, we have explicitly shown the behaviour of the normalised Rényi entropy in 3+1 D de Sitter space in the absence ( f p = 0) and the presence ( f p = 10 −7 ) of the axionic source with respect to the mass parameter ν 2 where we consider both small and large mass limit situation along with the limit case q → ∞, which produces the largest eigenvalue of the density matrix. Here it is important to note that the feature obtained for the q → ∞ case looks similar to the result obtained for normalised 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 the ν 2 > 0 region. It also establishes that the results obtained for q → ∞ and q → 1 are not same. Here also in Fig. 10a, b, for ν = 1/2 (black dashed line), ν = 3/2 (orange dashed line) and ν = 7/2 (magenta dashed line) we get the same maxima of the oscillations at S q,intr /S q,ν=1/2 max = 1.

Conclusion
To summarise, 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 the entanglement entropy in de Sitter space for a bipartite quantum field theory of the 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. • Furthermore, we have discussed about the basic setup for computing all possible UV-divergent and UV-finite contribution of the 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 the axion effective potential and its role in the present context. • Furthermore, we have derived the expression for the wave function of the 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 the presence of the 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 a Bogoliubov transformation for the Bunch-Davies vacuum state using which we construct the expression for the density matrix by tracing over the exterior region of the prescribed axionic quantum field theory in de Sitter space. • Furthermore, using the expression for the density matrix and also using the Von Neumann measure of entropy we have derived the new formula for the entanglement entropy in 3 + 1 dimensional de Sitter space in the presence of the 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 limit range we have provided the exact analytical expression for the entanglement entropy. But to analyse the correctness of our derived result and to compare with Ref. [15] we have further used numerical techniques to study the behaviour of the 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 the presence of the axion source in 3+1 dimensional de Sitter space using which after taking q → 1 we have found that in the presence of the axion linear contribution it is not possible to get the exact analytical formula for the entanglement entropy using Von Neumann entropy measure. Once we switch off the contribution from the linear term of the axion in the effective potential we can get back the exact formula for the entanglement entropy as appearing in Ref. [15] by taking the q → 1 limit. Our analysis also clearly implies that the definition of the Rényi entropy is not universal for any arbitrary structures of the effective potential. Only for a free theory where no such linear contribution appears, the definition is universal. Here the large mass limit range we have provided the exact analytical expression for the Rényi entropy. But to analyse the correctness of our result and to compare with Ref. [15] we have further used numerical techniques to study the behaviour of the entanglement entropy with mass parameter ν 2 for Case I and Case II. Additionally, by using numerical techniques and taking q → ∞ we have studied the behaviour of largest eigenvalue of the density matrix with respect to the mass parameter ν 2 for both cases in the presence of the 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) originating from quantum fluctuations during the initial inflationary era as it acts as a theoretical tool using which one might be able to lift 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 generalise this study of quantum entanglement in de Sitter space in the presence of the axion using α vacua [91]. • 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 the Bunch-Davies vacuum and the α vacua. • The conditions 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 the presence of the axion will also play an important role. We have not addressed all these issues in this paper.

B Derivation of the Rényi entropy
For 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| Consequently the expression for the entanglement entropy in terms of the complementary and particular part of the obtained solution can be expressed as (B.4)