Numerical techniques for solving the quantum constraint equation of generic lattice-refined models in loop quantum cosmology

To avoid instabilities in the continuum semiclassical limit of loop quantum cosmology models, refinement of the underlying lattice is necessary. This lattice refinement leads to new dynamical difference equations which, in general, do not have a uniform step size, implying complications in their analysis and solutions. We propose a numerical method based on Taylor expansions, which can give us the necessary information to calculate the wave function at any given lattice point. The method we proposed can be applied in any lattice-refined model, while in addition the accuracy of the method can be estimated. Moreover, we confirm numerically the stability criterion which was earlier found following a von Neumann analysis. Finally, the “motion” of the wave function due to the underlying discreteness of the space-time is investigated, for both a constant lattice, as well as lattice refinement models.

Figure 1

This figure demonstrates the problem of evolving a difference equation that has been defined on a varying lattice. Given the value of $\Psi$ defined on two adjacent lattice points (solid circles) one can use the difference equation to evaluate $\Psi$ on the subsequent lattice site (open circle). However, when one attempts to evaluate $\Psi$ at the next lattice point, one finds that one does not have the necessary data point (square).

Figure 2

The evolution of a wave function using the Taylor expansion scheme. The circles are the approximated values calculated on the varying lattice, while the line is the exact wave function. We use ${\mu }_0=1.0$, $A=-0.5$, ${\mathcal{H}}_{\varphi }=1.0\times {10}^{-4}$, and ${\mu }_{\mathrm{initial}}=1.0$, and initial conditions ${\Psi }_1=0.0$ and ${\Psi }_2=0.5$.

Figure 3

The difference between the exact approximate wave functions is less than 1%. This error is almost entirely due to the slight difference in the initial conditions between the two grids, i.e., the difference between the initial separation of ${\Psi }_1$ and ${\Psi }_2$. We use ${\mu }_0=1.0$, $A=-0.5$, ${\mathcal{H}}_{\varphi }=1.0\times {10}^{-4}$, and ${\mu }_{\mathrm{initial}}=1.0$, and initial conditions ${\Psi }_1=0.0$ and ${\Psi }_2=0.5$.

Figure 4

The value of the next order term in the Taylor expansion. We use ${\mu }_0=1.0$, $A=-0.5$, ${\mathcal{H}}_{\varphi }=1.0\times {10}^{-4}$, and ${\mu }_{\mathrm{initial}}=1.0$, and initial conditions ${\Psi }_1=0.0$ and ${\Psi }_2=0.5$.

Figure 5

(a) For the fixed lattice case the two-dimensional wave function can be calculated, given suitable initial conditions (solid circles). (b) In the case of a refining lattice, the data needed to calculate the value of the wave function at a particular lattice site (open square) are not given by previous iterations (solid squares). This is due to the fact that ${\delta }_{\mu }({\mu }_i,{\tau }_i)-{\delta }_{\mu }({\mu }_{i+1},{\tau }_i)\ne 0$, where ${\mu }_{i+1}={\mu }_i+4{\delta }_{\mu }({\mu }_i,{\tau }_i)$. To improve the clarity of the diagram, the explicit $\mu$ and $\tau$ dependence of the step sizes is suppressed, and the notation ${\delta }_{\overline{\mu }}\equiv {\delta }_{\mu }({\mu }_{i+1},{\tau }_i)$ and ${\delta }_{\overline{\tau }}\equiv {\delta }_{\tau }({\mu }_{i+1},{\tau }_i)$ is used.

Figure 6

Given the value of a function evaluated at three positions, $(x_1,y_1)$, $(x_2,y_2)$, and $(x_3,y_3)$, a Taylor expansion can be used to approximate the value at a fourth point, $(x_4,y_4)$.

Figure 7

