Scale-Dependent Gravity in Quantum Cosmology I: Singularity-Avoiding Sectors and Semiclassical Dynamics

1. Introduction

One of the central questions in quantum cosmology concerns whether the initial singularity of the classical Friedmann picture is an unavoidable feature of cosmological evolution, or whether quantum effects are capable of resolving it. In the standard classical treatment, generic matter sources drive the scale factor to zero, where the energy density and spacetime curvature diverge and the classical description breaks down [1,2]. This singular behavior has long served as a principal motivation for developing cosmological models within a quantum gravitational framework [3].

A broad class of quantum cosmological models produces a suppressed small-a regime together with a large-a oscillatory regime admitting a semiclassical interpretation [4,5,6,7]. In many instances, however, this outcome depends on the introduction of spatial curvature, exotic matter components [8], special boundary conditions [4,5], or nonstandard quantization prescriptions [9,10]. Here we ask a narrower question: how does the quantum-corrected Newtonian cosmology change when the gravitational coupling itself is allowed to run with the cosmological scale factor?

In this paper, we address this question by studying a cosmological model in which the gravitational coupling depends on the scale factor, $G\left(a\right)$ [11], while the remaining cosmological structure is kept as simple as possible. Starting from a quantum-corrected Newtonian cosmology [12,13], we derive a modified Friedmann equation in which the running of the gravitational coupling enters through effective couplings $G^{*}\left(a\right)$, $G_1^{*2}\left(a\right)$, and $G_2^{*2}\left(a\right)$, defined as the first integrals of the equation of motion. A scale-dependent G does more than rescale the standard Friedmann equation: after the first integral is taken, it changes the effective powers and coefficients that appear in the matter, classical-correction, and quantum-correction terms.

In the power-law case $G\left(a\right)=G_0{a}^n$, the small-a behavior is governed entirely by the scaling exponent n, so that the transition between a crunching and a bouncing cosmology is determined by the scaling exponent rather than fine-tuned parameter cancellations. For the flat dust model with positive $\mathsf{\Xi }$, the values $n=3/2$ and $n=2$ mark distinct analytic boundaries: $n=3/2$ is tied to the sign of the quantum-correction coefficient, while $n=2$ marks the degeneracy of the relevant powers of a. A genuine bounce occurs only when the root of $H^2\left(a\right)$ separates a forbidden small-a region from an allowed large-a branch. The global fate of the universe is therefore controlled by the scaling behavior of gravity itself, without modification of the matter content or spatial geometry.

We quantize this model in minisuperspace [14] to obtain the corresponding Wheeler–DeWitt equation [3]. For selected parameter choices, the theory admits exact analytic solutions of both Bessel and Airy type, as well as numerical solutions in physically relevant cases [15]. Recently, model-specific and nonstandard realizations of the Wheeler–DeWitt equation have appeared in the literature [16]. These solutions demonstrate that the wave function can be strongly suppressed near $a=0$, satisfying the DeWitt boundary condition [3] and indicating singularity avoidance, while in other sectors the wave function remains regular but nonvanishing at the origin [4]. The model therefore separates two cases that should not be conflated: a wave function that vanishes at $a=0$, and a wave function that is regular but nonzero there.

We also carry out a WKB analysis and show that, for appropriate parameter choices, the wave function exhibits the structure of a nonsingular quantum cosmology: a classically forbidden region at small scale factor, a turning point, and an oscillatory semiclassical regime at large scale factor. The conserved minisuperspace probability current provides a means to distinguish expanding and contracting branches, giving the oscillatory wave function a dynamical interpretation. Finally, we introduce a homogeneous scalar field as an internal clock [17], which allows wave packets to be constructed and followed in relational time. In the case studied here, the packet peak tracks the corresponding classical Friedmann trajectory at large scale factor.

The paper is organized as follows. In Section 2, we derive the modified Friedmann equation for a running gravitational coupling and define the effective couplings that control the dynamics. In Section 3, we quantize the model and present exact solutions of the Wheeler–DeWitt equation. In Section 4, we develop the WKB analysis and probability-current interpretation. In Section 5, we introduce the scalar-field clock and examine the quantum-to-classical correspondence in relational time. The final sections summarize the physical implications of the model.

2. Quantum–Corrected Newtonian Cosmology

2.1. Newtonian Cosmological Setup

We begin with the standard Newtonian cosmology construction. Assuming homogeneity and isotropy, consider a test particle of negligible mass located on the boundary of a spherical region of physical radius $R\left(t\right)$, filled with a uniform mass density $\rho \left(t\right)$. By symmetry, only the mass enclosed within radius R contributes to the gravitational force on the particle.

The physical radius of the shell is written as

$$R\left(t\right)=a\left(t\right)R_0,$$

(1)

where $a\left(t\right)$ is the cosmological scale factor and $R_0$ is a comoving coordinate label identifying the shell. The quantity $R_0$ is not observable and reflects a freedom in the normalization of comoving coordinates. The enclosed mass is

$$M\left(R\right)={\displaystyle \frac{4\pi }{3}}{R}^3\rho \left(t\right),$$

(2)

and the Hubble parameter is defined kinematically as

$$H\left(t\right)\equiv {\displaystyle \frac{\dot{a}}{a}}={\displaystyle \frac{\dot{R}}{R}}.$$

(3)

We assume that the cosmic fluid obeys a barotropic equation of state,

$$p=(\gamma -1)\rho c^2,$$

(4)

which implies, through the continuity equation,

$$\rho \left(t\right)={\rho }_0{a}^{-3\gamma },$$

(5)

where ${\rho }_0$ denotes the density at the present epoch.

2.2. Quantum–Corrected Newtonian Force Law

To incorporate classical and quantum corrections to Newtonian gravity, we adopt an effective radial force law motivated by quantum corrections to the gravitational potential [12,13,18,19,20]. The acceleration of the test particle is taken to be

$$\ddot{R}=-{\displaystyle \frac{G\left(R\right)M}{R^2}}-2\lambda {\displaystyle \frac{G{\left(R\right)}^2{M}^2}{c^2{R}^3}}+3\xi {\displaystyle \frac{\hslash G{\left(R\right)}^{2}M}{c^3{R}^4}}.$$

(6)

The first term is the standard Newtonian attraction, the second is the leading classical correction, and the third is the leading quantum correction proportional to ℏ. The parameters $\lambda$ and $\xi$ encode the relative strengths of the higher-order classical and quantum contributions.

In the present work we allow the gravitational coupling to depend on scale. Thus G is not treated as a fixed constant, but as a function of the cosmological scale through $G=G\left(R\right)$, or equivalently $G=G\left(a\right)$.

2.3. First Integral of the Equation of Motion

When G varies with scale, the force is not strictly conservative in the ordinary Newtonian sense, and a conserved mechanical energy should not be assumed at the outset. Instead, a first integral can be obtained from Newton’s second law. Multiplying Equation (6) by $\dot{R}$ and integrating with respect to time gives

$$\int \ddot{R}\dot{R}dt=\int f(R,G)dR,$$

(7)

where $f(R,G)$ denotes the right-hand side of Equation (6). Integration gives

$$\begin{array}{cc}\hfill {\displaystyle \frac{1}{2}}{\dot{R}}^2& ={\displaystyle \frac{M}{R}}\left[G-R\int {\displaystyle \frac{dG}{R}}\right]+\lambda {\displaystyle \frac{M^2}{c^2{R}^2}}\left[G^2-2R^2\int {\displaystyle \frac{G}{R^2}}dG\right]\hfill \\ \hfill & -\xi {\displaystyle \frac{M\hslash }{c^3{R}^3}}\left[G^2-2R^3\int {\displaystyle \frac{G}{R^3}}dG\right]+C,\hfill \end{array}$$

(8)

where C is an integration constant.

It is convenient to define the effective gravitational couplings

$$\begin{array}{cc}\hfill G^{*}& \equiv G-R\int {\displaystyle \frac{dG}{R}},\hfill \end{array}$$

(9)

$$\begin{array}{cc}\hfill G_1^{*2}& \equiv G^2-2R^2\int {\displaystyle \frac{G}{R^2}}dG,\hfill \end{array}$$

(10)

$$\begin{array}{cc}\hfill G_2^{*2}& \equiv G^2-2R^3\int {\displaystyle \frac{G}{R^3}}dG.\hfill \end{array}$$

(11)

In terms of these quantities, Equation (8) becomes

$${\displaystyle \frac{1}{2}}{\dot{R}}^2={\displaystyle \frac{G^{*}M}{R}}+\lambda {\displaystyle \frac{G_1^{*2}{M}^2}{c^2{R}^2}}-\xi {\displaystyle \frac{\hslash G_2^{*2}M}{c^3{R}^3}}+C.$$

(12)

The constant C is the Newtonian integration constant and will be identified with the spatial curvature term. The quantities $G^{*}$, $G_1^{*2}$, and $G_2^{*2}$ are not new independent constants. They are effective couplings generated by the scale dependence of G. In the constant-G limit, the integral contributions vanish and the usual Newtonian/Friedmann structure is recovered. For a genuinely running coupling, however, the starred quantities need not scale in the same way as $G\left(R\right)$ itself. They therefore summarize how the running of gravity modifies the matter, classical-correction, and quantum-correction sectors of the cosmological dynamics.

2.4. Friedmann–Type Equation

Dividing Equation (12) by $R^2$ and using $H=\dot{R}/R$, we find

$$H^2={\displaystyle \frac{2G^{*}M}{R^3}}+\lambda {\displaystyle \frac{2G_1^{*2}{M}^2}{c^2{R}^4}}-\xi {\displaystyle \frac{2\hslash G_2^{*2}M}{c^3{R}^5}}+{\displaystyle \frac{2C}{R^2}}.$$

(13)

Substituting

$$M={\displaystyle \frac{4\pi }{3}}{R}^3\rho$$

(14)

gives

$$H^2={\displaystyle \frac{8\pi }{3}}{G}^{*}\rho +\lambda {\displaystyle \frac{32{\pi }^2}{9c^2}}{G}_1^{*2}{R}^2{\rho }^2-\xi {\displaystyle \frac{8\pi }{3}}{\displaystyle \frac{\hslash }{c^3}}{G}_2^{*2}{\displaystyle \frac{\rho }{R^2}}+{\displaystyle \frac{2C}{R^2}}.$$

(15)

Using $R=a{R}_0$, the last term becomes

$${\displaystyle \frac{2C}{R^2}}={\displaystyle \frac{2C}{a^2{R}_0^2}}.$$

(16)

We identify the curvature parameter by

$$k\equiv -{\displaystyle \frac{2C}{c^2{R}_0^2}},$$

(17)

so that

$${\displaystyle \frac{2C}{R^2}}=-{\displaystyle \frac{k{c}^2}{a^2}}.$$

(18)

The comoving scale $R_0$ is arbitrary and unobservable. Any residual appearance of $R_0$ in the higher-order terms may therefore be absorbed into the redefinitions

$$\mathsf{\Lambda }\equiv \lambda R_0^2,\mathsf{\Xi }\equiv \xi R_0^{-2}.$$

(19)

This removes the fiducial comoving length from the final cosmological equation.

Substituting $\rho ={\rho }_0{a}^{-3\gamma }$, the quantum-corrected Friedmann equation becomes

$$\begin{array}{cc}\hfill H^2\left(a\right)& ={\displaystyle \frac{8\pi }{3}}{G}^{*}\left(a\right){\rho }_0{a}^{-3\gamma }+\mathsf{\Lambda }{\displaystyle \frac{32{\pi }^2}{9c^2}}{G}_1^{*2}\left(a\right){\rho }_0^2{a}^{2-6\gamma }\hfill \\ \hfill & -\mathsf{\Xi }{\displaystyle \frac{8\pi }{3}}{\displaystyle \frac{\hslash }{c^3}}{G}_2^{*2}\left(a\right){\rho }_0{a}^{-(3\gamma +2)}-{\displaystyle \frac{k{c}^2}{a^2}}.\hfill \end{array}$$

(20)

This is the basic dynamical equation used throughout the remainder of the paper. It reduces to the standard Friedmann equation when the higher-order corrections are removed and $G^{*}\to G_0$. At small scale factor, however, the correction terms can become significant and may alter the approach to the classical singularity.

2.5. Normalized Density Formulation

It is often convenient to rewrite Equation (20) in terms of dimensionless density parameters. We introduce the present-day critical density

$${\rho }_{c,0}\equiv {\displaystyle \frac{3H_0^2}{8\pi G_0}},$$

(21)

where $H_0$ and $G_0$ are the present values of the Hubble parameter and gravitational coupling, respectively. The corresponding density parameter is

