Quantum nature of Wigner function for inflationary tensor perturbations

We study the Wigner function for the inflationary tensor perturbation defined in the real phase space. We compute explicitly the Wigner function including the contributions from the cubic self-interaction Hamiltonian of tensor perturbations. Then we argue that it is no longer an appropriate description for the probability distribution in the sense that quantum nature allows negativity around vanishing phase variables. This comes from the non-Gaussian wavefunction in the mixed state as a result of the non-linear interaction between super- and sub-horizon modes. We also show that this is related to the explicit infrared divergence in the Wigner function, in contrast to the trace of the density matrix.


Introduction
Quantum mechanics has been established as a framework for describing nature down to the microscopic scale, yet there still remain many conceptual questions to be answered. In particular, a large gap between the classical and quantum description has brought about the attempts to explain the quantum-to-classical transition dynamically. Remarkably, cosmology is expected to provide a natural testing ground for this issue (for a recent review, see [1]). That comes from the fact that, whereas the primordial fluctuations as we observe from the cosmic microwave background are classical, it is likely to be originated from the quantum fluctuations [2][3][4][5][6] as supported by the inflationary cosmology [7][8][9].
In order to find out the observable that reveals the quantum origin of the primordial fluctuations, we need to make clear the meaning of the quantum-to-classical transition in the cosmological context. Since the typical size of the interaction between gravitons or curvature perturbation is given by H/m Pl 1 or even more suppressed by the factor of the slow-roll parameter, the dominant loss of the quantum nature can be discussed in terms of the quadratic action only. More concretely, the quadratic action encodes the incessant interaction between the quantum fluctuations and gravitational background as well. Then a set of quantum fluctuations behaves coherently, forming the two-mode squeezed state about which the quantum effects coming from the non-commutativity of field operators are suppressed [10][11][12][13][14]. Here the two-mode means that the quantum fluctuation with threemomentum k and that with −k appear pairwise. But still, such a squeezing is a unitary evolution which transforms the pure state into the pure state. Hence we need an additional JHEP03(2020)060 mechanism that converts the pure state into the mixed state, in which the (diagonalized) density matrix elements correspond to the quantum analog of classical probability distribution. This is achieved by decoherence, in which the system and the environment interact with each other through higher-order non-linear interaction (for reviews, see [15,16]). In the inflationary universe, the horizon of size 1/H becomes a natural cutoff distinguishing the system from the environment: an effective field theory below the energy scale H for the super-horizon modes is an open system that interacts with the sub-horizon modes as an environment [17][18][19][20][21][22][23][24]. Then the evolution of the super-horizon modes is not unitary, as characterized by the Lindblad terms [25]. They describe the transition from the pure to the mixed state in decoherence through the time evolution of the density matrix [26]. One of natural choices of the basis for the density matrix is a set of states obtained from the unitary time evolution of different particle number states. Since the unitary time evolution includes squeezing, each of basis states is the coherent superposition of multiparticle states, but still forms an orthogonal set of basis. Then we have the classical probability that the universe has evolved from different initial states, in which some of quantum fluctuations have excited [24].
With the mechanism described above, one may naively expect that cosmological perturbations as we observe are completely described in terms of classical physics. However, unlike the name stands, such a "quantum-to-classical transition" still preserves the quantum nature, from which we may confirm the quantum origin of cosmological perturbations. For example, whereas the density matrix is a quantum analog of the classical probability distribution, they are represented in terms of the quantum states we observe. 1 On the other hand, in statistical mechanics, the classical probability distribution is defined in the phase space: it is a function of canonical variables which are eigenvalues of non-commuting quantum operators. The Wigner (distribution) function [27] is devised to construct the classical probability distribution of such type from the density matrix (see [28][29][30] for early discussions in the context of cosmology). However, only limited system allows the Wigner function to be interpreted as the probability distribution since its positivity is not guaranteed [31]. The Wigner function is positive definite only when the states comprising the density matrix is represented by the Gaussian wavefunction. In the case of the harmonic oscillator, for example, since the Gaussian wavefunction appears only in the ground state, the Wigner function is positive definite only if we have the pure ground state. Regarding the squeezed state, the squeezing of the vacuum state is represented by the Gaussian wavefunction but when the squeezed states of the initial quantum excitations are mixed, we can find the region in the phase space in which the Wigner function becomes negative. This shows that the probability distribution in the phase space is not always well-defined, reflecting the quantum nature of the primordial fluctuations even after decoherence. Therefore, the negativity of the Wigner function is a signal of quantum origin of the system. We also note that the positivity of the Wigner function does not guarantee the absence of the 1 Since the quantum interference effects disappear through decoherence, one may attempt to explain the collapse of wavefunction after measurement in terms of it. However, the probabilistic nature still remains, so we need to employ another interpretation scheme in addition, such as the many-world interpretation. For more discussion, see, e.g. [16] and references therein.

