Similarity solutions for the Wheeler-DeWitt equation in $f\left( R\right) $-cosmology

In the case of a spatially flat Friedmann--Lema\^{\i}tre--Robertson--Walker universe in $f\left( R\right) $-gravity we write the Wheeler-DeWitt equation of quantum cosmology. The equation depends on the functional form of $f\left( R\right) $. We select to work with four specific functions of $f\left( R\right) $, in which the field equations for the classical models are integrable and solvable through quadratures. For these models, we determine similarity solutions for the Wheeler-DeWitt equation by determine Lie-B\"{a}cklund transformations. In addition we show how the classical limit is recovered by the similarity solutions of the Wheeler-DeWitt equation.


INTRODUCTION
Modified theories of gravity [1,2] are an alternate approach to the dark energy models to explain recent observational phenomena [3][4][5][6]. The common characteristic of the modified theories of gravity is the modification of the Einstein-Hilbert Action by adding new invariant terms to the gravitational Action. The novelty of that approach is that new geometrodynamical components are introduced into the field equations which drive the dynamics to explain the observations.
In the literature there have been proposed a plethora of different modified theories of gravity. A specific class of models which have drawn the attention are the so-called f −theories. In f −theories of gravity a function f (Q) is introduced into the Einstein-Hilbert Action, where Q is an invariant function. Some theories which belong to that class of models are the: f (R)-gravity in the metric formalism [7], f R -gravity in the Palatini formalism [8], f (G)-Gauss Bonnet theory [9], while in the Teleparallel formalism of gravity the f (T ) −theory has been widely studied the last decade [10][11][12][13]. For other modified theories which belong to that class we referee the reader to [14][15][16][17][18][19][20][21][22][23][24][25][26][27][28] and references therein.
In this work we are interested in f (R)-gravity in the metric formalism [29], where Action Integral in a fourdimensional manifold is given by the expression S = dx 4 √ −gf (R). In this theory, variable R corresponds to the Ricci scalar of the underlying geometry with line element g µν ; consequently, General Relativity with or without the cosmological constant is fully recovered when f (R) is a linear function. Various specific functional forms of f (R) −theory have been proposed in the literature in order to describe the various phases of the universe. The quadratic model f (R) = R + αR 2 can describe well the inflationary era of our universe [30,31]. The natural extension of the latter inflationary model is the f (R) = R + αR n model [36] which provides power-law attractors. For other f (R)-models with applications in the late acceleration phase of the universe see [32][33][34][35] and references therein.
f (R) −theory is a fourth-order theory and it is dynamical equivalent to the Brans-Dicke theory with zero value for the Brans-Dicke parameter. The scalar field attributes the extra degree of freedom such that the theory is written as a second-order theory but with extra dependent variables so that the total degrees of freedom are the same. The theory is nonlinear and there are few exact solutions, either for spacetimes with one free function such as the Friedmann-Lemaître-Robertson-Walker metric (FLRW) which is usually applied in modern cosmology. Indeed, in the case of a spatially flat FLRW spacetime the de Sitter solution, R = R 0 , is recovered when there exists a solution to the algebraic equation R 0 f ′ (R 0 ) = 2f (R 0 ) [31]. In addition, power law solutions, which describe an ideal gas with constant equation parameter, are recovered when f (R) = f 0 R n , n = 0, 1, 2. However, the latter exact solutions do not describe the generic analytic solution for the corresponding field equations because they are valid only for specific initial conditions. Some analytic solutions have been found by searching for conservation laws for the field equations and making a conclusion about the integrability of the gravitational model by writing the analytic solution with the use of closed-form functions or making use of theorems from the theory of Analytic Mechanics, for instance see [38][39][40][41].
We focus on the determination of exact solutions of the Wheeler-DeWitt (WdW) equation [42] in f (R)-cosmology. The WdW equation is mainly applied in quantum cosmology. Recall that in modern cosmology we assume that the spacetime is described by the FLRW metric with zero spatial curvature. WdW is an equation of Klein-Gordon type, where the dependent variable is denoted to describe the wavefunction of the universe and the independent variables are the dynamical variables of the classical system. There are various issues such that there is not a unique way for one to define probability [43,44]. Also there is the so-called problem of time, because time is involved in the wavefunction through the dynamical variables [45][46][47][48].
A previous analysis of the exact solutions of WdW in f (R)-cosmology was published in [49][50][51]. Specifically in [49] there was found that for the special power-law theory f (R) = R 3 2 , the classical solution can be recovered from the solution of the WdW equation. The case f (R) = R 3 2 describes an integrable cosmological model which admits a conservation law linear in the momentum. That approach has been extended and applied in other gravitational models, such as anisotropic universes [44], static spherical symmetric spacetimes [52][53][54], inhomogeneous spacetimes [55] and electromagnetic three-dimensional pp-wave spacetimes [56].
In our consideration we determine a family of Lie-Bäcklund transformations for the WdW equation for some specific models of f (R)-cosmology. The models of f (R)-cosmology that we study form integrable dynamical systems where the conservation laws which ensure the integrability are constructed by point transformations which leave the variational integral invariant. The plan of the paper is as follows.
In Section 2 we present the basic equations of f (R)-cosmology. The main mathematical materials necessary for the analysis of the present work are given in Section 3. Specifically, we show how Lie-Bäcklund transformations can be constructed for the conformally invariant Klein-Gordon equation by using the point symmetries of the classical Hamiltonian system. In addition we show how the Lie-Bäcklund operators are applied in order to determine similarity solutions for the WdW equation. The context of the one-dimensional optimal system is discussed. Section 4 includes the main material of our analysis. For four integrable classical models of f (R)-cosmology we write the WdW equation and we determine the infinitesimal generators of the point transformations where the WdW equation is invariant. From the infinitesimal generators we construct the Lie-Bäcklund operators and we find the similarity solutions. In order our results to be completed the one-dimensional optimal system is determined for each model. For the models of our study we observe that the classical limit is always recovered. Finally in Section 5, we discuss our results and we draw our conclusions.