$${\mathsf{\Omega }}_0\equiv {\displaystyle \frac{{\rho }_0}{{\rho }_{c,0}}}={\displaystyle \frac{8\pi G_0{\rho }_0}{3H_0^2}}.$$

(22)

Thus

$$\begin{array}{cc}\hfill {\displaystyle \frac{8\pi }{3}}{\rho }_0& ={\displaystyle \frac{{\mathsf{\Omega }}_0{H}_0^2}{G_0}},\hfill \end{array}$$

(23)

$$\begin{array}{cc}\hfill {\displaystyle \frac{32{\pi }^2}{9}}{\rho }_0^2& ={\displaystyle \frac{{\mathsf{\Omega }}_0^2{H}_0^4}{2G_0^2}}.\hfill \end{array}$$

(24)

Substitution into Equation (20) gives

$$\begin{array}{cc}\hfill H^2\left(a\right)& =\left({\displaystyle \frac{{\mathsf{\Omega }}_0{H}_0^2}{G_0}}\right)G^{*}\left(a\right)a^{-3\gamma }+\mathsf{\Lambda }\left({\displaystyle \frac{{\mathsf{\Omega }}_0^2{H}_0^4}{2G_0^2{c}^2}}\right)G_1^{*2}\left(a\right)a^{2-6\gamma }\hfill \\ \hfill & -\mathsf{\Xi }\left({\displaystyle \frac{{\mathsf{\Omega }}_0{H}_0^2}{G_0}}\right){\displaystyle \frac{\hslash }{c^3}}{G}_2^{*2}\left(a\right)a^{-(3\gamma +2)}-{\displaystyle \frac{k{c}^2}{a^2}}.\hfill \end{array}$$

(25)

Dividing by $H_0^2$, one obtains the form

$$\begin{array}{cc}\hfill {\displaystyle \frac{H^2\left(a\right)}{H_0^2}}& ={\mathsf{\Omega }}_0{\displaystyle \frac{G^{*}\left(a\right)}{G_0}}{a}^{-3\gamma }+\mathsf{\Lambda }{\mathsf{\Omega }}_0^2{\displaystyle \frac{H_0^2}{2G_0^2{c}^2}}{G}_1^{*2}\left(a\right)a^{2-6\gamma }\hfill \\ \hfill & -\mathsf{\Xi }{\mathsf{\Omega }}_0{\displaystyle \frac{\hslash }{G_0{c}^3}}{G}_2^{*2}\left(a\right)a^{-(3\gamma +2)}-{\displaystyle \frac{k{c}^2}{a^2{H}_0^2}}.\hfill \end{array}$$

(26)

This form is especially advantageous for numerical work, because the cosmological evolution is expressed in terms of the density parameters and the scale-dependent effective couplings.

2.6. Worked Examples and Physical Interpretation

To illustrate the physical meaning of Equation (20), we now consider two simple examples in order to show how the higher-order correction terms can change the small-scale dynamics even in a homogeneous and isotropic Newtonian cosmology.

Throughout this subsection we restrict attention to dust matter, $\gamma =1$, and spatial flatness, $k=0$, so that the effects of the classical and quantum corrections can be isolated.

2.6.1. Constant Gravitational Coupling

First consider the constant-G limit,

$$G\left(a\right)=G_0.$$

(27)

Then

$$G^{*}\left(a\right)=G_0,G_1^{*2}\left(a\right)=G_0^2,G_2^{*2}\left(a\right)=G_0^2.$$

(28)

For dust and flat geometry, Equation (20) reduces to

$$H^2\left(a\right)=\mathcal{A}{a}^{-3}+\mathcal{B}{a}^{-4}-\mathcal{C}{a}^{-5},$$

(29)

where

$$\begin{array}{cc}\hfill \mathcal{A}& \equiv {\displaystyle \frac{8\pi }{3}}{G}_0{\rho }_0,\hfill \end{array}$$

(30)

$$\begin{array}{cc}\hfill \mathcal{B}& \equiv \mathsf{\Lambda }{\displaystyle \frac{32{\pi }^2}{9c^2}}{G}_0^2{\rho }_0^2,\hfill \end{array}$$

(31)

$$\begin{array}{cc}\hfill \mathcal{C}& \equiv \mathsf{\Xi }{\displaystyle \frac{8\pi }{3}}{\displaystyle \frac{\hslash }{c^3}}{G}_0^2{\rho }_0.\hfill \end{array}$$

(32)

As $a\to 0$, the ordinary dust contribution diverges as $a^{-3}$, the classical correction as $a^{-4}$, and the quantum correction as $a^{-5}$. Thus the quantum term is the most singular contribution in the constant-G case.

A bounce occurs when $H^2$ vanishes at a finite nonzero scale factor $a_b$, with a forbidden region at smaller a. From Equation (29), the condition $H^2\left(a_b\right)=0$ gives

$$\mathcal{A}{a}_b^2+\mathcal{B}{a}_b-\mathcal{C}=0.$$

(33)

The physically relevant root is

$$a_b={\displaystyle \frac{-\mathcal{B}+\sqrt{{\mathcal{B}}^2+4\mathcal{A}\mathcal{C}}}{2\mathcal{A}}}.$$

(34)

Thus, even without a varying gravitational coupling, the quantum correction can prevent the scale factor from reaching zero. A contracting solution reaches a minimum size and then re-expands, giving a nonsingular cosmological bounce [13,21].

2.6.2. Power-Law Running of the Gravitational Coupling

We now consider the scale-dependent ansatz

$$G\left(a\right)=G_0{a}^n.$$

(35)

This power law is simple enough to be solved in closed-form while capturing the effect of a running gravitational coupling [22,23]. Using Eqs. (10)–(12), one finds

$$G^{*}\left(a\right)={\displaystyle \frac{G_0}{1-n}}{a}^n,G_1^{*2}\left(a\right)={\displaystyle \frac{G_0^2}{1-n}}{a}^{2n},G_2^{*2}\left(a\right)=-{\displaystyle \frac{3G_0^2}{2n-3}}{a}^{2n},$$

(36)

valid for $n\ne 1$ and $n\ne 3/2$. Substitution into Eq. (27) gives

$${\displaystyle \frac{H^2\left(a\right)}{H_0^2}}={\displaystyle \frac{{\mathsf{\Omega }}_0}{1-n}}{a}^{n-3}+\mathsf{\Lambda }{\mathsf{\Omega }}_0^2{\displaystyle \frac{H_0^2}{2c^2(1-n)}}{a}^{2n-4}+3\mathsf{\Xi }{\mathsf{\Omega }}_0{\displaystyle \frac{\hslash G_0}{c^3}}{\displaystyle \frac{1}{2n-3}}{a}^{2n-5}-{\displaystyle \frac{k{c}^2}{a^2{H}_0^2}}.$$

(37)

This expression shows why the exponent n controls the small-a behavior of the Friedmann equation rather than rescaling it. For dust, the matter, classical-correction, and quantum-correction terms scale as

$$a^{n-3},a^{2n-4},a^{2n-5},$$

respectively. Two exponents fix the small-scale structure. The quantum term is the most singular contribution whenever

$$2n-5<n-3⟹n<2,$$

(38)

while the coefficient of the quantum contribution changes sign at

$$n=\frac{3}{2}.$$

(39)

The values $n=3/2$ and $n=2$ therefore mark the natural analytic boundaries of the small-scale dynamics: the former is the sign change of the quantum-correction coefficient, while the latter is the power-degeneracy at which the matter and quantum terms scale with the same power of a.

For the simplified flat model with $\mathsf{\Lambda }=0$, the relevant two-term equation is

$${\displaystyle \frac{H^2\left(a\right)}{H_0^2}}={\displaystyle \frac{{\mathsf{\Omega }}_0}{1-n}}{a}^{n-3}+3\mathsf{\Xi }{\mathsf{\Omega }}_0{\displaystyle \frac{\hslash G_0}{c^3}}{\displaystyle \frac{1}{2n-3}}{a}^{2n-5}.$$

(40)

These two boundaries divide the equation into physically distinct regimes, and it is essential to keep separate two questions that are easily conflated: which term dominates the small-a behavior, and whether that behavior corresponds to a bounce.

In the interval

$$\frac{3}{2}<n<2,$$

(41)

the quantum term is the most singular contribution and therefore controls the small-a dynamics. Its coefficient is positive throughout this window, since $2n-3>0$, while the running renders the matter coefficient negative, since $1-n<0$. The dominant small-a term is thus positive, so $H^2\left(a\right)>0$ as $a\to 0$ and the small-scale region remains classically allowed. Any root of $H^2\left(a\right)$ in this interval is consequently an upper turning point, or recollapse scale, rather than a singularity-avoiding bounce. Dominance of the quantum correction is not, by itself, a bounce.

A genuine bounce requires the opposite sign structure. For $n>2$ the matter term becomes the most singular contribution, and because the running carries the factor $1/(1-n)$ its coefficient is negative; the small-a region is then classically forbidden, with $H^2\left(a\right)<0$, while the quantum term dominates at larger a and restores an allowed expanding branch. The root that separates these regions is a true bounce. The boundary case $n=2$ is degenerate: both terms scale with the same power of a, so no isolated turning point exists, and the bounce branch opens up only for $n>2$.

This distinction is applied in the numerical survey below, where a root of $H^2\left(a\right)$ is classified as a bounce only when it separates a forbidden small-a region from an allowed expanding branch. The condition $H^2\left(a_b\right)=0$ alone is not sufficient: a bounce requires both the vanishing of $H^2$ and the correct sign structure of $H^2\left(a\right)$ on the two sides of the root.

The small-a behavior is therefore controlled by the scaling behavior of gravity itself. Even with dust matter and flat spatial geometry, a running gravitational coupling can make the effective Friedmann equation behave like a multi-fluid system, and the transition between re-collapsing and bouncing behavior can be achieved by the power-law form of $G\left(a\right)$, without the need to add exotic matter by hand.

The quantum correction parameter $\mathsf{\Xi }$ is not an independent phenomenological input. It depends on the coefficient $\xi$ in the quantum-corrected force law (6), which is fixed by the one-loop gravitational scattering amplitude. Specifically, the leading quantum correction to the Newtonian potential takes the form [12,13,18]

$$V\left(r\right)=-{\displaystyle \frac{GM}{r}}\left[1+{\displaystyle \frac{3G(M+m)}{r{c}^2}}+{\displaystyle \frac{122G\hslash }{15\pi r^2{c}^3}}+\cdots \right],$$

(42)

where the coefficient $122/\left(15\pi \right)$ multiplying the $\hslash /r^2$ term is a definite prediction of perturbative quantum gravity, independent of any unknown UV completion [12]. After the rescaling $\mathsf{\Xi }\equiv \xi R_0^{-2}$ that removes the fiducial comoving length (Section 2.4), this value is carried through as the representative value of ${\mathsf{\Xi }}_0=122/\left(15\pi \right)$. As shown in Figure 1B, the bounce scale is insensitive to order-of-magnitude variations in $\mathsf{\Xi }$ for $n>2$, so the qualitative results are not fine-tuned to this choice.

2.7. Possible Origins of Scale-Dependent Power-Law Running

The ansatz $G\left(a\right)=G_0{a}^n$ is phenomenological, but it appears naturally in several extensions of gravity. In asymptotically safe quantum gravity, the running Newton constant near a UV fixed point scales as $G\left(k\right)\sim k^{-2}$, which translates to $G\left(a\right)\sim a^2$ in a cosmological setting where the renormalization group scale is identified with the Hubble parameter or inverse scale factor [24,25]. In scalar-tensor theories, a coupling function $F\left(\varphi \right)\sim e^{\alpha \varphi }$ after a field redefinition can yield an effective power-law $G\left(a\right)$ when $\varphi$ evolves logarithmically with a [26]. Kaluza–Klein compactifications with a time-varying modulus also produce $G_{\mathrm{eff}}\left(a\right)\sim a^n$ in the four-dimensional effective theory [27]. We do not derive $G\left(a\right)$ from these frameworks here, but note that each provides a physical rationale for the phenomenological ansatz.

While the strict power-law ansatz $G\left(a\right)=G_0{a}^n$ serves as a simple, analytically tractable description of the deep ultraviolet quantum regime ($a\to 0$), it must undergo an infrared (IR) freeze-out in the late universe to comply with tight observational bounds on $\dot{G}/G$ from Big Bang Nucleosynthesis (BBN) and solar system tests. Globally, this can be accommodated by an interpolating profile such as:

