{'fraction_non_alphanumeric': 0.03632679850469423, 'fraction_numerical': 0.011157720514796108, 'mean_word_length': 4.608834827144686, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 7, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 5, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Quantum Phase Space Description of a Cosmological Minimal Massive Bigravity Model

Bimetric gravity theories describes gravitational interactions in the presence of an extra spin-2 field. The Hassan-Rosen (HR) nonlinear massive minimal bigravity theory is a ghost-free bimetric theory formulated with respect a flat, dynamical reference metric. In this work the deformation quantization formalism is applied to a HR cosmological model in the minisuperspace. The quantization procedure is performed explicitly for quantum cosmology in the minisuperspace. The Friedmann-Lema\^itre-Robertson-Walker (FLRW) model with flat metrics is worked out and the computation of the Wigner functions for theHartle-Hawking, Vilenkin and Linde wavefunctions are done numerically and, in the Hartle-Hawking case, also analytically. From the stability analysis in the quantum minisuper phase space it is found an interpretation of the mass of graviton as an emergent cosmological constant and as a measure of the deviation of classical behavior of the Wigner functions. Also, from the Hartle-Hawking case, an interesting relation between the curvature and the mass of graviton in a cusp catastrophe surface is discussed


Introduction
One of the most successful theories we have in physics is the Einstein's General Relativity Theory (GR) [1]. We can observe that the majority of the phenomena related with gravity can be explained with the Einstein's Field Equations (EFE), from the classical Newtonian gravity force description to the light bend and gravitational waves observations: all of this can be simulated and reproduced in a proper limit with a very high precision. Nevertheless, in standard cosmology there is a important problem related with the present acceleration of Universe. The usual explanation is consider the addition of a cosmological constant Λ which must be adjusted to fit experimental data. An alternative is to look for modifications to GR in order to get weaker interaction at cosmological scales but that remains good to describe short distance behavior, and this is a important issue if we are interested in study the cosmology in the firsts instants of Big Bang. Many proposals for extend EFE has been developed over the years, where string theory is a well-know candidate for a quantum theory of gravity. GR is expressed as a non-linear self-interacting massless spin-2 particle field theory; from this point of view, a modification to gravity could be achieved through a massive spin-2 particle. Fierz and Pauli [2] linearized a massive spin-2 field fluctuation, and Boulware and Deser showed that any nonlinear extension to Fierz-Pauli theory would exhibit a ghost inestability [3,4], but one of the most significant advance was the construction of a theory of ghost-free nonlinear massive gravity by de Rham, Gabadadze and Tolley (dRGT) [5]. Later, Hassan and Rosen (HR) demonstrates that in a Hamiltonian constraint analysis, the Boulware-Deser ghost is abscent [6].
dRGT theory develop the dynamics of a massive spin-2 field in a flat space in which two metric tensor components are playing roles, one as dynamical g µν and the other as a non-dynamical reference metric f µν . As a result, open Friedmann-Lemaître-Robertson-Walker (FLRW) universes are allowed [7] but problems with strong coupling and ghostlike inestabilities arise [8]. Attemping to avoid the problems with flat FLRW solutions, HR [9] extended the massive gravity theory beyond dRGT to a theory with two dynamical symmetric tensors g µν and f µν , taked as foreground and background metrics, respectively, and having a symmetric and intertwined role. They called this a massive bigravity theory, and this shows a ghost-free description bimetric model containing nonlinear interactions of a massless and a massive spin-2 field in a dynamical background. In order to incorporate bimetric gravity models to cosmology, Rosen [

Canonical Quantization of a Minimal Bigravity Theory
The present section is an overview of the general procedure of a minimal massive gravity model. Only the principal results for the next work are presented, and detailed information can be found in Refs. [32,9,6,33]

Hamiltonian formalism
To set the Arnowitt-Deser-Misner (ADM) decomposition, is necessary consider a pseudo-Riemannian manifold (Ω, g) which is globally hyperbolic. Associated to this pseudo-Riemannian metric tensor g there is another pseudo-Riemannian metric tensor f , the reference metric. In general dRGT theory, f is necessary to ensure a ghost-free non-linear massive gravity theory. Thus, spacetime Ω admits the decomposition Ω = R × Σ, where Σ is a hypersurface on which the foliations ds 2 Here γ ij and ξ ij are the extrinsic metrics on Σ, N and M the lapse functions, and N i and M i the shift vectors related to g and f , respectively.
The space of all Riemannian 3-metrics and scalar matter configurations φ defined on Σ is called the x ∈ Σ} denotes the infinite-dimensional space of Riemannian metrics on Σ, then the moduli space M is defined as the superspace modulo the group of diffeomorphisms Diff(Σ), M = Riem(Σ)/ Diff(Σ) or for pure (bi)gravity, M = Met(Σ)/ Diff(Σ).
The phase space for pure (bi)gravity (Wheeler's phase superspace) Γ * ∼ = T * Met(Σ), the cotangent bundle of Met(Σ), is given by the pairs Γ * = {(γ ij , π ij ), (ξ ij , ij )}, where π ij = ∂L ∂γij , ij = ∂L ∂ξij , and L is the HR bigravity massive Lagrangian. All this fields are supposed to be defined at The dynamics of minimal bigravity theory is encoded in the HR bigravity massive action which has the form where g µν and f µν are the two dynamical metric tensors components and R (g) and R (f ) are their respective scalar curvatures in four dimensions. The matter action S m = Ω d 4 x √ det gL m (g, φ) is the matter source containing a scalar field φ, and for the case subject to a potential V(φ), has the form The parameter M eff is an effective mass defined as the harmonic mean of squared field masses M g and M f associated to g and f , respectively: and m is the mass of graviton. The meaning of the tensor g −1 f is such that and the trace χ = tr χ µ ν serves to write the Lovelock-like Lagrangians [9,34,35] L 0 (χ) = 1 Using Eq. (5) and the HR minimal deformation of determinant structure expansion ( [36], Eq.(3.3)) it is obtained the non-trivial minimal action [33] This action is highly non-linear in the lapse and shift variables. However, according with [9], is possible get a linearized version with the introduction of three new variables n i which are functions of the lapse and shift variables and spatial metrics The matrix D = D i j is given by where Q i j = αδ i j + n i n m ξ mj (10) and The matrix D is independent of the lapses N and M , and the shifts N i and M i . With the replacement of N i in favor of n i , the minimal action then becomes linear in the lapses N and M , and the shift M i : where R 0(f ) = det ξ R 3(f ) + 1 det ξ π φ = ∂L ∂φ and , j denotes the derivative ∂ j . The canonical Hamiltonian can be written as and the corresponding equations of motion for N and M yields the g, f -Hamiltonian constraint and M i yield the momentum constraint: and where are the g and f DeWitt metrics, respectively.
In this way, the Poisson brackets between canonical conjugate fields γ ij and π ij in one hand and ξ ij and ij by the other, are given by which are important structures in the process of quantization.

Canonical quantization
In order to apply the canonical quantization formalism in a minimal bigravity theory, is required that the promotion of canonical coordinates γ ij (x), ξ ij (x), φ(x), π ij , ij and π φ be done in the following form: This operators satisfy the commutation relations The constraints (18) and (19) have to be imposed at the quantum level in the form In the coordinate representation, it follows that where This constraint is called the Wheeler-DeWitt equation and in general is not amenable to extract physical information of the system. it is a common practice reduce the number of degrees of freedom. In the following sections its considered a reduction of the model to one degree of freedom working in the so called minisuperspace.

Deformation Quantization and Wigner Function of a Minimal Bigravity Model
There are several methods to obtain the quantization of a classical system. As already made it in the previous section, the first choice is starting with the canonical quantization developed by Heisenberg, Schrödinger, Dirac and many others [37]. By other way, with the path integral formulation established by Feynman [38], another set of technical and conceptual advances has been acquired in formal quantization. Additional to this paradigms, another techniques has been constructed in the phase space, as the formalism of geometric quantization and deformation quantization, which emerges to deal with some problems of the previous methods [39]. The Weyl-Wigner-Groenewold-Moyal (WWGM) formalism is an approach to the quantization of phase space [40,41,42]. Under this formalism, the Wigner function is the most important object; this function contains all the quantum information of the system. The key proposal of WWGM formalism is the connection between classical and quantum dynamics through the formal identification of usual canonical quantum observables with their classical counterparts.
With the aim of introduce the Weyl-Wigner-Groenewold-Moyal (WWGM) formalism for gravitational systems, first note that this is possible in a direct way since the minisuperspace in the minimal bigravity model is flat, and this ensures the existence of the Fourier transform, a key prerequisite. In this section a detailed description of WWGM formalism is left aside; the deformation quantization of gravity in ADM formalism and constrained systems, and in other cosmological models can be found in more detail in Refs. [43,44,30,31].

The Stratonovich-Weyl Quantizer
Let F [γ ij , π ij ; ξ ij , ij ; φ, π φ ] be a functional on the Wheeler phase superspace Γ * associated to hypersurface Σ and letF [µ ij , ν ij ; ζ ij , θ ij ; κ, λ] be its Fourier transform. By a pure analogy with the quantum structure, the Weyl quantization rule is defined aŝ where the functional measures are given by and so on, and the family of unitary operators 6 are expressed bŷ withγ ij ,π ij ,ξ ij ,ˆ ij ,φ andπ φ field operators defined bŷ As it can be expected, this states form a complete basis while operators satisfies the usual canonical commutation relations. Using Campbell-Baker-Hausdorff theorem [45], the completeness relations and standard formulas, Eq. (29) can be written in the form where denotes the deformation parameter which for non-null values resembles the deformation quantization. Following properties emerge from (31): and where trÂ is the trace of operatorÂ in some convenient representation.
Using the previous results, the Weyl quantization rule can be written aŝ where the operatorΩ 7 is the called Stratonovich-Weyl quantizer, which satisfies the relationŝ The Stratonovich-Weyl quantizer is useful to get the inverse map of W: and using (31) on (35), a most convenient form of Stratonovich-Weyl quantizer it is obtained:

Moyal Star Product
The introduction of Moyal -product obey the need of a (noncommutative) product of usual functionals instead operators.
be two functionals defined on Γ * with associated field operatorsF andĜ, respectively: The functional related with operator productFĜ is denoted by (F G)[γ ij , π ij ; ξ ij , ij ; φ, π φ ] and have the realization where the bidiferential functional operator ← → P is defined by The Eq. (39) is know as the Moyal -product for fields.

The Wigner Functional
The fundamental concept in deformation quantization formalism is the Wigner function, so the Wigner functional takes the same role in a straightforward generalization for infinite degrees of freedom. Ifρ is the density operator associated to a quantum state, then the Wigner functional ρ[γ ij , π ij ; ξ ij , ij ; φ, π φ ] corresponding to ρ reads as In case of a pure state of the systemρ = |Ψ Ψ|, Eq. (41) implies where γ ij , ξ ij , φ|Ψ = Ψ[γ ij , ξ ij , φ] is the wavefuntion of universe. The expected value for a general operatorF can be expressed throughρ as Follow the prescription to construct the equations of evolution of the quantum system, it is possible to write the equivalent of Hamiltonian constraint in terms of Wigner functional and Moyal -product as Eq. (44) is called the Wheeler-DeWitt-Moyal equation, and is the deformation quantization version of Wheeler-DeWitt-Schrödinger equation. The first point take into account is that the dynamics of the system is entirely determined by this equation and their explicit form depends of the particularity of the problem.

Quantum Cosmology of a Minimal Bigravity Model in the Minisuperspace
Although the general construction exposed in Secs. 2 and 3 in terms of functional integrals defined in the whole superspace of three-metrics is necessary in the more general case, the restriction to work in minisuperspace is natural due to its simplicity and because has been widely studied in the literature by other methods. The principal goal is to obtain the Wigner function for a minimal bigravity model and to motivate a further study using deformation quantization. This first step is necessary to gain some experience to eventually extend the deformation quantization to curved phase spaces. It is known that deformation quantization formalism admits a natural extension to these scenarios and allows suitable generalizations.

Minimal Bigravity Cosmology in the Minisuperspace
For the development of a standard cosmological theory, it is need to obtain the point-like Hamiltonian of this theory assuming two Friedmann-Lemaître-Robertson-Walker-like (FLRW) dynamical line elements such that they are SO(4) invariant metrics in a Ω = R × S 3 topology in a homogeneus and isotropic spacetime: 9 and where N (t) and M (t) are the lapse functions, a(t) and b(t) are the scale factors of metrics g µν and f µν , respectively, κ is the space curvature, assumed to be the same for both metrics, and l p = 2/3L p , where L p denotes the Planck length. So, using the expressions of R 3(g) = 6κ/a 2 , R 3(f ) = 6κ/b 2 and the metrics above on the field equations, and using the assumption M 2 g = M 2 f = M 2 eff /2, it is obtained a point-like gravitational Hamiltonian in the minisuperspace {N, a, M, b}: where µ := ml p and H m is the matter Hamiltonian. As mentioned in Sec. 2.1, in ADM formalism the Hamiltonian constraint can be written as H 0 = 0. In canonical quantization, this constraint in the form of a operatorĤ is the Wheeler-DeWitt equation, where the wavefunction of the universe Ψ belongs to kerĤ. On the other hand, it can be noted that the matter term in the action is independent of the modifications of the bigravity model by mass term. So, the matter part of the Hamiltonian can simply added to the geometric one to obtain the total Hamiltonian expression. In order to model the matter part Hamiltonian, the equation of state for a perfect fluid matter it is used: In the Schutz representation [46,47] for a perfect fluid, the matter part Hamiltonian can be expressed as where the conjugate momentum p T associated to dynamical variable T (the collective thermodynamic parameters of perfect fluid) shows a convenient interpretation as a time-related magnitude. So, including this as dynamical parameter, the Hamiltonian constraint is expressed as As it is seen from construction, only the matter is coupled with the metric g µν . Taking the Hamilton canonical classical equationṡ it is possible to calculate the equations of motion corresponding to the minisuperspace variables N and M : or defining the Hubble-like parameters H =ȧ/N a and K =ḃ/b, There is an interesting discussion about Eqs. (55) and (56). A comparision with the standard Friedmann equation reveals that the reduced mass of graviton µ 2 plays the role of a comological constant Λ, provided b a takes a constant value, which can be positive (negative) for b > a and negative (positive) if b < a. This emerging parameter can be interpret the dynamical competition between the two scale factors as a phase transition from deceleration to acceleration era. So, to taking in account this point, is natural to make the assumption σ := b/a = const and, with the use of the momentum constraint (19), write and then, p b = σp a . With all this, Eq. (50) then can be written as where p := p a for simplicity.

Canonical Quantization in the Minisuperspace
In the canonical quantization formalism, the Wheeler-DeWitt equationĤ 0 Ψ = 0 is the operator version of the Hamiltonian constraint H 0 = 0; the solution of this equation is the wavefunction of universe Ψ, which encodes all the information of the evolution of universe. So, after the replacements p →p and p T →p T over Eq. (58), the Wheeler-DeWitt reads aŝ Eq. (59) admits the separation Ψ(a, T ) = e iET ψ(a), with E a constant. With this, Wheeler-DeWitt equation takes the form As it can see, the coupling between the matter source and the gravitation g µν states the identification of two metrics g µν and f µν and this reduces from two-variable system to one-variable system in favor just the scale factor a, getting a formal reduction Ψ(a, b) ≡ Ψ(a, σ) ≡ ψ(a). To avoid further confusion, Eq. (61) is referred as Wheeler-DeWitt-Schrödinger (WDWS) equation.

Exact Solution
As mentioned by Darabi and Mousavi [33], in this representation the conjugate momentum associated with T appears linearly in the Hamiltonian of the model, so without loss of generality T can take the role of timevariable in the Wheeler-DeWitt-Schrödinger equation (61). The emergency of a Schrödinger-like equation in this model give us useful information about T -parameter system evolution. On the other hand, through T the factor order operator for the operatorsâ andp a can be considered by means of Using Eq. (62), the next equation it is obtained: Taking ω = −1 and Taking the curvature κ = 0 in agreement with the current observations [48], the solution takes the form α are the Hankel functions of first and second kind, respectively. Its important to note that the definition of satisfy the constraint 3E < 2µ 2 (3 − 4σ + σ 4 ) [33], which ensures a (bounded) oscillatory wavefunction.
The parameter λ do not modify the wavefunction in semiclassical regime. Without the presence of λ, i.e. for the value λ = 0, Eq.(65) resembles a one-dimensional Schrödinger equation for a particle modeled by a a(t) coordinate with zero energy moving in a potential V (a) = − 2 a 4 , and this is equivalent to the asymptotic approximation [49] for values of a 0. As is usual in quantum cosmology [20], is precise get ψ(a) in the form ψ(a) = c 1 ψ + (a) + c 2 ψ − (a), where ψ ± (a) are the WKB solutions for (65).

WKB Semiclassical Approximation
In order to get the WKB description, a suitable potential function is needed. However, the potential V (a) = − 2 a 4 has no barrier; so, in this case, to get an expression to compare with Eq.(65) is necessary perform the separation ψ(a) = a − λ 2 ϕ(a). So, it is obtained the equation with the effective potential and the turning point at a 0 = − λ 2 2 1 − λ 2 1 6 (as can be seen in Fig. 1). In the classical allowed region a > a 0 the WKB solutions to Eq.(68) is given by where Taking this into account, it is possible express WKB solutions ψ ± (a) for a > a 0 as Using the asymptotical expansions for the Hankel functions [50], it is obtained by comparation between Eq.(74) and Eq.(67) that the exact solution in terms of expanding and contracting wavefunctions is where J α is the first kind Bessel function. For the characteristic value α = 1 3 , it follows where Ai and Bi are the Airy function of first and second kind, respectively.

Boundary Conditions
The wavefunction (76) need to satisfy specific boundary conditions prescribed by an particular physical environment. In quantum cosmology, there is a boundary condition setting problem [21], motivated by the different evolution paradigms in cosmology. For this reason, this paper studies three of the most employed boundary conditions: Hartle-Hawking [18], Vilenkin [19,20] and Linde [22] proposals. With this, the aim is analyze some quantum properties of this conditions and their relation with the classical behavior. This three boundary conditions in quantum cosmology are Hartle-Hawking: Vilenkin: Linde: where τ is the turning point. As discussed above, the semiclassical approximation is not fullfilled, but τ can be setted as τ = 0, so Eqs. (78)-(80) can be written for a > 0 as and taking Eqs. (51) and (52) in account, the Eq.(84) gives the interpretation for ψ ± (a) that negative values of p a correspond to an expanding universe. Therefore, the Vilenkin wavefunction represents tunneling and includes only a outgoing term corresponding to a expanding universe while Hartle-Hawking and Linde wavefunctions includes expanding and contracting universes with equal weight. In fact, the explicit form of this wavefunctions in WKB approximation are given by and where S(a) is given by Eq. (73). On the other hand, its exact forms are given by and As before, for α = 1 3 , and Its important to note that there is an open discussion about which could be the correct boundary condition in quantum cosmology. An exploration of this problem in the context of deformation quantization, a formalism beyond canonical quantization, will be developed in next section.

Deformation Quantization in the Minisuperspace
The next step in deformation quantization is apply the non-commutative Moyal -product in the construction of Wheeler-DeWitt-Moyal equation (44) in order to get the Wigner functional in a minimal bigravity model. The exposition in previous section put the focus on functional integrals defined in the whole superspace of three-metrics; however, it is unnecessary take this approach because the initial restriction to work at the minisuperspace.

Exact Solution: Hartle-Hawking Case
With the restriction of work in minisuperspace it is possible to search the Wigner function and rewrite the Hamiltonian constraint (50) in terms of the Moyal -product and the Wigner function as where the corresponding Moyal -product is given by or recalling the assumptions made it to get Eq. (65), aside the change of variables p := p a , q := 2 a 4 , Using the Bopp shift relation [51], results in If a new variable z is defined as z := p 2 − q, then the imaginary part of Eq.(98) is identically equals to zero and the real part takes the form whose solution is Eq. (99) indeed admits another solution corresponding to Bi 2 2 3 (p 2 − 2 a 4 ) , but in this case, the Airy function of second kind Bi(x) is divergent for x > 0 and this part cannot be included as a suitable Wigner function. With the purpose of obtaining the next result, it is necessary the assumption of the existence of the Wigner transform. As is pointed in Refs. [30,52], there is exist a problem with the range of integration over all real line for the variable q (and then the scale factor a). So, in order to avoid this and get a solution to (44), it is convenient use the following integral representation of Wigner function: Using the convolution theorem and the Fourier transform for the Airy function is possible get the Wigner function for the Hartle-Hawking wavefunction ρ W (a, p) = 3 √ 2π Eq.(102) is the principal object in deformation quantization: essentially, provides all the physical information of the system, in this case, the dynamics of a universe in a minimal bigravity in minisuperspace phase space. To describe the coincident classical trajectories and compare with classical behavior [64,53], is necessary to consider the Hamilton equations for a and p associated to classical Hamiltonian constraint (58) (for p T → E, κ = 0 and ω = −1):ȧ Eqs. (103) and (104) defines the stationary vector field in phase space (ȧ,ṗ) = (2p, 4 2 a 3 ), which is showed in Fig. 3. The locus [54,55] of extrema of Wigner function (102) is the deformed version of the classical Hamiltonian given by where a n is the n-th zero of the derivative of Airy function, d/dx Ai(−x)| x=a n = 0 (see Table 1). Eq. (105) presents the most probable cosmological solutions. Before continue with Vilenkin and Linde cases, some remarks are necessary about Hartle-Hawking Wigner function. As is pointed in [56], Wigner functions may exhibit a cusp catastrophe behavior in nonclassical regions of phase space. Catastrophes of ρ W occurs at points (a, p) where the stationary phase approximation diverges, and resembles the caustics of families of trajectories in configuration or momentum space, although, due to Liouville theorem, there are no caustics in phase space. If the complete potential (66) is extended with a control term V c = λa, it obtains the prototypical cusp catastrophe potential In Fig. 4 is showed the cusp catastrophe plot, a surface where a satisfying dV (a)/da = 0 for parameters α = 4κ/ 2 and β = λ/ 2 , drawn for parameter β continuously varied and for several values of parameter α. Outside the cusp locus of bifurcations, for each point (α, β) in parameter space there is only one extreme value of a. Inside the cusp, there are two different values of a giving local minima of V (a) for each (α, β), separated Table 1: Zeroes a n of Ai(−x) and zeroes a n of d dx Ai(−x) for n = 1, 2, 3, 4 n a n a n by a value of a giving a local maximum. From these parameters, it is possible get a qualitative bound for κ and in terms of stability criteria. A qualitative dimensional comparation between (106) and (105) suggest that V c can be related with 2 2 3 a n and then λ = 2 2 3 . On the other hand, the trichotomy of parameter α = 4κ/ 2 (α < 0 for two unstable solutions, α = 0 for bifurcation point and α > 0 for one stable solution) shows that in this cosmological model, definite-positive values for curvature κ implies stable solutions, giving another argument for the election of flat geometries in a minimal bigravity theory.

Numerical Analysis of Wigner Function
Eq. (102) is obtained integrating out from minus to plus infinity in variable a, even though the wavefunction only is defined for positive values of the scale factor. Because the exponential decay of Airy function Ai(x) for x > 0, the wavefunction and its associated Wigner function are very similar. As is mentioned in [52], the restriction to positive values for the calculation of an analytical expression from Eq. (101) is very difficult even using the WKB approximation, so it is necessary perform a numerical integration. Through an implementation in Python language the Wigner functions for Hartle-Hawking, Vilenkin, and Linde wavefunctions are calculated and the results depicted in Figs. 5, 7, and 9. Their respective density-contour plots are showed in Figs. 6, 8, and 10, respectively.
After obtaining the Wigner functions for the different boundary conditions, some remarks are necessary. All this plots can be described by the WKB approximation for the wavefunctions, as it is pointed out in Refs. [27,30,31].
For the Hartle-Hawking Wigner function an oscillatory symmetric form around zero momentum is observed (see Fig. 5), where its largest peak corresponds to the classical trajectory given by which are represented by thick lines in the density plots (see Figs. 6, 8, and 10) for the characteristic value = 1. It is observed that the extrema (maxima and minima) of Wigner function decrease with the distance to the classical trajectory. Besides the Wigner function has negative values which are interpreted by some authors [57,58] as nonclassical effects (similar characteristics are presented in the other cases).
The Vilenkin Wigner function (see Fig. 7), is constructed with the wavefunction corresponding to an expanding universe. In this case there are not interference terms, as the presented in Hartle-Hawking case. This behavior corresponds to the prototypical tunneling process. For this reason the highest peaks of the Wigner function are near to the classical trajectory associated with negative values of the momenta (expanding universe; see Fig. 8).
The Linde Wigner function (see Fig. 9) has a similar behavior as the Hartle-Hawking case, but the classical trajectory is not on the highest peaks (in fact, seems that is on the minima of Wigner function), as it is showed in Fig. 10 For Hartle-Hawking and Linde cases the Wigner function have interference terms corresponding to a contracting and expanding universe; the difference between them is an opposite sign in the interference term of the wavefunctions (78) and (80).
From point of view of quantum cosmology, this results suggest that the classical limit is difficult to achieve since there exist oscillations of the Wigner function for the three cases considered. Until the classical trajectory is on the highest peaks in Hartle-Hawking Wigner function, the decoherence of the Vilenkin function could be easier to obtain because it has fewer oscillations to the other cases due the interference terms not being present since there is only an expanding universe around p = 0. Besides, with the presence of a tunneling scenario, it is possible avoid the classical initial singularity by means the probability for creation from nothing at a small nonvanish scale factor. By extremizing the marginal it is found that the size of initial scale factor corresponding to the maximum contribution of |ψ(a)| 2 is determinated by the frequency term and, as it can seen from Eq. (64), by the mass of graviton m. Explicitly, where G is the gravitational constant, m P is the Plank mass and a 1 is the first zero of Airy function Ai(−x) (see Table 1). In the context of minimal bigravity, the mass of graviton besides play a role of cosmological constant [33] (taking a constant σ = b a = 1 in account). So, in concordance with [27], the deviation of classical behavior of Wigner functions is in function of the cosmological constant and therefore, in function of m. In addition, Sasaki et al. [59] have shown that at early universe for Hartle-Hawking boundary conditions, dRGT massive gravity give a dynamical character to the mass of graviton. They also shown that massive gravity sector improve the probability of a large number of e-foldings for a sufficiently large graviton mass comparable with the Hubble parameter during inflationary stage of order m ≈ 10 12 GeV.

Conclusions
The phase space quantum description of a physical system is a construction where the coordinates and momenta coexists simultaneously and in where some physical properties can be appreciated in a more complete form with respect to a particular representation. The key piece in this approach is the Wigner function which contains all quantum information in the phase space. The Wigner quasi-probabilistic distribution and their properties are useful in the description of semiclassical effects and to obtain the classical limit.
In this paper the quantum phase space description was employed to analyze the physical behavior of a nonlinear massive minimal bigravity model. In sake of simplification, Friedmann-Lemaître-Robertson-Walkerlike metrics with common curvature κ were supposed. Due the presence of interaction between two scale 19 factors in the massive bigravity action and by means of momentum constraint, it was possible the reduction of the two-variable minisuperspace into a one-variable minisuperspace including one scale factor, a. With the particular selection of perfect fluid ω = −1 and flat curvature κ = 0, corresponding to an early universe, the Wheeler-DeWitt-Schrödinger and Wheeler-DeWitt-Moyal equations were obtained. For the Wheeler-DeWitt-Schrödinger an exact and a WKB solution were found, meanwhile by numerical analysis the Wigner function associated to a three distinct types of boundary conditions, Hartle-Hawking, Vilenkin and Linde were calculated, and also an exact analytical Wigner function was found for the Hartle-Hawking boundary case.
The behavior of the Hartle-Hawking Wigner function has an oscillatory symmetric form around zero momentum and presents high frequency interference terms between expanding and contracting universe, with the highest peak corresponding with the classical trajectory.
The Linde Wigner function shows a very similar behavior with the Hartle-Hawking case, where there are also expanding and contracting terms of the wavefunction. As it can seen, there are fewer oscillations outside the region of phase space bounded by the classical trajectory, interchanging maxima with minima. This behavior is due to the sign difference in the interference term, and it has as a consequence that the classical trajectory is not on the highest peak of the Wigner function.
For the Vilenkin tunneling boundary condition its associated Wigner function shows oscillations that rapidly damp since there is only one term of expanding universe, that can be explained because there are non-interference terms. Besides, the classical trajectory is closer to the maxima of Wigner function and the process of decoherence seems to be feasible to obtain because the damping of interference terms.
These results indicates that in a minimal bigravity model of universe the classical limit is difficult to obtain for the three boundary conditions as a consequence of the oscillations in Wigner function, despite the Vilenkin Wigner function could be simpler to obtain because there is only an expanding universe. In addition, a relation between this oscillations with frequency and the mass of graviton m through the cosmological constant interpretation indicates the role of the mass in quantum cosmological bigravity models. Also, a relation between curvature κ and stability of Wigner function through catastrophe cusp plot found, giving an argument to discriminate possible geometries.
The phase space quantum analysis developed in this work drive the study to other minimal bigravity models where the curvature are not flat or where the topology of space is nontrivial [60,61]. For example, it can be applied to analyze the properties in changes of topology [62,63,21] by means of the Wigner function, where the case of the different boundary conditions can be worked with more detail.
Besides, deformation quantization formalism is up to now the most general method to quantize arbitrary systems [39], and could guide to obtaining analytical expressions for the Wigner function of a wide class of cosmological objects. The previous approach could allow new results and some would be analyzed in the classical limit in a more systematic way. Both the nontrivial topology spaces cosmology quantization an the issues corresponding to the classical limit are the subject to current research. In future works eventually it will discussed the application of deformation quantization formalism to curved phase spaces and other spaces without trivial topology, and an research of this work under the Wigner flow analysis [64,53,54]  as a natural extension.