f (R)-COSMOLOGY
For a spatially flat FLRW background space with line element and Ricciscalar the gravitational field equations of f (R)-gravity are calculated to be [7] 3f where H (t) is the Hubble function, H (t) =ȧ (t) a(t) , dot indicates derivative with respect to the independent variable "t"; and a prime denotes derivative with respect to the Ricciscalar, that is, f ′ (R) = df (R) dR . The latter field equations can be written in an equivalently form as follows [7] where now G µ ν , is the Einstein tensor, k ef f is a varying "Einstein-constant" defined as k ef f = 1 f ′ (R) , and T µ f ν is the effective energy momentum tensor which attributes the geometrodynamical degrees of freedom of the higher-order of gravity. Indeed, the energy-momentum tensor T µ f ν is defined as where the energy density ρ f and pressure term p f are defined as [7] ρ Hence, the field equations are while the equation of state parameter for the effective fluid Note that the latter expression for f (R) = R − 2Λ gives w f = −1 which means that the theory of General Relativity with the cosmological constant is recovered.

Minisuperspace approach
The gravitational field equations (2), (3) and (4) can be derived by a variation principle of the Action integral A = L N, a,ȧ, R,Ṙ dadRdN (10) where L N, a,ȧ, R,Ṙ is defined as [39] L N, a,ȧ, where N (t) is a generic lapse function for the FLRW metric, such that the Hubble function is defined H (t) = 1 Nȧ a . We note that Lagrangian (11) is a singular Lagrangian since ∂L ∂Ṅ = 0. Lagrangian (11) defines a constraint system, with constraint equation ∂L ∂N = 0. The two-second order order differential equations follow by the variation with respect to the variables {a, R}, that is, d dt ∂L ∂ȧ − ∂L ∂a = 0 and d dt ∂L where The second-rank tensor G AB defines the space where the dynamical variables {a, R} evolve. We can define a canonical momenta for the variable {a, R}, and write the point-like Lagrangian (11) as a Hamiltonian system. It follows that the two momentum p a = ∂L ∂a , and p R = ∂L ∂Ṙ are so the Hamiltonian function is or equivalently where P A = (p a , p R ) is the canonical momentum. Hence, the constraint equation provides while the rest of the field equations are given by the Hamilton equationṡ that is 1 and 1