$$G\left(a\right)=G_0{\left[1+{\left({\displaystyle \frac{a_x}{a}}\right)}^n\right]}^{-1},$$

(43)

where $a_x$ represents a quantum transition scale. For the bounce scales identified in Section 2.8, the turning point satisfies $a_t\ll 1$, so the quantum regime $a\lesssim a_x$ must itself lie well below the BBN epoch ($a_{\mathrm{BBN}}\sim {10}^{-10}$); any transition scale $a_x$ satisfying this bound leaves the standard thermal history entirely unaffected, and the power-law scaling is then a self-consistent description of gravity in the pre-BBN ultraviolet regime without further constraint from late-time observations. The early-time limit $a\ll a_x$ yields $G\left(a\right)\simeq G_0{(a/a_x)}^n$, reproducing the same power-law scaling up to an overall normalization. For $a\gg a_x$, the coupling approaches $G_0$. In this paper we use the pure power law as a local early-time scaling model; a complete phenomenological treatment would require a specified transition scale and a comparison with observational bounds. The observational consequences of allowing this running to persist into later epochs — including the BBN era, the epoch of structure formation, and recombination — are examined in [37].

2.8. Bounce-Scale Parameter Survey

The power-law model also allows a simple analytic and numerical survey of the turning-point scale. To isolate the role of the running exponent, we consider the dust, flat, and purely quantum-corrected case,

$$\gamma =1,k=0,\mathsf{\Lambda }=0,$$

with

$$G\left(a\right)=G_0{a}^n.$$

In this limit the Friedmann equation reduces to

$${\displaystyle \frac{H^2\left(a\right)}{H_0^2}}={\displaystyle \frac{{\mathsf{\Omega }}_0}{1-n}}{a}^{n-3}+3\mathsf{\Xi }{\mathsf{\Omega }}_0{\displaystyle \frac{\hslash G_0}{c^3}}{\displaystyle \frac{1}{2n-3}}{a}^{2n-5}.$$

(44)

For the normalized survey shown in Figure 1, we set $c=\hslash =G_0=H_0=1$, so that the quantum correction is controlled by the parameter $\mathsf{\Xi }$. The specific value $\mathsf{\Xi }=122/\left(15\pi \right)$ is chosen as a representative parameter: it makes the quantum correction term comparable to the matter term at the bounce scale without fully dominating it. Over the range $\mathsf{\Xi }\in [0.1,10]$, the bounce scale shifts by less than a factor of 3 for $n=3$ (see Figure 1), so the qualitative results are not fine-tuned.

A turning point occurs when $H^2\left(a_t\right)=0$. From Equation (44), this gives

$$a_t(n,\mathsf{\Xi })={\left[{\displaystyle \frac{2n-3}{3\mathsf{\Xi }(n-1)}}\right]}^{1/(n-2)},n\ne 1,n\ne {\displaystyle \frac{3}{2}},n\ne 2.$$

(45)

Classification of turning points.

A root $a_t$ of $H^2\left(a_t\right)=0$ can be of three types:

• Type I (Bounce): $H^2\left(a\right)<0$ for $a<a_t$ (forbidden), $H^2\left(a\right)>0$ for $a>a_t$ (allowed). The universe contracts, stops, and re-expands.
• Type II (Recollapse): $H^2\left(a\right)>0$ for $a<a_t$, $H^2\left(a\right)<0$ for $a>a_t$. The universe expands, stops, and contracts.
• Type III (Degenerate): No sign change across $a_t$.

Only Type I corresponds to a nonsingular bounce that replaces the classical singularity.

However, not every turning point is a true bounce. A bounce requires more than $H^2\left(a_t\right)=0$: the region $a<a_t$ must be classically forbidden, while the region $a>a_t$ must be classically allowed. In other words, the root must separate an inaccessible small-scale region from an expanding cosmological branch.

For $3/2<n<2$, the equation may possess a positive root, but the small-a region remains allowed. Such a solution is better interpreted as an upper turning point or recollapse scale, not as a singularity-avoiding bounce. By contrast, for $n>2$, the leading small-a behavior makes $H^2\left(a\right)<0$ near the origin, while the larger-a branch becomes allowed. Thus, in the parameter range emphasized here, the true bounce branch begins at $n=2$. The point $n=2$ itself is degenerate because the two terms in Equation (44) scale with the same power of a, so there is no isolated turning point.

Figure 1 summarizes the structure of the bounce branch. Panel (A) shows that, for the reference value ${\mathsf{\Xi }}_0=122/\left(15\pi \right)$, the true bounce scale grows from very small values just above $n=2$ and approaches order unity as n becomes large. For example, at $n=3$,

$$a_t={\displaystyle \frac{1}{2{\mathsf{\Xi }}_0}}\simeq 0.193.$$

(46)

Panel (B) shows that the dependence on $\mathsf{\Xi }$ is highly sensitive near the degenerate value $n=2$, but becomes much milder for larger n. Panel (C) displays the same information as a phase diagram in the $(n,\mathsf{\Xi })$ plane, and panel (D) translates the bounce scale into the corresponding bounce redshift.

The existence and location of the bounce are controlled primarily by the scaling exponent n. The parameter $\mathsf{\Xi }$ shifts the location of the bounce, but the qualitative transition between a true bounce and a non-bouncing solution is governed by the power-law behavior of the gravitational coupling itself.

2.9. Connection to the Wheeler–DeWitt Analysis

The modified Friedmann equation derived in this section provides the starting point for the minisuperspace quantization in Section 3. In the quantum theory, the same scale dependence of $G\left(a\right)$ enters the Wheeler–DeWitt equation through the effective minisuperspace potential. Consequently, the classical turning-point structure and the quantum behavior of the wave function near $a=0$ are governed by the same underlying running of the gravitational coupling.

Both classical and quantum descriptions share the same scale dependence. Classically, the quantum-corrected Friedmann equation determines whether the universe reaches a singularity, turns around, or undergoes a nonsingular bounce. Quantum mechanically, the corresponding Wheeler–DeWitt potential determines whether the wave function is suppressed near the classical singularity, remains regular, or becomes oscillatory. The analysis that follows therefore treats the Friedmann dynamics and the Wheeler–DeWitt equation as two sides of the same scale-dependent gravitational model.

3. Minisuperspace Quantization and the Wheeler–DeWitt Equation

The Wheeler–DeWitt framework and its minisuperspace realization have a long history in quantum cosmology [3,14,28]. In this section we quantize the scale-factor degree of freedom and derive a Wheeler–DeWitt-type of equation associated with the modified Friedmann dynamics obtained in Section 2. We quantize an effective minisuperspace of a one-dimensional cosmological model where the Hamiltonian constraint is chosen so that its classical limit reproduces the modified Friedmann equation derived above. This is sufficient for the present purpose, which is to determine how the running gravitational coupling modifies the quantum behavior of the scale factor.

The scale factor $a\left(t\right)$ is taken to be dimensionless, with the physical radius written as $R\left(t\right)=a\left(t\right)R_0$. The constants c, ℏ, and $G_0$ are kept throughout. This convention makes it clear where the usual minisuperspace structure is being modified by the scale dependence of the gravitational coupling. However, for several numerical examples and parameter surveys we adopt units in which $c=\hslash =G_0=1$.

3.1. Reduced Minisuperspace Dynamics

We begin with the effective gravitational kinetic term

$$L_{\mathrm{kin}}(a,\dot{a})=-{\displaystyle \frac{3\pi c^3}{4G\left(a\right)}}a{\dot{a}}^2,$$

(47)

and write the reduced Lagrangian as

$$L(a,\dot{a})=L_{\mathrm{kin}}(a,\dot{a})-V\left(a\right),$$

(48)

where $V\left(a\right)$ denotes the effective minisuperspace potential. In the present model $V\left(a\right)$ contains the matter contribution, the curvature term, the classical correction proportional to $\mathsf{\Lambda }$, and the quantum correction proportional to $\mathsf{\Xi }$.

The canonical momentum conjugate to a is

$$p_a\equiv {\displaystyle \frac{\partial L}{\partial \dot{a}}}=-{\displaystyle \frac{3\pi c^3}{2G\left(a\right)}}a\dot{a},$$

(49)

so that

$$\dot{a}=-{\displaystyle \frac{2G\left(a\right)}{3\pi c^3}}{\displaystyle \frac{p_a}{a}}.$$

(50)

The corresponding Hamiltonian is

$$\mathcal{H}=p_a\dot{a}-L=-{\displaystyle \frac{G\left(a\right)}{3\pi c^{3}a}}{p}_a^2+V\left(a\right).$$

(51)

The classical dynamics are obtained from the Hamiltonian constraint

$$\mathcal{H}=0.$$

(52)

Using Equation (49), the constraint implies

$$H^2\left(a\right)\equiv {\left({\displaystyle \frac{\dot{a}}{a}}\right)}^2={\displaystyle \frac{4G\left(a\right)}{3\pi c^3}}{\displaystyle \frac{V\left(a\right)}{a^3}}.$$

(53)

Rearranging,

$$V\left(a\right)={\displaystyle \frac{3\pi c^3}{4G\left(a\right)}}{a}^3{H}^2\left(a\right).$$

(54)

Thus the effective potential is not introduced independently; it is fixed by the modified Friedmann function $H^2\left(a\right)$.

For the scale-dependent gravity model considered here,

$$\begin{array}{cc}\hfill H^2\left(a\right)=& {\displaystyle \frac{8\pi }{3}}{G}^{*}\left(a\right){\rho }_0{a}^{-3\gamma }+\mathsf{\Lambda }{\displaystyle \frac{32{\pi }^2}{9c^2}}{G}_1^{*2}\left(a\right){\rho }_0^2{a}^{2-6\gamma }\hfill \\ \hfill & -\mathsf{\Xi }{\displaystyle \frac{8\pi }{3}}{\displaystyle \frac{\hslash }{c^3}}{G}_2^{*2}\left(a\right){\rho }_0{a}^{-(3\gamma +2)}-{\displaystyle \frac{k{c}^2}{a^2}}.\hfill \end{array}$$

(55)

This is the same modified Friedmann equation derived in Section 2. The same function that determines the classical turning-point structure will also determine the effective potential in the quantum theory.

3.2. Factor Ordering and Wheeler–DeWitt Equation

Upon quantization, we promote

$$p_a\to -i\hslash {\displaystyle \frac{d}{da}}.$$

(56)

Because the scale factor is a minisuperspace variable, the ordering of $p_a^2$ is not unique. We adopt the standard one-parameter factor-ordering prescription

$$p_a^2\psi \left(a\right)⟶-{\hslash }^2\left({\displaystyle \frac{d^2\psi }{d{a}^2}}+{\displaystyle \frac{q}{a}}{\displaystyle \frac{d\psi }{da}}\right),$$

(57)

where q is the factor-ordering parameter. This convention includes commonly used choices in minisuperspace quantum cosmology, while making explicit that different orderings can shift the subleading behavior of the wave function near $a=0$ [3,14,28]. The ordering $q=1$ corresponds to the Laplace–Beltrami operator on the minisuperspace metric $\mathrm{d}{s}^2=-a\mathrm{d}{a}^2$, and is the choice most directly motivated by the DeWitt quantization prescription [3,14]. The value $q=3/4$ is used in Section 3.4.1 as an example that simultaneously satisfies the DeWitt boundary condition $\psi \left(0\right)=0$ and yields a damped oscillatory envelope at large a; it is not singled out as physically preferred over $q=1$, but rather chosen to exhibit both properties in a single analytic example. The model therefore does not derive singularity avoidance from the choice of ordering: the same scale-dependent potential structure that produces a classically forbidden small-a region in the Friedmann equation also suppresses the wave function near the origin, independent of the factor-ordering parameter.

The Wheeler–DeWitt equation $\widehat{\mathcal{H}}\psi =0$ is therefore

$${\displaystyle \frac{{\hslash }^{2}G\left(a\right)}{3\pi c^{3}a}}\left({\psi }^{″}\left(a\right)+{\displaystyle \frac{q}{a}}{\psi }^{\prime }\left(a\right)\right)+V\left(a\right)\psi \left(a\right)=0.$$

(58)

Eliminating $V\left(a\right)$ with Equation (54), we obtain the compact form

$${\psi }^{″}\left(a\right)+{\displaystyle \frac{q}{a}}{\psi }^{\prime }\left(a\right)+{\displaystyle \frac{9{\pi }^2{c}^6}{4{\hslash }^2}}{\displaystyle \frac{a^4}{G{\left(a\right)}^2}}{H}^2\left(a\right)\psi \left(a\right)=0.$$

(59)