JHEP03(2020)060
quantum effect 2 [32,33]. For instance, we may find out a Bell inequality violation even if the Wigner function is non-negative. Checking the positivity of the Wigner function is just a basic test for the quantum nature of the system.
In this article, we study the appearance of such negativity of the Wigner function in the context of cosmological perturbations. Our conclusion is that the higher-order nonlinear interaction gives rise to the possibility for our universe to have the non-Gaussian wavefunction or equivalently, the squeezed states of the initial quantum excitation through the super-and sub-horizon interaction, hence the absolute positivity of the Wigner function is not guaranteed. Moreover, the negative Wigner function is related to its infrared divergence. This is in contrast to the unit value of the trace of the density matrix, which is the result of the cancellation between the infrared divergence in the no-graviton-excitation state and that in the soft graviton states. Such an explicit infrared divergence implies that the Wigner function is no longer a meaningful description for the probability distribution of the cosmological perturbations. For this purpose, we consider the tensor perturbations, or the graviton fluctuations, for more quantitative analysis. Whereas the cosmological perturbations also contain the curvature perturbation, reflecting a spontaneous breaking of the de Sitter (dS) isometry in quasi-dS spacetime [34][35][36][37], the mechanism for the negativity of the Wigner function is essentially the same as that of the tensor perturbations: the only qualitative difference is a factor of the slow-roll parameter which is an order parameter for the breaking of the dS isometry. In section 2, we briefly review the squeezing of tensor perturbations, which is useful in our discussions. From this, in section 3, we present the Wigner function explicitly to show the breaking of its positivity as a result of nonlinear interaction which was intensively studied in [24]. Implication of the Wigner function is visited in section 4 after which we conclude. We also provide appendix sections to give some details not presented in the main text.

Quantum state for tensor perturbations
We consider the spatial metric as where dτ = dt/a is the conformal time and h ij is the pure tensor perturbations with Since there are two physical degrees of freedom for h ij , we introduce the polarization tensor e ij (λ), with λ being the polarization index, so that h ij = 2 λ=1 h λ e ij (λ). For canonical normalization, we introduce then there are two copies of the identical action of a canonically normalized scalar field v λ for each polarization state so from now on we just consider only one polarization state and drop the subscript λ. At quadratic level, adopting the Heisenberg picture, the conjugate JHEP03(2020)060 pair for tensor perturbations can be written aŝ where the time-dependence is given, not to the mode functions, but to the creation and annihilation operators. Here, the initial mode functions u k and v k are found to be What is important here is that we cannot isolate the evolution of the pure sub-system of the mode with k: when we excite the k mode, −k mode is at the same time excited as well. This becomes clear if, from the above expression, we express the creation and annihilation operators in terms of the canonical variables as where it is very important to note that a † k (τ ) is given by the canonical pair of not k but −k. This is obviously because the particle creation always occurs in pair with the opposite momenta. To isolate the mode with k from that with −k, we "define" the positionq k and momentump k which are purely written in terms of the operators with the corresponding momentum k, not involving the opposite momentum −k as Now we make two steps further. First, we "discretize" the momentum as d 3 k/(2π) 3 → L −3 k where L 3 is the volume under consideration. With a given volume L 3 , we may now isolate the volume from the canonical creation and annihilation operators a k and a † k in such a way that, since a k has a mass dimension of −3/2, we define and likewise for a † k . That means, the new dimensionless operatorsâ k andâ † k now satisfy the following commutation relation: Second, we note thatq k andp k are of mass dimension −2 and −1 respectively. To apply directly the wisdom of the standard quantum harmonic oscillators, we rescale the dimensionful factors to introduce the dimensionless positionq k and momentump k aŝ