Wheeler-DeWitt equation
Constraint equation (17) yields the WdW equationĤΨ(q) = 0, whereĤ is the Hamiltonian operator under canonical quantization, in which ∆ L is the conformal Laplace operator defined as where R is the Ricciscalar of the minisuperspace G AB and n = dim G AB and ∆ is the Laplace operator, that is, For the second-rank tensor (13) we calculate n = 2, which means that ∆ L = ∆. The conformal Laplace operator ∆ L has the property that it is invariant under conformal transformations, such a requirement it is necessary in the case of quantum cosmology since the theory should be conformal invariant because of the arbitance of the lapse function N (t). While in the case where n = 2 operator ∆ L follows from the canonical quantization P A ≃ 1 √ G ∂ ∂q A , for higher-dimensional spaces, n ≥ 3, the conformal Laplace operator ∆ L follows from the canonical quantization only for conformally flat spaces, and in general the term n−2 4(n−1) R should be added by hand. However, which quantization process which provides the operator ∆ L for n ≥ 3 from a point-Lie Hamiltonian function is still an open problem.
In general and in terms of the 1+3 decomposition notation of GR the WdW equation it follows from the Hamiltonian constraint where G ijkl is defined as is the the metric of superspace, the space of all 3-geometries with metric h ij and Ricci scalar R, and the matter configuration. We note that in general the WdW equation (26) is a hyperbolic functional differential equation on superspace, where in the case of the minisuperspace approximation it is reduced to a single equation for all the points of the superspace.
As far as our model of f (R)-cosmology is concerned, with constraint equation (15), the WdW equation in the minisuperspace approach is found to be For the latter linear second-order partial differential equation we shall determine exact solutions for specific forms of f (R) function. In the following section we present the main mathematical tools which will be applied in order to determine solutions for equation (28).

CONSTRUCTING SIMILARITY SOLUTIONS
Consider the partial differential equation H q A , Ψ, Ψ ,A , Ψ ,AB , ... ≡ 0 where q A = q 1 , q 2 , ..., q n denotes the n−independent variables and Ψ = Ψ q A is the dependent variable.
When the differential equation H q A , Ψ, Ψ ,A , Ψ ,AB , ... is invariant under the infinitesimal one-parameter point transformation q n → q n + ε, then we shall say that the differential equation admits the Lie point symmetry X = ∂ q n and vice versa.
In general, the differential equation H q A , Ψ, Ψ ,A , Ψ ,AB , ... is invariant under the infinitesimal one-parameter point transformation if and only if there exists a function λ q B , Ψ such thatξ where The Lie symmetries for the conformal invariant Klein-Gordon equation (23) have been studied before in [60]. Specifically it has been found that the generic Lie symmetry has the form in which a 0 is a constant, b q A is a solution of the original equation (23) and represents the infinity number of solutions, since the equation is linear, and ξ A q A is a conformal vector field for the minisuperspace G AB q C , with conformal factor ψ q B , that is, L ξ denotes the Lie derivative with respect to the vector field ξ.
In addition, the conformal vector field and the potential U(q C ) satisfy the constraint condition By definition, if X = ξ A q B , Ψ ∂ A + η q B , Ψ ∂ Ψ is a Lie point symmetry for the differential equation H q A , Ψ, Ψ ,A , Ψ ,AB , ... , the symmetry vectorX = η q B , Ψ − ξ A q B , Ψ Ψ ,A ∂ Ψ is a Lie-Bäcklund symmetry. Vector fieldX is the canonical form of the vector field X.
A Lie-Bäcklund symmetry preserves the set of solutions for the differential equation, that is, Condition (33) provides a constraint equation which will be used in the following to solve the WdW equation (28). The symmetry condition (32) where ξ A q B is a conformal vector field of the minisuperspace has been found before in [38] for the variational symmetries for singular Lagrangians of the form of (11). Indeed, for every variational symmetry of (11) a Lie point symmetry and consequently a Lie-Bäcklund can be constructed for the WdW equation (28). However, that it is not the only relation between variational symmetries of classical Lagrangians and Lie symmetries of the WdW equation.
If we consider that the lapse-function N (t) in the Lagrangian (11) is fixed, then we can apply the results for the variational symmetries of regular Lagrangians [39]. As we shall see for the f (R)-theory we recover the results of [61] while also we found new Lie-Bäcklund operators which will be used to determine new similarity solutions for equation (28).
However, the Lie point symmetries of regular systems can be time-dependent, something which is not true for the WdW equation. Below we show two cases of important interest where we show how Lie-Bäcklund operators are constructed by using the time-dependent symmetries of the regular Lagrangian.

