A fully quantum model of Big Bang

In the paper the closed Friedmann-Robertson-Walker model with quantization in the presence of the positive cosmological constant and radiation is studied. For analysis of tunneling probability for birth of an asymptotically deSitter, inflationary Universe as a function of the radiation energy a new definition of a"free"wave propagating inside strong fields is proposed. On such a basis, tunneling boundary condition is corrected, penetrability and reflection concerning to the barrier are calculated in fully quantum stationary approach. For the first time non-zero interference between the incident and reflected waves has been taken into account which turns out to play important role inside cosmological potentials and could be explained by non-locality of barriers in quantum mechanics. Inside whole region of energy of radiation the tunneling probability for the birth of the inflationary Universe is found to be close to its value obtained in semiclassical approach. The reflection from the barrier is determined for the first time (which is essentially differs on 1 at the energy of radiation close to the barrier height). The proposed method could be easily generalized on the cosmological models with the barriers of arbitrary shape, that has been demonstrated for the FRW-model with included Chaplygin gas. Result is stable for variations of the studied barriers, accuracy are found to be 11--18 digits for all coefficients and energies below the barrier height.


I. INTRODUCTION
In order to understand what really happens in the formation of the Universe, many people came to the point of view that a quantum consideration of this process is necessary. After the publication of the first paper on the quantum description of Universe formation [14,49], a lot of other papers appeared in this topic (for example, see Refs. [4,19,22,37,[43][44][45]52] and some discussions in Refs. [38,46] with references therein).
Today, among all variety of models one can select two approaches which are the prevailing ones: these are the Feynman formalism of path integrals in multidimensional spacetime, developed by the Cambridge group and other researchers, called the "Hartle-Hawking method" (for example, see Ref. [19]), and a method based on direct consideration of tunneling in 4-dimensional Euclidian spacetime, called the "Vilenkin method" (for example, see Refs. [43][44][45][46]). Here, according to Ref. [47], in the quantum approach we have the following picture of the Universe creation: a closed Universe with a small size is formed from "nothing" (vacuum), where by the word "nothing" one refers to a quantum state without classical space and time. A wave function is used for a probabilistic description of the creation of the Universe and such a process is connected with transition of a wave through an effective barrier. Determination of penetrability of this barrier is a key point in the estimation of duration of the formation of the Universe, and the dynamics of its expansion in the first stage.
However, in the majority of these models, with the exception of some exactly solvable models, tunneling is mainly studied in details in the semiclassical approximation (see Refs. [38,46]). An attractive side of such an approach is its simplicity in the construction of decreasing and increasing partial solutions for the wave function in the tunneling region, the outgoing wave function in the external region, and the possibility to define and to estimate in an enough simply way the penetrability of the barrier, which can be used for obtaining the duration of the nucleation of the Universe. The tunneling boundary condition [46,47] could seem to be the most natural and clear description, where the wave function should represent an outgoing wave only in the enough large value of the scale factor a. However, is really such a wave free in the asymptotic region? In order to draw attention on the increase of the modulus of the potential with increasing value of the scale factor a and increasing magnitude of the gradient of such a potential, acting on this wave "through the barrier", then one come to a serious contradiction: the influence of the potential on this wave increases strongly with a! Now a new question has appeared: what should the wave represent in general in the cosmological problem? This problem connects with another and more general one in quantum physics -the real importance to define a "free" wave inside strong fields. To this aim we need a mathematical stable tool to study it. It is unclear whether a connection between exact solutions for the wave function at turning point and "free" wave defined in the asymptotic region is correct.
Note that the semiclassical formula of the penetrability of the barrier is constructed on the basis of wave which is defined concerning zero potential at infinity, i.e. this wave should be free outgoing in the asymptotic region. But in the cosmological problem we have opposite case, when the force acting on the wave increases up to infinity in the asymptotic region. At the same time, deformations of the shape of the potential outside the barrier cannot change the penetrability calculated in the framework of the semiclassical approach (up to the second order). An answer to such problem can be found in non-locality of definition of the penetrability in quantum mechanics, which is reduced to minimum in the semiclassical approach (i. e. this is so called "error" of the cosmological semiclassical approach).
The problem of the correct definition of the wave in cosmology is reinforced else more, if one wants to calculate the incident and reflected waves in the internal region. Even with the known exact solution for the wave function there is uncertainty in determination of these waves! But, namely, the standard definition of the coefficients of penetrability and reflection is based on them. In particular, we have not found papers where the coefficient of reflection is defined and estimated in this problem (which differs essentially from unity at the energy of radiation close to the height of the barrier and, therefore, such a characteristics could be interesting from a physical point of view). Note that the semiclassical approximation put serious limits to the possibility of its definition at all [21].
Thus, in order to estimate probability of the formation of the Universe as accurately as possible, we need a fully quantum definition of the wave. Note that the non-semiclassical penetrability of the barrier in the cosmological problems has not been studied in detail and, therefore, a development of fully quantum methods for its estimation is a perspective task.
Researches in this direction exist [2], and in these papers was estimated the penetrability on the basis of tunneling of wave packet through the barrier. However, a stationary boundary condition has uncertainty that could lead to different results in calculations of the penetrability. The stationary approach could allow to clarify this issue. It is able to give stable solutions for the wave function (and results in Ref. [29] have confirmed this at zero energy of radiation), using the standard definition of the coefficients of the penetrability and reflection, is more accurate to their estimation.
Aims of this Chapter are: (1) to define the wave in the quantum cosmological problem; (2) to construct the fully quantum (non-semiclassical) stationary method of determination of the coefficients of penetrability of the barriers and reflection from them on the basis of such a definition of the wave; (3) to estimate how much the semiclassical approach differs in the estimation of the penetrability from the fully quantum one. In order to achieve this goal, we need to construct tools for calculation of partial solutions of the wave function. In order to resolve the questions pointed out above, we shall restrict ourselves to a simple cosmological model, where the potential has a barrier and internal above-barrier region.