JHEP03(2020)060
This new pair is obviously Hermitian, i.e.q † k =q k andp † k =p k , and satisfies the canonical commutation relation: q k ,p q = i â k ,â † q = iδ kq . (2.14) We note that typical cosmological observables are correlators of field operators with respect to the initial vacuum, from which we can obtain information on (2.12) and (2.13) indirectly. However, the Hermiticity of these operators implies that there must in principle be a clever way to "measure" them. Given the canonical variables (2.12) and (2.13), we can now make use of the position basis vector |q k , q −k for the two-mode state with k and −k. From noting that H = k H k with we can write the (free) evolution operator as a product of U k which arises solely from H k : Then the state |Ψ that has evolved from the initial vacuum |0 at τ 0 , using the complete basis |q k , q −k , can be written where we have defined the wavefunction Ψ(q k , q −k ) in the position basis |q k , q −k as the representation of |Ψ . The evolution operator U k (τ, τ 0 ) can be written as the product of the "rotation" operator R k and the "squeezing" operator S k , i.e. U k = S k R k . Further, while the only role of R k is to change the phase, S k is as the name stands solely responsible for the two-mode squeezing of the initial vacuum state. More concretely, in the rotation operator the exponentsâ † kâ k andâ † −kâ −k are the number operators counting the number of k and −k mode excitations, respectively. Then when R k is applied to |0 , it just provides an overall phase e −iθ k . In the density matrix we will obtain in section 3, the phase does not play any role in the matrix element for In the same way, the phase does not contribute to the diagonal density matrix elements like U kâ † k |0 0|â k U † k as well. Thus, from now on we drop the irrelevant rotation operator and consider only S k in the wavefunction. The phase will be written explicitly when the JHEP03(2020)060 off-diagonal density matrix element is discussed [see (3.20)]. The wavefunction we are interested in is given by where the n-particle excited state, with the momenta k and −k each, is Then it can be shown that Ψ(q k , q −k ) is given by Here, r k and ϕ k are the parameters of the Bogoliubov transformation given by (B.4) and (B.5). The detail of the derivation of (2.22) is given in [38] (see also [39]).

Wigner function for tensor perturbations
Now we consider explicitly the reduced density matrix ρ red of the inflationary tensor perturbations on super-horizon scales under the influence of cubic interaction with sub-horizon modes. We obtain the Wigner function from ρ red . If we ignore cubic or higher nonlinear interactions, super-and sub-horizon modes behave independently, then the two-mode squeezed state comprised of super-horizon modes form the pure state. Since the wavefunction for the squeezed vacuum state in (2.22) is Gaussian, the Wigner function for this state is also Gaussian and positive definite. In this case, it has been shown that the two-point correlators in terms of the quantum and the classical stochastic formalism are not distinguishable regardless of squeezing [30]. On the other hand, studies in this section show that when we take into account the non-linear interaction, the inevitable interaction between super-and sub-horizon modes results in the violation of the positivity of the Wigner function. This indicates that the Wigner function cannot be used as a conceivable probability distribution under decoherence, invalidating correlators defined upon the Wigner function as physical quantities. Under the interaction between the super-("system" S) and sub-horizon modes ("environment" E), tracing out the environmental sub-horizon modes introduces the non-unitary terms in the Lindblad equation:

JHEP03(2020)060
where dots indicate the irrelevant unitary evolution terms and the Lindblad operators L 1 and L 2 are defined in terms of the free, unitary evolution operator for the environment U 0,E with the subscript 0 standing for the free evolution and the interaction Hamiltonian H int by [22,24] Here the interaction Hamiltonian with subscripts S and I correspond to those written in the Schrödinger and interaction picture, respectively. Since the interaction Hamiltonian is dominated by cubic term for H/m Pl 1, L † 1 L 2 or L † 2 L 1 in (3.1) provide at most six creation operators for super-horizon modes depending on how many sub-horizon excitations are traced out in (3.2). Since these excitations act on the initial vacuum, they evolve in time by the free unitary operator U 0,S dominantly. From now on, we omit the subscript S as obviously we only consider the evolution of super-horizon modes. Hence, to the lowest order in H/m Pl , ρ red for each super-horizon mode k is written in the schematic form where the evolution operator, creation and annihilation operators and the vacuum state are those for the super-horizon modes. Thus we can see that ρ red contains inherently a set of basis and using this basis can be written in the matrix form, with the subscript ab denoting the element at a-th row and b-th column in the matrix, where ρ 20 = ρ * 02 . The result is written to the leading order, O(H 2 /m 2 Pl ). If we take the quartic or higher interactions into account, the components for the reduced density matrix will be filled more. For the detailed calculation, see [24]. Note that only for the 00 component the evolution operator acts upon the vacuum state, it is the only one for which we can apply directly the wavefunction (2.22) which, as can be seen from (2.20), has evolved from the vacuum state. Other components such as 11 are from the excited initial states, so require more care. 00 component. First we consider the simplest case -the 00 component:

JHEP03(2020)060
where the coefficient is given by [24] with the infrared cutoff in the Planck unit ε 1. Since we are dealing with a two-mode excited state, we need to expand the simplest Wigner function (A.1) to include the other, opposite momentum −k as [1] As we have noted, concentrating only on S k for the evolution operator with the subscript 1 and 2 denoting respectively k and −k, the contribution of ρ red | 00 to the Wigner function becomes (3.9) Note that as mentioned before we can write the Wigner function in terms of the wavefunction only for the 00 component. After some calculations, with the detail given in appendix C, we find The exponent for the exponential is a four-variable function (p k , p −k , q k , q −k ) and is thus hard to visualize. But nevertheless what is easy to note is that the whole exponential factor is always positive definite. Thus we conclude that the Wigner function for the 00 component (3.10) is always positive.
11 component. Now we move to the first non-trivial part, the 11 component of ρ red . Schematically, from (3.3) where the coefficient is given by [24] (with typo corrected) where q a = −q b ≡ q andq ≡ q/H <

JHEP03(2020)060
Thus among all k-modes contained in the Hamiltonian, only the mode which has the same momentum as the external one survives. This can be more directly understood by noting that in (3.11), the evolution operator U 0 is operational on the state a † q |0 ≡ |1 q , (3.14) the one-particle excited state with the momentum q that belongs to the system sector, q < H. As the evolution operator is a free one, only the mode with q is responding to |1 q and all the other modes simply disappear. From now on, without losing generality we identify the super-horizon external momentum q as k. Now let us proceed the practical calculations for the Wigner function for the 11 component. Rewriting (3.11) gives where we have dropped the rotation operator R k in the evolution operator U k = R k S k .
Here, for the annihilation operator it should be a −k due to the delta function in ρ 11 as shown in (3.12), but as the squeezing operator S k will excite both k and −k modes, we do not distinguish k and −k assigned to the operators. Since S k is operating on not |0 but a † k |0 , we cannot directly apply the wavefunction (2.20) as we did for the 00 component. Instead, we have to calculate how S k evolves the one-particle excited state a † k |0 . The detail of the calculations is given in appendix C, and we find This is a complicated four-variable function of (p k , p −k , q k , q −k ) so not easy to visualize. But nevertheless we can collect the relevant parts. For this purpose, we note that the exponential factor is the same as W 00 and is always positive definite. Thus, the only part we need to see if negative at all is the polynomial terms inside the curly brackets. They are, for any moderate value of the squeezing parameter r k , exponentially positive. Nevertheless, they become negative for very small values of (p k , p −k , q k , q −k ) all around zero. We will return to this issue in section 4.
where the coefficient is given by ρ 20 = ρ 11 e 2ikτ . Again, we have dropped the rotation operator R k in the evolution operator U k = R k S k . As before, since S k is operating on not |0 but