Higher-order Lie-Bäcklund operators
We show, for two models of special interest, how to construct Lie-Bäcklund operators for the conformal Laplace equation (23) by using time-dependent point symmetries of regular Lagrangians.

Oscillator
Consider the point-like regular Lagrangian where h AB y C and F y C are arbitrary functions. Lagrangian (34) admits the variational symmetries, ∂ t and e ±µt ∂ x , with respective gauge function f t, x, y A = µe ± x. Consequently, from the two-latter variational symmetries for the dynamical system with Lagrangian (34) we can construct the time-dependent conservation laws It is easy to show that the combined integral I 0 = I + I − is time independent and equals The corresponding conformal invariant Klein-Gordon equation is Equation (37) does not admit any Lie point symmetry for general h AB , F y C while R y C is the Ricciscalar for the metric h AB .
We observe that equation (37) is separable with respect to x. Indeed the solution can be written in the form Ψ x, y A = w (x) S y A . This implies that the operator

Ermakov-Pinney system
The second case we consider is that of the Ermakov-Pinney system. Let us assume the generic regular Lagrangian function where h AB y C and F y C are arbitrary functions. The dynamical system described by the Lagrangian (40) admits the time-dependent conservation laws where h is the value for the integral of motion described by the Hamiltonian for Lagrangian (40). In a similar way as before we construct the autonomous first integral [62] Φ 0 = h 2 − I + I − , which equals This is the well known Ermakov invariant, also known as Lewis invariant. Consider now the conformal invariant Klein-Gordon equation where R y C is the Ricciscalar of the metric h AB y C . The latter equation does not have any Lie point symmetries. However, the latter Klein-Gordon equation is separable, in the sense that Ψ r, y C = w (r) S y C . Then we shall say that the operatorΦ satisfies the equationΦΨ = 0 which means that the Klein Gordon equation (45) admits the Lie-Bäcklund symmetry with generatorX 3.2. One-dimensional optimal system However, Lie symmetries are used to find new similarity solutions for other similarity solutions by applying the adjoint representation of the admitted Lie group for the given differential equations. Hence, it is important to determine the one-dimensional optimal system for the admitted Lie algebra for the equation of our study. In that case we will determine all the unique similarity solutions which can not derived by adjoint transformation. In the following we give the definition of the adjoint operator as also when two Lie point symmetries are connected through the adjoint representation. Let a given differential equation H q A , Ψ, Ψ ,A , Ψ ,AB , ... to admit a n-dimensional Lie algebra G n with elements X 1 , X 2 , ... X n . Then we shall say that the two vector fields [58,59] are equivalent if and only if W = lim n j=i Ad (exp (ε i X i )) Z or bi ai = c , c = const, in which Ad (exp (ε i X i )) is the adjoint operator defined as