A. Dynamics of Universe in the Friedmann-Robertson-Walker metric
Let us consider a simple model of the homogeneous and isotropic Universe in Friedmann-Robertson-Walker (FRW) metric (see Ref. [48], p. 438; also see Refs. [9,23,40,41]): where t and r, θ, φ are time and space spherical coordinates, the signature of the metric is (−, +, +, +) as in Ref. [41] (see p. 4), a(t) is an unknown function of time and k is a constant, the value of which equals +1, 0 or −1, with appropriate choice of units for r. Further, we shall use the following system of units:h = c = 1. For k = −1, 0 the space is infinite (Universe of open type), and for k = +1 the space is finite (Universe of closed type). For k = 1 one can describe the space as a sphere with radius a(t) embedded in a 4-dimensional Euclidian space. The function a(t) is referred to as the "radius of the Universe" and is called the cosmic scale factor. This function contains information of the dynamics of the expansion of the Universe, and therefore its determination is an actual task. One can find the function a(t) using the Einstein equations taking into account the cosmological constant Λ in this metric (we use the signs according to the chosen signature, as in Ref. [41] p. 8; the Greek symbols µ and ν denote any of the four coordinates t, r, θ and φ): where R µν is the Ricci tensor, R is the scalar curvature, T µν is the energy-momentum tensor, and G is Newton's constant. From (1) we find the Ricci tensor R µν and the scalar curvature R: The energy-momentum tensor has the form (see [41], p. 8): T µν = (ρ + p) U µ U ν + p g µν , where ρ and p are energy density and pressure. Here, one needs to use the normalized vector of 4-velocity U t = 1, U r = U θ = U φ = 0. Substituting the previously calculated components (2) of the Ricci tensor R µν , the scalar curvature (4), the components of the energy-momentum tensor T µν and including the component ρ rad (a), describing the radiation in the initial stage (equation of state for radiation: p(a) = ρ rad (a)/3), into the Einstein's equation (2) at µ = ν = 0), we obtain the Friedmann equation with the cosmological constant (see p. 8 in Ref. [41]; p. 3 in Ref. [9]; p. 2 in Ref. [47]): whereȧ is derivative a at time coordinate. From here, we write a general expression for the energy density:

B. Action, lagrangian and quantization
We define the action as in Ref. ([47], see (1), p. 2): Substituting the scalar curvature (4), then integrating item atä by parts with respect to variable t, we obtain the lagrangian (see Ref. [47], (11), p. 4): Considering the variables a andȧ as generalized coordinate and velocity respectively, we find a generalized momentum conjugate to a: and then hamiltonian: The passage to the quantum description of the evolution of the Universe is obtained by the standard procedure of canonical quantization in the Dirac formalism for systems with constraints. In result, we obtain the Wheeler-De Witt (WDW) equation (see Ref. [47], (16)- (17), in p. 4, [14,39,49]), which can be written as where ϕ(a) is wave function of Universe. This equation looks similar to the one-dimensional stationary Schrödinger equation on a semiaxis (of the variable a) at energy E rad with potential V (a). It is convenient to use system of units where 8π G ≡ M −2 p = 1, and to rewrite V (a) in a generalized form as In particular, for the Universe of the closed type (k = 1) we obtain A = 36, B = 12 Λ (this potential coincides with Ref. [2]).
C. Potential close to the turning points: non-zero energy case In order to find the wave function we need to know the shape of the potential close to the turning points. Let us find the turning points a tp, in and a tp, out concerning the potential (12) at energy E rad : Let us expand the potential V (a) (13) in powers of q out = a − a tp (where the point a tp, in or a tp, out is used as a tp . Expansion is calculated at these points), where (for small q) we restrict ourselves to the liner term: where the coefficients V 0 and V 1 are: Now eq. (15) transforms into a new form at variable q with potential V (q):

III. TUNNELING BOUNDARY CONDITION IN COSMOLOGY
A. A problem of definition of "free" wave in cosmology and correction of the boundary condition Which boundary condition should be used to obtain a wave function that describes how the wave function leaves the barrier accurately? A little variation of the boundary condition leads to change of the fluxes concerning the barrier and, as result, it changes the coefficients of penetrability and reflection. So, a proper choice of the boundary condition is extremely important. However before, let us analyze how much the choice of the boundary condition is natural in the asymptotic region.
• In description of collisions and decay in nuclear and atomic physics potentials of interactions tend to zero asymptotically. So, in these calculations, the boundary conditions are imposed on the wave function at infinity. In cosmology we deal with another, different type of potential: its modulus increases with increasing of the scale factor a. The gradient of the potential also increases. Therefore, here there is nothing common to the free propagation of the wave in the asymptotic region. Thus, a direct passage of the application of the boundary condition in the asymptotic region into cosmological problems looks questionable.
• The results in Ref. [29], which show that when the scale factor a increases the region containing solutions for the wave function enlarges (and its two partial solutions), reinforce the seriousness of this problem. According to Ref. [29], the scale factor a in the external region is larger, the period of oscillations of each partial solution for the wave function is smaller. One has to decrease the time-step and as a consequence increase the calculation time. This increases errors in computer calculations of the wave function close the barrier (if it is previously fixed by the boundary condition in the asymptotic region). From here a natural conclusion follows on the impossibility to use practically the boundary condition at infinity for calculation of the wave (in supposition if we know it maximally accurately in the asymptotic region), if we like to pass from the semiclassical approach to the fully quantum one. Another case exists in problems of decay in nuclear and atomic physics where calculations of the wave in the asymptotic region are the most stable and accurate.
• One can add a fact that it has not been known yet whether the Universe expands at extremely large value of the scale factor a. Just the contrary, it would like to clarify this from a solution of the problem, however imposing a condition that the Universe expands in the initial stage.
On such a basis, we shall introduce the following definition of the boundary condition: The boundary condition should be imposed on the wave function at such value of the scale factor a, where the potential minimally acts on the wave, determined by this wave function.
The propagation of the wave defined in such a way is close to free one for the potential and at used value of the scale factor a (we call such a wave conditionally "free"). However, when we want to give a mathematical formulation of this definition we have to answer two questions: 1. What should the free wave represent in a field of a cosmological potential of arbitrary shape? How could it be defined in a correct way close to an arbitrary selected point?