The initial data (solid circles) are given on two adjacent $\mu$ rows and the lattice points along the left-hand (small $\tau$) diagonal. This allows us to evaluate the wave function on lattice points (open circles) within a triangular region of base ${\mu }_{\max }-{\mu }_{\mathrm{initial}}$ and height ${\tau }_{\max }-{\tau }_{\mathrm{initial}}$.

Figure 8

The wave function is calculated by iterating the difference equation using first-order Taylor expansion to evaluate the necessary data points at each site. The initial conditions are a Gaussian on ${\tau }_{\mathrm{initial}}=5000$, centered between ${\mu }_{\mathrm{initial}}=100$ and ${\mu }_{\max }=735$, of width $\sigma =25$. The lattice refinement model used is ${\delta }_{\mu }(\mu ,\tau )={\mu }^{-1/2}$, ${\delta }_{\tau }(\mu ,\tau )={\tau }^{-1/2}$. Only the relevant section of the wave function is plotted and to improve clarity, only every 20th $\tau$ lattice point is shown. In (a) the propagation of the wave function is shown, while in (b) the value of the second-order terms is plotted. The second-order corrections are typically of the order of ${10}^{-2}\%$ over this range.

Figure 9

A perturbation of $\Psi ={10}^{-6}$ was put on an otherwise empty initial $\mu$ row at $\mu =230$, for the lattice refinement model ${\delta }_{\mu }(\mu ,\tau )={\mu }^{-1}$, ${\delta }_{\tau }(\mu ,\tau )={\tau }^{-1}$. The maximum amplitude of the resulting wave function is plotted as a function of $\tau$, for the initial perturbation set at $\tau =113$, $\tau =114$, and $\tau =116$. Analytically, the region of stability is given for $\tau >\mu /2$, i.e., the amplitude of the perturbation grows exponentially for $\tau <115$ and oscillates for $\tau >115$. Here we confirm this numerically. Note the different scale on the graph starting within the stable region.

Figure 10

For the lattice refinement models ${\delta }_{\mu }={\mu }^{-A_{\mu }}$ and ${\delta }_{\tau }={\tau }^{-A_{\tau }}$, the condition of stability has been found. Here we plot $\mu /\tau$ along the critical line, i.e., the line in which the system is just becoming unstable, as a function of the parameter $A$.

Figure 11

The maximum amplitude for a Gaussian wave packet centered on $\mu =1060$ of width $\sigma =5.25$, using the lattice refinement model ${\delta }_{\mu }={\mu }_0\sqrt{\tau }{\mu }^{-1}$ and ${\delta }_{\tau }={\tau }_0{\tau }^{-1/2}$. For small $\tau$ the solution does not vary significantly on scales of the order of the lattice size and hence the solution is stable; however, as soon as the solution begins to vary this stability is lost.

Figure 12

The instability present in the difference equation occurs for large $\mu$ and small $\tau$. If one is interested in reaching the value of the wave function at large $\tau$, one needs to have a large initial range of $\mu$ to have sufficient data. This can result in a portion of the domain being in the unstable region, which makes the code very sensitive to numerical inaccuracies. To avoid this, one can simply use a restricted $\mu$ range and evaluate the wave function as far as possible. This wave function can then be used to create the initial conditions for a subsequent evaluation, provided the wave function has not reached the edge of the lattice.

Figure 13

The twist induced on the wave function can be shown to disappear on large scales by considering all the lattice points being zero except at one point on the initial $\mu$ row. From this the resulting wave function can be calculated; however, it is only in the large-scale limit that there is no deviation in the $\mu$ ($\Theta$) direction. When the discreteness cannot be ignored, there is a net motion toward larger values of $\mu$.

Figure 14

(a) A Gaussian centered on $(\mu ,\tau )=(800,200)$ twists to larger values of $\mu$ as $\tau$ is increased. (b) As we move into the continuum limit this effect disappears. In this case the Gaussian is centered on $(\mu ,\tau )=(800,5000)$. Here the constant lattice case is shown; however, this effect persists when lattice refinement is modeled.