SIMILARITY SOLUTIONS OF THE WHEELER-DEWITT EQUATION
As we discussed before, in order to solve the WdW equation (28) we will construct differential operators by using the variational symmetries for the classical system. Such an analysis was performed before in [39] where the unknown function f (R) which defines the theory is constrained by the requirement the field equations in f (R)-cosmology to admit conservation laws generated by point symmetries.
In Lagrangian (11) we assume that N (t) = 1. Therefore, for arbitrary function f (R) the dynamical system is autonomous and admits the point symmetry ∂ t . However ,the latter symmetry does not provide any differential operator for the WdW equation (28).
In addition, there are four specific functions of f (R)-function where Lagrangian (11) is transformed such that the variation of the Action Integral (10) to be invariant. Specifically, the cases we shall study are (A) f (R) = R Indeed, model D is the Λ 1 3 2 CDM while model E corresponds to the Λ 1 7 8 CDM. At this point it is important to mention that because the WdW equation (28) is a linear second-order partial differential equation, it admits for arbitrary function f (R) the two symmetry vectors X Ψ = Ψ∂ Ψ and X b = b (a, R) ∂ Ψ , in which b (a, R) is a solution of (28). Symmetry X b denotes the infinity number of solutions for the partial differential equations. However, X b plays no role in the derivation of similarity solutions and for that reason we will omit it.
However, under the change of coordinates {a, R} → {z, w} with the relation a = 9 2 − 1 3 √ z , R = w 2 z the point-like Lagrangian (50) is simplified as follows, Consequently, the field equations in the Hamiltonian formalism become [ The latter system can be easily integrated and the exact solution is presented in [39]. From the Hamiltonian (52) results the WdW equation which admit the Lie point symmetries or in canonical form the Lie-Bäcklund operatorŝ The commutators and the Adjoint representation of the admitted Lie algebra are presented in tables I and II. Therefore, from the Adjoint representation we determine the one-dimensional optimal system Hence, we shall determine seven invariant solutions for the WdW equations (54) which are not related through adjoint transformation By using {X 1 } and X 2 we infer that Ψ (z, w) = 0, which is a trivial solution. On the other hand by using {X 3 } we find where I 0 (x) , K 0 (x) are the modified Bessel functions and Ψ 0 3(1) , Ψ 0 3(2) are constants. In addition, from the symmetry vector {X 1 + γX 2 } we calculate the travel-wave like wavefunction In a similar way, the rest of the similarity solutions are determined to be and We observe that solutionsΨ 1 (z, w) , Ψ 2 (z, w) and Ψ 12 (z, w) are equivalent, hence we have found in total five independent similarity solutions. Because the WdW equation is linear the generic similarity solution by point transformations is written as where the sum is on all the free parameters of the solutions. Recall that no boundary conditions have been applied to constrain the similarity solutions. The boundary conditions in quantum cosmology is still an open problem. However, for the classical system and specifically from (52) the Hamilton-Jacobi equation follows ∂S ∂z ∂S ∂w − w 3 9 = 0 with the constraint equation ∂S ∂z = S 0 , which is nothing else than the conservation lawṗ z = 0. Consequently, the generic solution of the Hamilton-Jacobi equation is which is nothing else than the exponent function of the similarity solutionΨ 1 (z, w). Therefore, we can infer that solutionΨ 1 (z, w) is the one which recovers the classical solution where parameter δ is related with the conservation law p z = S 0 . In addition, we observe that the solution of the Hamilton-Jacobi equation is included in solutionΨ 4 (z, w), but not in the rest of the solutions, namely Ψ 3 (z, w) andΨ 3 (z, w). In Fig. 1 we give the qualitative evolution of the wavefunction Im Ψ 1 (a, R) for δ = i 10 , that is, Im Ψ 1 (a, R) ∼ sin (S (a, R)) where S (a, R) is the solution of the Hamilton-Jacobi equation.
where the point-like Lagrangian takes the simple form The latter Lagrangian describes the two-dimensional Ermakov-Pinney system without the oscillatory term, while constant V 0 has the value V 0 = − 1 42 . The Hamiltonian constraint is is the Ermakov-Pinney invariant, also known as Lewis invariant. As far as the WdW equation (28) is concerned it is calculated to be The later partial differential equation is invariant under the one-parameter point transformations with generators the vector fields where in the canonical forms arê In tables III and IV we present the commutators and the adjoint representation of the admitted point symmetries. From table IV we find that the one-dimensional optimal system to be From the one-dimensional algebras {Y 2 } and {Y 3 } we find the trivial solutions Ψ (ρ, σ) = 0. From the other onedimensional algebras it follows √ γ e 6ζ , ζ = y + 1 6 ln γρ 12 − 1 ρ 6 (71) while from {Y 3 + δY 4 } we get the solution of Ψ 23 (ρ, σ) multiplied by the function exp δ 12γ ρ −6 e 6y . However as we discussed in the previous section for the Ermakov-Pinney system the Lewis invariant (68) can be used to construct the Lie-Bäcklund operator By using the later constraint we find the wavefunction where a 1 , a 2 , b 1 and b 2 are integration constants and J i (x) , Y i (x) are the Bessel functions. We observe that solutions Ψ 1 (ρ, σ) andΨ 1 (ρ, σ) are included in the latter generic solution. In total we have found four different solutions, hence, the generic wavefunction is expressed as In order to relate any quantum solution with the classical universe, we should solve the Hamilton-Jacobi equation (67) with the use of the constraint (68) where p ρ = ∂S ∂ρ and p σ = ∂S ∂σ . We find that In the new coordinates, the Hamiltonian constraint is written while the field equations becomesρ Finally, the field equations admit the Lewis invariant which is written as The WdW equation (28) is written as follows and has no other point symmetries except the trivial ones. However, as we discussed in Section 3 from the Lewisinvariant we construct the differential operator Hence from (97) and (98), withΦΨ = 0 we find the similarity solution where, I a (x) , J a (x), K a (x) and Y a (x) are the Bessel functions. We observe that in order the wavefunction to be total periodic, then Λ < 0. As far as concerns the classical limit, in a similar approach with Model B, that is recovered in the limit where e 6σ → +∞. The qualitative evolution of Ψ LB (ρ, σ) for Ψ 0 LB(2) = Ψ 0 LB(3) = Ψ 0 LB(4) = 0 and for Λ < 0 is presented in Fig. 4.

CONCLUSIONS
In this work we focused on the determination of similarity solutions for the WdW equation in quantum cosmology and more specifically in f (R)-gravity in a spatially flat FLRW universe. The WdW equation is a linear equation of Klein-Gordon class which by definition is conformal invariant. For the cosmological of our consideration the WdW equation provides the solution of the wavefunction Ψ in terms of the two indepedent variables of the theory, the scale factor a (t) and the Ricciscalar R (t).
We recall, that f (R) is a fourth-order theory and the Ricciscalar R (t) has been added as a Lagrangian multiplier in order to attribute the higher-order derivatives, such that the f (R) −gravity to be of second-order but with more degrees of freedom. Because of the latter property the theory is dynamical equivalent with scalar-tensor theories while a point-like Lagrangian description is possible, which is necessary for our approach on the problem.
For the function form of f (R) which specifies the theory, we considered four models which were found before and are integrable by one-parameter point transformations. Two of the models are power-law while the other two models belong to the family of Λ bc CDM. For these specific models we write the WdW equation and we determine the infinitesimal generators of the one-parameter transformations where the WdW equations are invariant. We use the infinitesimal generators to define Lie-Bäcklund operators which are used as constraint equations to solve the WdW equation. These solutions are called similarity solutions.
A novel observation for the solutions that were found by that approach is that in the classical limit, that is, in the WKB approximation, the solution of the Hamilton-Jacobi equation for the classical system is recovered, consequently the classical limit is recovered. We can say that the similarity solutions which provide the classical limit are preferred. Indeed there are not initial and boundary conditions to constrain the solutions of the WdW equation, however by the requirement the similarity solution to provide the classical limit we can construct a family of boundary conditions. Because a similarity solution is invariant under the infinitesimal transformations which have been applied for the determination, the boundary conditions should be also invariant under the same infinitesimal transformations [64,65].
The similarity solutions can be used to define probability, or calculate the quantum potential of Bohmiam mechanics. However such applications is not the scope of the present work and such analysis will be published elsewhere.