2.
Where should the boundary condition be imposed?
To start with, let us consider the second question namely where we must impose the boundary condition on the wave function. One can suppose that this coordinate could be (1) a turning point (where the potential coincides with energy of radiation), or (2) a point where a gradient from the potential (having a sense of force of interaction) becomes zero, or (3) a point where the potential becomes zero. But the clear condition of free propagation of the wave is the minimal influence of the potential on this wave. So, we define these coordinate and force in the following way: The point in which we impose the boundary condition is the coordinate where the force acting on the wave is minimal. The force is defined as minus the gradient of the potential.
It turns out that according to such a (local) definition the force is minimal at the external turning point a tp, out . Also, the force, acting on the wave incident on the barrier in the internal region and on the wave reflected from it, has a minimum at the internal turning point a tp, in . Thus, we have just obtain the internal and external turning points where we should impose the boundary conditions in order to determine the waves. B. Boundary condition at a = 0: stationary approach versus non-stationary one A choice of the proper boundary condition imposed on the wave function is directly connected with the question: could the wave function be defined at a = 0, and which value should it be equal to at this point in such a case? The wave function is constructed on the basis of its two partial solutions which should be linearly independent. In particular, these two partial solutions can be real (not complex), without any decrease of accuracy in determination of the total wave function. For any desirable boundary condition imposed on the total wave function, such methods should work. In order to achieve the maximal linear independence between two partial solutions, we choose one solution to be increasing in the region of tunneling and another one to be decreasing in this tunneling region. For the increasing partial solution we use as starting point a x the internal turning point a tp, in at E rad = 0 or zero point a x = 0 at E rad = 0. For the second decreasing partial solution the starting point a x is chosen as the external turning point a tp, out . Such a choice of starting points turns out to give us higher accuracy in calculations of the total wave function than starting calculations of both partial solutions from zero or from only one turning point.
In order to obtain the total wave function, we need to connect two partial solutions using one boundary condition, which should be obtained from physical motivations. According to analysis in Introduction and previous section, it is natural not to define the wave function at zero (or at infinity), but to find outgoing wave at finite value of a in the external region, where this wave corresponds to observed Universe at present time. But, in practical calculations, we shall define such a wave at point where forces minimally act on it. This is an initial condition imposed on the outgoing wave in the external region 1 .
Let us analyze a question: which value has the wave function at a = 0? In the paper the following ideas are used: • the wave function should be continuous in the whole spatial region of its definition, • we have outgoing non-zero flux in the asymptotic region defined on the basis of the total wave function, • we consider the case when this flux is constant.
The non-zero outgoing flux defined at arbitrary point requires the wave function to be complex and non-zero. The condition of continuity of this flux in the whole region of definition of the wave function requires this wave function to be complex and non-zero in the entire region. If we include point a = 0 into the studied region, then we should obtain the non-zero and complex wave function also at such point. If we use the above ideas, then we cannot obtain zero wave function at a = 0. One can use notions of nuclear physics, field in which the study of such questions and their possible solutions have longer history then in quantum cosmology. As example, one can consider elastic scattering of particles on nucleus (where we have zero radial wave function at r = 0, and we have no divergences), and alpha decay of nucleus (where we cannot obtain zero wave function at r = 0). A possible divergence of the radial wave function at zero in nuclear decay problem could be explained by existence of source at a point which creates the outgoing flux in the asymptotic region (and is the source of this flux). Now the picture becomes clearer: any quantum decay could be connected with source at zero. This is why the vanishing of the total wave function at a = 0, after introduction of the wall at this point (like in Ref. [2]), is not obvious and is only one of the possibilities.
If we wanted to study physics at zero a = 0, we should come to two cases: • If we include the zero point into the region of consideration, we shall use to quantum mechanics with included sources. In such a case, the condition of constant flux is broken. But a more general integral formula of nonstationary dependence of the fluxes on probability can include possible sources and put them into frameworks of the standard quantum mechanics also (see eq. (19.5) in Ref. [21], p. 80). Perhaps, black hole could be another example of quantum mechanics with sources and sinks.
• We can consider only quantum mechanics with constant fluxes and without sources. Then we should eliminate the zero point a = 0 from the region of our consideration. In this way, the formalism proposed in this paper works and is able to calculate the penetrability and reflection coefficients without any lost of accuracy.
This could be a stationary picture of interdependence between the wave function at zero and the outgoing flux in the asymptotic region. In order to study the non-stationary case, then we need initial conditions which should define also the evolution of the Universe. In such a case, after defining the initial state (through set of parameters) it is possible to connect zero value of wave packet at a = 0 (i. e. without singularity at such a point) with non-zero outgoing flux in the asymptotic region. In such direction different proposals have been made in frameworks of semiclassical models in order to describe inflation, to introduce time or to analyze dynamics of studied quantum system (for example, see [17,42]).