This equation is the central quantum equation used in the remainder of the paper. It shows explicitly how the running gravitational coupling enters the Wheeler–DeWitt potential through both $G\left(a\right)$ and the modified Friedmann function $H^2\left(a\right)$.

3.3. Power-Law Running

We now specialize to the power-law running

$$G\left(a\right)=G_0{a}^n.$$

(60)

Using the starred couplings derived in Section 2,

$$G^{*}\left(a\right)={\displaystyle \frac{G_0}{1-n}}{a}^n,G_1^{*2}\left(a\right)={\displaystyle \frac{G_0^2}{1-n}}{a}^{2n},G_2^{*2}\left(a\right)=-{\displaystyle \frac{3G_0^2}{2n-3}}{a}^{2n},$$

(61)

valid for $n\ne 1$ and $n\ne 3/2$, Equation (59) becomes

$${\psi }^{″}\left(a\right)+{\displaystyle \frac{q}{a}}{\psi }^{\prime }\left(a\right)+{\displaystyle \frac{9{\pi }^2{c}^6}{4{\hslash }^2{G}_0^2}}{a}^{4-2n}{H}^2\left(a\right)\psi \left(a\right)=0.$$

(62)

The corresponding Friedmann function is

$$\begin{array}{cc}\hfill H^2\left(a\right)=& {\displaystyle \frac{8\pi }{3}}{\displaystyle \frac{G_0}{1-n}}{\rho }_0{a}^{n-3\gamma }+\mathsf{\Lambda }{\displaystyle \frac{32{\pi }^2}{9c^2}}{\displaystyle \frac{G_0^2}{1-n}}{\rho }_0^2{a}^{2n+2-6\gamma }\hfill \\ \hfill & +\mathsf{\Xi }{\displaystyle \frac{8\pi \hslash }{c^3}}{\displaystyle \frac{G_0^2}{2n-3}}{\rho }_0{a}^{2n-(3\gamma +2)}-{\displaystyle \frac{k{c}^2}{a^2}}.\hfill \end{array}$$

(63)

Substitution gives the expanded Wheeler–DeWitt equation

$$\begin{array}{cc}\hfill {\psi }^{″}\left(a\right)+{\displaystyle \frac{q}{a}}{\psi }^{\prime }\left(a\right)+[& {\displaystyle \frac{6{\pi }^3{c}^6}{{\hslash }^2}}{\displaystyle \frac{{\rho }_0}{(1-n)G_0}}{a}^{4-n-3\gamma }+{\displaystyle \frac{8{\pi }^4{c}^4}{{\hslash }^2}}{\displaystyle \frac{\mathsf{\Lambda }{\rho }_0^2}{1-n}}{a}^{6-6\gamma }\hfill \\ \hfill & +\mathsf{\Xi }{\displaystyle \frac{18{\pi }^3{c}^3}{\hslash }}{\displaystyle \frac{{\rho }_0}{2n-3}}{a}^{2-3\gamma }-{\displaystyle \frac{9{\pi }^{2}k{c}^8}{4{\hslash }^2{G}_0^2}}{a}^{2-2n}]\psi \left(a\right)=0.\hfill \end{array}$$

(64)

The running exponent n changes the powers of a appearing in the Wheeler–DeWitt potential. Thus the small-scale behavior of the wave function is controlled by the matter equation of state, the factor ordering, and by the scaling behavior of the gravitational coupling itself.

3.4. Representative Analytic and Numerical Solutions

We now examine representative solutions of Equation (64) where we wish to display the main types of behavior admitted by the model. We focus on flat models with

$$k=0,\mathsf{\Lambda }=0,$$

so that the competition between the ordinary matter term and the $\mathsf{\Xi }$-induced quantum correction can be highlighted. Exact and semiclassical solutions of the Wheeler–DeWitt equation are known in special minisuperspace settings [15,28], and the present model provides additional tractable examples in which the effect of the running gravitational coupling can be seen explicitly.

For $n=2$, the ordinary matter contribution and the $\mathsf{\Xi }$-dependent quantum correction in Equation (64) scale with the same power of a. The Wheeler–DeWitt equation then reduces to

$${\psi }^{″}\left(a\right)+{\displaystyle \frac{q}{a}}{\psi }^{\prime }\left(a\right)+{\lambda }_{\mathrm{eff}}{a}^{2-3\gamma }\psi \left(a\right)=0,$$

(65)

where

$${\lambda }_{\mathrm{eff}}=-{\displaystyle \frac{6{\pi }^3{c}^6{\rho }_0}{{\hslash }^2{G}_0}}+{\displaystyle \frac{18{\pi }^3\mathsf{\Xi }{c}^3{\rho }_0}{\hslash }}.$$

(66)

For the representative value of

$$\mathsf{\Xi }={\displaystyle \frac{122}{15\pi }},$$

(67)

this becomes

$${\lambda }_{\mathrm{eff}}=-{\displaystyle \frac{6{\pi }^3{c}^6{\rho }_0}{{\hslash }^2{G}_0}}+{\displaystyle \frac{732}{5}}{\displaystyle \frac{{\pi }^2{c}^3{\rho }_0}{\hslash }}.$$

(68)

In the examples below we take ${\lambda }_{\mathrm{eff}}>0$, so that the solutions are oscillatory at sufficiently large scale factor. If ${\lambda }_{\mathrm{eff}}<0$, the corresponding solutions are expressed in terms of modified Bessel functions, and the oscillatory interpretation must be modified.

3.4.1. Dust-Era Bessel Solution Satisfying the DeWitt Condition

For a dust-dominated universe, $\gamma =1$, Equation (65) becomes

$${\psi }^{″}\left(a\right)+{\displaystyle \frac{q}{a}}{\psi }^{\prime }\left(a\right)+{\lambda }_{\mathrm{eff}}{a}^{-1}\psi \left(a\right)=0.$$

(69)

The general solution is

$$\psi \left(a\right)=a^{(1-q)/2}\left[C_1{J}_{1-q}\left(2\sqrt{{\lambda }_{\mathrm{eff}}a}\right)+C_2{Y}_{1-q}\left(2\sqrt{{\lambda }_{\mathrm{eff}}a}\right)\right].$$

(70)

The $Y_{1-q}$ branch is singular at the origin and is discarded if regularity at $a=0$ is required. The regular branch is therefore

$$\psi \left(a\right)=C{a}^{(1-q)/2}{J}_{1-q}\left(2\sqrt{{\lambda }_{\mathrm{eff}}a}\right).$$

(71)

Near the origin,

$$\psi \left(a\right)\sim a^{1-q}.$$

(72)

Thus this branch satisfies the DeWitt boundary condition,

$$\psi \left(0\right)=0,$$

(73)

whenever $q<1$. The DeWitt condition is commonly interpreted as suppression of the classical singular geometry in the quantum theory [3].

At large a, the Bessel asymptotics give

$$\psi \left(a\right)\sim a^{(1-2q)/4}cos\left[2\sqrt{{\lambda }_{\mathrm{eff}}a}-{\displaystyle \frac{\pi }{2}}(1-q)-{\displaystyle \frac{\pi }{4}}\right].$$

(74)

The envelope decays for $1/2<q<1$, is asymptotically constant for $q=1/2$, and grows for $q<1/2$.

An illustrative parameterization is

$$q={\displaystyle \frac{3}{4}},$$

(75)

for which

$$\psi \left(a\right)=C{a}^{1/8}{J}_{1/4}\left(2\sqrt{{\lambda }_{\mathrm{eff}}a}\right).$$

(76)

This solution vanishes at $a=0$, satisfies the DeWitt criterion, and exhibits damped oscillations at large scale factor. It therefore provides a simple analytic example of a state that is suppressed near the classical singularity and becomes semiclassical in the large-a regime.

3.4.2. Regular Dust-Era Solution Without DeWitt Suppression

The case $q=1$ results in a valuable comparison. Equation (70) becomes

$$\psi \left(a\right)=C_1{J}_0\left(2\sqrt{{\lambda }_{\mathrm{eff}}a}\right)+C_2{Y}_0\left(2\sqrt{{\lambda }_{\mathrm{eff}}a}\right).$$

(77)

Since $Y_0\left(x\right)$ diverges logarithmically as $x\to 0$, regularity at the origin selects

$$\psi \left(a\right)=C{J}_0\left(2\sqrt{{\lambda }_{\mathrm{eff}}a}\right).$$

(78)

However,

$$J_0\left(0\right)=1,$$

(79)

so this solution is finite but nonzero at $a=0$. It therefore satisfies a regularity condition, but not the stronger DeWitt condition. A regular nonzero wave function at the origin is more naturally interpreted as a regular quantum origin than as strict DeWitt suppression of the singular geometry. This is closer in spirit to the Hartle–Hawking regularity viewpoint than to imposing $\psi \left(0\right)=0$ as a boundary condition [4].

Figure 2. Dust-case Wheeler–DeWitt wave function for $n=2$, $\gamma =1$, $q=\frac{3}{4}$, $k=\mathsf{\Lambda }=0$, and $\mathsf{\Xi }=122/\left(15\pi \right)$. The exact solution $\psi \left(a\right)\propto a^{1/8}{J}_{1/4}\left(2\sqrt{{\lambda }_{\mathrm{eff}}a}\right)$ vanishes at $a=0$, satisfying the DeWitt boundary condition, and exhibits damped oscillations for increasing scale factor.

Figure 2. Dust-case Wheeler–DeWitt wave function for $n=2$, $\gamma =1$, $q=\frac{3}{4}$, $k=\mathsf{\Lambda }=0$, and $\mathsf{\Xi }=122/\left(15\pi \right)$. The exact solution $\psi \left(a\right)\propto a^{1/8}{J}_{1/4}\left(2\sqrt{{\lambda }_{\mathrm{eff}}a}\right)$ vanishes at $a=0$, satisfying the DeWitt boundary condition, and exhibits damped oscillations for increasing scale factor.

For large a,

$$\psi \left(a\right)\sim a^{-1/4}cos\left(2\sqrt{{\lambda }_{\mathrm{eff}}a}-{\displaystyle \frac{\pi }{4}}\right),$$

(80)

so the solution is again oscillatory with a slowly decaying envelope. The factor-ordering parameter therefore controls whether the origin is regular or suppressed.

Figure 3. Dust-case Wheeler–DeWitt wave function for $n=2$, $\gamma =1$, $q=1$, $k=\mathsf{\Lambda }=0$, and $\mathsf{\Xi }=122/\left(15\pi \right)$. The regular solution $\psi \left(a\right)\propto J_0\left(2\sqrt{{\lambda }_{\mathrm{eff}}a}\right)$ remains finite and nonzero at $a=0$. It therefore satisfies regularity at the origin, but not the DeWitt boundary condition.

Figure 3. Dust-case Wheeler–DeWitt wave function for $n=2$, $\gamma =1$, $q=1$, $k=\mathsf{\Lambda }=0$, and $\mathsf{\Xi }=122/\left(15\pi \right)$. The regular solution $\psi \left(a\right)\propto J_0\left(2\sqrt{{\lambda }_{\mathrm{eff}}a}\right)$ remains finite and nonzero at $a=0$. It therefore satisfies regularity at the origin, but not the DeWitt boundary condition.

3.4.3. A Non-Bessel Airy Example

The model also admits exact analytic sectors that are not of Bessel type. For

$$n=2,\gamma ={\displaystyle \frac{1}{3}},q=0,k=\mathsf{\Lambda }=0,$$

Equation (65) becomes

$${\psi }^{″}\left(a\right)+{\lambda }_{\mathrm{eff}}a\psi \left(a\right)=0.$$

(81)

Defining

$$x=-{\lambda }_{\mathrm{eff}}^{1/3}a,$$

(82)

this becomes the standard Airy equation,

$${\displaystyle \frac{d^2\psi }{d{x}^2}}-x\psi =0.$$

(83)

The general solution is

$$\psi \left(a\right)=C_1\mathrm{Ai}\left(-{\lambda }_{\mathrm{eff}}^{1/3}a\right)+C_2\mathrm{Bi}\left(-{\lambda }_{\mathrm{eff}}^{1/3}a\right).$$

(84)

A DeWitt-compatible linear combination is obtained by imposing $\psi \left(0\right)=0$. Since

$${\displaystyle \frac{\mathrm{Ai}\left(0\right)}{\mathrm{Bi}\left(0\right)}}={\displaystyle \frac{1}{\sqrt{3}}},$$

(85)

one obtains

$$\psi \left(a\right)=C\left[\mathrm{Ai}\left(-{\lambda }_{\mathrm{eff}}^{1/3}a\right)-{\displaystyle \frac{1}{\sqrt{3}}}\mathrm{Bi}\left(-{\lambda }_{\mathrm{eff}}^{1/3}a\right)\right].$$