JHEP03(2020)060
we cannot directly apply the wavefunction (2.20) as we did for the 00 component. Instead, we have to calculate how S k evolves the two-particle excited state a † −k a † k |0 . The detail is given in appendix C, and we find W 20 as where the factor e −3iθ k at the end comes from the rotation operator (2.19) acting on |2 k . Note that with the coefficient of the 02 component of ρ red being given by ρ 02 = ρ * 20 , computing W 02 explicitly gives Total Wigner function. Now we consider the total Wigner function for the reduced density matrix. To begin with, we note that the Wigner function elements (3.10), (3.16) and (3.20) are all for a single mode with k. The "total" Wigner function of the system of our interest should include all the modes. To see how the contributions from other momentum modes are included, let us consider simply two modes, k 1 and k 2 : where |n i corresponds to the n-particle excited state with the momentum k i , e.g. |1 1 = |1 k 1 , and ρ ijkl is an arbitrary coefficient. Now, to compute the "total" Wigner function of this system, since there are two momenta, 3 we introduce two dummy integra-3 Precisely speaking, there are two excitations each, with positive and negative momenta. Further, as there are infinitely many discrete momenta, the whole state can be written as so that the states given in (3.22) [as well as in (3.14) and (3.19)] are in fact and so on.

JHEP03(2020)060
tion variables s 1 and s 2 such that where we have written in such a way that the Wigner function of each mode is more manifest. We should first note that when we are exciting a certain mode, the others are in their vacuum states. Thus the terms inside the square brackets, multiplied by the vacuumvacuum component of the Wigner function for the k 2 mode is the total Wigner function only for the k 1 mode, except that the component for which all the modes are in the vacuum states are equally distributed, in the current case divided by two as there are two modes. This is because, as the case for which all the modes are in the vacuum is unique, this should equally contribute to the Wigner function of each mode. This suggests that for a huge number of discrete modes we should divide (3.10) by the volume of the super-horizon momentum space with the ultraviolet cutoff given by H, V SH = 4πH 3 /3 × 1/2, with the factor 1/2 being due to the pairwise appearance of k and −k in the two-mode squeezed state. Then, (3.23) is written as Expanding ρ 20 = ρ 11 e 2iqτ , we may just sum W 11 , W 20 and W 02 to find (3.25) Thus, (3.24) can be written as

JHEP03(2020)060
This can be immediately extended to a large number of discrete momenta to give (3.27) Now using the de Sitter solutions for r k and ϕ k given respectively by (B.4) and (B.5), we find sinh r k = 1 2kτ , (3.28) then it is trivial to find the explicit time dependence. Especially, in the late-time limit kτ → 0, we find This shows that for non-zero value of (q k , q −k , p k , p −k ), W 11 + W 20 + W 02 becomes large and positive at late time, −kτ → 0, as the first term in the curly bracket is dominant. On the other hand, as all the value of (q k , q −k , p k , p −k ) approach zero, only the second term, −1 in the curly bracket remains, which results in the negative Wigner function.