A. Wave function of Universe: calculations and analysis
The wave function is known to oscillate above the barrier and increase (or decrease) under the barrier without any oscillations. So, in order to provide an effective linear independence between two partial solutions for the wave function, we look for a first partial solution increasing in the region of tunneling and a second one decreasing in this tunneling region. To start with, we define each partial solution and its derivative at a selected starting point, and then we calculate them in the region close enough to this point using the method of beginning of the solution presented in Subsection IV D 1. Here, for the partial solution which increases in the barrier region, as starting point we use the internal turning point a tp, in at non-zero energy E rad or equals to zero a = 0 at null energy E rad , and for the second partial solution, which decreases in the barrier region, we choose the starting point to be equal to the external turning point a tp, out . Then we calculate both partial solutions and their derivatives in the whole required range of a using the method of continuation of the solution presented in Subsection IV D 2, which is improvement of the Numerov method with constant step. In this way, we obtain two partial solutions for the wave function and their derivatives in the whole studied region.
In order to clarify how the proposed approach gives convergent (stable) solutions, we compare our results with the paper of [2]. Let us consider the behavior of the wave function. The first partial solution for the wave function and its derivative in my calculation are presented in Fig. 1, which increase in the tunneling region and have been obtained at different values of the energy of radiation E rad . From these figures one can see that the wave function 15   satisfies the rules satisfied by the wave function inside the sub-barrier and in above-barrier regions [51]. Starting from very small a, the wave function has oscillations and its maxima increase monotonously with increasing of a. This corresponds to the behavior of the wave function in the internal region before the barrier (this becomes more obvious after essential increasing of scale, see left panel in Fig. 2). Moreover, for larger values of a, the wave function increases monotonously without any oscillation, that points out on transition into the tunneling region (one can see this in a logarithmic presentation of the wave function, see central panel in Fig. 2). A boundary of such a transformation in behavior of the wave function must be the point of penetration of the wave into the barrier, i. e. the internal turning point a tp, in . Further, with increasing of a the oscillations are appeared in the wave function, which could be possible inside the above barrier region only (in the right panel of Fig. 2 one can see that such a transition is extremely smooth that characterizes the accuracy of the method positively). The boundary of such a transformation in the behavior of the wave function should be the external turning point a tp, out . Like Ref. [29], but at arbitrary non-zero energy E rad we obtain monotonous increasing of maximums of the derivative of the wave function and smooth decreasing of this wave function in the external region. One can see that the derivative is larger than the wave function. At large values of a we obtain the smooth continuous solutions up to a = 100 (in Ref. [2] the maximal presented limit is a = 30). In Fig. 3, it is presented the second partial solution of the wave function and its derivative at different values of the energy of radiation E rad According to the analysis, this solution close to the turning points, in the tunneling region, in the sub-barrier and above-barrier regions looks like the first partial solution, but with the difference that now the maxima of the wave function and their derivatives are larger essentially in the external region in a comparison with the internal region, and amplitudes in the tunneling region decrease monotonously.
Comparing the previous pictures of the wave function with the results of Ref. [2], one can see that the wave function, in this approach, is essentially more continuous, has no divergencies and its behavior is everywhere clear. From here we conclude that the developed method for the determination of the wave function and its derivative at arbitrary energy of radiation is essentially more quick, more stable and accurate in comparison with the non-stationary quantum approach  in Ref. [2]. Note that: • With increasing a, the period of the oscillations, both for the wave function and its derivative, decreases uniformly in the external region and increases uniformly in the internal region (this result was partially obtained earlier in Ref. [29] at E rad = 0).
• At larger distance from the barrier (i. e. for increasing values of a, in the external region, and at decreasing value of a, in the internal region) it becomes more difficult to get the convergent continuous solutions for the wave function and its derivative (this result was partially obtained earlier in Ref. [29] at E rad = 0).