(86)

This solution vanishes at $a=0$ by construction and becomes oscillatory for large positive a. Although $\gamma =1/3$ corresponds to an exotic equation of state rather than ordinary dust or radiation, the example is revealing because it shows that the scale-dependent model supports exact non-Bessel wave-function sectors.

Figure 4. Airy-type Wheeler–DeWitt wave function for $n=2$, $\gamma =\frac{1}{3}$, $q=0$, $k=\mathsf{\Lambda }=0$, and $\mathsf{\Xi }=122/\left(15\pi \right)$. The linear combination $\psi \left(a\right)\propto \mathrm{Ai}(-{\lambda }_{\mathrm{eff}}^{1/3}a)-\mathrm{Bi}(-{\lambda }_{\mathrm{eff}}^{1/3}a)/\sqrt{3}$ is chosen so that $\psi \left(0\right)=0$, thereby satisfying the DeWitt condition.

Figure 4. Airy-type Wheeler–DeWitt wave function for $n=2$, $\gamma =\frac{1}{3}$, $q=0$, $k=\mathsf{\Lambda }=0$, and $\mathsf{\Xi }=122/\left(15\pi \right)$. The linear combination $\psi \left(a\right)\propto \mathrm{Ai}(-{\lambda }_{\mathrm{eff}}^{1/3}a)-\mathrm{Bi}(-{\lambda }_{\mathrm{eff}}^{1/3}a)/\sqrt{3}$ is chosen so that $\psi \left(0\right)=0$, thereby satisfying the DeWitt condition.

3.4.4. Numerical Radiation-Era Example

As a radiation-era example, consider a radiation-dominated universe with

$$\gamma ={\displaystyle \frac{4}{3}},n=-1,q=0,k=\mathsf{\Lambda }=0,\mathsf{\Xi }={\displaystyle \frac{122}{15\pi }}.$$

In this case Equation (64) has the schematic form

$${\psi }^{″}\left(a\right)+\left(Aa-{\displaystyle \frac{B}{a^2}}\right)\psi \left(a\right)=0,A>0,B>0,$$

(87)

after collecting constants. The inverse-square term dominates as $a\to 0$ and produces a suppressed regular branch. At larger scale factor the positive linear term dominates, leading to an oscillatory regime. The numerical solution therefore transitions from suppression near $a=0$ to oscillatory behavior at larger a.

3.5. Local Airy form Near a Turning Point

Before turning to the full WKB analysis, we examine the local behavior near a simple turning point. Removing the first-derivative term by writing

$$\psi \left(a\right)=a^{-q/2}u\left(a\right),$$

(88)

the Wheeler–DeWitt equation takes the reduced form

$$u^{″}\left(a\right)+U_{\mathrm{eff}}\left(a\right)u\left(a\right)=0.$$

(89)

Suppose that the effective potential has a simple turning point $a_t$, so that

$$U_{\mathrm{eff}}\left(a_t\right)=0,U_{\mathrm{eff}}^{\prime }\left(a_t\right)\ne 0.$$

(90)

Expanding linearly,

$$U_{\mathrm{eff}}\left(a\right)\simeq U_{\mathrm{eff}}^{\prime }\left(a_t\right)(a-a_t),$$

(91)

Equation (89) reduces to Airy’s equation. The local solution may be written as

$$u\left(a\right)=c_1\mathrm{Ai}\left({\displaystyle \frac{a_t-a}{s}}\right)+c_2\mathrm{Bi}\left({\displaystyle \frac{a_t-a}{s}}\right),s={\left[U_{\mathrm{eff}}^{\prime }\left(a_t\right)\right]}^{-1/3},$$

(92)

up to the sign convention associated with $U_{\mathrm{eff}}^{\prime }\left(a_t\right)$. Selecting the $\mathrm{Ai}$ branch gives exponential suppression on the forbidden side of the turning point and oscillatory behavior on the allowed side. This is the standard local structure underlying WKB matching in quantum cosmology [4,5].

Figure 5. Numerical Wheeler–DeWitt wave function for the radiation-dominated case $\gamma =\frac{4}{3}$, with $n=-1$, $q=0$, $k=0$, $\mathsf{\Lambda }=0$, and $\mathsf{\Xi }=122/\left(15\pi \right)$. The regular branch is suppressed near $a=0$, satisfying the DeWitt condition, and becomes oscillatory at larger scale factor.

Figure 5. Numerical Wheeler–DeWitt wave function for the radiation-dominated case $\gamma =\frac{4}{3}$, with $n=-1$, $q=0$, $k=0$, $\mathsf{\Lambda }=0$, and $\mathsf{\Xi }=122/\left(15\pi \right)$. The regular branch is suppressed near $a=0$, satisfying the DeWitt condition, and becomes oscillatory at larger scale factor.

Figure 6. Illustrative Airy-type semiclassical Wheeler–DeWitt wave function near a single turning point $a_t$. The local form $\psi \left(a\right)=a^{-q/2}\mathrm{Ai}((a_t-a)/s)$, shown here for $q=1$, is exponentially suppressed on the forbidden side of the turning point and oscillatory on the allowed side. This figure illustrates the universal Airy matching behavior that appears near a simple Wheeler–DeWitt turning point.

Figure 6. Illustrative Airy-type semiclassical Wheeler–DeWitt wave function near a single turning point $a_t$. The local form $\psi \left(a\right)=a^{-q/2}\mathrm{Ai}((a_t-a)/s)$, shown here for $q=1$, is exponentially suppressed on the forbidden side of the turning point and oscillatory on the allowed side. This figure illustrates the universal Airy matching behavior that appears near a simple Wheeler–DeWitt turning point.

3.6. Summary of the Wheeler–DeWitt Sectors

The examples above show that the scale-dependent model admits several distinct quantum behaviors. The dust solution with $q=\frac{3}{4}$ is revealing: it is analytic, physically transparent, satisfies the DeWitt condition, and has a damped oscillatory large-a limit. The $q=1$ dust solution shows that regularity at $a=0$ is not the same as DeWitt suppression. The Airy example illustrates that the model is not restricted to Bessel-type sectors, while the radiation-era numerical solution provides a more familiar matter source in which the small-a region is suppressed and the large-a regime becomes oscillatory.

These results prepare us for the WKB analysis in Section 4. There, the same effective potential structure will be used to identify forbidden and allowed regions, turning points, and the branch interpretation of the minisuperspace wave function.

4. WKB Dynamics, Turning Points, and Probability Current

The exact solutions discussed in Section 3 show explicitly that the scale-dependent model admits wave functions that are either suppressed at the classical singularity or regular at the origin. We now turn to the WKB regime, where the connection between the Wheeler–DeWitt wave function and classical cosmological dynamics becomes explicit. The purpose of this section is to show how a forbidden small-a region, a turning point, and an oscillatory large-a branch arise in a representative case, and how the probability current distinguishes expanding and contracting semiclassical sectors.

We consider the flat dust-dominated model

$$\gamma =1,n=3,q=1,k=0,\mathsf{\Lambda }=0,$$

(93)

with the quantum-correction parameter

$$\mathsf{\Xi }={\displaystyle \frac{122}{15\pi }}.$$

(94)

This example is instructive because it reduces the Wheeler–DeWitt equation to a single-turning-point problem while retaining the effect of the running gravitational coupling. The resulting WKB structure is the familiar one: exponential behavior in a forbidden region, Airy matching near the turning point, and oscillatory behavior in the allowed region [4,5,29].

4.1. Reduced Wheeler–DeWitt Equation for a Chosen Model

We begin with the power-law Wheeler–DeWitt equation derived in Section 3,

$${\psi }^{″}\left(a\right)+{\displaystyle \frac{q}{a}}{\psi }^{\prime }\left(a\right)+{\displaystyle \frac{9{\pi }^2{c}^6}{4{\hslash }^2{G}_0^2}}{a}^{4-2n}{H}^2\left(a\right)\psi \left(a\right)=0.$$

(95)

For $q=1$ and $n=3$, this becomes

$${\psi }^{″}\left(a\right)+{\displaystyle \frac{1}{a}}{\psi }^{\prime }\left(a\right)+{\displaystyle \frac{9{\pi }^2{c}^6}{4{\hslash }^2{G}_0^2}}{a}^{-2}{H}^2\left(a\right)\psi \left(a\right)=0.$$

(96)

For dust matter, $\gamma =1$, with $k=0$ and $\mathsf{\Lambda }=0$, the modified Friedmann function reduces to

$$H^2\left(a\right)={\displaystyle \frac{8\pi }{3}}{\displaystyle \frac{G_0}{1-3}}{\rho }_0{a}^{3-3}+\mathsf{\Xi }{\displaystyle \frac{8\pi \hslash }{c^3}}{\displaystyle \frac{G_0^2}{2·3-3}}{\rho }_0{a}^{6-(3+2)}.$$

(97)

Hence

$$H^2\left(a\right)=-{\displaystyle \frac{4\pi }{3}}{G}_0{\rho }_0+\mathsf{\Xi }{\displaystyle \frac{8\pi \hslash }{3c^3}}{G}_0^2{\rho }_{0}a.$$

(98)

For the representative value $\mathsf{\Xi }=122/\left(15\pi \right)$, this becomes

$$H^2\left(a\right)=-{\displaystyle \frac{4\pi }{3}}{G}_0{\rho }_0+{\displaystyle \frac{976}{45}}{\displaystyle \frac{\hslash G_0^2{\rho }_0}{c^3}}a.$$

(99)

Substituting Equation (98) into Equation (96) gives

$${\psi }^{″}\left(a\right)+{\displaystyle \frac{1}{a}}{\psi }^{\prime }\left(a\right)+\left({\displaystyle \frac{\beta }{a}}-{\displaystyle \frac{\mu }{a^2}}\right)\psi \left(a\right)=0,$$

(100)

where

$$\mu ={\displaystyle \frac{3{\pi }^3{c}^6{\rho }_0}{{\hslash }^2{G}_0}},\beta ={\displaystyle \frac{6{\pi }^3\mathsf{\Xi }{c}^3{\rho }_0}{\hslash }}.$$

(101)

For $\mathsf{\Xi }=122/\left(15\pi \right)$,

$$\beta ={\displaystyle \frac{244}{5}}{\pi }^2{\displaystyle \frac{c^3{\rho }_0}{\hslash }}.$$

(102)

To remove the first-derivative term, we define

$$\psi \left(a\right)=a^{-1/2}\varphi \left(a\right).$$

(103)

Then Equation (100) becomes

$${\varphi }^{″}\left(a\right)+Q\left(a\right)\varphi \left(a\right)=0,$$

(104)

with

$$Q\left(a\right)={\displaystyle \frac{\beta }{a}}-{\displaystyle \frac{\alpha }{a^2}},\alpha =\mu -{\displaystyle \frac{1}{4}}={\displaystyle \frac{3{\pi }^3{c}^6{\rho }_0}{{\hslash }^2{G}_0}}-{\displaystyle \frac{1}{4}}.$$

(105)

Equivalently,

$$Q\left(a\right)={\displaystyle \frac{\beta a-\alpha }{a^2}}.$$

(106)

4.2. Turning Point and WKB Regions

The turning point is determined by

$$Q\left(a_t\right)=0.$$

(107)

For $\alpha >0$ and $\beta >0$, this gives the single positive turning point

$$a_t={\displaystyle \frac{\alpha }{\beta }}.$$

(108)

The sign of $Q\left(a\right)$ then separates minisuperspace into two regions:

$$Q\left(a\right)<0,0<a<a_t,$$

(109)

and

$$Q\left(a\right)>0,a>a_t.$$

(110)

Thus the small-a region is WKB-forbidden, while the large-a region is WKB-allowed. This yields a forbidden small-a region and an allowed oscillatory large-a region.

For this choice of parameters,

$$G_0={\rho }_0=\hslash =c=1,\mathsf{\Xi }={\displaystyle \frac{122}{15\pi }},$$

(111)

one finds

$$\alpha \simeq 92.769,\beta \simeq 481.637,a_t\simeq 0.193.$$

(112)

This is the turning point shown in Figure 7.

4.3. Forbidden-Region and Allowed-Region Solutions

In the forbidden region $0<a<a_t$, define

$$\kappa \left(a\right)=\sqrt{\left|Q\right(a\left)\right|}={\displaystyle \frac{\sqrt{\alpha -\beta a}}{a}}.$$

(113)

The WKB solution for $\varphi \left(a\right)$ is