Discussions
In the previous section, we have computed the Wigner function contributions from the nonlinearly evolved elements in the reduced density matrix of tensor perturbations. While the exponential factor, which the 00 component also contains, remains positive definite, the factor −1 in the polynomial terms in (3.25) can make the Wigner function negative for small enough values of (p k , p −k , q k , q −k ). One may then argue naively as follows. These variables are written essentially in terms of the mode function and its time derivative. In any reasonable model of inflation, typically the mode function becomes frozen on superhorizon scales and accordingly its time derivative almost vanishes. Further, the mode function amplitude contains the factor H/m Pl 1. Thus, very naively, it seems that we are naturally led to have vanishingly small (p k , p −k , q k , q −k ). Interestingly, Wigner function becomes negative around this region, which signals the quantum nature of the cosmological perturbation still remains even after the horizon excape.
Before getting into detail, we make one point clear. For this purpose, let us consider the Bogoliubov transformation (B.1), where the time-and momentum-dependent coefficients α k (τ ) and β k (τ ) are subject to the constraint (B.2). The well-known "mode function" JHEP03(2020)060 v k (τ ) is given by the coefficient of the initial creation and annihilation operators: from (2.4) and (B.1), we can writê Thus, while both α k (τ ) and β k (τ ) contain the plane-wave e ±ikτ , the mode function v k (τ ) as we know is given by a specific combination of them. Meanwhile, from (2.12) and (2.13) we can read that 1) the position and momentum operatorsq k andp q have no specific combination of α k (τ ) and β k (τ ) that gives the mode function solution v k (τ ) and thus they cannot be written in terms of v k (τ ), and 2)q k andp q contain both initial annihilation and creation operators for both k and −k modes, and thus are superpositions of ±k modes.
Then what about the state |q k , q −k ? As in the case of textbook quantum harmonic oscillator, the state upon which the annihilation and creation opeartors act is the particle number state: the state that denotes 0 particle with a specific momentum k (a k |0 = 0), one-particle (a † k |0 = |1 k ), and so on. Thus the position state |q k and the particle number state |n k -the collection of states that can be written in terms of the vacuum |0 and the creation operator a † k -are independent. What we can do is to write one state in terms of the other as, given the competeness n |n k n k | = 1, where Ψ n (q k ) is the "wavefunction", such as (A.12) for one-dimensional harmonic oscillator. Now, we can address what is the nature of q k and so on that appear in (3.25). Inserting (2.21) |n, k; n, −k n, k; n, −k|q k , q −k =Ψ * n (q k ,q −k ) = infinite sum of |n−1; n , |n; n+1 , |n+1; n , |n; n−1 . Thus we conclude we find a collection of infinitely many excitations, and unlike the naive expectation we have seen at the beginning of this section, no direct relation to the mode function v k (τ ) can be found.
Since the wavefunction of the squeezed vacuum (2.22) is Gaussian while those of the squeezed excitations are not [see (C.25) and (C.35)], the Wigner function is not positive definite. Even worse, the Wigner function is not free of the infrared divergence. These show that Wigner function as the classical probability distribution in the phase space is not well-defined contrary to the density matrix, which is positive definite and infrared JHEP03(2020)060 finite. To be more quantitative on the infrared finiteness, consider the Wigner function at q k = q −k = 0 and p k = p −k = 0, then we have Comparing this to the trace of the density matrix which is taken over the (q k , q −k ) space [we denote them by (x, y) for comparison], we find that the crucial difference is the additional negative sign in the ket state. For the density matrix to be well-defined description for the mixed state probability distribution, the trace of the density matrix is infrared finite by the cancellation between the infrared divergent parts of ρ 00 and ρ ii with i ≥ 1. This comes from the fact that the extremely soft excitations are not distinguishable to the vacuum state, which is originally used to argue the infrared finiteness of the quantum transition amplitude [40]. Then we can normalize the trace by Tr[ρ] = 1 reflecting the total probability is given by 1. Explicit calculation in [24] shows that at least to the order of H 2 /m 2 Pl , the infrared finiteness of the trace, where d 3 q/(2π) 3 ρ 11 (q) = ρ 00 , is well satisfied from the cancellation between infrared divergences in ρ 00 and ρ 11 . Moving to W (0, 0; 0, 0), since the first (second) excitation state is odd (even) under x → −x and y → −y, i.e. −x, −y|1 = − x, y|1 and −x, −y|2 = x, y|2 , we find easily This can be checked explicitly from (3.10), (3.16) and (3.20): putting q k = q −k = 0 and p k = p −k = 0, we obtain respectively W 00 (0, 0; 0, 0) = 4 1 − ρ 00 , W 11 (0, 0; 0, 0) = −4ρ 11 , and W 20 (0, 0; 0, 0) = 0, which is consistent with the above. The relative minus sign between two terms shows that the infrared divergence in ρ 11 is no longer cancelled with that in ρ 00 hence W (0, 0; 0, 0) cannot be infrared finite.
We can also see that if ρ 00 > 1/2, W (0, 0; 0, 0) becomes negative. From (3.10), this is possible either for large value of H 2 /m 2 Pl or small value of ε. For large H 2 /m 2 Pl , however, we should expect the next-to-leading corrections in terms of the interaction Hamiltonian of O(H 4 /m 4 Pl ) should be significant so that our perturbative results should be modified. Meanwhile, the term containing the infrared cutoff ε may call for different perturbative expansion in terms of (H 2 /m 2 Pl ) log ε and even resummed to some analytic function as JHEP03(2020)060 can be seen in QCD. No matter how ε behaves, the point is that the infrared cutoff ε is introduced to regulate the infrared divergence, which is eventually cancelled by the soft graviton contribution in ρ 11 to make Tr[ρ] = 1 while not in the Wigner function. This indicates that the standard interpretation of the Wigner function as the classical probability in the phase space is invalid because of not only being negative, but also exhibiting "naked" infrared divergence.

