{'fraction_non_alphanumeric': 0.03721107745151826, 'fraction_numerical': 0.011926627388307135, 'mean_word_length': 4.3277417715084505, 'pattern_counts': {'":': 0, '<': 6, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Tunneling in quantum cosmology: numerical study of particle creation

We consider a minisuperspace model for a closed universe with small and positive cosmological constant, filled with a massive scalar field conformally coupled to gravity. In the quantum version of this model, the universe may undergo a tunneling transition through an effective barrier between regions of small and large scale factor. We solve numerically the minisuperspace Wheeler--De Witt equation with tunneling boundary conditions for the wave function of the universe, and find that tunneling in quantum cosmology is quite different from that in quantum mechanics. Namely, the matter degree of freedom gets excited under the barrier, provided its interaction with the scale factor is not too weak, and makes a strong back reaction onto tunneling. In the semiclassical limit of small values of cosmological constant, the matter energy behind the barrier is close to the height of the barrier: the system ``climbs up'' the barrier, and then evolves classically from its top. These features are even more pronounced for inhomogeneous modes of matter field. Extrapolating to field theory we thus argue that high momentum particles are copiously created in the tunneling process. Nevertheless, we find empirical evidence for the semiclassical-type scaling with the cosmological constant of the wave function under and behind the barrier.


Introduction and summary.
A traditional formalism in quantum cosmology is based on the Wheeler-De Witt equation [1,2], often truncated to minisuperspace. Though this equation is not applicable at length scales near and below the Planck scale, it should correctly describe quantum gravity phenomena at larger distances. Among these phenomena, penetration through classically forbidden regions is of particular interest. Indeed, various proposals for "quantum creation of the Universe" [3,4,6,7,8] invoke this possible quantum gravity effect in one or another way (for discussion and references see, e.g., Refs. [9,10]).
A prototype example is an empty closed Universe with small and positive cosmological constant Λ. In the minisuperspace approximation, the only dynamical degree of freedom is the radius of the Universe a. In Planck units and with appropriate rescaling, the gravitational action has the following form (see, e.g., Refs. [9,10] and references therein), where η is the time variable (in N = 1 gauge, η has the meaning of conformal time), π a is momentum conjugate to a, N is the lapse function and is the (conformal) Hamiltonian for the scale factor. Correspondingly, the minisuperspace Wheeler-De Witt equation iŝ where Ψ is the wave function of the Universe, andâ andπ a are treated as operators with standard commutation relations 1 . In the coordinate representation one has Ψ = Ψ(a) andπ a = −i∂/∂a. One may modify eq. (1) slightly by adding a spatially homogeneous mode of a massless conformal scalar field. Then the Wheeler-De Witt equation becomesĤ where ǫ is the conformal energy of the scalar mode. There are various proposals for the wave function of the Universe [3,5,6,11] which, in the language of the Wheeler-De Witt equation, correspond to different boundary conditions for Ψ(a). In this paper we consider tunneling boundary conditions [3,7,8]. These resemble closely the boundary conditions corresponding to tunneling transitions in conventional quantum mechanics. Namely, eq. (2) has the form of the stationary Schrödinger equation at energy ǫ in a potential see Fig. 1. At small Λ (in Planck units) the system is semiclassical: it is straightforward to see that Λ plays the role ofh. At sufficiently large ǫ (again in Planck units), but ǫ < V max = 1/(16Λ), there are classically allowed regions on the two sides of the barrier. In conventional quantum mechanics, tunneling from small a to large a would be described by the wave function which contains only outgoing waves as a → +∞; if there are degrees of freedom other than a, then the incoming component of the wave function has to be specified on the left of the barrier. In this paper we consider the same boundary conditions in the context of quantum cosmology.
At small Λ, the decaying part of the wave function in the classically forbidden region is so the wave function behind the barrier is suppressed by the standard semiclassical exponential (assuming that Ψ ∼ 1 in the classically allowed region on the left of the barrier), where a 1 and a 2 are the two turning points. The analogy to the stationary Schrödinger equation is no longer perfect if matter degrees of freedom (or gravitons) are added. Consider, as an example, a massive scalar field conformally coupled to gravity, whose action is In the minisuperspace approximation, the scalar field is spatially homogeneous, Φ = φ(η) a(η) , and its action, with appropriate rescaling, is We note in passing that if one blindly considered an inhomogeneous mode of conformal momentum k, then the scalar Hamiltonian would be With the scalar field added, the Wheeler-De Witt equation becomes (Ĥ a +Ĥ φ )Ψ + ǫΨ = 0 .
Unlike in conventional quantum mechanics, this equation contains a negative "kinetic term" for the scale factor, (−π 2 a /2), and the potential V (a) enters with negative sign also. This feature, which at first sight appears technical (and is often ignored in literature, see, however, Ref. [12]), in fact reflects the difference between tunneling in quantum cosmology and in quantum mechanics. Consider, as an example, tunneling from the ground state of the φ-oscillator. In quantum mechanics, the excitation of a fluctuating degree of freedom (similar to the field φ) would lower the energy of the tunneling degree of freedom (the scale factor a), and thus suppress tunneling even further. Hence, the fluctuating degree of freedom does not get excited much and the semiclassical exponent (4) does not change (the fluctuating degree of freedom affects the pre-exponential factor only). In quantum cosmology, on the contrary, the excitation of the field φ makes the tunneling of the Universe exponentially easier [7]: crudely speaking, the effective value of ǫ gets increased. Hence, the φ-quanta tend to be copiously created. Because of their back reaction, the semiclassical expansion about the solution (3) may break down, and the entire picture of tunneling in quantum cosmology may get substantially modified when matter degrees of freedom or gravitons are added.
It is worth pointing out, that at given ǫ there exists only one real Euclidean solution in the theory with the action (S a + S φ ). This solution, the instanton (periodic instanton for ǫ = 0), has φ = 0 and coincides with the instanton in a theory truncated to a single variable a. The above arguments imply that this instanton may become completely irrelevant if matter degrees of freedom and/or gravitons are kept.
In this paper we substantiate our qualitative observations by a numerical study of quantum systems of the Wheeler-De Witt type. Namely, we consider systems with two degrees of freedom whose own Hamiltonians are of opposite signs. To obtain a clear picture, we modify eq. (7) in such a way that the potential V (a) is well localized in a, and the interaction between φ and a is also localized in a. Then the tunneling boundary conditions are well defined even for a finite interval of a (which one has to work with in a numerical study), while the peculiarities of the Wheeler-De Witt equation remain there. To verify that the main results of our study are not sensitive to the details of the a-φ interaction, we modify the interaction in different ways, and consider three different models in subsequent sections. Our particular systems are described in section 2.1.
If the interaction between φ and a is weak enough, the matter degree of freedom may be treated as perturbation (provided it is in its ground state on the left of the barrier, otherwise the parameter ǫ should be redefined). To the zeroth approximation, the wave function in the classically forbidden region is then given by eq. (3). In section 2.2 we present a perturbative treatment of our systems about the solution (3) and get a heuristic understanding of the range of parameters in which the perturbation theory breaks down. Unlike in conventional quantum mechanics, this range includes part of the would-be semiclassical region of small Λ. We will be mostly interested in this range of parameters.
Our main purpose is to study our systems at the quantum level. We have solved the Wheeler-De Witt equation for these systems numerically, without making use of any approximation (except for approximations inevitable in numerical studies, such as discretization, truncation of the matter Hamiltonian to finite number of levels, etc.; these approximations are safely under control). We describe our results in section 3. We begin with the numerical analysis of the classical evolution (section 3.1) to see that, in a wide range of parameters, the potential barrier cannot be overcome classically. We then proceed in section 3.2 to the quantum treatment. We find that unless the interaction between a and φ is weak, the matter degree of freedom does indeed get strongly excited under the barrier, and the system behind the barrier contains a large number of φ-quanta. Furthermore, in the region of the parameter space which in conventional quantum mechanics would be adequately described semiclassically, with the wave function (3), the energy of the matter degree of freedom is close, behind the barrier, to the height of the barrier V max ∼ 1/Λ. The system "climbs up" the potential barrier, and the classically allowed region effectively starts again just right of the point of maximum, a max . These features are even more pronounced for a larger a-independent term in the oscillator frequency, k 2 in eq. (6), again in contrast to quantum mechanics. Even if the interaction between matter and gravity is weak (i.e., the parameter m in eq. (6) is small), the φ-quanta are copiously created due to tunneling, provided that their frequency is high (k is large enough). Extrapolating to field theory in the tunneling Universe, we thus argue that high momentum modes get strongly excited in the tunneling process.
In section 3.3 we concentrate on scaling properties in Λ of our quantum solutions in the region of small Λ. We find, somewhat surprisingly, that the wave functions in fact scale in a semiclassical way, ln Ψ = F /Λ, where F depends on appropriate combinations of other parameters and Λ. This empirical evidence suggests that tunneling in quantum cosmology may still be described semiclassically, but the relevant classical solutions are entirely different from the instanton existing in this model. This conjecture is further supported by our study of the dependence on the a-φ coupling parameter (m in eq. (5)) at small but fixed Λ. As we already mentioned, at small m the matter degree of freedom is not excited due to tunneling, while at larger m it is, and its energy behind the barrier is close to V max . We find numerically that the transition between the two regimes is very sharp in m and becomes close to step function as Λ gets smaller. This suggests that the wave function is a combination of two semiclassical exponentials, where the first exponential is due to tunneling along φ = 0 and coincides with eq. (3), while the second exponential comes from an entirely different classical path. In spite of considerable effort, we were unable to identify the classical solution responsible for tunneling with excitation of the matter degree of freedom, which would be of considerable interest to quantum cosmology. Overall, we have found that tunneling in quantum cosmology is quite different from that in quantum mechanics. This observation may have interesting implications to the discussion of initial conditions for the classical stage of the cosmological evolution.
Appendix A contains details of our numerical procedure. We describe various checks of this procedure in Appendix B. These checks ensure that our results are reliable, even though we are dealing with very small numbers typical for tunneling. The calculations have been performed at the supercomputer SGI Origin 2000 of the Center for Computational Science of Boston University.

The model and perturbative solution.
It is convenient to rescale the variable a by introducing b = √ Λa and write the Wheeler-De Witt equation (7) in the following form, where the coordinate representation in b is chosen, while the representation of operatorsπ φ andφ is not yet fixed. HereĤ φ is the Hamiltonian of an oscillator with b-dependent frequency,Ĥ where and we introduced a parameter ω for generality 2 . The parametersǫ and M are related to original parameters as follows, It is worth noting that for arbitrary U(b) and f (b), eq. (9) corresponds to the action where N is the lapse function.
In the original Wheeler-De Witt equation, one has and f (b) = b 2 . In our numerical study, however, we found it convenient to choose U(b) and f (b) is such a way that they are well localized.

Modified models.
We note that, since we are imposing free boundary conditions at b = 0 (as U(0) = 0, f (0) = 0), the range of b can be extended to (−∞; +∞). It is then convenient to shift b in such a way that the maximum of the potential U(b) occurs at b = 0. To be able to impose the boundary conditions in a finite interval of b, we will work with the following potential and interaction functions, and We have chosen three sets of parameters α 1 , α 2 , α 3 and b 1 , b 2 , b 3 , summarized in Table 1. Since not all α i have the same sign, the functions f 1 , f 2 and f 3 all decay exponentially as |b| → ∞. We have chosen values of c such that the maximum values of these functions are close to 1. Then the strength of b-φ interaction is determined by the parameter M in eq. (11).
The potential U(b) and the three interaction functions are shown in Fig. 2. The shapes of the interaction functions are quite different: while f 1 is substantial in front of the barrier (and hence the matter degree of reedom φ may be excited there), f 2 and f 3 vanish in front of the barrier, with f 3 being non-zero in the classically forbidden region only. Thus, these three interaction functions constitute a representative set. We will see that the main features of tunneling do not depend on the shape of the interaction function. Table 1: Parameters of the interaction functions. For fixed choice of f (b), there are four parameters in the model, namely, Λ,ǫ, ω and M. Naively, one would expect either from eq. (9) or eq. (13) that as Λ → 0 with ǫ, ω and M fixed, tunneling through the barrier U(b) is described semiclassically, and, furthermore, that the matter degree of freedom makes a small effect (provided that tunneling occurs from a low-lying state of the φ-oscillator). Let us consider the corresponding perturbation theory to see whether this is indeed the case.

Perturbative treatment.
Let us consider the case Λ ≪ 1 withǫ, ω and M fixed, and treatĤ φ in eq. (9) as perturbation. We take the φ-oscillator at left of the barrier in its ground state. NeglectingĤ φ altogther, we obtain that in the classically forbidden region, except near the turning points, the dominant part of the wave function is where and b 1 is the left turning point. To include the effect ofĤ φ perturbatively, we write and obtain, to the first non-trivial order, that Ψ φ (b) obeys where we introduced a variable τ , instead of b, such that and τ = 0 corresponds to the left turning point. Equation (16) is a "wrong sign" Euclidean Schrödinger equation; the "wrong sign" (as compared to conventional quantum mechanics) originates from the "wrong sign" of the Hamiltonian for the scale factor [7].
To deal with eq. (16), let us decompose Ψ φ (τ ) into eigenstates of the "instantaneous" Hamiltonian (10), i.e., the eigenstates Ψ (Ω) n (τ ) of an oscillator of frequency We find that (16) translates into a set of equaions for the coefficients c n , This set is still hard to solve. To get an idea of its solutions, we make use of the adiabatic approximation, and treat the last two terms in eq. (17) perturbatively.
Since we consider an oscillator in its ground state on the left of the barrier, the zeroth order solution is To the first order we obtain that the only non-vanishing coefficient is The first factor in the integrand, accounts for the creation of a pair of φ-quanta at point b ′ in the forbidden region; it is small in the adiabatic regime, which occurs at large ω or small M. However, there is a growing exponential factor in the integrand of eq. (19), which reflects the property that the creation of φ-quanta enhances tunneling. The perturbation theory about the solution (15) breaks down when this factor becomes dominant; in that case one expects that matter has strong back reaction on the process of tunneling. Our perturbative treatment shows that at given ω and small M, tunneling in quantum cosmology is indeed described by the wave function (15) and effects of the matter degree of freedom are small. However, as M increases, the tunneling regime changes, and matter effects become important. A similar phenomenon occurs if M is kept fixed and small, while ω increases instead: high energy quanta are copiously created even if their interaction to gravity is weak 3 . This is in sharp contrast to conventional quantum mechanics.
The system with the action (13) may either tunnel through the potential barrier or overcome it classically, depending on the parameters. Before studying tunneling, it is instructive to discuss briefly the classical solutions and find the range of parameters in which classical transitions over the barrier do not occur. Let the classical system begin its evolution from large negative b towards the barrier. In terms of the variables b and ϕ = √ Λφ the classical equations of motion read Here we fixed the gauge by introducing the standard should also be satisfied. Note that the parameter Λ does not enter the classical equations; as we pointed out above, it is analogous toh and appears in the quantum problem only.
At large negative b the variables decouple. One can freely choose the value of b in this region at initial time. Then there are two independent initial data, which we choose as the φ-oscillator energy E ϕ (t = 0) ≡ E ϕ,0 and its initial phase θ 0 . The "energy" associated with the scale factor, which we define as 3 According to eq. (12) this may happen even if the physical masses and momenta are small.
is then determined at the initial moment of time from the constraint (21), For given initial b the latter relation determinesḃ(t = 0) ≡ −π b (t = 0), so the initial data are completely defined. Conversely, one may choose the initial value of π b , and hence E b,0 , arbitrarily, and findǫ from eq. (22). Note that sinceǫ is non-negative. For every pair of parameters (M, ω) there exists a region in the plane (E ϕ,0 ; E b,0 ) in which at some values of the phase θ 0 classical over-barrier transitions occur. The rest of this plane may be called "the no-transition region" of initial data: irrespectively of the oscillator phase, the barrier cannot be overcome classically for initial data belonging to this region. Let us note that at E ϕ,0 = 0 the φ-oscillator remains in its classical ground state ϕ = 0 during the entire evolution, so the initial data (E ϕ,0 = 0; E b,0 > U max /Λ ≡ 1/Λ) give rise to classically allowed over-barrier transitions, while the interval (E ϕ,0 = 0; 0 < E b,0 < 1/Λ) belongs to the no-transition region of initial data.
For E ϕ,0 = 0 classical over-barrier transitions may occur even if E b,0 is lower than the barrier height 1/Λ. Indeed the parametric amplification of ϕ-oscillations may take place and give rise to over-barrier transitions for E b,0 < 1/Λ. Hence, the boundary of the no-transition region of initial data is not a straight line; this boundary depends on parameters M and ω.
We have found the regions in the plane of initial data (E ϕ,0 ; E b,0 ) in which classical over-barrier transitions occur at some values of the initial oscillator phase θ 0 , and those where these transitions do not happen for any θ 0 , by solving numerically the system (20), (21). Details of the numerical procedure are given in Appendix A. We have found that even at relatively large interaction strengh M, classical transitions do not occur for E ϕ,0 ≪ 1/Λ and E b,0 < 1/Λ. An example is shown in Fig. 3. In what follows we always choose the parameters in such a way that the classically allowed over-barrier transitions are impossible.

Solutions to the Wheeler-De Witt equation.
Our main purpose is to study solutions to the Wheeler-de Witt equation (9) for our model systems numerically, without making use of any approximation. The tunneling boundary conditions are as follows. At large |b| the variables b and φ decouple, and the frequency of the φ-oscillator is equal to ω. Hence, the asymptotic wave function may be decomposed into the oscillator eigenstates, At b → ±∞ the asymptotic solutions are, respectively, These conditions are direct analogs of the tunneling boundary conditions in quantum mechanics.
Unless explicitly stated, we will consider tunneling from the ground state of the φ-oscillator, n 0 = 0 .
We will discuss the dependence on n 0 towards the end of this section. We describe our numerical procedure in Appendix A. Various checks of this procedure are summarized in Appendix B. Here we present the results of our calculations.
To describe our solutions at finite b, including classically forbidden region, it is convenient to work in the basis of the "instantaneous" Hamiltonian of the φoscillator. The eigenfunctions obeŷ The wave function is expanded We will be interested in the behavior of the occupation numbers |C n (b)| 2 as functions of the scale factor b and excitation number n. We will see that at small Λ, which is of primary interest to us, the occupation numbers are exponential, To discuss these plots, we recall first that the potential (14) is centered at b = 0 and its width equals √ 2, so the potential is small at b = ±2, the first and last points displayed in Fig. 4. In Fig. 5 we show also the occupation numbers at larger b where the interaction function is large and additional creation of φ-quanta occurs, now in the classicaly allowed region. We also note that we consider quite small values of Λ, so the occupation numbers are actually very small (because of tunneling suppression) in the forbidden region and behind the barrier. Nevertheless we are confident that the results are reliable, partially due to checks of our numerical procedure summarized in Appendix B. It is clear from Figs. 4 and 5 that on the left of the barrier, the φ-oscillator is mostly in its ground state (since we have chosen n 0 = 0), while in the classically forbidden region its wave function gradually becomes peaked at the value of n such that the total matter energy, (E (n) φ +ǫ/Λ), is close to the height of the barrier, V max = 1/Λ. Comparing Figs. 4 and 5, we see that this property is essentially independent of the shape of the interaction function, as well as on parameters of the model (provided that M is sufficiently large, see below). We have checked that this is indeed the case by studying the solutions for numerous sets of parameters and also for the interaction function f 3 .
The mechanism due to which tunneling is accompanied by the creation of φquanta is clear. The φ-b interaction results in weak creation of these quanta (either in front of the barrier, as is the case for the interaction function f 1 , or in the forbidden region, as is the case for f 2 and f 3 ). The low modes get strongly suppressed in the forbidden region, while the modes of energy E (n) φ ∼ 1/Λ survive. The latter effect occurs because the total energy of matter is close to the barrier height. Hence, there is strong back reaction of created φ-quanta on tunneling.
One remark is in order. Even though we impose the boundary condition (25) with n 0 = 0, the wave function on the left of the barrier still contains a small admixture of the excited states of the φ-oscillator. This is due to the presence of the reflected wave.
Another way of studying the properties of our solutions is to consider the average energy of the φ-oscillator. For finite b we define it as follows, In the asymptotic region b → +∞, a more appropriate definition would involve the flux factor, as in conventional quantum mechanics, Since |C n | are sharply peaked at a certain n, the two definitions do not differ much: . We will use the definition (27) for general b, and the definition (28) when discussing the asymptotics b → +∞. As one anticipates from Figs. 4 and 5, the total matter energy at b → +∞ should be close to V max = 1/Λ, i.e., This is shown in Fig. 6. In fact, the relation (29) should not be precise, because the relevant quantity is the oscillator energy near the top of the barrier rather than its asymptotic value. For given n the former is equal to (n + 1/2)Ω(b ≈ 0) = (n + 1/2) ω 2 + M 2 f (b ≈ 0) which is larger than the asymptotic value (n + 1/2)ω. Hence, it is more appropriate to study the oscillator energy (27) as function of b for different values of parameters. The corresponding plots are shown in Fig. 7, which clearly displays how the system climbs up the potential barrier, and then stays at constant |C n | 2 , with energy decreasing due to the effect just described. as the exponential factor compensates for small M at larger b. In fact, the properties we discuss should be absent for very small M, and the perturbative regime of section 2.2 should set in instead. We will come back to this point later on.
It is instructive to see what happens at different Λ. We still consider tunneling from the oscillator ground state, n 0 = 0. At large Λ the barrier is effectively absent, so the φ-oscillator remains in its ground state with energy At small Λ the energy of the φ-oscillator increases substantially, to The overall behavior of the asymptotic value of the oscillator energy is shown as function of Λ in Fig. 8 in log-log scale. One indeed observes the asymptotics (31) at small Λ and (30) at large Λ. The oscillatory behavior of E φ asymp as function of Λ, shown in detail in Fig. 9, may be understood as follows. For finite Λ, the number of "under-barrier" oscillator levels (which have energy below (1 −ǫ)/Λ) is finite, and as Λ decreases, some over-barrier levels become under-barrier. Since the system tends to climb on the top of the barrier, but not much above the top, the structure of levels near the top matters. This structure changes with Λ in an oscillatory way, hence the oscillatory behavior shown in Fig. 9.
A similar dependence on Λ is inherent in the magnitude of the wave function where we again consider tunneling from the ground state of the φ-oscillator. We find numerically that at small Λ the exponent scales as i.e., it is independent of Λ. This is shown in Fig. 10. Note that in the regime we consider, Γ/Λ is much larger than the naive tunneling exponent entering eq. (4). Thus, creation of φ-quanta has dramatic effect on the magnitude of the wave function behind the barrier. In the regime of relatively large interaction parameter M, the dominant part of the wave function behind the barrier corresponds to φ-oscillator excited to a particular energy, such that the total matter energy near the top of the barrier is close to V max = 1/Λ. Hence, this state is essentially independent of the state of the φ-oscillator on the left of the barrier, i.e., of n 0 , up to the overall penetration factor (the latter depends on n 0 strongly). This is illustrated in Fig. 11.
Let us now discuss the case of small M, when the couping between φ and b is weak. According to the perturbative analysis of section 2.2 in that case the φ-oscillator does not get excited due to tunneling, provided that ω is small. However, one may expect that at large enough ω the situation is different and that the system is back in the regime of strong particle creation. Hence, as ω increases, the asymptotic energy of the oscillator should change from a small value to   Fig. 12. Of course, the value of ω at which the transition between the two regimes occurs depends on M and other parameters, but in any case the scalar quanta with large enough ω are copiously created due to tunneling. In the context of quantum cosmology this means that the Universe after tunneling is filled with particles, possibly of large mass and high spatial momenta.

Scaling with Λ.
We have seen in section 3.2 that at small Λ and sufficiently large M and/or ω, the φ-oscillator gets excited in such a way that the total energy of matter is close to the height of the barrier. This means that the wave function of the oscillator is peaked near a certain excitation number n ∝ 1/Λ. One may consider this property as an indication that the tunneling transition is of semiclassical type. Further support of this conjecture comes from the semiclassical-type scaling of the overall magnitude of the wave function behind the barrier, eq. (33). Let us describe further empirical evidence in favor of the semiclassical nature of the tunneling transitions in the regime where the φ-oscillator gets strongly excited.
One observation is that not only the overall magnitude, but also the entire wave function has the semiclassical form both in the classically forbidden region and behind the barrier. As an example, we show in Figs. 13 and 14 real and imaginary parts of F ≡ Λ ln C n behind the barrier as functions of ΛE   Figure 14: The same as in Fig. 13, but for imaginary part.
More evidence comes from the study of the transition from the perturbative regime, occuring at small values of the coupling parameter M, to the regime of strong excitation of the φ-oscillator at larger M. As shown in Fig. 15, this transition becomes more sharp as Λ decreases, and becomes essentially step function at small Λ. We consider this as an indication that the wave function behind the barrier and in the forbidden region is a combination of two exponentials of semiclassical type, eq. (8).
What classical solutions describe the tunneling transitions with strong excitation of the φ-oscillator, remains an open problem.

A.1 Classical solutions.
To obtain the solutions to eqs. (20) and (21), we chose the initial value b = −15, well outside the barrier, and imposed the constraint (21) at the initial moment of time. We found the classical solution by making use of the fourth order Runge-Kutta method with adaptive step size (see, e.g., Refs. [13]) for the first order Hamiltonian equations of motion. We checked that the constraint (21) was satisfied during entire evolution with relative precision 10 −9 . For given initial oscillator phase θ 0 , the boundary of the no-transition region of initial data was found with precision 10 −6 by "division by 2" search in E b,0 . The phase θ 0 was discretized into 100 points, and the minimum in θ 0 of the boundary of the no-transition regions was found.

A.2 The discretized Wheeler-De Witt equation.
We solved the Wheeler-De Witt equation (9) by expanding the wave function in eigenstates of the oscillator of fixed frequency ω, where Ψ (ω) n are independent of b. 5 Asymptotically,C n coincide with the coefficients entering eq. (24), so the boundary conditions are precisely (25). The Wheeler-De Witt equation takes the form of an ordinary differential equation for the vector C, where the matrix A is defined as (n + 1)(n + 2)δ n ′ , n+2 + n(n − 1)δ n ′ , n−2 .
To solve eq. (36) numerically, we truncated the system to a finite number of levels, i.e., set 0 ≤ n ≤ (N 0 − 1). We considered the system in a finite interval of b and discretized the variable b by introducing a lattice of spacing ∆. The sites were chosen at points b i = i∆ , One way to discretize eq. (36) would be to make use of the symmetric approximation for the second derivative However, we made use of the Numerov-Cowling method which has higher precision. Namely, the Taylor expansion ofC(b) gives Making use of eq. (36), one writes the second term on the right hand side as where eq. (37) was used for the derivative of the vector (AC). In this way we obtain the discretized version of eq. (36) with error O(∆ 6 ), The unknowns of this finite difference equation are numbersC i,n with i = −N b , . . . N b , n = 0, . . . (N 0 − 1). So, the number of unknowns is N 0 (2N b + 1). The system (39) is a set of N 0 (2N b − 1) equations. The boundary conditions (25) after discretization become 2N 0 linear equations Hence, the total number of equations coincides with the number of unknowns, sõ C i,n are uniquely determined. As we will see below, the solution of the discretized system well approximates the solution of the original one only if the numbers of lattice points and oscillator levels are quite large, N b ∼ 20000, N 0 ∼ 150. Numerical methods for generic system of equations would fail in this case, since the total number of complex equations is more than 3 million. However, following Ref. [14] we exploited the special form of eq. (39) to derive the solution by an iterative procedure. Indeed, by inverting numerically a set of N 0 × N 0 matrices 6 , one writes eq. (39) in the following form, where L i and R i are again N 0 ×N 0 matrices. Equations (41) can now be used to eliminate a set ofC i variables with given i, leading to new equations, which, with some matrix algebra, can be cast again into the form ( It is worth pointing out that the procedure in which intermediate variables are excluded involves operations with real matrices only. Note also that the interaction between differentC n in eq. (36) occurs only if they have the same parity. Hence, it is sufficient to consider only even (or only odd) modes.
The algorithm described above lends itself to parallel implementation: the vari-ablesC i at odd sites are iliminated in parallel, and the procedure is repeated in a recursive manner. Our code scaled well with the number of processors N proc , and most of the results presented in this paper were obtained with N proc = 64.
One must of course be careful that the effects of discretization and truncation of the original Wheeler-De Witt equation are kept small. We used a value for the cut-off of the number of modes N 0 such that the truncation energy exceeds the height of the barrier by a factor of 2, The lattice spacing was chosen equal to 1/12 of the minumum wavelength along b, We used the lattice of the size so that the potential is very small near the end points. We performed various checks, summarized in Appendix B, which show that this choice of parameters is appropriate.
A.3 Iterative improvement of precision.
In our calculation we used a double precision representation of real numbers, which entails round-off errors of order 10 −15 to 10 −16 . The exponential behavior in b and n gives origin to numbersC i,n varying over a range that extends well beyond the limits of double precision. What makes it still possible for us to use double precision is that only coefficientsC i,n with neighboring values of i and n and comparable magnitude are coupled by the equations that must be solved. Thus one never encounters additions or subtractions of variables with wildly differing exponents. The elimination procedure used for solving the equations can nevertheless propagate round-off errors, so that special care is needed to ensure that the values of the smallerC i,n are accurate. We accomplished this by an iterative improvement of precision, which we briefly describe here. LetC i,n be the initial solution obtained by the elimination procedure. The propagation of round-off errors will give origin to absolute errors of order ǫ = 10 −13 to 10 −14 in the coefficientsC i,n . Thus those coefficients whose true value is smaller than ǫ will turn out to be totally wrong and the equations coupling them will have errors of the same order of magnitude. It is convenient to denote these equations symbolically byÂ C = 0. Substituting into the equations the numbers C (0) we obtain insteadÂ C (0) = B (0) , where B (0) ∼ ǫ. We introduce therefore a correction δC (0) , which we find by solvingÂ δC (0) = −B (0) . Note that the equations for δC (0) only involve right-hand side terms of order ǫ, so that the propagation of roundoff errors will now give origin to errors of much smaller absolute magnitude ∼ ǫ 2 . With δC (0) we form the first iterate of the solution C (1) = C (0) + δC (0) . While the numbers C (0) and δC (0) will typically not be smaller than ǫ, the first iterates C (1) corresponding to coefficients whose true value is less than ǫ will accordingly take smaller values, being now affected by errors ∼ ǫ 2 . Inserting the first iterate C (1) into the equations we find now errors B (1) =Â C (1) ∼ ǫ 2 and the procedure can be repeated until the values of the smallest coefficients become stable.
We show in Fig. 17 the relative error as function of n after different numbers of iterations. This quantity is defined as the ratio of the modulus of the left hand side of eq. (41) to the sum of moduli of all terms in this equation. The absolute value ofC i,n at b = 0 at different steps of the iteration procedure is presented in Fig. 16. We conclude that the iteration procedure indeed enables one to obtain the solutions with relative, not absolute, error of order 10 −13 .

A.4 Total computational time.
Our main emphasis in this paper is the study of the solutions to the Wheeler-De Witt equation at small Λ. Let us estimate the dependence of the time of calculations on this parameter. At each step of the iteration procedure of section A.3, N b · N 3 0 floating point operations are performed. The number of iterations is proportional to logarithm of the minimum value of the wave function. The latter scales as N iter ∝ 1/Λ, see eq. (35). Thus, the time of calculations is proportional to  This estimate shows that reaching lower values of Λ occurs with dramatic increase of processor time.
In practice, for Λ = 0.0225 and ω = 0.6, the time of calculations with 64 processors was equal to 30 to 40 minutes. At larger values of ω the time of calculations was somewhat smaller.

A.5 Bogoliubov transformation.
Numerical calculations were conveniently performed in the basis of the eigenfunctions of the oscillator of b-independent frequency ω. The interpretation of the results is done, however, in the basis of the "instantaneous" Hamiltonian, describing the oscillator of frequency Ω(b). So, for every b we performed the Bogolubov transformation and obtained the coefficients C n starting fromC n .
The key element of this procedure is a set of transformation coefficients 0 Ω |k ω k ω |α n Ω |m ω ; n, m < N 0 .
Here α Ω denotes the annihilation operator in the basis of the "instantaneous" Hamiltonian. The first factor on the right hand side is given explicitly by while the matrix elements k ω |α n Ω |m ω can be calculated numerically, using the harmonic oscillator algebra.
B Checks of numerical procedure.
Since we were dealing with very small values of the coefficientsC i,n , we had to perform a series of checks of the entire procedure. To make sure that round-off errors do not spoil our calculations, we repeated the calculations with quadruple precision (32 decimal places), and found that the results for the wave function coincided with those obtained with double precision, with relative error 10 −13 . This check was performed for the interaction function f 3 and the following values of parameters: √ Λ = 0.3, 0.21 and 0.17; M = 0.2; ω = 0.6;ǫ = 0.27. We also checked that the effects of discretization of b and truncation of the oscillator energy levels are small, by varying the parameters N b , N 0 and ∆ by a factor of 2. The interaction function in this check was f 3 , the parameters were √ Λ = 0.21; M = 0.2; ω = 0.6;ǫ = 0.27. The relative difference between the solutions was always less than 10 −3 .
Another way to check our procedure is to make use of the current which is conserved in the continuum theory as a consequence of the Wheeler-De Witt equation, i.e., ∂ b J = 0 , while the discretized version does not have this property. By calculating the value of the current for every lattice site, we checked the discretization error. The precision of the calculation of J for given set of vectorsC i must be O(∆ 6 ), otherwise the conservation of current would hold with lower precision than the precision the solution was found. From the Taylor expansion we find Here we made use of eq. (36) to express higher derivatives ofC i . The first and second derivatives of the matrix A were found by straightforward differentiation of the expression (37). Making use of the calculated current J, we estimated the relative error of the calculation ofC i , where the exact value of the current was found in the asyptotic region of large positive b, Because of the exponential suppression of the wave function behind the barrier, eq. (44) is satisfied with better precision, as compared to the precision at which the current is calculated in the classically forbidden region. Equations (42), (43) and (44) give an estimate of the relative error in the calculation ofC i,n , which we found to be less than 10 −7 .
As one more check we changed the sign of the Hamiltonian of the scale factor, and obtained a conventional quantum mechanical system of two degrees of freedom. We checked for that system, that the sum of reflection and transmission coefficients was equal to 1 with precision 10 −10 .
Finally, we compared the numerical solutions at small M and Λ to the analytic solutions obtained within the adiabatic approximation of section 2.2. The first two modes of these two solutions coincided within 3 %. The latter, fairly large deviation we ascribe to the approximations involved in obtaining the analytic solutions.