Hunting for Gravitational Quantum Spikes

: We present the result of our examination of quantum structures called quantum spikes. The classical spikes that are known in gravitational systems, occur in the evolution of the inhomogeneous spacetimes. A different kind of spikes, which we name strange spikes, can be seen in the dynamics of the homogeneous sector of the Belinski–Khalatnikov–Lifshitz scenario. They can be made visible if the so-called inhomogeneous initial data are used. The question to be explored is whether the strange spikes may survive quantization. The answer is in the afﬁrmative. However, this is rather a subtle effect that needs further examination using sophisticated analytical and numerical tools. The spikes seem to be of fundamental importance, both at classical and quantum levels, as they may serve as seeds of real structures in the universe.


Introduction
The dynamics underlying the Belinskii-Khalatnikov-Lifshitz (BKL) scenario, which concerns a generic gravitational singularity [1,2], can be described by the nonlinear coupled system of ODEs for the three effective directional scale factors (see Part I of [3]). These dynamics have been recently quantized [4,5]. The quantum BKL scenario predicts that a gravitational singularity can be avoided by a quantum bounce, occurring in the unitary evolution of a given gravitational system.
A different approach to solve the problem of a singularity in the BKL scenario has been proposed by Ashtekar et al. [6]. It can be used, after the successful quantization, to tackle a generic gravitational singularity as well. Furthermore, even if one restricts to its homogeneous sector, the model can be explored from the perspective of another interesting issue, which is the emergence of gravitational structures known as spikes. The aim of this paper is to uncover such structures at the classical level and to investigate if they can survive the quantization. The spikes that we name here the "strange spikes" are different from the (transient or permanent) spikes observed in the dynamics of inhomogeneous spacetimes (see [7][8][9][10][11][12][13][14][15] and references therein). The latter have well-understood properties, whereas our spikes have been discovered and preliminarily examined in the context of quantum physics only recently, in [16]. Results of the latter paper suggest that quantum (strange) spikes do not exist. However, the issue of time has not been treated satisfactorily due to the fact that [16] deals mainly with the vacuum case. In the present paper, we couple the system to a massless scalar field, so that it can be used as a reference clock at both the classical and quantum levels. Moreover, [16] has included only simplified analyses of quantum observables of a spike. Our paper fills this gap as well.
Let us stress that we do not address the issue of possible resolution of a generic gravitational singularity, which is predicted by the BKL conjecture and follows from the To connect with the notation that is more common for affine algebras, we perform the partial redefinition of variables (C I , P J ) =: (C I , −2D J ), which leads to the affine Poisson brackets [16] {D I , D J } = 0 = {C I , C J } , {C J , D I } = δ I J C I .
An algebra with such brackets is called an affine Lie algebra. The dynamics of the system is specified by the equations (We have the extra factor 2 in Equation (3) and (4) that is missing in the corresponding equations of [16].) where D = D 1 + D 2 + D 3 and C = C 1 + C 2 + C 3 . There is no summation ∑ I in the rhs of (3) and (4). The solutions to (3)- (6) have to satisfy the Hamiltonian constraint where κ = ±1 defines two possible dynamics (two different signatures of the corresponding bilinear forms) with respect to the field φ. Unlike the traditional momentum, which serves to translate the canonical coordinate C I , the variable D I serves to dilate C I . The set of equations (3)- (7) incorporates the dynamics of all Bianchi-type A models. It presents a coupled system of nonlinear equations that have not been solved in the general case analytically yet. To get some insight into the local geometry of the space of solutions to these equations, we apply the dynamical systems method [20,21].
It is easy to see that space S of the critical points of the dynamics, defined by the vanishing of the right-hand sides of (3)- (6) and satisfying the constraint (7), reads S = {(C 1 , C 2 , C 3 , D 1 , D 2 , D 3 , π, φ) ∈ R 8 | (C I = 0 = π) ∧ (D 2 = 2 ∑ where I = 1, 2, 3 and D = D 1 + D 2 + D 3 . The Jacobian of the system (3)-(6) is easily found to be where Thus, the characteristic polynomial associated with J S reads so that the eigenvalues are the following: Since the real parts of all, but one, eigenvalues of the Jacobian J S are equal to zero, the fixed points defined by Equation (8) [20,21].). Thus, getting insight into the structure of the space of solutions to the dynamics near such points require an examination of the exact form of the dynamics. The information obtained from linearized set of equations is unable to reveal the nature of dynamics in the neighborhood of such fixed points.
In the next subsection, we present explicit, but approximate, solution to our dynamics characterizing the strange spike. It is obtained by solving the dynamics with the socalled inhomogeneous initial data (This definition and example explaining the idea of the inhomogeneous initial data are due to David Sloan.). The latter means that the initial data is not just a set of 3 Cs and 3 Ds per point in phase space, but the related data on some curve in this space. For instance, let us choose (C 1 , C 2 , C 3 ) := (x, 0.8, 0.4) and (D 1 , D 2 , D 3 ) := ( f (x), 2, 7), where f (x) is the value that solves the Hamiltonian constraint (7) for C 1 =x. Next, we allowx to vary from −0.1 to 0.1. Then, instead of solving one set of equations for each point of space, we solve a whole continuous (in practical calculations, discrete) family of them by taking the sequence ofx ∈ (−0.1, 0.1). The plot of the Cs and the Ds as functions ofx reveals a peculiar structure that emerges inx aroundx = 0 that we call the strange spike (We usex to denote the initial data in phase space sticking to the notation of Section 2 of Ref. [16].) It results from Equation (5) and (6) that φ is a monotonic function of time. Thus, it can be used as an evolution parameter of the dynamics. Dividing both sides of (3) and (4) bẏ φ = κπ, we obtain which defines the relative dynamics with respect to the variable φ.
In our case the equations of motion (13)- (14) and the constraint (7) are constructed from elements of the affine algebra (2). This algebra can be realized in terms of the Poisson algebra by the adjoint action where X and Y are linear combinations of basic elements C J and D I of affine algebra. Exponentiation of this algebra gives the operators representing the affine group. For fixed I (where I = 1, 2, 3) the elements of the affine group in "one direction" are represented by the following operations in the phase space where α I ∈ R and β I ∈ R + . The full group is composed of all three independent transformations g(α 1 , β 1 , α 2 , β 2 , α 3 , β 3 ) = g 1 (α 1 , β 1 )g 2 (α 2 , β 2 )g 3 (α 3 , β 3 ). This property allows, in most cases, to perform calculations for fixed I and to generalize the result to other I.
The action of this group on the basic elements of the affine algebra can be summarized as The manifold of the variables (C I , D I ) splits into orbits with respect to the affine group. The orbit of the element (C I0 , D I0 ) is defined as Using (17), it is easy to check that one gets two large orbits which we denote by Π I − and Π I + and continuum of orbits consisting of single points Π I D I In fact, due to the constraint (7) the sets of points in the orbits are reduced to some submanifolds.
Since C I = 0 is a critical point of the system (3)-(4), the sign of each C I along any dynamical trajectory is fixed by the initial conditions. The orbits Π I − and Π I + carry such solutions. Every single-point orbit Π I D I separates positive (C I > 0) and negative (C I < 0) parts of the kinematical trajectory. However, if the system enters the orbit Π I D I , it is not able to leave it. This property is very strong. It is even fulfilled not only for infinitesimal perturbation of the motion, but also for finite difference form of the equations of motion obtained by the Euler algorithm. Changing in (13)- (14) derivatives into finite differences, one gets where φ n+1 = φ n + ∆φ.
Assuming that for a given φ n the point (C I (φ n ) = 0, D I (φ n )) belongs to the orbit Π I D I , the resulting point (C I (φ n+1 ) = 0, D I (φ n+1 ) = D I (φ n )) also belongs to the same orbit Π I D I independently of ∆φ.
These single-point orbits separate classical trajectories into C I > 0 and C I < 0 regions. In fact, the region to which belong given trajectory depends on sign of C I (φ 0 ) of the initial conditions (C I (φ 0 ) = 0, D I (φ 0 )). This is shown in the next section.
Space S (see (8)) consists of the nonhyperbolic critical points so that the neighborhood of each such point includes rapidly changing trajectories. Thus, the trajectories approaching asymptotically the orbits Π I D I are very sensitive to the choice of the initial data for the dynamics. This can be seen in Figures 1 and 2 (for the case I = 1). These neighborhoods represent the strange spikes.
It is expected that quantization may smear out the regions around the single-point orbits so that the corresponding spikes may become more smooth.

Parametrization of Dynamics by a Scalar Field
In order to derive the spike solutions, within the dynamics parameterized by the scalar field φ, we follow the approach presented in Section 2 of Ref. [16].
Let us assume that the initial conditions for D I and C I at φ = φ 0 have the form: D 1 < D 2 < D 3 < 0 and 1 C I > 0. Then, it follows from (13)-(14) that C 2 and C 3 almost instantly vanish, while D 2 and D 3 turn out to be essentially constant. For later convenience we define D ± := D 2 + D 3 ± 2 √ D 2 D 3 . Therefore, the problem reduces to finding the evolution of C 1 and D 1 , which is governed by the equations (here we denote by prime the derivative with respect to φ): where the last one results from the constraint (7). Inserting the right-hand side of (26) into (25), we obtain an equation independent of C 1 , whose solution can be written as In the above expression the initial condition D 1 (φ 0 ) = D 10 has been replaced by due to the relation (26) for C 1 (φ 0 ) = C 10 . Furthermore, (27) and (26) give where "sgn" denotes the sign function (its value for C 10 = 0 is irrelevant since then C 1 (φ) = 0). One can verify that (29) together with (27) solve the equations (24) and (26).
Choosing the simple parametrization C 10 =x, we can now draw D 1 and C 1 as functions of both the evolution parameter φ and the initial conditionx, or as functions of onlyx, for different fixed values of φ. Figures 1 and 2 present the corresponding plots for the setting of other quantities: D 2 = −2, D 3 = −1, κ = 1, π = 1 and φ 0 = 0. One can see that D 1 (φ) and C 1 (φ) behave in the same way as P 1 (t) and C 1 (t) presented in Ref. [16], which is expected as the evolution parameter φ is a monotonic function of the evolution parameter t owing to Equations (5) and (6).

Parametrization of Dynamics by the Arc Length
The arc length of the curve r(x) ≡ (C 1 (x), D 1 (x)) is given by the integral wherex 0 is a certain chosen minimal value ofx. Calculating (30) we obtain where F denotes the elliptic integral of the first kind and E of the second kind. This allows us to express the curve r(x) as a function of s, which needs to be calculated numerically.
Introducing the (normalized) Frenet vectorŝ one can define the generalized curvature of r(s) as follows (see, e.g., [22]) In Figure 3 we depict the generalized curvature of the curve (C 1 (s), D 1 (s)) as a function of the normalized arclengths corresponding tox ∈ [−5, 5] (i.e., s divided by the maximal value s(x = 5), for a given φ) for different values of the evolution parameter φ. The values of κ, π, φ 0 and D 2 , D 3 are kept the same as in the previous subsection. Moreover, dots on the horizontal axis denote the value ofs(x = 0) for a given φ, which naturally coincides with the middle of the spike. The double peak corresponds to the two inflection points of the curve visible on the right plot in  Figure 3 shows that the spike is created at some moment in the evolution of the gravitational system and seems to be permanent. The shape of the spike depends on time and changes from a plateau to a singular structure.

Representation of the Affine Group
The quantum version of the Lie algebra (2) is defined by the algebraic quantization principle: C I →Ĉ I and D I →D I , such that (Throughout the paper we chooseh = 1 and use Planck's units except where otherwise stated.) where I, J = 1, 2, 3. The commutation relations (34) are the same as for the generators of the affine group [23]. The affine group Aff(R + ) I generated by the pairĈ I andD I has two inequivalent unitary representations U − (p, q) I and U + (p, q) I . They are constructed in two carrier spaces of square-integrable functions L 2 (R − , dν(x I )) and L 2 (R + , dν(x I )), dν(x I ) = dx I /|x I |, which correspond to the negative and positive spectrum of the position operatorĈ I , respectively. Because of physical interpretation we needs the full spectrum of the position operator. This requirement enforces using the reducible representation of the affine group in the carrier space K I := L 2 (R − , dν(x I )) ⊕ L 2 (R + , dν(x I )). The general form of the vector f ∈ K I can be written as a direct sum of the functions f ∓ ∈ L 2 (R ∓ , dν(x I ): The scalar product of such two vectors is the sum of the appropriate partial scalar products: The action of the affine group Aff(R + ) I in this carrier space K I can be written as where p ∈ R, q ∈ R + and This structure allows for extension of this affine action to the whole straight line. For this purpose it is enough to extend the appropriate functions from half-line to the full straight line: f − (x I ) = 0 for x I ≥ 0 and f + (x I ) = 0 for x I ≤ 0. Then, denoting by |x I ⊕ x I the "position" vector in the space K I , every function belonging to K I can be represented as: It is obvious that the space K I ⊂ L 2 (R, dν(x I )) and that the scalar product (36) can be rewritten as The action of the affine group Aff(R + ) I in this new carrier space, which we denote again by K I , can be written as The explicit representation of the generators of this group are given by the following operatorsD where I = 1, 2, 3. The corresponding unitary operators representing elements of the affine group are: Taking into account three variables x I (I = 1, 2, 3), the carrier space K for the representation of the algebra (34) can be defined to be and the scalar product is constructed according to prescription for tensor product of Hilbert spaces: The "total" affine group used in this paper is the direct product of the three affine groups Aff 0 = Aff(R + ) 1 3 . This realization of the affine group allows for the physical interpretation of quantized C I and D I variables.

Quantum Dynamics
The quantum dynamics of our system may be derived, to some extent, from the quantum version of the Hamiltonian constraint defined by Equation (7). In a standard approach, one maps the dynamical constraint into an operator defined in kinematical Hilbert space. Its kernel may be used to construct physical Hilbert space. However, such an approach leads to the problem of time at the quantum level.
The reason for having the scalar field in the Hamiltonian (7), is the hope that it may resolve the problem of time both at classical and quantum levels. Such an approach works in the classical case as it leads to the relative dynamics, defined by Equations (13)-(14), parameterized by the scalar field φ. However, an extension of this strategy to the quantum level faces serious difficulty. Namely, quantization of the scalar field algebra {φ, π} = 1 as followsπ so that [φ,π] = iI, leads to the inconsistency. In this case, according to the standard approach to quantum mechanics, φ represents an additional degree of freedom of our quantum system. The field φ is the variable involved in every required wave function and it cannot be considered as a parameter representing a reference clock we want to introduce. One needs to notice that every quantum amplitudes is independent of the variable φ because they are obtained by calculating the appropriate scalar product containing among others integration over φ.
To parameterize with φ the reference classical clock uncoupled to our quantum system one needs to construct a hybrid approximation of deterministic unitary quantum evolution: the field should evolve in a classical way and quantum states of the system should evolve according to a unitary prescription.
Let us treat the field φ as a classical field which value is considered as a parameter showing tics of a classical clock, that is φ is a parameter enumerating changes of our Hamiltonian system.
We propose to modify Schrödinger type unitary evolution operator to the form containing both: evolution of the classical field and evolution of the quantum system itself. This operator we denote by U (φ, φ 0 ). It is defined by a series of natural conditions: • First of all, the operator U (φ, φ 0 ) evolves the quantum state of our gravitational system from the "time" φ 0 and the state Ψ 1 to the "time" φ and the state Ψ 2 as follows where Ψ 1 Ψ 2 ∈ K, and where K is a Hilbert space. Thus U changes the state vector in the Hilbert state space and the time parameter by changing the field. Since the time parameter φ does not couple to the gravitational field in (7), we can factorize the evolution operator U (φ, φ 0 ) into two independent operations: (a) the unitary operator V K (φ, φ 0 ) acting on the spatial dependence of state vectors in the Hilbert space K while the field φ is changing, (b) the operation V π (φ, φ 0 ) acting on the parametric dependence of the state vectors of the field φ.
In what follows, we assume the dependence of the evolution operator on the difference τ = φ − φ 0 between the final and initial value of the field φ, i.e., we assume the translational invariance of the evolution operator with respect to the parameter φ. This means that U does not depend on the choice of the initial time φ 0 , but only on τ.
Thus, the full evolution operator can be written as • The evolution operator fulfils the standard conditions for quantum evolution: The first one represents the fact, that if there is no shift in time, the state vector stays the same. The second means that every evolution can be split into intermediate steps. These two conditions are expected to hold for both the classical end quantum evolution. The last line represents the unitarity condition which is related to the probabilistic interpretation of quantum mechanics.
To fulfil the last condition the parametric part of the evolution operator has to transform as the complex conjugation: Let us now consider a formal shift operation with respect to the field φ. For this purpose, we define a kind of adjoint action of the field φ and its canonically conjugate momentum π on the classical phase space. For an arbitrary function g(φ, π) on this phase space the adjoint action is defined to be where the Poisson bracket is given by One can directly check that where The powers of the adjoint action are understood as The isomorphic realization of the classical shift operation (56) with respect to the field φ in the state space is given by e τ{π,·} → e τ ∂ ∂φ .
As a consequence the formal shift of the state Ψ(φ, x) in respect to the time φ is given by where x := (x 1 , x 2 , x 3 ). The comparison of the operations (60) and (56) suggests that the shift generator ∂ ∂φ =:π, defined in the quantum state space, may play a similar role to the classical momentum π acting (by the adjoint action) in the phase space. Working in the quantum state space we postulate the replacement of the classical momentum π with the operationπ. The classical evolution in the phase space can be written in terms of the adjoint action as e τ{E(π),·} , where E(π) is a real function of the momentum π generating evolution of this free field. As a consequence of (59) the appropriate realization of this operation to be a part of the evolution operator is the replacement E(π) → −iE(π). The imaginary unit has to be added to fulfil the unitarity requirement (53). The classical part of the evolution operator is expected to be: where E(π) is some real function ofπ. Making use of (60)-(61) and the factorization (49), we rewrite (48) as follows Taking derivative of (62) with respect to τ, at τ = 0, leads to the local evolution equation:π IntroducingŴ we can rewrite (63) in the form Assuming enables rewriting (65) in the separable form which leads to the two eigenequations: andŴ We assume, according to (49), that the quantum evolution operator corresponding to the classical constraint consists of the quantized position-dependent part of (7) and the shifted parametric part of (7). This way we avoid quantization of the algebra {φ, π} = 1, and consequently quantization of the classical time variable φ. Both classical and quantum evolutions are now parameterized by a single variable φ that we call the time.
Since there are no products of C I and D I in (7), and due to (34), the mapping of H defined by (7) into a Hamiltonian operatorĤ is straightforward. We get which implies that 3.2.1. Solving the Eigenequation (68) Analytically Making use of (72), we get (68) in the form The solution to (73), for κλ = 1/2, is found to be whereas for κλ = 1/2 one has where A λ and B λ are arbitrary constants. In what follows we denote the solutions (74)-(75) as ω λ (A λ , B λ ; φ).

Solving the Eigenequation (69) by Variational Method
The eigenequation (69) can be solved numerically in terms of given finite basis of functions {ψ n } N n=0 , by taking the solution ψ λ in the form where c n are unknown coefficients to be determined. The functions ψ n should be consistent with the boundary conditions. It means, they should vanish sufficiently fast at zero and infinity to satisfy the condition The coefficients c n can be found by considering the following functional: It is clear that (78) vanishes identically if ψ N λ is an exact solution to the Equation (69). If this is not the case but R[ψ λ ] 1, then we have an approximate solution. The smaller R[ψ λ ], the better the approximation. The latter fact suggests a method of finding the numerical solution. Namely, one can minimize (78) with respect to all unknown coefficients, including the eigenvalue λ. This fixes all the parameters in Equation (76) and determines the error R[ψ λ ].
To start the procedure one should fix the basis {ψ n }. It is reasonable to incorporate the fact that the operatorŴ is invariant under S 3 group of permutations of the variables {x 1 , x 2 , x 3 }. Therefore, when looking for the basis, it is reasonable to consider functions sharing this symmetry, i.e., requiring they are symmetric with respect to the replacements x i ↔ x j . A convenient choice is provided by the following ansatz where α ≥ 1 2 , γ > 0,γ ∈ R, while ∑ n 1 +n 2 +n 3 ≤N stands for the sum over n 1 , n 2 , n 3 ∈ [0, N] such that n 1 + n 2 + n 3 ≤ N, i.e., the series (79) is terminated at the N-th order (N = n 1 + n 2 + n 3 ). The bracket (n 1 , n 2 , n 3 ) denotes ordering operation, e.g., c (023) = c 023 , c (203) = c 023 , etc. The operation guaranties that the function ψ N λ consists of symmetric terms with respect to the replacement x i ↔ x j . For instance, there are two second-order (N = 2) terms in (79): 1 2 c 011 (ln x 2 1 ln x 2 2 + ln x 2 1 ln x 2 3 + ln x 2 2 ln x 2 3 ) and 1 2 c 002 ( ln x 2 The additional weights 1/(n 1 + n 2 + n 3 )! are introduced for technical simplicity. They guarantee that the coefficients c n 1 n 2 n 3 , fixed by the minimization procedure, are of a similar order. The latter improve the minimization. Note that the number of c n 1 n 2 n 3 grows fast with order N.
The coefficients c n 1 n 2 n 3 are listed in Appendix A. As in the case of ψ λ 1 the point-like precision (81) gives E[ψ λ 2 ] 0.0058.

Solving the Eigenequation (69) by Spectral Method
We start with the ansatz where γ > 0, α ≥ 1/2. The eigenequation (69) can be rewritten aŝ With the ansatz (84), the problem of solving the eigenequation (69) reduces to the problem of solving the corresponding eigenequation forF operator, i.e., We solve Equation (87) by using the spectral methods [24]. This form is much more convenient from the numerical point of view because of the lack of two terms |x 1 x 2 x 3 | α and exp(−γ/2(|x 1 | + |x 2 | + |x 3 |)), causing additional numerical errors (More precisely, within the spectral method, these terms result in the combination of small and large numbers (components of the matrix representing an approximate form of the eigenequation at a lattice).).
It is convenient, for numerical treatment, to assume where N > 1 is the cut-off, while f n 1 n 2 n 3 stand for a fixed basis of functions. The standard procedure involves cosine function, however, it will be convenient to adopt a different choice. The solution to the eigenequation (87) is specified by fixing the unknown coefficients c n 1 ,n 2 ,n 3 . They can be determined by demanding the eigenequation to be satisfied at a lattice composed of fixed points. To illustrate this, let us restrict for simplicity to the one-dimensional case, rewriting the ansatz (88) as The eigenequation readsF Let {x n } N n=1 stands for the lattice. For instance, one could consider Tchebychev's nodes. These are defined as roots of the Tchebychev polynomial of the first kind of the degree n [24]. On the finite interval [−1, 1], they read This can be extended on [a, b] defining Here, it is important that the number of points should match the number of coefficients c n . At the lattice Equation (90) can be rewritten as where c = (c n ) = {c 1 , ..., c N } is a vector built out of unknown coefficients, while ( f nm ) and (F nm ) stand for N × N matrices defined as Solving Equation (93) one guaranties the combination (89) satisfies the eigenequation (90) at x = x n , n = 1, ..., N. Equation (93) is a generalized eigenequation: having specified matrices ( f nm ) and (F αγ nm ) one obtains unknown coefficients c n and the eigenvalue λ solving algebraic eigenequation (93). The coefficients specify approximate solution of the differential equation. The denser the grid {x n }, the better the precision. Now, we consider three-dimensional case. The first element of the construction is the functional basis f n 1 n 2 n 3 (x 1 , x 2 , x 3 ). It turns out, a convenient choice is the basis We can now search for solutions to the eigenequation (87). Finding the numerical solution f n 1 n 2 n 3 of Equation (87), one finds the solution to the eigenequation (69), given by Equation (84). In order to do so, consider a three-dimensional grid {x n } = {x As we are interested in covering both positive and negative x i ∈ R in (95), it is reasonable to allow negative n i (we exclude n i = 0 because of the form of the right hand side of Equation (95)). The sum (89) becomes ∑ −1 n=−N c n f n (x) + ∑ N n=1 c n f n (x) . Due to the presence of logarithmic function in Equation (95), this choice should respect the fact that terms in Equation (95) become highly oscillating in the limit x i → 0. Hence, it is reasonable to make the grid denser close to zero. However, this is not the case for the original Tchebychev's nodes (92). This can be achieved adopting the following, modified Tchebychev's nodes: where b ± stands for two real parameters, positive b + and negative b − . They provide respectively positive and negative nodes. Clearly, this holds for all three dimensions. Because of terms |x 1 x 2 x 3 | α exp(−γ/2(|x 1 | + |x 2 | + |x 3 |)), the function (84) vanishes close to zero, |x I | 1, and close to infinity, |x I | 1. Therefore, one gets a good approximation restricting to a relatively small finite number of nodes. This justifies the choice (96). Having specified the grid and functional basis, we are ready to find the solutions. Choosing (It turns out that takeing |b − | = |b + | significantly improves numerical precision.) and adopting the basis (95) with Here ψ s stands for a normalized function found by renormalization of the numerical solution ψ s → ψ s / ψ s 1/2 . We have obtained a fairly good precision despite considering a small grid. More precisely, taking N = 5 means we adopted ten points per dimension (the whole three-dimensional lattice is composed of 1000 points). The rationale for this is the following. First, the functions (95) provides a good basis in the sense that a combination involving a small number of terms results in a good approximation to the solution of the eigenequation (in the sense the numerical error turns out to be small). This is because of the presence of logarithmic function; something that has already been observed discussing variational method. Second, the eigenfunction vanishes fast for |x I | 1, and so we can restrict our analysis to covering a small, finite region In addition to the symmetric function (95), one can consider the antisymmetric one, adopting the basis Adopting the choices (97) Choosing, for instance, the first eigenvalue, one finds the antisymmetric solution ψ a and the corresponding error (81): The functions ψ s and ψ s are orthogonal and they both were constructed as normalized.

Imposition of the Dynamical Constraint
Equation (65) is the Schrödinger-like equation corresponding to the classical dynamics defined by Equation (13)- (14). However, the latter is constrained by the condition H = 0, with H given by (7). The Dirac quantization scheme applied in this paper consists in mapping the classical constraint to the quantum constraintĤ = 0, which according to Equation (70) reads:Ĥ Therefore, not all solutions to (65) are physical but only the ones satisfying (104). It turns out, however, that the solution to (104), in the form (66) with ω λ defined by (74)-(75), can only be the trivial one Ψ(x) = 0. To address this difficulty, we propose to impose, instead of (104), the weak form of the Dirac condition: which has to be satisfied by a given linear combination of the products of eigenfunctions where ω λ are defined by (74)-(75) (up to arbitrary constants A λ and B λ ), and ψ λ is determined numerically via (76) and (84). The symbol ∑ λ denotes summation or integration depending on the solutions to the eigenequations (68)-(69).
For λ < 0, one has Equation (108) leads to the condition where A λ B λ = 0. Assuming the orthonormality condition ψ λ |ψ λ = δ(λ , λ), we get Equations (110)-(112) lead to the condition Let us consider a special solution including only two eigenvalues λ = λ . In such a case Equations (111)-(113) give and The solution to (114)-(115) reads Therefore, one of the possible solutions to the constraint (105) has the form which is defined by the specification of any pair of λ = λ .

Quantum Spikes
The expectation values of our basic observables (42) in a state described by the wavefunction (106) (satisfying the weak Dirac condition (105)) read and The coefficients A λ , B λ occurring here can be fixed by imposing the initial conditions at a certain value of φ = φ 0 so that for the case I = 1 we have (Since the classical spike has been derived for the case I = 1, we stick to this case at the quantum level as well. However, the system is symmetric with respect to the choice of I so that the same is true for two other cases.): withx C 1 ,x D 1 ∈ R. In principle, φ 0 can be arbitrary but we choose φ 0 = 0 below, as we also did for classical variables in Figures 1-3.
In the next two subsections, we apply our numerical results for the wavefunction Ψ to calculate (118) and (119).

Using the Results of the Variational Method
Let us first consider the solutions (79)-(80) and (82)-(83). We calculate where (C 1 ) ij := ψ λ i |Ĉ 1 |ψ λ j , and we get The full wavefunction Ψ(φ, x) is given by Equation (117) with λ = λ 1 and λ = λ 2 . Due to the condition A λ B λ = 0 (cf. (111)), we may assume that e.g., A λ 2 = B λ 1 = 0. Then, A λ 1 and B λ 2 remain two independent complex parameters. Let us parameterize them as The absolute values |A λ 1 |, |B λ 2 | are fixed by Equations (116): The considered numerical solutions (80) and (83) correspond, respectively, to the values λ 1 −0.0821 and λ 2 −0.0957. Since we are interested in the oscillatory case, we choose κ = 1, so that (74) gives us complex functions The final form of the wavefunction reads Calculation of the expectation value (118) for the state (127) leads to the result with the following parameters The values of β, δβ and χ for the case of (123) and λ 1 , λ 2 mentioned below (125) are Meanwhile, the parameter ∆ϕ can be eliminated from Equation (128) by imposing on the latter the initial condition (120), which gives where we simplified the notation by takingx ≡x C 1 . Equation (133) allows to express ∆ϕ as a function ofx. Assuming that φ 0 = 0, we find two different solutions ∆ϕ (±) = ∆ϕ (±) (x): where atan2(.,.) stands for two-argument arctangent function. Substituting (134) into Equation (128), we ultimately obtain In order to detect quantum spikes, we now examine Ĉ 1 (±) for a fixed φ but as a function ofx. The dependence onx is trivial for φ = 0. For φ = 0, the domain is restricted to the intervalx ∈ [− β 2 + δβ 2 , β 2 + δβ 2 ] and the function Ĉ 1 (±) (x) has the non-trivial derivative: In particular, atx = 0 we have which vanishes if The second derivative of Equation (135) reads Equations (137)-(139) show that, depending on the values of χφ and δβ, Ĉ 1 (±) reach a local maximum or minimum atx = 0. In particular, taking (132) and k = 0 in (138), one finds that this happens for φ 0.72. We consider such a phenomenon to be the quantum analogue of a classical spike. In general, these quantum spikes occur only at specific moments of time φ, belonging to a periodic discrete set with the period ∆φ = π/χ 1.45, determined by Equation (138).
Analogously to (128), we obtain Here, D 1 = const because there are no off-diagonal components in the matrix (140). More precisely, the off-diagonal components are non-zero, but are small of order 10 −20 . This result is almost unaffected by change of the order N of the numerical approximation and off-diagonal components are actually becoming smaller with growing (For instance, for N = 8 one finds them to be equal 7.7 · 10 −20 , while for N = 10 one gets 2.2 · 10 −20 .) N. In conclusion, the numerical results indicate that D 1 does not evolve with time φ.
It is also worth stressing that, at least in the case of restriction to a superposition of two eigenstates (117), the presence of a spike-like structure, associated with the observablê C 1 is unaffected by the choice of a pair of numerical solutions, i.e., one symmetric and one antisymmetric wavefunction. Adopting different ones modifies the values of coefficients β, δβ, χ but one can still find the value of φ corresponding to the quantum spike. In the case ofĈ 1 the latter is given by Equation (138); this is a simple function of two eigenvalues λ 1 , λ 2 . In fact, the crucial requirement for the occurrence of spikes in the presence of non-zero components in the matrix (C 1 ).

Using the Results of the Spectral Method
One can now perform the analogous analysis for numerical solutions obtained via the spectral method. Taking ψ 1 = ψ s , λ 1 = λ s given by Equation (100) and ψ 2 = ψ a , λ 2 = λ a given by Equation (103), we obtain where c 0 = −1.28 · 10 −5 , while (D 1 ) ij = 0, ∀i, j. The matrices (C 1 ) and (D 1 ) are defined as before in Equations (122) and (140) but in contrast to the former case, the numerical solutions do not contain imaginary terms. For this reason, we have (D 1 ) = 0, as well as δc = 0. On the other hand, we can achieve the much better precision. Following the same steps as described in the previous subsection, one can again express the expectation value Ĉ 1 as a simple function of φ (cf. (128)): where ∆ϕ, β and χ are given by Equations (129)-(131). The numerical values of the latter two constants for the considered case of (100) and (103) are β 1.28 · 10 −5 , χ 6.05 .
Similarly to the Equation (133), we eliminate the parameter ∆ϕ by imposing the boundary condition (120), which now becomes Solving this equation for ∆ϕ, one finds two solutions (as in (134) before) Substitution of ∆ϕ = ∆ϕ (±) into Equation (144) finally gives The two obtained solutions are depicted in Figure 5. Quantum spikes (represented by solid lines) occur in both cases at time φ 0.26 and are periodic in φ, with the period ∆φ = π/χ 0.52. The existence of spikes is ensured by the presence of non-zero elements in the matrix (C 1 ). Starting with a different pair of symmetric and antisymmetric functions ψ s and ψ a , corresponding to different eigenvalues λ 1 and λ 2 , one also expects to find spikes (unless (C 1 ) ij = 0, ∀i, j).

Conclusions
In this paper, we have attempted to uncover the existence of quantum strange spikes (in short, quantum spikes), i.e., certain distinctive features in the quantum evolution of the considered gravitational system. Such features are expected to be analogs of steep structures, called by us strange spikes, that arise in the corresponding classical evolution. Figure 1 shows that the classical spikes presented in Figure 2 (and previously in the paper [16]) are not apparent effects of the projection of a 3D plot into a 2D plot, but real structures. Furthermore, it seems that an extremum (maximum, minimum) or an inflection point, occurring locally at the classical level, may turn into a similar structure at the quantum level. Quantization does not need to suppress the classical spikes as it was preliminarily concluded in [16].
Let us discuss the latter claim in more detail. On the basis of numerical results, we conjecture that a quantum spike is an extremum in the evolution of the expectation value of a given quantum observable. The inflection-type classical spike C 1 (x, φ) presented in Figures 1 and 2 becomes the extremum-type quantum spike presented in Figures 4 or 5. Another difference is that our classical spike is a monotonic function of time, whereas the corresponding candidate for a quantum spike is periodic. To be specific, we have restricted our analyses to just one quantum observableĈ 1 (the results forD 1 have insufficient accuracy).
Compared to the classical spikes, quantum spikes differ in two ways. Namely, they are rather mild and periodic in time. At the quantum level, classical structures become smoother and of a slightly different type, being specific only to discrete moments in time. Nevertheless, it can be conjectured that classical spikes survive quantization in this sense.
Constructing solutions of the quantum evolution, we have applied the variational and spectral methods, which are quite different. These methods are not only different conceptually but also use different bases of auxiliary functions. Still, the obtained quantum spikes are quite similar. In particular, they are periodic in time. The two completely dissimilar methods lead to similar structures, which may serve as the robustness test of our results.
For simplicity, we have identified quantum spikes by making use of only two classes (obtained by two different numerical methods) of solutions to the quantum dynamics. Many other classes of solutions are possible and hence other types of quantum spikes could exist. Some of them might look similar to classical spikes and be monotonic functions of time. However, increasing the number of solutions in a wave packet makes the construction much more technically involved. This is the reason we have restricted ourselves to a pair of numerical solutions, but it was sufficient to test the method.
Another issue that we need to stress is our application of a classical massless scalar field in the role of a clock at both the classical and quantum levels. The scalar field is not coupled to the gravitational degrees of freedom in the Hamiltonian (7) and hence such treatment is justified. In fact, this allows us to avoid the inconsistency that is commonly ignored in literature: the situation when time is a parameter in the classical theory but a quantized variable in the quantum theory. In both cases, it should be just the same evolution parameter. In this paper, we did not wish to consider the quantization of time.
The implementation of the dynamical constraint at the quantum level has been performed in the weak sense. Such a way of imposing constraints is practiced in other branches of quantum physics and quantum chemistry, especially in variational methods (see, e.g., [25][26][27] and references therein).
A study of the corresponding issues in the case of inhomogeneous spacetimes would be highly interesting since they naturally favor structure formation. Different to our strange spikes, the spikes found in such spacetimes (see [7][8][9][10][11][12][13][14][15] and references therein) have never been quantized. Let us also stress that we do not investigate the possible relation between the latter "inhomogeneous" spikes and our "homogeneous" ones, as it is beyond the scope of the present paper.
When we think about the real spikes, which might occur in the early observed universe, we rather think in terms of possible structures in spacetime. In contrast, the spikes that we study in this paper arise in the phase space of the Hamiltonian framework. The path from the dynamics in phase space to the dynamics in spacetime is complicated due to the Hamiltonian constraint. Apart from this, the truncation of the full system to the homogeneous sector considered in [6], which underlies our paper, introduces additional complexity. Thus, the interpretation of our spikes in terms of the spikes in spacetime is rather difficult. These difficulties are enhanced by the procedure of quantization. We postpone the examination of quantum spikes that are described directly in spacetime to our future work on the quantization of spikes known in the context of Gowdy space [7]. Our paper is about the possible existence of quantum spikes. The theoretical framework has been established but much more effort is necessary to prove that quantum spikes are a generic feature of the quantum gravitational systems. In particular, the preliminary results presented here could be extended by the further examination of the eigenequation problem (69). New classes of its solutions could lead to new types of quantum spikes. This extension of our research definitely requires making use of sophisticated analytical and numerical tools so that is far from being complete. More activity in this direction is needed.
Author Contributions: Authors contributed equally to this manuscript. All authors have read and agreed to the published version of the manuscript.
Funding: This research received no external funding.