Conclusions
In this article, we have tested whether the Wigner function can be a conceivable physical quantity measuring the probability distribution on the real phase space of the primordial tensor perturbation. We first observe that even our universe has started from the Bunch-Davies vacuum, it would eventually become the mixed state containing the squeezed states of the initial excitations as a result of the non-linear interaction. The wavefunctions of squeezed states are not Gaussian unless the Bunch-Davies vacuum is squeezed, indicating that the positivity of the Wigner function is not guaranteed. We find that around the vanishing values of the real phase space variables (q k , q −k , p k , p −k ), the negativity of the Wigner function is evident. This in fact has to do with the cancellation of the infrared divergence to achieve Tr[ρ] = 1. Such a cancellation does not take place in the Wigner function, which not only leads the Wigner function to the negative value, but also makes it uncontrollably divergent. These two enable us to conclude that the Wigner function is not a good description to the probability distribution of the cosmological perturbations, indicating that the quantum nature is not completely lost even after decoherence.

JHEP03(2020)060
given by Ψ(x) = x|Ψ , Then, clearly W (q, p) is real. If integrated over either position or momentum, we find Given a Hermitian operator A, we define the so-called "Weyl transformation" A(q, p) as, similar to the Wigner function, Integrating the product of A and the Wigner function W over q and p gives We now consider the quantum harmonic oscillator with m = ω = 1. The Hamiltonian operator is written as where the pair (q,p) is canonical: We can define the annihilation and creation operators respectively bŷ which satisfy the canonical commutation relation: Then the normalized wavefunction in the position space Ψ(x) = x|Ψ is found to be where n is an integer and H n (x) is the Hermite polynomial.
With the wavefunction Ψ n (x), we can compute straightly the Wigner function for each state. In table 1 we show the results for n = 0, 1 and 2 explicitly. As we consider higher excited states, the corresponding Wigner functions become more complicated. But they all exhibit negative values in certain regime, not always positive definite. Table 1. The Wigner functions W n (q, p) of a quantum harmonic oscillator and the region where W n ≤ 0 for the states n = 0, 1 and 2.

B Bogoliubov transformation
The general solutions for the Hamiltonian equations for the creation and annihilation operators are given by the linear combination of the initial ones: which is the so-called Bogoliubov transformation. Then standard commutation relations lead to the constraint α k (τ ) and β k (τ ) are always subject to: Then we can parametrize them in terms of the hyperbolic functions as α k (τ ) = e −iθ k cosh r k , β k (τ ) = e i(θ k +2ϕ k ) sinh r k .

(B.3)
Assuming a perfect de Sitter background so that a /a = −1/τ , we can find analytically the solutions as [13,14] r k = sinh −1 1 2kτ , which are used to derive the well-known mode function solution v k (τ ).