$${\varphi }_{\mathrm{forb}}\left(a\right)\approx {\left|Q\left(a\right)\right|}^{-1/4}\left[A{e}^{-I\left(a\right)}+B{e}^{+I\left(a\right)}\right],$$

(114)

where

$$I\left(a\right)={\int }_a^{a_t}\kappa \left(a^{\prime }\right)d{a}^{\prime }.$$

(115)

Since $\psi \left(a\right)=a^{-1/2}\varphi \left(a\right)$, the corresponding WKB form for the original wave function is

$${\psi }_{\mathrm{forb}}\left(a\right)\approx {(\alpha -\beta a)}^{-1/4}\left[A{e}^{-I\left(a\right)}+B{e}^{+I\left(a\right)}\right].$$

(116)

The branch proportional to $e^{-I\left(a\right)}$ is suppressed as $a\to 0$. Selecting this branch gives a wave function that is exponentially small in the classically forbidden small-a region.

In the allowed region $a>a_t$, define

$$p\left(a\right)=\sqrt{Q\left(a\right)}={\displaystyle \frac{\sqrt{\beta a-\alpha }}{a}}.$$

(117)

The WKB solution for $\varphi \left(a\right)$ is

$${\varphi }_{\mathrm{all}}\left(a\right)\approx Q{\left(a\right)}^{-1/4}\left[C{e}^{i\mathsf{\Theta }\left(a\right)}+D{e}^{-i\mathsf{\Theta }\left(a\right)}\right],$$

(118)

where

$$\mathsf{\Theta }\left(a\right)={\int }_{a_t}^{a}p\left(a^{\prime }\right)d{a}^{\prime }.$$

(119)

Thus the original wave function is

$${\psi }_{\mathrm{all}}\left(a\right)\approx {(\beta a-\alpha )}^{-1/4}\left[C{e}^{i\mathsf{\Theta }\left(a\right)}+D{e}^{-i\mathsf{\Theta }\left(a\right)}\right].$$

(120)

4.4. Airy Matching Near the Turning Point

The connection between the forbidden and allowed WKB regions is obtained by matching to the local Airy solution near $a=a_t$. Expanding $Q\left(a\right)$ about the turning point gives

$$Q\left(a\right)\simeq Q^{\prime }\left(a_t\right)(a-a_t),$$

(121)

where, using $Q\left(a\right)=\beta /a-\alpha /a^2$,

$$Q^{\prime }\left(a_t\right)={\displaystyle \frac{{\beta }^3}{{\alpha }^2}}>0.$$

(122)

Defining the local Airy scale

$$s={\left[Q^{\prime }\left(a_t\right)\right]}^{-1/3}={\left({\displaystyle \frac{{\alpha }^2}{{\beta }^3}}\right)}^{1/3},$$

(123)

the reduced equation takes the standard Airy form near the turning point. The usual connection formula (see, e.g., [29,30]) gives the real matched solution in the allowed region as

$${\psi }_{\mathrm{all}}\left(a\right)\approx {(\beta a-\alpha )}^{-1/4}cos\left(\mathsf{\Theta }\left(a\right)-{\displaystyle \frac{\pi }{4}}\right).$$

The $-\pi /4$ phase shift is the standard result for matching a decaying exponential on the forbidden side to an oscillatory solution on the allowed side.

This standing-wave form is the standard WKB continuation of a real decaying solution in the forbidden region. A purely outgoing or purely incoming branch may instead be selected by imposing an additional boundary condition in the allowed region.

4.5. Closed-Form WKB Actions and Asymptotic Behavior

For this case, the WKB actions can be evaluated analytically. In the forbidden region,

$$I\left(a\right)=2\sqrt{\alpha }\left[artanh\sqrt{1-{\displaystyle \frac{a}{a_t}}}-\sqrt{1-{\displaystyle \frac{a}{a_t}}}\right].$$

(124)

This expression vanishes at the turning point and diverges logarithmically as $a\to 0$. Consequently, the decaying branch in Equation (116) is strongly suppressed near the classical singularity. In fact,

$$\psi \left(a\right)\sim a^{\sqrt{\alpha }},a\to 0,$$

(125)

up to an overall constant.

In the allowed region,

$$\mathsf{\Theta }\left(a\right)=2\sqrt{\alpha }\left[\sqrt{{\displaystyle \frac{a}{a_t}}-1}-arctan\sqrt{{\displaystyle \frac{a}{a_t}}-1}\right].$$

(126)

For large a,

$$Q\left(a\right)\sim {\displaystyle \frac{\beta }{a}},$$

(127)

and the phase behaves as

$$\mathsf{\Theta }\left(a\right)\sim 2\sqrt{\beta a}-\pi \sqrt{\alpha },a\to \infty .$$

(128)

Therefore the large-a wave function has the asymptotic form

$$\psi \left(a\right)\sim {\left(\beta a\right)}^{-1/4}cos\left(2\sqrt{\beta a}-\pi \sqrt{\alpha }-{\displaystyle \frac{\pi }{4}}\right),a\to \infty .$$

(129)

The large-scale wave function is therefore oscillatory with a slowly decaying envelope. This is the usual WKB signature of a semiclassical regime.

The result is a WKB realization of the usual tunneling-type structure: a suppressed small-a region, a turning point, and an oscillatory large-a semiclassical branch. In the present model, this structure is generated by the scale dependence of the gravitational coupling and the associated $\mathsf{\Xi }$-dependent correction.

4.6. Probability Current and Branch Interpretation

The WKB wave function becomes physically clearer when supplemented by the minisuperspace probability current. For the factor-ordered Wheeler–DeWitt equation

$${\psi }^{″}\left(a\right)+{\displaystyle \frac{q}{a}}{\psi }^{\prime }\left(a\right)+\mathsf{\Omega }\left(a\right)\psi \left(a\right)=0,$$

(130)

with real $\mathsf{\Omega }\left(a\right)$, the conserved current is

$$J_q\left(a\right)={\displaystyle \frac{1}{2i}}{a}^q\left({\psi }^{*}{\psi }^{\prime }-\psi {\psi }^{*\prime }\right).$$

(131)

For the present model $q=1$, this becomes

$$J\left(a\right)={\displaystyle \frac{1}{2i}}a\left({\psi }^{*}{\psi }^{\prime }-\psi {\psi }^{*\prime }\right).$$

(132)

In the allowed region $a>a_t$, the outgoing and incoming WKB branches are

$${\psi }_{\mathrm{out}}\left(a\right)\approx {(\beta a-\alpha )}^{-1/4}{e}^{i\mathsf{\Theta }\left(a\right)},{\psi }_{\mathrm{in}}\left(a\right)\approx {(\beta a-\alpha )}^{-1/4}{e}^{-i\mathsf{\Theta }\left(a\right)}.$$

(133)

Substitution into Equation (132) gives

$$J_{\mathrm{out}}=+1,J_{\mathrm{in}}=-1,$$

(134)

up to the overall normalization of the wave function. Thus the outgoing branch carries positive flux, while the incoming branch carries negative flux. In this sense, the sign of the conserved current distinguishes expanding and contracting semiclassical branches [5].

By contrast, the real standing-wave combination

$${\psi }_{\mathrm{stand}}\left(a\right)\approx 2{(\beta a-\alpha )}^{-1/4}cos\left(\mathsf{\Theta }\left(a\right)-{\displaystyle \frac{\pi }{4}}\right)$$

(135)

has

$$J_{\mathrm{stand}}=0.$$

(136)

Similarly, a purely real decaying solution in the forbidden region carries no current. A nonzero conserved flux therefore appears only after selecting a complex traveling-wave branch in the WKB-allowed region. This is why the current is an invaluable diagnostic: oscillatory behavior alone does not specify a direction of semiclassical evolution, while the current does.

For the illustrative choice of values in Equation (112), one finds numerically that the current is constant across the allowed region:

$$\langle J_{\mathrm{out}}\rangle =1,\langle J_{\mathrm{in}}\rangle =-1,\langle J_{\mathrm{stand}}\rangle =0.$$

(137)

This confirms the branch interpretation of the WKB solutions.

4.7. Hamilton–Jacobi Limit and Classical Recovery

The WKB phase also explains how the classical modified Friedmann dynamics are recovered. Starting from the compact Wheeler–DeWitt equation,

$${\psi }^{″}\left(a\right)+{\displaystyle \frac{q}{a}}{\psi }^{\prime }\left(a\right)+{\displaystyle \frac{9{\pi }^2{c}^6}{4{\hslash }^2}}{\displaystyle \frac{a^4}{G{\left(a\right)}^2}}{H}^2\left(a\right)\psi \left(a\right)=0,$$

(138)

we use the WKB ansatz

$$\psi \left(a\right)=A\left(a\right)exp\left({\displaystyle \frac{i}{\hslash }}S\left(a\right)\right).$$

(139)

At leading order in the WKB expansion, the phase satisfies

$${\left({\displaystyle \frac{dS}{da}}\right)}^2={\displaystyle \frac{9{\pi }^2{c}^6}{4G{\left(a\right)}^2}}{a}^4{H}^2\left(a\right).$$

(140)

Identifying

$$p_a={\displaystyle \frac{dS}{da}},$$

(141)

and using the minisuperspace momentum

$$p_a=-{\displaystyle \frac{3\pi c^3}{2G\left(a\right)}}a\dot{a},$$

(142)

one immediately obtains

$${\left({\displaystyle \frac{dS}{da}}\right)}^2={\displaystyle \frac{9{\pi }^2{c}^6}{4G{\left(a\right)}^2}}{a}^2{\dot{a}}^2={\displaystyle \frac{9{\pi }^2{c}^6}{4G{\left(a\right)}^2}}{a}^4{H}^2\left(a\right).$$

(143)

Thus the WKB Hamilton–Jacobi equation is exactly the classical modified Friedmann dynamics written in canonical form. The phase of the Wheeler–DeWitt wave function therefore encodes the classical cosmological trajectory [14,17].

The factor-ordering term affects the transport equation for the WKB amplitude, but not the leading Hamilton–Jacobi relation. The $\mathsf{\Xi }$-dependent term should also be interpreted carefully: it originates from the quantum-corrected Friedmann equation before minisuperspace quantization. Thus the WKB expansion is an effective semiclassical expansion of an already quantum-corrected cosmological model, rather than a derivation of the $\mathsf{\Xi }$ term from the Wheeler–DeWitt quantization itself.

The physical picture is then consistent across the classical and quantum descriptions. The roots of $H^2\left(a\right)$ determine the classical turning points. The same roots enter the Wheeler–DeWitt potential and determine where the WKB wave function changes from exponential to oscillatory behavior. In the oscillatory regime, the probability current separates expanding and contracting branches. In this way the WKB analysis connects the modified Friedmann dynamics, the minisuperspace wave function, and the emergence of semiclassical cosmological evolution.

5. Wheeler–DeWitt Equation with a Scalar-Field Clock

The Wheeler–DeWitt equation derived above is a constraint equation, so it does not contain an external time variable. To give the quantum theory a relational interpretation, we enlarge the minisuperspace from the single variable a to $(a,\varphi )$, where $\varphi$ is a homogeneous free scalar field. Since the scalar field momentum is conserved, $\varphi$ can be used as an internal clock on any classical branch for which $p_{\varphi }\ne 0$ [17,31].

We begin with the gravitational kinetic term used in the minisuperspace construction,

$$L_{g,\mathrm{kin}}=-{\displaystyle \frac{3\pi c^3}{4G\left(a\right)}}{\displaystyle \frac{a{\dot{a}}^2}{N}},$$

(144)

where the lapse function $N\left(t\right)$ has been restored. For a homogeneous free scalar field, $V\left(\varphi \right)=0$, the corresponding kinetic term is

$$L_{\varphi ,\mathrm{kin}}={\displaystyle \frac{{\pi }^2{a}^3}{N}}{\dot{\varphi }}^2.$$

(145)

The canonical momenta are therefore

$$p_a=-{\displaystyle \frac{3\pi c^3}{2G\left(a\right)}}{\displaystyle \frac{a\dot{a}}{N}},p_{\varphi }={\displaystyle \frac{2{\pi }^2{a}^3}{N}}\dot{\varphi }.$$

(146)

Solving for the velocities and performing the Legendre transform gives the kinetic part of the Hamiltonian constraint,

$${\mathcal{H}}_{\mathrm{kin}}=-{\displaystyle \frac{G\left(a\right)}{3\pi c^{3}a}}{p}_a^2+{\displaystyle \frac{1}{4{\pi }^2{a}^3}}{p}_{\varphi }^2.$$

(147)