• A number of oscillations of the wave function in the internal region increases with increasing of the energy of radiation E rad (this is a new result).
B. Definition of the wave minimally interacting with the potential Now we shall be looking for a form of the wave function in the external region, which describes accurately the wave, whose propagation is the closest to the "free" one in the external region at the turning point a tp, out and is directed outside. Let us return back to eq. (16) where the variable q = a − a tp, out has been introduced. Changing this variable to this equation is transformed into From quantum mechanics we know two linearly independent exact solutions for the function ϕ(ξ) in this equation -these are the Airy functions Ai (ξ) and Bi (ξ) (see Ref. Furthermore, we shall be interested in the solution ϕ(ξ) which describes the outgoing wave in the range of a close to the a tp point. However, it is not clear what the wave represents in general near the point a tp , and which linear combination of the Ai (ξ) and Bi (ξ) functions defines it in the most accurate way.
The clearest and most natural understanding of the outgoing wave is given by the semiclassical consideration of the tunneling process. However, at the given potential the semiclassical approach allows us to define the outgoing wave in the asymptotic region only (while we can join solutions in the proximity of a tp by the Airy functions). But it is not clear whether the wave, defined in the asymptotic region, remains outgoing near the a tp . During the whole path of propagation outside the barrier the wave interacts with the potential, and this must inevitably lead to a deformation of its shape (like to appearance of a phase shift in the scattering of a wave by a radial potential caused by interaction in scattering theory). Does the cosmological potentials deform the wave more than the potentials used for description of nuclear collisions in scattering theory? Moreover, for the given potential there is a problem in obtaining the convergence in the calculation of the partial solutions for the wave function in the asymptotic region. According to our calculations, a small change of the range of the definition of the wave in the asymptotic region leads to a significant increase of errors, which requires one to increase the accuracy of the calculations. Therefore, we shall be looking for a way of defining the outgoing wave not in the asymptotic region, but in the closest vicinity of the point of escape, a tp . In a search of solutions close to the point a tp , i. e. at small enough |ξ|, the validity of the semiclassical method breaks down as |ξ| approaches zero. Therefore, we shall not use the semiclassical approach in this paper.
Assuming the potential V (a) to have an arbitrary form, we define the wave at the point a tp in the following way.
Definition 1 (strict definition of the wave). The wave is a linear combination of two partial solutions of the wave function such that the change of the modulus ρ of this wave function is approximately constant under variation of a: According to this definition, the real and imaginary parts of the total wave function have the closest behaviors under the same variation of a, and the difference between possible maximums and minimums of the modulus of the total wave function is the smallest. For some types of potentials (in particular, for a rectangular barrier) it is more convenient to define the wave less strongly.

Definition 2 (weak definition of wave):
The wave is a linear combination of two partial solutions of wave function such that the modulus ρ changes minimally under variation of a: According to this definition, the change of the wave function caused by variation of a is characterized mainly by its phase (which can characterize the interaction between the wave and the potential). Subject to this requirement, we shall look for a solution in the following form: where where T is an unknown normalization factor, f (ξ) is an unknown continuous function satisfying f (ξ) → const at ξ → 0, and u max is the unknown upper limit of integration. In such a solution, the real part of the function f (ξ) gives a contribution to the phase of the integrand function, while the imaginary part of f (ξ) deforms its modulus. Let us find the first and second derivatives of the function Ψ(ξ) (a prime denotes a derivative with respect to ξ): From this we obtain: Considering the solutions at small enough values of |ξ|, we represent f (ξ) in the form of a power series: where f n are constant coefficients. The first and second derivatives of f (ξ) are Substituting these solutions into eq. (24), we obtain Considering this expression at small |ξ|, we use the following approximation: Then from eq. (18) we obtain the following condition for the unknown f n : Requiring that this condition is satisfied for different ξ and with different powers n, we obtain the following system: Assuming that the coefficients f 0 and f 1 are known, we find the following solutions for the unknown f 2 , f 3 and f n : where the following notations for the integrals have been introduced: Thus, we see that the solution (22) taking into account eq. (23) for the function ϕ (ξ) has arbitrariness in the choice of the unknown coefficients f 0 , f 1 and the upper limit of integration, u max . However, the solutions found, eqs. (32), define the function f (ξ) so as to ensure that the equality (22) is exactly satisfied in the region of a close to the escape point a tp . This proves that the function ϕ (ξ) in the form (22) For such a choice of the coefficients f 0 and f 1 , the integrand function in the solution (23) (up to ξ 2 ) has a constant modulus and a varying phase. Therefore, one can expect that the solution (22) at the turning point a tp describes the wave accurately.

C. Total wave function
Having obtained two linearly independent partial solutions ϕ 1 (a) and ϕ 2 (a), we can write the general solution (a prime is for the derivative with respect to a) as: where T is a normalization factor, C 1 and C 2 are complex constants found from the boundary condition introduced above: the ϕ (a) function should represent an outgoing wave at turning point a tp, out . Fig. 4 plots the total wave function calculated in this way for the potential (12) with parameters A = 36, B = 12 Λ at Λ = 0.01 at different values of the energy of radiation E rad . One can see that the number of oscillations of the wave function in the internal region increases with increasing of the energy of radiation. Another interesting property are the larger maxima of the wave function in the internal region at smaller distances to the barrier for arbitrary energy (result found for the first time). In Fig. 5 it has been shown how the modulus of this wave function changes at selected values of the energy of radiation. From these figures it becomes clear why the coefficient of penetrability of the barrier is extremely small (up to the energy E rad = 2000). In order to estimate, how effective is the boundary condition introduced above in building up the wave on the basis of the total wave function close to the external turning point a tp, out , it is useful to see how the modulus of this wave function changes close to this point. In Fig. 6 we plot the modulus of the found wave function close to the turning points at the energy of radiation E rad = 2000 is shown. Here, one can see that the modulus at a tp, out is practically constant (see left panel in Fig. 6). It is interesting to note that the modulus of the wave function, previously defined, does not change close to the internal turning point a tp, in , and is close to maximum (see right panel in Fig. 6). Here, we look for the regular partial solution of the wave function close to an arbitrary selected point a x . Let us write the wave function in the form: and rewrite the potential through the variableā: where Substituting the wave function (39), its second derivative and the potential (40) into Schrödinger equation, we obtain recurrent relations for unknown b n :  FIG. 6: The behavior of the modulus of the total wave function at the energy of radiation E rad = 2000, close to the turning points (for atp, in = 8.58, atp, out = 15.04, see also Table 1): (a) the modulus decreases monotonously in the tunneling region, with increasing of a. It shows maxima and holes connected with the oscillations of the wave function in the external region, but the modulus is not equal to zero (thispoints out the existence of a non-zero flux); (b) when a increases, the modulus reaches a minimum close to the external turning point atp, out (this demonstrates the practical fulfillment of the definition for the wave at such a point); (c) transition close to atp, in is shown, where at increasing of a the modulus with maximums and holes is transformed rapidly into a monotonously decreasing function without maximums and holes. This is connected with transition to the region of tunneling.
Given the values of b 0 and b 1 and using eqs. (42)-(44) one can calculate all b n needed. At limit E rad → 0 and at a x = 0 all found solutions for b i transform into the corresponding solutions (40), early obtained in [29] at E rad = 0. Using c 2 = 1, from eqs. (39) we find: So, on the basis of the coefficients b 0 and b 1 one can obtain the values of the wave function and its derivative at point a x . Imposing two different boundary conditions via b 0 and b 1 , we obtain two linear independent partial solutions ϕ 1 (a) and ϕ 2 (a) for the wave function. Using the internal turning point a tp, in as the starting point, we calculate the first partial solution which increases in the barrier region (we choose: b 0 = 0.1, b 1 = 1), and using the external turning point a tp, out as the starting point, we calculate the second partial solution which decreases in the barrier region (we choose: b 0 = 1, b 1 = −0.1). Such a choice provides effectively a linear independence between two partial solutions.

Method of continuation of the solution
Let us rewrite equation (18) in such a form 2 : Let a n be a set of equidistant points a n = a 0 + nh. Denoting the values of the wave function ϕ (a) at points a n as ϕ n , we have constructed an algorithm of the ninth order to determine ϕ n+1 and ϕ ′ n when ϕ n and ϕ n−1 are known: ϕ n+1 = ϕ n−1 g 11 + g 01 g 01 − g 11 + ϕ n g 01 g 10 − g 00 g 11 where A local error of these formulas at point a n equals to:

E. The penetrability and reflection in the fully quantum approach
Let us analyze whether a known wave function in the whole region of its definition allows us to determine uniquely the coefficients of penetrability and reflection.

Problem of interference between the incident and reflected waves
Rewriting the wave function ϕ total in the internal region through a summation of incident ϕ inc wave and reflected ϕ ref wave: we consider the total flux: where The j mixed component describes interference between the incident and reflected waves in the internal region (let us call it mixed component of the total flux or simply flux of mixing). From the constancy of the total flux j total we find the flux j tr for the wave transmitted through the barrier, and: Now one can see that the mixed flux introduces ambiguity in the determination of the penetrability and reflection for the same known wave function.

F. Determination of the penetrability, reflection and interference coefficients
In quantum mechanics the coefficients of penetrability and reflection are defined considering the potential as a whole, including asymptotic regions. However, in the radial calculation of quantum decay such a consideration depends on how the incident and reflected waves are defined inside finite internal region from the left of the barrier. The question is: does the location of such a region influence the penetrability and reflection? In order to obtain these coefficients, we shall include into definitions coordinates where the fluxes are defined (denote them as x left and x right ): So, the T and R coefficients determine the probability of transmission (or tunneling) and reflection of the wave relatively the region of the potential with arbitrary selected boundaries x left , x right . When x right tends to the asymptotic limit, the coefficient defined before should transform into standard ones. From eqs. (53) and (54) we obtain (j tr and j ref are directed in opposite directions, j inc and j tr -in the same directions): Now we see that the condition |T | + |R| = 1 has sense in quantum mechanics only if there is no interference between incident and reflected waves, and for this is enough that: A new question appears: does this condition allow to separate the total wave function into the incident and reflected components in a unique way? It turns out that the choice of the incident and reflected waves has essential influence on the barrier penetrability, and different forms of the incident ϕ inc and reflected ϕ ref waves can give zero flux j mix . Going from the rectangular internal well to the fully quantum treatment of the problem would become more complicated.

G. Wave incident on the barrier and wave reflected from it in the internal region
One can define the incident wave to be proportional to the function Ψ (+) and the reflected wave to be proportional to the function Ψ (−) : where I and R are new constants found from continuity condition of the total wave function ϕ total and its derivative at the internal turning point a tp, int : On the basis of these solutions we obtain at the internal turning point a tp, int the flux incident on the barrier, the flux reflected from it and the flux of mixing. The flux transmitted through the barrier was calculated at the external turning point a tp, ext .

H. Penetrability and reflection: fully quantum approach versus semiclassical one
Now we shall estimate through the method described above the coefficients of penetrability and reflection for the potential barrier with parameters A = 36, B = 12 Λ, Λ = 0.01 at different values of the energy of radiation E rad . We shall compare the coefficient of penetrability obtained with the values given by the semiclassical method. In the semiclassical approach we shall consider two definitions of this coefficient: where  From this table one can see that inside the entire range of energy, the fully quantum approach gives value for the coefficient of penetrability enough close to its value obtained by the semiclassical approach. This differs essentially from results in the non-stationary approach [2]. This difference could be explained by difference in a choice of the boundary condition, which is used in construction of the stationary solution of the wave function.

I. The penetrability in the FRW-model with the Chaplygin gas
In order to connect universe with dust and its accelerating stage, in Ref. [20] a new scenario with the Chaplygin gas was proposed. A quantum FRW-model with the Chaplygin gas has been constructed on the basis of equation of state instead of p (a) = ρ rad (a)/3 (where p (a) is pressure) by the following (see also Refs. [5,6]): where A is positive constant and 0 < α ≤ 1. In particular, for the standard Chaplygin gas we have α = 1. Solution of equation of state (62) gives the following dependence of density on the scale factor: where B is a new constant of integration. Using parameter α, this model describes transition between stage, when Universe is filled with dust-like matter, and its accelerating expanding stage (through scenario of Chaplygin gas applied to cosmology, for details, see Refs. [7,8,20], also historical paper [13]). Let us combine expression for density which includes previous forms of matter and the Chaplygin gas in addition. At limit α → 0 eq. (63) transforms into the ρ dust component plus the ρ Λ component. From such limit we find and obtain the following generalized density: Now we have:ȧ After quantization we obtain the Wheeler-De Witt equation where For the Universe of closed type (at k = 1) at 8π G ≡ M −2 p = 1 we have (see eqs. (6)-(7) in Ref. [7]): Let us expand the potential (69) close to arbitrary selected pointā by powers of q = a −ā and restrict ourselves to linear terms: For coefficients V 0 and V 1 we find: and eq. (67) has the form: After the change of variable eq. (72) becomes: After the new change From eqs. (73) and (75) we have: Using such corrections after inclusion of the density component of the Chaplygin gas, we have calculated the wave function and on its basis the coefficients of penetrability, reflection and mixing by the formalism presented above. Now following the method of Sec. III A, we have defined the incident and reflected waves relatively to a new boundary which is located in the minimum of the hole in the internal region. Results are presented in Tabl. 3. One can see that penetrability changes up to 100 times, in such a coordinate, in dependence on the location of the boundary or in the internal turning point (for the same barrier shape and energy E rad )! This confirms that the coordinate where incident and reflected waves are defined has essential influence on estimation of the coefficients of penetrability and reflection. This result shows that the method proposed in the present paper has physical sense. In the next Tabl. 4, we demonstrate the fulfillment of the property (55) inside the entire energy range, which is calculated on the basis of the coefficients of penetrability, reflection and mixing obtained before.

A. Passage to non-stationary WDW equation: motivations
Tunneling is a pure quantum phenomenon characterized by the fact that a particle crosses through a classicallyforbidden region of the barrier. By such a reason, the process of incidence of the particle on the barrier and its further tunneling and reflection are connected by unite cause-effect relation. So, the dynamical consideration of the tunneling process through cosmological barriers is a natural one. The rejection of the dynamical consideration of tunneling from quantum cosmology limits the possible connection between initial stage, when the wave is incident on the barrier, and next propagation of this wave. This leads to uncertainties in determination of penetrability and rates. According to quantum mechanics, a particle is a quantum object having properties both particle and wave. In the classically forbidden regions the wave properties of the studied object are evident. So, the wave description of tunneling is natural.
So, we define a non-stationary generalization of WDW equation as where τ is a new variable describing dynamics of evolution of the wave function being analog of time. According to quantum mechanics, the penetrability and reflection are stationary characteristics, and such characteristics, obtained in the following, are independent on the parameter τ . Note that all these characteristics are solutions of stationary WDW equation, while non-stationary consideration of multiple packets moving along barrier gives clear understanding of the process.
In order to give a basis to readers to estimate ability of the approach developed in this paper, let us consider results in [35] (see eq. (19)). Here was studied the non-stationary WDW equation with the potential for the closed FRW model with the included generalized Chaplygin gas.
After change of variable a new = a old √ 12 the non-stationary eq. (79) transforms into our eq. (78) simce the V eff potential is independent on the τ variable (such a choice allows a correspondence between energy levels, convenient in comparative analysis). The potential (79) after such a transformation is shown in figs. 8. We shall analyze the Now let us come to another more difficult problem, namely that a packet penetrating through the radial barrier of arbitrary shape in a cosmological problem. In order to apply the idea of multiple internal refections for study the packet tunneling through the real barrier, we have to generalize the formalism of the multiple internal reflections presented above. We shall assume that the total potential has successfully been approximated by finite number N of rectangular steps: where V i are constants (i = 1 . . . N ). Let us assume that the packet starts to propagate outside inside the region with some arbitrary number M (for simplicity, we denote its left boundary a M−1 as a start ) from the left of the barrier. We are interested in solutions for energies above that of the barrier while the solution for tunneling could be obtained after by change i ξ i → k i . A general solution of the wave function (up to its normalization) has the following form: at a min ≤ a ≤ a 1 (region 1), where α j and β j are unknown amplitudes, A T and A R are unknown amplitudes of transmission and reflection, k i = 1 h 2m(E − V i ) are complex wave numbers. We have fixed the normalization so that the modulus of the starting wave e ikM a equals to one. We look for a solution of such a problem by the approach of the multiple internal reflections.
Let us consider the initial stage when the packet starts to propagate to the right in the region with number M . According to the method of the multiple internal reflections, propagation of the packet through the barrier is considered by steps of its propagation relatively to each boundary (see [10,25,27], for details). Each next step in such a consideration of propagation of the packet will be similar to the first 2N − 1 steps. From analysis of these steps recurrent relations are found for calculation of all unknown amplitudes A (n) and β (n) j for arbitrary step n (for region with number j), summation of these amplitudes are calculated.
We shall look for the unknown amplitudes, requiring the wave function and its derivative to be continuous at each boundary. We shall consider the coefficients T ± 1 , T ± 2 . . . and R ± 1 , R ± 2 . . . as additional factors to amplitudes e ±i k a . Here, the bottom index denotes the number of the region, upper (top) signs "+" and "−" denote directions of the wave to the right or to the left, correspondingly. To begin with, we calculate T ± 1 , T ± 2 . . . T ± N −1 and R ± 1 , R ± 2 . . . R ± N −1 : Analyzing all possible "paths" of the propagations of all possible packets inside the barrier and internal well, we obtain: Choosing as starting points, the following: i and R ± i . Finally, we determine coefficients α j and β j : the amplitudes of transmission and reflection: and coefficients T and R describing penetration of the packet from the internal region outside and its reflection from the barrier Choosing a min = 0, we assume full propagation of the packet through such a boundary (with no possible reflection) and we have R − 0 = −1 (it could be interesting to analyze results with varying R − 0 ). We use the test: Now if energy of the packet is located below then height of one step with number m, then the following change should be used for description of transition of this packet through such a barrier with its tunneling. In the case of a barrier consisting from two rectangular steps of arbitrary heights and widths we have already obtained coincidence between amplitudes calculated by method of MIR and the corresponding amplitudes found by standard approach of quantum mechanics up to first 15 digits. Even increasing the number of steps up to some thousands has the right accuracy to fulfill the property (90).
In particular, we reconstruct completely the pictures of the probability and reflection presented in figs. 9 (a) and (b), figs. 10 (a) and (b), figs. 11 (b), but using such a standard technique. So, the result concerning the oscillating dependence of the penetrability on the position of the starting point a start in such figures is independent on the fully quantum method chosen for calculations.
This is an important test which confirms reliability of the method MIR. So, we have obtained full coincidence between all amplitudes, calculated by method MIR and by standard approach of quantum mechanics. This is why we generalize the method MIR for description of tunneling of the packet through potential, consisting from arbitrary number of rectangular barriers and wells of arbitrary sizes.

C. Results
We have applied the above method to analyze the behavior of the packet tunneling through the barrier (80) (we used a new → √ 12 a old ). The first interesting result is a visible change of the penetrability on the displacement of the starting point a min ≤ a ≤ a 1 , where we put the packet. Using the possibility of decreasing the width of intervals up to an enough small value (and choosing, for convenience, the width of each interval to be the same), we choose a min as starting point (and denote it as a start ), from where the packet begins to propagate outside. We have analyzed how the position of such a point influences the penetrability. In fig. 9 (a) one can see that the penetrability strongly changes in dependence of a start for arbitrary values of energy of radiation E rad : it has oscillating behavior. Difference between its minimums and maximums is minimal at a start in the center of the well (i. e. its change tends to zero in the center of the well), this difference increases with increasing value of a start and achieves the maximum close to the turning point. With this result, we may conclude that exists a dependence of penetrability on the starting point a start of the packet. The coefficients of reflection, oscillations and penetration on the position of the starting point a start are presented in next figs. 9 (b), (c), (d) and have similar behavior.
Usually, in cosmological quantum models the penetrability is determined by the barrier shape. In the non-stationary approach one can find papers where the role of the initial condition is analyzed in calculations of rates, penetrability etc. 3 But, the stationary limit does not give us any choice on which to work. We conclude: (a) the penetrability 20   should be connected with the initial condition (not only in non-stationary consideration, but also in the stationary one), which determines position (coordinate) of maximum of the packet which begins to propagate outside (at initial time moment t = 0). (b) Even in the stationary consideration, the penetrability of the barrier should be determined in dependence on the initial condition. The first question is how much these results are reliable. In particular, how stable will such results be if we shift the external boundary outside? The results of such calculations are presented in fig. 10, where it is shown how the penetrability changes with a max (for clearness sake, we have fixed the starting point: a start = 10). One can see that all calculations are well convergent, that confirms efficiency of the method of the multiple internal reflections. On the basis of such results we choose a max = 70 for further calculations. However, one can see that inclusion of the external region can change the coefficients of penetrability and penetration up to 2 times for the chosen energy level.
The second question is how strong this affects the calculations of the penetrability. If it was small than, the semiclassical approaches would have enough good approximation. From figs. 9 it follows that the penetrability is not strongly changed in dependence on shift of the starting point. However, such small variations are connected with relatively small height of the barrier and depth of the well, while they would be not small at another choice of parameters (the coefficient of oscillation and penetration turn out to change at some definite energies of radiation, see below). So, this effect is supposed to be larger at increasing height of the barrier and depth of the well, and also for near-barrier energies (i. e. for energies comparable with the barrier height) and above-barrier energies of radiation).
We have analyzed how these characteristics change in dependence on the energy of radiation. We did not expect the results that we got (see figs. 11). The coefficient of penetration has oscillations with peaks clearly shown. These peaks are separated by similar distances and could be considered as resonances in energy scale. So, by using the fully quantum approach we observed for the first time clear pictures of resonances which could be connected with some early unknown quasi-stationary states. At increasing energy of radiation the penetrability changes monotonously and determines a general tendency of change of the coefficient of penetration, while the coefficient of oscillations introduces the peaks. Now the reason of the presence of resonances has become clearer: oscillations of the packet inside the internal well produce them, while the possibility of the packet to penetrate through the barrier (described by the penetrability of the barrier) has no influence on them. In general, we observe 44   In the last fig. 12 one can see that we have achieved |T bar + R bar − 1| < 10 −15 inside whole region of changes of a start and a max (such data were used in the previous figs. 9 and 10). This is the accuracy of the method of the multiple internal reflections in obtaining T bar and R bar . D. The fully quantum penetrability versus semiclassical one in cosmology: a quick comparison Does the penetrability, determined according to the semiclassical theory by a shape of the barrier between two turning points, give exhaustive answers and the best estimations of rates of evolution of universe? If we look at figs. 9 (a), we shall see that this is not the case. The penetrability is depended on the position (coordinate) of maximum of the packet which begins to propagate outside at time moment t = 0. So, the penetrability should be a function of some parameters of the packet at its start. For the first time, it has been demonstrated the difference between the fully quantum approach and the semiclassical However, let us perform a general analysis.
(1) If we wanted to check the semiclassical approach, we should miss some of the parameters. One can use test of T + R = 1 (where T and R are the penetrability through the barrier and reflection from it). But, note that the semiclassical approximation neglects the reflected waves in quantum mechanics (see [21], eq. (46.10), p. 205, p. 221-222). Therefore, we cannot use the test above for checking T in the semiclassical theory.
(2) If we would like to determine the reflection coefficient, then we should find a more accurate semiclassical approximation (in order to take into account both decreasing and increasing components of the wave function in the tunneling region). In such a case, we shall face another problem, namely the presence of a non-zero interference between the incident and reflected waves. Now the relation T + R = 1 cannot be used as test, and one needs to take the third component M of interference into account (see [31]   wave function in the incident and reflected waves 4 , the interference component should increase without limit. In such a case, the penetrability and reflection can freely exceed unit and increase without limit. What is now the general meaning of the penetrability? (3) We shall give only some examples from quantum mechanics. (i) If we consider two-dimensional penetration of the packet through the simplest rectangle barrier (with finite size), we shall see that the penetrability is directly dependent on direction of tunneling of the packet. So, the penetrability is not a single value but a function. (ii) If we consider one-dimensional tunneling of the packet through the simplest rectangular barrier, we shall obtain "interference picture" of its amplitude in the transmitted region, which is dependent on time and space coordinates and is an exact analytical solution. Of course, the stationary part of such a result exactly coincides with well known stationary solutions [27].
(4) A tunneling boundary condition [46,47] seems to be natural and clear, where the wave function should represent an outgoing wave at large scale factor a. However, is such a wave free? In contrast to problems of quantum atomic and nuclear physics, in cosmology we deal with potentials, which modules increase with increasing the scale factor a (their gradients increase, which have sense of force acting on the wave). Therefore, in quantum cosmology we should define the boundary condition on the basis of the waves propagating inside strong fields (see [31]).
These points destroy the semiclassical basis of the cosmological models. Now the statement concerning reliability of the semiclassical approach become a question of " faith" (note that this is widespread [31,32]). The semiclassical approach could be compared with "black box", where deeper and more detailed information about the dynamics of the universe is hidden. of oscillation. The formula found, seems to be the fully quantum analogue of the semiclassical formula of Γ width of decay in quasistationary state proposed in Ref. [18]. Here, the coefficient of oscillations is the fully quantum analogue for the semiclassical F factor of formation and the coefficient of penetration is analogue for the semiclassical Γ width. 5. The penetrability of the barrier visibly changes in dependence of the position of the starting point R start inside the internal well, where the packet begins to propagate (see figs. 9). We note the following peculiarities: the penetrability has oscillating behavior, difference between its minimums and maximums is minimal at R start in the center of the well, with increasing R start this difference increases achieving to maximum near the turning point. The coefficients of reflection, oscillations and penetration have similar behavior. We achieve coincidence (up to the first 15 digits) between the amplitudes of the wave function obtained by such a method, and the corresponding amplitudes obtained by the standard approach of quantum mechanics (see Appendix B in [33] where solutions for amplitudes were calculated in general quantum decay problem). This confirms that this result does not depend on a choice of the fully quantum method applied for calculations. Such a peculiarity is shown in the fully quantum considerations and it is hidden after imposing the semiclassical restrictions.
6. In the non-stationary and stationary considerations the penetrability of the barrier should be connected with the initial condition. We suggest that a possible introduction of the initial condition into the known stationary semiclassical models can change the obtained results.
7. If one takes into account the external tail of the barrier, the penetrability is visibly changed. For example, the penetrability is changed up to 2 times (see figs. 10) for the barrier (8)

9.
A dependence of the penetrability on the starting point has maxima and minima. This allows to predict some definite initial values of the scale factor, when the universe begins to expand. Such initial data is direct result of quantization of the cosmological model. 10. The modulus of the wave function in the internal and external regions has minima and maxima which were clearly established in [29,31]. This indicates, in terms of values of the scale factor, where the probable "appearance" of the universe is the maximal or minimal. So, the radius of the universe during its expansion changes not continuously, but consequently passes through definite discrete values connected with these maxima. It follows that space-time of universe on the first stage after quantization seems to be rather discrete than continuous. According to results [29,31], difference between maxima and minima is slowly smoothed with increasing of the scale factor a. In this way, we obtain the continuous structure of the space-time at latter times. The discontinuity of space-time is direct result of quantization of cosmological model. This new phenomenon is the most strongly shown on the first stage of expansion and disappears after imposition of the semiclassical approximations.