JHEP03(2020)060
C Explicit calculations for the Wigner functions

C.1 00 component
We first consider the explicit calculations of W 00 . From (3.9), using the explicit form of Ψ given by (2.22), where for the second equality we have used the fact that q k is Hermitian and have defined the coefficient |C| 2 as which is positive definite. Further, since the exponential factor in the integral contains the mixing terms of the integration variables x and y as dxdy exp(· · · + xy + · · · ), the integration of the exponential factor is not in the simplest form. One typical way to separate the variables is to introduce new variables as then the integration measure transforms as Finally, we note that the exponent of (2.22) can be written as thus we can identify A = −(α + β)/2 and B = α − β so that α = 1 + e 2iϕ k tanh(r k /2) 2 1 − e 2iϕ k tanh(r k /2) = 1 + i sin(2ϕ k ) sinh r k 2 cosh r k − cos(2ϕ k ) sinh r k , (C.6) β = 1 − e 2iϕ k tanh(r k /2) 2 1 + e 2iϕ k tanh(r k /2) = 1 − i sin(2ϕ k ) sinh r k 2 cosh r k + cos(2ϕ k ) sinh r k . (C.7)

JHEP03(2020)060
Collecting all these, we can write W 00 analytically to find which can be integrated analytically to give (3.10).

C.2 11 component
We begin with the original expression for S k as Operating uponâ † k |0 first gives 10) and the commutator reads S k ,â † k = ABC,â † k = A,â † k BC + A B,â † k C + AB C,â † k . (C.11) Three simplifications are ahead. First, as we can see from (C.9), A only containsâ † −kâ † k , and thus the commutator withâ † k simply vanishes: A,â † k = 0 . (C.12) Next, if we expand the exponential in C, C = 1 + e −2iϕ k tanh r k 2 â −kâk + 1 2! e −2iϕ k tanh r k 2 2 â −kâk 2 + · · · , (C. 13) which only containsâ −kâk . Thus, in C,â † k , while one annihilation operator brings withâ † k a delta function via the canonical commutation relation, still more than one annilation operators remain which eliminate the vacuum state |0 multiplied to the right. Thus, we have C,â † k |0 = 0 . (C.14) Finally, if C is acting directly on the vacuum state, only the first term in (C.13) survives and the original vacuum state remains identical:

JHEP03(2020)060
Thus, we only have S k ,â † k |0 = A B,â † k |0 . (C. 16) To proceed further, we consider B more closely, we first single outâ † kâ k and expand in terms of logarithm to write Multiplying the vacuum |0 to the right and separately considering n = 0, Sinceâ † k andâ † −k commute,â † k may go in front of theâ † −kâ −k term. Expandingâ † −kâ −k similarly gives Thus, finally, (C.23) Since now we have exchanged the position ofâ † k and S k such that the squeezing operator S k is directly acting on the vacuum state from the left, we are in the half way of making use of the explicit form of the wavefunction Ψ(q k , q −k ) = q k , q −k S k (τ, τ 0 ) 0 0 given by (2.20). But still we have a † k lurking around, which we should work out. A hint on how to deal withâ † k comes from the explicit functional form of the wavefunction Ψ(q k , q −k ). It is, as explicitly shown in (2.20), a function of q k and q −k . This could well be anticipated as it is the representation of the squeezed state evolved from the vacuum S k |0 in the JHEP03(2020)060 Therefore, we finally find q k , q −k S kâ † −kâ † k 0 = 1 2 cosh 2 (r k /2) Now we have all the ingredients to write the Wigner function W 20 . We find explicitly By changing the variables as x = (s + t)/2 and y = (s − t)/2 with A = −(α + β)/2 and B = α − β, we can perform the integrations to find (3.20).
Open Access. This article is distributed under the terms of the Creative Commons Attribution License (CC-BY 4.0), which permits any use, distribution and reproduction in any medium, provided the original author(s) and source are credited.