The modified Friedmann equation can be written in the compact form

$$H^2\left(a\right)=F\left(a\right),H={\displaystyle \frac{\dot{a}}{a}},$$

(148)

where $F\left(a\right)$ denotes the full right-hand side of the scale-dependent Friedmann equation. The Hamiltonian constraint consistent with this classical limit is

$$\mathcal{H}=-{\displaystyle \frac{G\left(a\right)}{3\pi c^{3}a}}{p}_a^2+{\displaystyle \frac{1}{4{\pi }^2{a}^3}}{p}_{\varphi }^2+{\displaystyle \frac{3\pi c^3}{4G\left(a\right)}}{a}^{3}F\left(a\right)=0.$$

(149)

Hence, if $p_{\varphi }=0$ and $N=1$, then $p_a=-(3\pi c^3/2G)a\dot{a}$, and Equation (149) reduces to ${\dot{a}}^2=a^{2}F\left(a\right)$, as required.

Upon quantization we take

$$p_a\to -i\hslash {\partial }_a,p_{\varphi }\to -i\hslash {\partial }_{\varphi },$$

(150)

and use the same factor ordering as above,

$$p_a^2\to -{\hslash }^2\left({\partial }_a^2+{\displaystyle \frac{q}{a}}{\partial }_a\right).$$

(151)

The Wheeler–DeWitt equation in the two-dimensional minisuperspace is then

$$\left[{\displaystyle \frac{G\left(a\right){\hslash }^2}{3\pi c^{3}a}}\left({\partial }_a^2+{\displaystyle \frac{q}{a}}{\partial }_a\right)-{\displaystyle \frac{{\hslash }^2}{4{\pi }^2{a}^3}}{\partial }_{\varphi }^2+{\displaystyle \frac{3\pi c^3}{4G\left(a\right)}}{a}^{3}F\left(a\right)\right]\Psi (a,\varphi )=0.$$

(152)

Multiplying by $3\pi c^{3}a/\left[G\left(a\right){\hslash }^2\right]$, we obtain the more convenient form

$$\left[{\partial }_a^2+{\displaystyle \frac{q}{a}}{\partial }_a-{\displaystyle \frac{3c^3}{4\pi G\left(a\right)a^2}}{\partial }_{\varphi }^2+{\displaystyle \frac{9{\pi }^2{c}^6}{4{\hslash }^{2}G{\left(a\right)}^2}}{a}^{4}F\left(a\right)\right]\Psi (a,\varphi )=0.$$

(153)

This equation shows the structural change caused by the scalar clock: the scalar-field kinetic term enters with a coefficient proportional to ${\left[G\left(a\right)a^2\right]}^{-1}$. Thus the same running gravitational coupling that modifies $F\left(a\right)$ also modifies the clock term itself.

For the power-law model

$$G\left(a\right)=G_0{a}^n,$$

(154)

the scalar-clock contribution scales as

$$-{\displaystyle \frac{3c^3}{4\pi G\left(a\right)a^2}}{\partial }_{\varphi }^2=-{\displaystyle \frac{3c^3}{4\pi G_0}}{a}^{-(n+2)}{\partial }_{\varphi }^2.$$

(155)

Hence the standard inverse-square clock term is recovered only when $n=0$. For $n>0$, the scalar-clock term becomes more singular as $a\to 0$, while for $n<0$ it is softened.

We now separate variables by writing

$$\Psi (a,\varphi )=e^{i\omega \varphi }{\chi }_{\omega }\left(a\right),$$

(156)

where $\omega$ is the conserved scalar-field frequency. Since ${\partial }_{\varphi }^2\Psi =-{\omega }^2\Psi$, Equation (153) becomes

$${\chi }_{\omega }^{″}\left(a\right)+{\displaystyle \frac{q}{a}}{\chi }_{\omega }^{\prime }\left(a\right)+\left[{\displaystyle \frac{3c^3{\omega }^2}{4\pi G\left(a\right)a^2}}+{\displaystyle \frac{9{\pi }^2{c}^6}{4{\hslash }^{2}G{\left(a\right)}^2}}{a}^{4}F\left(a\right)\right]{\chi }_{\omega }\left(a\right)=0.$$

(157)

For $G\left(a\right)=G_0{a}^n$, this is

$${\chi }_{\omega }^{″}\left(a\right)+{\displaystyle \frac{q}{a}}{\chi }_{\omega }^{\prime }\left(a\right)+\left[{\displaystyle \frac{3c^3{\omega }^2}{4\pi G_0}}{a}^{-(n+2)}+{\displaystyle \frac{9{\pi }^2{c}^6}{4{\hslash }^2{G}_0^2}}{a}^{4-2n}F\left(a\right)\right]{\chi }_{\omega }\left(a\right)=0.$$

(158)

Using the normalized-density Friedmann function,

$$\begin{array}{cc}\hfill F\left(a\right)& =\left({\displaystyle \frac{{\mathsf{\Omega }}_0{H}_0^2}{G_0}}\right)G^{*}\left(a\right)a^{-3\gamma }+\mathsf{\Lambda }\left({\displaystyle \frac{{\mathsf{\Omega }}_0^2{H}_0^4}{2G_0^2{c}^2}}\right)G_1^{*2}\left(a\right)a^{2-6\gamma }\hfill \\ \hfill & -\mathsf{\Xi }\left({\displaystyle \frac{{\mathsf{\Omega }}_0{H}_0^2}{G_0}}\right){\displaystyle \frac{\hslash }{c^3}}{G}_2^{*2}\left(a\right)a^{-(3\gamma +2)}-{\displaystyle \frac{k{c}^2}{a^2}},\hfill \end{array}$$

and the power-law effective couplings

$$G^{*}\left(a\right)={\displaystyle \frac{G_0}{1-n}}{a}^n,G_1^{*2}\left(a\right)={\displaystyle \frac{G_0^2}{1-n}}{a}^{2n},G_2^{*2}\left(a\right)=-{\displaystyle \frac{3G_0^2}{2n-3}}{a}^{2n},$$

(159)

with $n\ne 1$ and $n\ne 3/2$, Equation (158) may be written as

$$\begin{array}{cc}\hfill {\chi }_{\omega }^{″}\left(a\right)& +{\displaystyle \frac{q}{a}}{\chi }_{\omega }^{\prime }\left(a\right)+[{\displaystyle \frac{3c^3{\omega }^2}{4\pi G_0}}{a}^{-(n+2)}+A{a}^{4-n-3\gamma }+B{a}^{6-6\gamma }\hfill \\ \hfill & +C{a}^{2-3\gamma }-D{a}^{2-2n}]{\chi }_{\omega }\left(a\right)=0,\hfill \end{array}$$

where

$$\begin{array}{cccc}\hfill A& ={\displaystyle \frac{9{\pi }^2{c}^6}{4{\hslash }^2{G}_0^2}}{\displaystyle \frac{{\mathsf{\Omega }}_0{H}_0^2}{1-n}},\hfill & \hfill B& =\mathsf{\Lambda }{\displaystyle \frac{9{\pi }^2{c}^4}{8{\hslash }^2{G}_0^2}}{\displaystyle \frac{{\mathsf{\Omega }}_0^2{H}_0^4}{1-n}},\hfill \\ \hfill C& ={\displaystyle \frac{27{\pi }^2{c}^3}{4\hslash G_0}}{\displaystyle \frac{\mathsf{\Xi }{\mathsf{\Omega }}_0{H}_0^2}{2n-3}},\hfill & \hfill D& ={\displaystyle \frac{9{\pi }^2{c}^8}{4{\hslash }^2{G}_0^2}}k.\hfill \end{array}$$

(160)

The coefficient D carries the sign of k. Thus the curvature term should not be interpreted as positive-definite unless the curvature choice has already been fixed.

It is advantageous to remove the first-derivative term by defining

$${\chi }_{\omega }\left(a\right)=a^{-q/2}{u}_{\omega }\left(a\right).$$

(161)

The mode equation becomes

$$u_{\omega }^{″}\left(a\right)+V_{\mathrm{eff}}(a;\omega )u_{\omega }\left(a\right)=0,$$

(162)

with

$$\begin{array}{cc}\hfill V_{\mathrm{eff}}(a;\omega )& ={\displaystyle \frac{3c^3{\omega }^2}{4\pi G_0}}{a}^{-(n+2)}+A{a}^{4-n-3\gamma }+B{a}^{6-6\gamma }+C{a}^{2-3\gamma }-D{a}^{2-2n}-{\displaystyle \frac{q(q-2)}{4a^2}}.\hfill \end{array}$$

(163)

This is the scalar-clock version of the Schrödinger-like Wheeler–DeWitt equation. Note that the new term’s scaling, $a^{-(n+2)}$, is controlled by the same exponent n as the gravitational term. The clock is therefore not an external add-on to the model; it is tied to the scale-dependent gravitational dynamics.

The relational interpretation can also be seen at the classical level. For a free scalar field,

$$\dot{\varphi }={\displaystyle \frac{N}{2{\pi }^2{a}^3}}{p}_{\varphi },$$

(164)

and on a classical branch with $N=1$ and ${\dot{a}}^2=a^{2}F\left(a\right)$,

$${\left({\displaystyle \frac{da}{d\varphi }}\right)}^2={\displaystyle \frac{4{\pi }^4{a}^8}{p_{\varphi }^2}}F\left(a\right).$$

(165)

Thus the roots of $F\left(a\right)=0$ are also turning points of the relational evolution $a\left(\varphi \right)$. This is why the same structure that produces a classical bounce or forbidden small-a region also appears in the scalar-clock Wheeler–DeWitt problem.

5.1. Wave Packets and Classical Recovery

The separated functions ${\chi }_{\omega }\left(a\right)$ are mode functions. A semiclassical cosmological state should instead be built as a superposition of modes,

$$\Psi (a,\varphi )=\int d\omega A\left(\omega \right)e^{i\omega \varphi }{\chi }_{\omega }\left(a\right),$$

(166)

where $A\left(\omega \right)$ is sharply peaked around ${\omega }_0$. The Gaussian packet used below is not a new exact solution of the Wheeler–DeWitt equation, and it does not replace the Bessel, Airy, or WKB solutions discussed earlier. It is a narrow-packet approximation to a superposition of the Wheeler–DeWitt modes, used to describe how a localized state follows a classical trajectory in scalar-field time.

In the oscillatory regime, each mode may be written in WKB form as

$${\chi }_{\omega }\left(a\right)\simeq \mathcal{A}(a,\omega )e^{iS(a,\omega )}.$$

(167)

The phase of the wave packet is then $\omega \varphi +S(a,\omega )$. Stationary phase gives

$${\displaystyle \frac{\partial }{\partial \omega }}{\left[\omega \varphi +S(a,\omega )\right]}_{\omega ={\omega }_0}=0,$$

(168)

or

$$\varphi =-{\left.{\displaystyle \frac{\partial S(a,\omega )}{\partial \omega }}\right|}_{\omega ={\omega }_0}.$$

(169)

This relation defines the trajectory followed by the peak of the wave packet. In the semiclassical regime, that trajectory agrees with the classical relational solution determined by Equation (165).

The factor ordering also determines the natural minisuperspace measure on a constant-$\varphi$ slice,

$$d\mu \left(a\right)=a^{q}da.$$

(170)

Therefore expectation values are defined by

$$\langle f\left(a\right)\rangle \left(\varphi \right)={\displaystyle \frac{{\displaystyle {\int }_0^{\infty }da{a}^{q}f\left(a\right){\left|\Psi (a,\varphi )\right|}^2}}{{\displaystyle {\int }_0^{\infty }da{a}^q{\left|\Psi (a,\varphi )\right|}^2}}}.$$

(171)

In particular,

$$\langle a\rangle \left(\varphi \right)={\displaystyle \frac{{\displaystyle {\int }_0^{\infty }da{a}^{q+1}{\left|\Psi (a,\varphi )\right|}^2}}{{\displaystyle {\int }_0^{\infty }da{a}^q{\left|\Psi (a,\varphi )\right|}^2}}},{\left(\Delta a\right)}^2\left(\varphi \right)=\langle a^2\rangle -{\langle a\rangle }^2.$$

(172)

A state is semiclassical when

$${\displaystyle \frac{\Delta a}{\langle a\rangle }}\ll 1.$$

(173)

To see the classical limit explicitly, suppose the probability density is localized around the stationary-phase trajectory $\overline{a}\left(\varphi \right)$, with $\overline{a}\left(\varphi \right)=a_{\mathrm{cl}}\left(\varphi \right)$, and is well approximated near its peak by

$${\left|\Psi (a,\varphi )\right|}^2\simeq {\mathcal{N}}^2\left(\varphi \right)exp\left[-{\displaystyle \frac{{(a-\overline{a}\left(\varphi \right))}^2}{2{\sigma }_a^2\left(\varphi \right)}}\right],{\sigma }_a\left(\varphi \right)\ll \overline{a}\left(\varphi \right).$$

(174)

Since the packet is assumed to be localized away from $a=0$, the lower limit of the Gaussian integral may be extended to $-\infty$ at leading order. With $a=\overline{a}+x$, expansion of the measure gives

$$a^q={\overline{a}}^q\left[1+q{\displaystyle \frac{x}{\overline{a}}}+{\displaystyle \frac{q(q-1)}{2}}{\displaystyle \frac{x^2}{{\overline{a}}^2}}+\dots \right].$$

(175)

Using the standard Gaussian moments, one obtains

$$\langle a\rangle \left(\varphi \right)=\overline{a}\left(\varphi \right)\left[1+q{\displaystyle \frac{{\sigma }_a^2\left(\varphi \right)}{{\overline{a}}^2\left(\varphi \right)}}+\mathcal{O}\left({\displaystyle \frac{{\sigma }_a^4}{{\overline{a}}^4}}\right)\right],$$

(176)

and

$${\left(\Delta a\right)}^2\left(\varphi \right)={\sigma }_a^2\left(\varphi \right)+\mathcal{O}\left({\displaystyle \frac{{\sigma }_a^4}{{\overline{a}}^2}}\right).$$

(177)

Thus, to leading order,

$$\langle a\rangle \left(\varphi \right)\simeq a_{\mathrm{cl}}\left(\varphi \right),{\displaystyle \frac{\Delta a}{\langle a\rangle }}\simeq {\displaystyle \frac{{\sigma }_a}{\overline{a}}}.$$

(178)

For the Laplace–Beltrami ordering $q=1$, this reduces to

$$\langle a\rangle \simeq \overline{a}+{\displaystyle \frac{{\sigma }_a^2}{\overline{a}}},{\left(\Delta a\right)}^2\simeq {\sigma }_a^2.$$

(179)

The measure-induced shift of the packet center is therefore negligible whenever ${\sigma }_a/\overline{a}\ll 1$. Classical recovery fails near a turning point or near the small-a region, where the packet broadens and the quantum terms are no longer small.

5.2. Hamilton–Jacobi Limit and Probability Current

The same result follows from the Hamilton–Jacobi form of the Wheeler–DeWitt equation. Write

$$\Psi (a,\varphi )=R(a,\varphi )e^{iS(a,\varphi )/\hslash }.$$

(180)

Substitution into Equation (153) and taking the real part gives

$$-{\left({\partial }_{a}S\right)}^2+{\displaystyle \frac{3c^3}{4\pi G\left(a\right)a^2}}{\left({\partial }_{\varphi }S\right)}^2+{\displaystyle \frac{9{\pi }^2{c}^6}{4G{\left(a\right)}^2}}{a}^{4}F\left(a\right)+Q_{\hslash }(a,\varphi )=0,$$

(181)

where

$$Q_{\hslash }={\displaystyle \frac{{\hslash }^2}{R}}\left[{\partial }_a^{2}R+{\displaystyle \frac{q}{a}}{\partial }_{a}R-{\displaystyle \frac{3c^3}{4\pi G\left(a\right)a^2}}{\partial }_{\varphi }^{2}R\right].$$

(182)

When $Q_{\hslash }$ is negligible, Equation (181) reduces to the classical Hamilton–Jacobi equation associated with Equation (149). This is the precise sense in which the modified Friedmann dynamics are recovered from the Wheeler–DeWitt phase in the large-scale, oscillatory regime.

The two-variable Wheeler–DeWitt equation also admits a conserved minisuperspace current. Multiplying Equation (153) by $a^q$ and using the fact that the coefficients are real gives

$${\partial }_a{J}_a+{\partial }_{\varphi }{J}_{\varphi }=0,$$

(183)

with

$$J_a={\displaystyle \frac{a^q}{2i}}\left({\Psi }^{*}{\partial }_a\Psi -\Psi {\partial }_a{\Psi }^{*}\right),$$

(184)

and

$$J_{\varphi }=-{\displaystyle \frac{3c^3{a}^q}{8\pi G\left(a\right)a^{2}i}}\left({\Psi }^{*}{\partial }_{\varphi }\Psi -\Psi {\partial }_{\varphi }{\Psi }^{*}\right).$$

(185)

This generalizes the one-dimensional current used in Section 4. In the WKB regime the sign of $J_a$ separates outgoing and incoming branches, while $J_{\varphi }$ describes the flow with respect to the scalar-field clock [5].

5.3. Interpretation of the Scalar-Clock Figure

Figure 8 summarizes the scalar-clock interpretation. The dashed curve is the classical solution obtained from the same modified Friedmann function $F\left(a\right)$ used in the Wheeler–DeWitt equation. The solid curve and shaded band show the expectation value and spread of a narrow quantum packet. Near the turning point, the packet is broad and the classical trajectory is not a good description of the state. At larger scale factor, the packet narrows relative to its mean value, and $\langle a\rangle \left(\varphi \right)$ follows $a_{\mathrm{cl}}\left(\varphi \right)$. This illustrates the quantum-to-classical transition.

The scalar field does not create a new Friedmann model; it gives the Wheeler–DeWitt wave function a relational time parameter. When $F\left(a\right)$ contains a finite turning point, that same turning point appears in both the classical relational dynamics and the quantum mode structure. In the large-a oscillatory regime, the WKB phase, the conserved current, and the expectation value of a narrow packet all give the same semiclassical picture: the quantum state becomes peaked around the modified Friedmann trajectory. In the small-a region, the packet broadens and the quantum terms become critical, thus the classical description is no longer adequate.

6. Discussion

Allowing the gravitational coupling to depend on the scale factor modifies the first integral of the Newtonian equation of motion through the effective couplings $G^{*}\left(a\right)$, $G_1^{*2}\left(a\right)$, and $G_2^{*2}\left(a\right)$. Because these terms scale differently with a, the running reorganizes the small-scale structure of the modified Friedmann equation rather than simply rescaling Newton’s constant.

For the power-law ansatz $G\left(a\right)=G_0{a}^n$, this reorganization is controlled entirely by the exponent n. In the dust, flat, $\mathsf{\Lambda }=0$ case the matter and quantum-correction terms scale as $a^{n-3}$ and $a^{2n-5}$, respectively. These exponents coincide at $n=2$, where the two terms lock onto the same power of a and no isolated turning point can exist. For $n>2$ the matter term is the more singular contribution near $a=0$ and carries a negative coefficient, so $H^2\left(a\right)<0$ there, while the quantum-correction term restores $H^2>0$ at larger a. The root separating these regions is a genuine bounce. For $3/2<n<2$ the quantum term dominates but its coefficient is positive, so the small-a region remains classically allowed; any root in this window is a recollapse scale, not a bounce. The boundary $n=3/2$ marks the sign change of the quantum-correction coefficient.

A zero of $H^2\left(a\right)$ is therefore not sufficient for singularity avoidance. A true bounce requires the additional sign condition: $H^2<0$ below the root and $H^2>0$ above it. Without this check, a formal turning point may be misidentified. Every root in the parameter survey of Sec. 2.8 is classified on this basis.

The Wheeler–DeWitt analysis is the quantum counterpart of the same statement. The running of $G\left(a\right)$ enters the minisuperspace potential explicitly, and the resulting wave functions reflect the sign structure of the Friedmann equation. For the dust case with $n=2$, the equation reduces to a Bessel problem. With the ordering $q=3/4$, the regular solution satisfies the DeWitt condition $\psi \left(0\right)=0$; with $q=1$ it is proportional to $J_0$ and remains finite but nonzero at the origin. DeWitt suppression and regular finite behavior at $a=0$ are physically distinct. The factor-ordering parameter controls which case applies, while the scale-dependent potential determines the structure of the Wheeler–DeWitt equation itself.

The Airy and radiation-era examples extend this picture beyond Bessel sectors. Their significance is not that they add new qualitative behavior but that they confirm the same three possibilities—DeWitt suppression, regularity, oscillatory large-a behavior—within different matter and running-exponent combinations, all from the same potential structure.

The WKB analysis sharpens the semiclassical reading. For the values of $\gamma =1$, $n=3$, $q=1$, $k=\mathsf{\Lambda }=0$, the reduced Wheeler–DeWitt equation has a single turning point at $a_t\simeq 0.193$ (in the arbitrary units of Eq. (112)). The WKB solution is exponentially suppressed for $a<a_t$ and oscillatory for $a>a_t$. The conserved probability current then distinguishes the two oscillatory branches: the outgoing branch carries positive flux, the incoming branch negative flux, and a real standing-wave carries none. Oscillatory behavior alone does not select a direction of semiclassical evolution; the current does.

The scalar-field clock ties these threads together relationally. In the two-dimensional Wheeler–DeWitt equation (Eq. (153)), the clock term scales as $a^{-(n+2)}$, so the same exponent governing the running of $G\left(a\right)$ also controls the small-a behavior of the relational dynamics.

The scope of this model is limited. It is an effective, homogeneous, isotropic minisuperspace model derived from a quantum-corrected Newtonian cosmology. This construction should therefore be interpreted as an effective minisuperspace quantization motivated by the modified Friedmann dynamics rather than as a derivation from the full Einstein–Hilbert action. It does not address anisotropies, inhomogeneous perturbations, or the full constraint algebra of general relativity. The coefficient $\mathsf{\Xi }$ is treated phenomenologically. The power law $G\left(a\right)=G_0{a}^n$ is a local early-time scaling ansatz; any complete model must account for the transition to the observed, effectively constant G at late times, as sketched in Sec. 2.7. The model instead provides a simplified minisuperspace framework in which the consequences of a running gravitational coupling can be tracked through the modified Friedmann equation to the Wheeler–DeWitt potential and semiclassical limit.

7. Conclusions

Starting from a quantum-corrected Newtonian cosmology with a scale-dependent gravitational coupling, we derived a modified Friedmann equation and corresponding Wheeler–DeWitt equation in which the running of $G\left(a\right)$ changes the small-a structure of both the classical and quantum dynamics. For the power-law form $G\left(a\right)=G_0{a}^n$, the dust, flat, $\mathsf{\Lambda }=0$ model develops two natural analytic boundaries at $n=3/2$ and $n=2$, separating regimes with qualitatively different turning-point and singularity structures.

Minisuperspace quantization leads to a Wheeler–DeWitt equation whose effective potential is determined by the same scale-dependent structure. Exact Bessel and Airy solutions, together with a numerical radiation-era example, exhibit three distinct behaviors near $a=0$: DeWitt suppression ($\psi \left(0\right)=0$), regularity without suppression ($\psi \left(0\right)\ne 0$), and oscillatory semiclassical behavior at large a. The factor-ordering parameter controls which of the first two applies; the model separates these cases, rather than conflating regularity with singularity avoidance.

The WKB analysis for the representative model ($n=3$, $\gamma =1$, $q=1$) gives a forbidden small-a region, a single turning point, and an oscillatory large-a branch. The conserved probability current distinguishes expanding and contracting semiclassical branches, while the Hamilton–Jacobi limit recovers the modified Friedmann dynamics as the leading-order phase of the wave function. Introducing a homogeneous scalar field as an internal clock, wave packets can be constructed and followed in relational time: broad near the turning point, they narrow onto the classical trajectory at larger a. The clock term inherits the same running exponent n, coupling the timekeeping structure directly to the scale-dependent gravitational dynamics.

In the early-universe scaling regime considered here, the running of $G\left(a\right)$ is not just a perturbative correction. This term changes the qualitative structure of both the classical and quantum cosmological equations, and the same scale dependence appears consistently across the Friedmann dynamics, Wheeler–DeWitt potential, WKB turning point, and scalar-clock semiclassical limit.

Several extensions are needed to determine whether this mechanism survives contact with observations. The full parameter space $(n,\gamma ,\mathsf{\Xi },\mathsf{\Lambda },q)$ should be mapped systematically, with the Type I/II/III classification applied throughout, and the pure power law replaced by a transition profile such as Eq. (43) and matched to observational bounds. The phenomenological consequences of the running — its effect on the BBN expansion rate, the growth of linear matter perturbations, and the CMB acoustic scale — are examined in the companion paper [37], which shows that the strong running required for quantum singularity avoidance must cut off well before nucleosynthesis, and that the power-law model is best understood as an early-universe scaling limit rather than a global description of gravity.