Constraints from Solar System tests on a covariant loop quantum black hole

The spherically symmetric exterior geometry of this effective LQG black hole is described as follows Alonso-Bardaji:2021yls ; Alonso-Bardaji:2022ear :

$\displaystyle\mathrm{d}s^{2}=-f\left(r\right)\mathrm{d}t^{2}+\frac{1}{g\left(r\right)f\left(r\right)}\mathrm{d}r^{2}+r^{2}\left(\mathrm{d}\theta^{2}+\sin^{2}\theta\thinspace\mathrm{d}\phi^{2}\right),$

(1)

$\displaystyle f\left(r\right)=1-\frac{2M}{r},\quad g\left(r\right)=1-\frac{r_{0}}{r}.$

(2)

A new length scale $r_{0}$ is introduced as a result of quantum gravity effects:

$\displaystyle r_{0}=2M\frac{{\lambda}^{2}}{1+{\lambda}^{2}}.$

(3)

Here, $\lambda$ is a dimensionless parameter inspired by holonomies, and without loss of generality, we can assume that $\lambda>0$. It is evident that the quantum parameter $r_{0}$ defines a minimum area gap ${r_{0}}^{2}$. $M$ represents the constant of motion, which is associated with the ADM mass as:

$\displaystyle\tilde{M}\equiv M_{\mathrm{ADM}}=M+\frac{r_{0}}{2}.$

(4)

The ADM mass will also be identified as the celestial mass. In the limit $\lambda\to 0$, yielding $r_{0}=0$, the effective LQG geometry described by Eq. (1) regresses to the conventional Schwarzschild geometry of GR.

Before proceeding, we would like to provide some insights into the interior geometry of this effective LQG black hole. Upon incorporating the LQG correction, the interior classical singularity is resolved by a minimal spacelike hypersurface at $r=r_{0}$, resulting in a connected region between a black hole and a white hole. This region is characterized by two external asymptotically flat areas of equal mass. In this scenario, the two-sphere bounce surface characterized by a minimal area of $4\pi r_{0}^{2}$ always hide inside the event horizon, i.e. $r_{0}<2M$. This region is also referred to as the black bounce over the global spacetime structure Moreira:2023cxy . Notice that as the limit $M\to 0$ is approached, resulting in $r_{0}\to 0$, the spacetime geometry reduces to a Minkowski configuration for any value of $\lambda$.

For the sake of convenient calculations throughout the paper, we will express the components of the metric (1) in terms of the ADM mass $\tilde{M}$ and the dimensionless parameter $\tilde{r}_{0}$, which is defined as $\tilde{r}_{0}\equiv{r_{0}}/\tilde{M}={2\lambda^{2}}/\left(1+2\lambda^{2}\right)$:

$\displaystyle f\left(r\right)=1-\frac{2{\tilde{M}}}{r}\left(1-\frac{\tilde{r}_{0}}{2}\right),\quad g\left(r\right)=1-\frac{{\tilde{M}}{\tilde{r}_{0}}}{r}.$

(5)

To be able to reduce to the Newtonian limit, the gravitational constant in this effective LQG black hole can be related to the Newtonian gravitational constant by will2015gravity

$\displaystyle{G}_{\mathrm{N}}=G\left(1-\frac{\tilde{r}_{0}}{2}\right).$

(6)

From now on, we will set ${G}_{\mathrm{N}}=1$ instead of $G=1$, which is more convenient for subsequent calculations. Furthermore, for the remainder of the paper, we will eliminate the use of the tilde for the sake of simplicity.

We commence by considering the Lagrangian governing the motion of a test particle over the effective LQG black hole spacetime:

$\displaystyle\mathcal{L}=\frac{1}{2}mg_{\mu\nu}\dot{x}^{\mu}\dot{x}^{\nu}=\frac{1}{2}mg_{\mu\nu}\frac{\mathrm{d}x^{\mu}}{\mathrm{d}\tau}\frac{\mathrm{d}x^{\nu}}{\mathrm{d}\tau}\ .$

(7)

Here, $m$ is the mass of the test particle, while $\tau$ can be chosen as the proper time or affine parameter for massive or massless particles along geodesics, with the overdot indicating derivative with respect to $\tau$. Then the canonical momentum of the particle can be worked out as:

$\displaystyle p_{\mu}=\frac{\partial\mathcal{L}}{\partial\dot{x}^{\mu}}=mg_{\mu\nu}\dot{x}^{\nu}.$

(8)

Specially, we can explicitly express the four components of the canonical momentums as follows:

	$\displaystyle p_{t}=m\left(-f\left(r\right)\right)\dot{t},$		(9)
	$\displaystyle p_{r}=m\frac{1}{g\left(r\right)f\left(r\right)}\dot{r},$		(10)
	$\displaystyle p_{\theta}=mr^{2}\dot{\theta},$		(11)
	$\displaystyle p_{\phi}=mr^{2}\sin^{2}\theta\thinspace\dot{\phi}.$		(12)

Given that the Lagrangian does not depend on the variables $t$ and $\phi$, namely, ${\partial\mathcal{L}}/{\partial t}=0$ and ${\partial\mathcal{L}}/{\partial\phi}=0$, we have two Killing vectors, $\xi^{a}=\left({\partial}/{\partial t}\right)^{a}$ and $\eta^{a}=\left({\partial}/{\partial\phi}\right)^{a}$, which are associated with the energy ${E}$ and angular momentum ${l}$ of the test particle’s motion, respectively. These quantities are determined by Eqs. (9) and (12):

	$\displaystyle{E}$	$\displaystyle=$	$\displaystyle{-p_{t}}\thinspace{\xi^{t}}=mf\left(r\right)\dot{t},$		(13)
	$\displaystyle{l}$	$\displaystyle=$	$\displaystyle{p_{\phi}}\thinspace{\eta^{\phi}}=mr^{2}\sin^{2}\theta\thinspace\dot{\phi}.$		(14)

By employing the Legendre transformation, we obtain the Hamiltonian $\mathcal{H}$ as follows:

$\displaystyle\mathcal{H}$

$\displaystyle=$

$\displaystyle p_{\mu}\dot{x}^{\mu}-\mathcal{L}=\frac{1}{2m}\thinspace g^{\mu\nu}p_{\mu}p_{\nu}.$

(15)

Then, we can explicitly derive the EOMs for the system. These equations are determined by evaluating the Poisson brackets between the canonical phase space variables and the Hamiltonian:

	$\displaystyle\dot{t}=\left\{t,\mathcal{H}\right\}=-\frac{1}{mf\left(r\right)}\thinspace p_{t},$		(16)
	$\displaystyle\dot{p}_{t}=\left\{p_{t},\mathcal{H}\right\}=0,$		(17)
	$\displaystyle\dot{r}=\left\{r,\mathcal{H}\right\}=\frac{g\left(r\right)f\left(r\right)}{m}\thinspace p_{r},$		(18)
	$\displaystyle\dot{p}_{r}=\left\{p_{r},\mathcal{H}\right\}=-\frac{f^{\prime}\left(r\right)}{2mf^{2}\left(r\right)}\thinspace p_{t}^{2}-\frac{g^{\prime}\left(r\right)f\left(r\right)+g\left(r\right)f^{\prime}\left(r\right)}{2m}\thinspace p_{r}^{2}+\frac{p_{\theta}^{2}}{mr^{3}}+\frac{p_{\phi}^{2}}{mr^{3}\sin^{2}\theta},$		(19)
	$\displaystyle\dot{\theta}=\left\{\theta,\mathcal{H}\right\}=\frac{1}{mr^{2}}\thinspace p_{\theta},$		(20)
	$\displaystyle\dot{p}_{\theta}=\left\{p_{\theta},\mathcal{H}\right\}=\frac{\cos\theta}{mr^{2}\sin^{3}\theta}\thinspace p_{\phi}^{2},$		(21)
	$\displaystyle\dot{\phi}=\{\phi,\mathcal{H}\}=\frac{1}{mr^{2}\sin^{2}\theta}\thinspace p_{\phi},$		(22)
	$\displaystyle\dot{p}_{\phi}=\left\{p_{\phi},\mathcal{H}\right\}=0.$		(23)

Equations (17) and (23) also indicate the conservation of energy and angular momentum for the test particle.

Utilizing these two constants of motion, we can express the reduced Hamiltonian as follows:

$\displaystyle\mathcal{H}=\frac{1}{2m}\left(-\frac{{E}^{2}}{f\left(r\right)}+g\left(r\right)f\left(r\right)\thinspace p_{r}^{2}+\frac{p_{\theta}^{2}}{r^{2}}+\frac{{l}^{2}}{r^{2}\sin^{2}\theta}\right)=\frac{\sigma}{2m}.$

(24)

In the above equation, we have utilized the normalization condition for four-velocity, denoted as $g^{\mu\nu}p_{\mu}p_{\nu}={\sigma}$, where $\sigma$ takes different values depending on the nature of the particles involved. To be more specific, for massless particles, $\sigma=0$, while for massive particles, $\sigma=-m^{2}$.

By employing the variable separation technique, we can identify a third constant of motion as follows:

$\displaystyle\mathcal{K}$

$\displaystyle=$

$\displaystyle p_{\theta}^{2}+\frac{l^{2}}{\sin^{2}\theta}=r^{2}\thinspace\sigma+\frac{r^{2}E^{2}}{f\left(r\right)}-r^{2}g\left(r\right)f\left(r\right)\thinspace p_{r}^{2}.$

(25)

This separation constant $\mathcal{K}$ is commonly referred to as the Carter constant Carter:1968rr .

Next, we derive the equations for $r$ and $\theta$, which are expressed in terms of the aforementioned three constants of motion:

	$\displaystyle\left(\frac{\mathrm{d}\theta}{\mathrm{d}\tau}\right)^{2}=\frac{\tilde{\mathcal{K}}}{\thinspace r^{4}}-\frac{\tilde{l}^{2}}{r^{4}\sin^{2}\theta},$		(26)
	$\displaystyle\left(\frac{\mathrm{d}r}{\mathrm{d}\tau}\right)^{2}=-\frac{\tilde{\mathcal{K}}\thinspace g\left(r\right)f\left(r\right)}{r^{2}}+\tilde{\sigma}\thinspace g\left(r\right)f\left(r\right)+\tilde{E}^{2}\thinspace g\left(r\right).$		(27)

In the above calculations, we have introduced the variables per unit mass as: $\tilde{E}=E/m$, $\tilde{l}=l/m$, $\tilde{\mathcal{K}}=\mathcal{K}/m^{2}$, $\tilde{\sigma}=\sigma/m^{2}$. Once again, the tilde is also dropped for notational simplicity in the following paper. Thus, the geodesic equations for test particles then yield:

	$\displaystyle\frac{\mathrm{d}t}{\mathrm{d}\tau}$	$\displaystyle=$	$\displaystyle\frac{{E}}{f\left(r\right)},$		(28)
	$\displaystyle\left(\frac{\mathrm{d}r}{\mathrm{d}\tau}\right)^{2}$	$\displaystyle=$	$\displaystyle-\frac{{\mathcal{K}}g\left(r\right)f\left(r\right)}{r^{2}}+\sigma g\left(r\right)f\left(r\right)+{E}^{2}g\left(r\right),$		(29)
	$\displaystyle\left(\frac{\mathrm{d}\theta}{\mathrm{d}\tau}\right)^{2}$	$\displaystyle=$	$\displaystyle\frac{{\mathcal{K}}}{r^{4}}-\frac{{l}^{2}}{r^{4}\sin^{2}\theta},$		(30)
	$\displaystyle\frac{\mathrm{d}\phi}{\mathrm{d}\tau}$	$\displaystyle=$	$\displaystyle\frac{{l}}{r^{2}\sin^{2}\theta}.$		(31)

Once we have the above equations available, we can delve deeper into exploring constraints on the quantum parameter through solar system experiments.

Without loss of generality, we can solely focus on the evolution of motion is in the equatorial plane, where $\theta={\pi}/{2}$, then $\dot{\theta}=0$. As a result, the Carter constant in Eq. (25) simplifies to $\mathcal{K}=l^{2}$.

Let us consider a scenario where a light ray originates from infinity, gets deflected by the sun, and then escapes back to infinity. By utilizing Eqs. (29) and (31), we can determine the trajectory of the light ray as follows:

$\displaystyle\frac{\mathrm{d}\phi}{\mathrm{d}r}=\pm\left(\frac{E^{2}r^{4}g\left(r\right)}{l^{2}}-g\left(r\right)f\left(r\right)\thinspace r^{2}\right)^{-\frac{1}{2}}.$

(32)

In this equation, the minus sign corresponds to photons moving inward with decreasing $r$, while the plus sign indicates outward-moving photons with increasing $r$. Subsequently, we will introduce the impact parameter $b$, which represents the perpendicular distance from the straight line of motion to the centerline of the sun that is parallel to it. When the light ray at infinity, i.e. $r\to\infty$ and $\phi\to 0$, we get $\sin{\phi}\approx\phi={b}/{r}$. By solving the Eq. (32) under this scenario, we approximately have:

$\displaystyle r=\pm\frac{l}{E\phi}.$

(33)

It’s easy to infer that $b=l/E$. When the light ray approaches the sun and is affected by the gravitational field, it naturally reaches a turning point known as the closest approach. This point is located at a distance of $s$ from the center of the sun, at which the Eq. (32) vanishes; in other words, $\left({\mathrm{d}r}/{\mathrm{d}\phi}\right)_{r=s}=0$. This leads to the following relationship:

$$b=\left(\frac{s^{2}}{f\left(s\right)}\right)^{\frac{1}{2}}.$$

(34)

Particularly, the magnitude of the total change in the coordinate interval $\phi$ is just twice the change in angle from the turning point $r=s$ to infinity. Therefore, we can express it in terms of $b$ as follows:

$\displaystyle\phi=2\int_{s}^{\infty}\left(\frac{r^{4}g\left(r\right)}{b^{2}}-g\left(r\right)f\left(r\right)r^{2}\right)^{-\frac{1}{2}}\mathrm{d}r.$

(35)

In the absence of the sun, the light ray propagates in a straight line, where $\phi=\pi$. Consequently, the deflection angle $\Delta{\phi}$ is related to $\phi$ as: $\Delta\phi=\phi-\pi$. To proceed, we introduce a coordinate transformation $r={s}/{x}$ and define $\epsilon\equiv M/s$. And then, we can expand the equation of $\Delta{\phi}$ in the weak field approximation as follows:

$\displaystyle\Delta\phi=2\int_{0}^{1}\left\{\frac{1}{\sqrt{1-x^{2}}}+\frac{\left[-r_{0}+2\left(1+x+x^{2}\right)\right]\epsilon}{\left(2-r_{0}\right)\left(1+x\right)\sqrt{1-x^{2}}}+\mathcal{O}\left(\epsilon^{2}\right)\right\}\mathrm{d}x-\pi.$

(36)

Therefore, we can obtain the following approximate expression for the light deflection angle with the quantum-corrected term:

$\displaystyle\Delta\phi\approx\frac{4M}{s}\left(1+\frac{r_{0}}{4-2r_{0}}\right)=\Delta\phi_{\rm{GR}}\left(1+\frac{r_{0}}{4-2r_{0}}\right).$

(37)

Here, $\Delta\phi_{\mathrm{GR}}$ represents the deflection value in GR. Notice that in the above calculation, we have reverted $\epsilon$ back to $M/s$.

When the light ray just grazes the sun, we assume that the closest approach $s$ is equal to the radius of the sun $R_{\odot}$, and $M$ is the solar mass $M_{\odot}$. In this scenario, the parameterized post-Newtonian (PPN) formalism equation for light deflection is given as follows will2018theory :

$\displaystyle\Delta\phi\simeq 1.75^{\prime\prime}\left(\frac{1+\gamma}{2}\right),$

(38)

where $\gamma$ is the PPN deflection parameter robertson1962space ; eddington1923mathematical . According to the astrometric observation measuring $\gamma$ by the Very Long Baseline Array (VLBA) Fomalont:2009zg , we compare Eq. (37) with Eq. (38). Consequently, the constraint on the quantum-corrected parameter $r_{0}$ can be immediately determined as follows:

$\displaystyle 0<r_{0}<2.0\times 10^{-4}.$

(39)

This leads to corresponding constraints on $\lambda$ where

$\displaystyle 0<\lambda<1.0\times 10^{-2}.$

(40)

In this subsection, we will analyze the constraints on the LQG-corrected parameter through the study of the Shapiro time delay. Specially, we will consider a simplified scenario in which a radar signal is transmitted from a transmitter located on earth, denoted as $r=A$, then it passes through the closest approach to the sun at the turning point $r=s$, and finally returns to earth by reflection from a spacecraft-mounted reflector, denoted as $r=B$.

By combining Eqs. (28) and (29), we can derive the differential equation for massless particles between $t$ and $r$:

$\displaystyle\frac{\mathrm{d}t}{\mathrm{d}r}=\pm\frac{r}{f\left(r\right)\sqrt{g\left(r\right)\left(r^{2}-b^{2}f\left(r\right)\right)}},$

(41)

where the plus sign and the minus sign correspond to the outgoing and incoming radar waves, respectively.

Then, we can determine the travel time of the radar wave propagating from the transmitter at point $A$ to the turning point at point $s$:

$\displaystyle\Delta t_{A}=-\int_{A}^{s}\frac{r}{f\left(r\right)\sqrt{g\left(r\right)\left(r^{2}-b^{2}f\left(r\right)\right)}}\mathrm{d}r.$

(42)

To evaluate the integral $\Delta t_{A}$, we once again use the coordinate transformation $r=s/x$ and expand the integral using the small quantity $\epsilon\equiv M/s$. By integrating to the subleading order, we obtain:

$\displaystyle\Delta t_{A}$

$\displaystyle\approx$

$\displaystyle\sqrt{A^{2}-s^{2}}+M\sqrt{\frac{A-s}{A+s}}+M\left(\frac{4-r_{0}}{2-r_{0}}\right)\tanh^{-1}\left(\sqrt{1-\frac{s^{2}}{A^{2}}}\right).$

(43)

Similarly, the time it takes for the radar wave to travel between the turning point at $s$ and the reflector at $B$ can be determined in the same manner, and we will refer to it as $\Delta t_{B}$.

Based on the position of the spacecraft carrying the reflector, we usually classify it into two scenarios: inferior conjunction and superior conjunction when calculating the gravitational time delay. In the case of inferior conjunction, the spacecraft is situated between the earth and the sun. In the absence of gravitational effects, the total roundtrip time of the radar signal can be expressed as $\Delta t_{\mathrm{I-S}}=2\left(\sqrt{A^{2}-s^{2}}-\sqrt{B^{2}-s^{2}}\right)$. When considering gravitational effects with LQG corrections, the roundtrip time delay is calculated as follows:

$\displaystyle\Delta t_{\mathrm{I}-r_{0}}\approx 4M\left(\frac{4-r_{0}}{4-2r_{0}}\right)\ln\left(\frac{A}{B}\right)=\Delta t_{\mathrm{I-GR}}\left(\frac{4-r_{0}}{4-2r_{0}}\right),$

(44)

where $\Delta t_{\rm{I-GR}}$ represents the Shapiro time delay in GR. It is evident that the roundtrip time delay receives the LQG corrections.

In the case of superior conjunction, the spacecraft is located at the opposite side of the earth with respect to the sun. Therefore, in the absence of gravitational effects, the total roundtrip time for the radar wave to travel is $\Delta t_{\mathrm{S-S}}=2\left(\sqrt{A^{2}-s^{2}}+\sqrt{B^{2}-s^{2}}\right)$. Then, the roundtrip time delay with LQG corrections is calculated as follows:

$\displaystyle\Delta t_{\mathrm{S}-r_{0}}\approx 4M\left[1+\left(1+\frac{r_{0}}{4-2r_{0}}\right)\ln\left(\frac{4AB}{s^{2}}\right)\right].$

(45)

In fact, the above equation is reformulated in the PPN-like formalism of the Shapiro time delay, which is will2018theory ; Weinberg:1972kfs :

$\displaystyle\Delta t\simeq 4M\left[1+\left(\frac{1+\gamma}{2}\right)\ln\left(\frac{4AB}{s^{2}}\right)\right].$

(46)

It is equivalent to the relation between $r_{0}$ and $\gamma$ given in Sec. III.1.

Next, we will use the Cassini solar conjunction mission in 2002 iess2003cassini ; bertotti2003test to constrain LQG-corrected parameter. By utilizing a multifrequency link in the X and Ka bands to minimize the influence of solar corona noise, significant improvements have been achieved, resulting in $\gamma=1+\left(2.1\pm 2.3\right)\times 10^{-5}$ bertotti2003test ; iess1999doppler ; Deng:2015sua ; Deng:2017hkj . Then we can give rise to the upper constraint on $r_{0}$ as:

$\displaystyle 0<r_{0}<8.80\times 10^{-5}.$

(47)

The corresponding constraint on $\lambda$ is as follows:

$\displaystyle 0<\lambda<6.63\times 10^{-3}.$

(48)

On the other hand, we can also constrain the LQG parameter using the Doppler tracking of the Cassini spacecraft bertotti1992relativistic ; bertotti1993doppler . In contrast to the Shapiro time delay, which measures the time delay, the Doppler tracking directly measures the relative frequency variation. To achieve this, we differentiate Eq. (45) with respect to time $t$, leading to the fractional frequency shift for the radar signal iess1999doppler ; Deng:2015sua :

$\displaystyle\frac{\mathrm{d}\Delta t_{\mathrm{S}-r_{0}}}{\mathrm{d}t}\equiv\delta\nu=\frac{\nu\left(t\right)-\nu_{0}}{\nu_{0}}\approx\left[-\frac{8M}{s}-\frac{4Mr_{0}}{\left(2-r_{0}\right)s}\right]\frac{\mathrm{d}s(t)}{\mathrm{d}t}\approx\delta\nu_{\mathrm{GR}}-\frac{4Mr_{0}}{\left(2-r_{0}\right)s}\nu_{\oplus},$

(49)

where $\nu_{0}$ and $\nu(t)$ are the emitted and received frequencies, respectively, and ${\mathrm{d}s(t)}/{\mathrm{d}t}$ is approximately equivalent to the average orbit velocity of earth, denoted as $\nu_{\oplus}$. Therefore, the frequency shift caused by the quantum correction parameter $r_{0}$ can be expressed as Deng:2017hkj ; Zhu:2020tcf ; Liu:2022qiz :

$\displaystyle\left|\delta\nu_{\mathrm{S}-r_{0}}\right|\approx\frac{4Mr_{0}}{\left(2-r_{0}\right)s}\nu_{\oplus}=\frac{4M_{\odot}r_{0}}{\left(2-r_{0}\right)R_{\odot}}\frac{16}{27}\nu_{\oplus}<10^{-14},$

(50)

where $M_{\odot}$ and $R_{\odot}$ respectively denote the mass and radius of the sun. By taking these conditions into account, we obtain an upper constraint on $r_{0}$ within the range of $0<r_{0}<4.0\times 10^{-5}$, along with corresponding constraints on $\lambda$ falling in the range of $0<\lambda<4.47\times 10^{-3}$. It is evident that both the Shapiro time delay and the Doppler tracking of the Cassini spacecraft yield consistent results.

In this subsection, we investigate the constraints on the LQG-corrected parameter through the study of perihelion precession. To do so, we analyze the motion of a massive particle $(\sigma=-1)$ orbiting the sun.

First, by combining Eq. (29) and Eq. (31), we can derive the equation describing the orbit precession of the massive particle as follows:

$\displaystyle\left(\frac{\mathrm{d}x}{\mathrm{d}\phi}\right)^{2}=\left[\frac{s^{2}}{l^{2}}\left(1-f\left(\frac{s}{x}\right)\right)+\frac{s^{2}}{b^{2}}-f\left(\frac{s}{x}\right)x^{2}\right]g\left(\frac{s}{x}\right),$

(51)

where we have also introduced the coordinate transformation $r=s/x$ as previously utilized in Sec. III.1. In addition, the impact parameter relates both $l$ and $E$ as $b=l/\sqrt{E^{2}-1}$ in the case of massive particle.

When we differentiate the equation above with respect to $\phi$ and then expand the formula within the weak field limit, defined in terms of the small parameter $\epsilon$ where $\epsilon\equiv M/s$, we obtain an approximate description of the revolution of the orbits, given by:

$\displaystyle\frac{\mathrm{d}^{2}x}{\mathrm{d}\phi^{2}}+x-\frac{M^{2}}{l^{2}\epsilon}\approx 3\epsilon x^{2}+\frac{r_{0}\left[b^{2}\epsilon^{2}x^{2}l^{2}\left(3-8\epsilon x\right)-M^{2}\left(l^{2}+4b^{2}\epsilon x\right)\right]}{b^{2}l^{2}\left(2-r_{0}\right)\epsilon}.$

(52)

We can clearly observe that the LQG effect is distinctly manifested in the second term on the right-hand side of Eq. (52). In what follows, we will solve the above differential equation using the perturbation methods.

To do this, we begin by expressing $x\left(\phi\right)$ as $x\left(\phi\right)=x_{0}\left(\phi\right)+x_{1}\left(\phi\right)$, with the condition that $x_{1}\left(\phi\right)\ll x_{0}\left(\phi\right)$. When the left-hand side of Eq. (52) equals zero, the situation reverts to the Newtonian gravity theory. In this case, the solution reads

$\displaystyle x_{0}\left(\phi\right)=\frac{M^{2}}{l^{2}\epsilon}\left(1+e\cos\phi\right),$

(53)

which is the unperturbed part and is commonly known as the conic section formula with eccentricity $e$ involved.

Next, we will determine the perturbation part $x_{1}\left(\phi\right)$. To this end, we substitute the expression $x\left(\phi\right)=x_{0}\left(\phi\right)+x_{1}\left(\phi\right)$, with $x_{0}\left(x\right)$ obtained in Eq. (53), into Eq. (52), while considering the initial conditions $x_{1}\left(0\right)=0$, ${\mathrm{d}x_{1}\left(0\right)}/{\mathrm{d}\phi}=0$. This yields the following equation:

$\displaystyle\frac{\mathrm{d}^{2}x_{1}}{\mathrm{d}\phi^{2}}+x_{1}=\sum_{i=0}^{3}\chi_{i}\cos^{i}\phi,$

(54)

where

	$\displaystyle\chi_{0}$	$\displaystyle=$	$\displaystyle\frac{b^{2}\left[-8M^{6}r_{0}+2l^{2}M^{4}(3-2r_{0})\right]-l^{6}M^{2}r_{0}}{b^{2}l^{6}(2-r_{0})\epsilon},$		(55)
	$\displaystyle\chi_{1}$	$\displaystyle=$	$\displaystyle\frac{4M^{4}e\left[l^{2}(3-r_{0})-6M^{2}r_{0}\right]}{l^{6}(2-r_{0})\epsilon},$		(56)
	$\displaystyle\chi_{2}$	$\displaystyle=$	$\displaystyle\frac{6M^{4}e^{2}\left(l^{2}-4M^{2}r_{0}\right)}{l^{6}(2-r_{0})\epsilon},$		(57)
	$\displaystyle\chi_{3}$	$\displaystyle=$	$\displaystyle-\frac{8M^{6}e^{3}r_{0}}{l^{6}(2-r_{0})\epsilon}.$		(58)

We can therefore obtain the solution:

	$\displaystyle x_{1}(\phi)$	$\displaystyle=$	$\displaystyle\chi_{0}+\frac{\chi_{2}}{2}-\chi_{0}\cos\phi-\frac{\chi_{2}}{3}\cos\phi+\frac{\chi_{3}}{32}\cos\phi-\frac{\chi_{2}}{6}\cos\left(2\phi\right)-\frac{\chi_{3}}{32}\cos\left(3\phi\right)$		(59)
		$\displaystyle+$	$\displaystyle\frac{\chi_{1}}{2}\phi\sin\phi+\frac{3\chi_{3}}{8}\phi\sin\phi.$

Regarding the perihelion precession, when the terms involving $\phi\sin\phi$ in Eq. (59) are absent, the test particle remains on a closed orbit without deflection. As time progresses, the cumulative effect makes the perihelion precession of planetary orbits observable. Therefore, in this scenario, the remaining terms in Eq. (59) can be omitted. Finally, we obtain the approximate solution to Eq. (52) as follows:

$\displaystyle x\left(\phi\right)\approx\frac{Ms}{l^{2}}\left(1+e\cos\phi\right)+\left(\frac{\chi_{1}}{2}+\frac{3\chi_{3}}{8}\right)\phi\sin\phi\approx\frac{Ms}{l^{2}}\left(1+e\cos\left(\phi-\phi_{0}\right)\right),$

(60)

where we have transformed $\epsilon$ back into $M/s$, and now we can relate the precession angle $\delta\phi_{0}$ using the expression $\phi_{0}=\left(\delta\phi_{0}/2\pi\right)\phi$ as follows:

$\displaystyle\delta{\phi_{0}}\approx\frac{4\pi M^{2}}{l^{2}}\left(\frac{3-r_{0}}{2-r_{0}}\right).$

(61)

In particular, the radial distance $r$ attains its minimum value at perihelia, where the condition $\phi-\phi_{0}=0$ yields the equation $s/r_{-}={Ms\left(1+e\right)}/{l^{2}}$. On the other hand, the radial distance achieves its maximum value at aphelia, where the condition $\phi-\phi_{0}=\pi$ results in the equation $s/r_{+}={Ms\left(1-e\right)}/{l^{2}}$. While for any bound orbit, we can determine the semimajor axis $a$ using the following formula:

$\displaystyle a=\frac{r_{-}+r_{+}}{2}=\frac{l^{2}}{M(1-e^{2})}.$

(62)

Combining Eqs. (61) and (62), we obtain the angle of perihelion precession per revolution deviated from the GR prediction:

$\displaystyle\Delta{\phi}=\delta{\phi_{0}}\approx\frac{6\pi M}{a\left(1-e^{2}\right)}\left(\frac{6-2r_{0}}{6-3r_{0}}\right)=\Delta{\phi}_{\rm{GR}}\left(1+\frac{r_{0}}{6-3r_{0}}\right).$

(63)

Additionally, we can express the $\Delta{\phi}_{\rm{GR}}$ in terms of the solar mass $M_{\odot}$ as $\Delta{\phi}_{\rm{GR}}={6\pi M_{\odot}}/{\left[a\left(1-e^{2}\right)\right]}$.

For the observation of Mercury’s anomalous perihelion advance, the MESSENGER mission provided highly accurate measurements park2017precession , yielding a value of $\Delta{\phi}=\left(42.9799\pm 0.0009\right)^{\prime\prime}$ per century. Using this measured data, we can establish upper bounds on the parameters $r_{0}$ and $\lambda$, resulting in the following constraints:

$\displaystyle 0<r_{0}<1.26\times 10^{-4},\quad 0<\lambda<7.93\times 10^{-3}.$

(64)

This result aligns with the expectations that the contribution from LQG effect is less than the observational error $0.0009^{\prime\prime}$ per century. In Appendix V, we provide further discussions on this topic using the PPN method.

It should be noted that Eq. (63) only considers the gravitoelectric perihelion shift resulting from the influence of the solar mass $M_{\odot}$, whereas the gravitomagnetic component, commonly referred to as the Lense-Thirring effect, is not taken into account within in this equation lense1918einfluss . Based on the summary provided in park2017precession , it is observed that the level of uncertainty associated with the total precession rate is comparatively smaller than the estimated contributions attributed to the Lense-Thirring effect over a century. Consequently, the measurement of the Lense-Thirring effect has not yet attained a commensurate level of precision. As a result, for the purposes of this paper, it is deemed appropriate to disregard this effect and consider the central object as nonrotating.

Hence, the LAGEOS satellites are taken into consideration due to their ability to yield precise outcomes through the measurement of the relativistic precession of LAGEOS II’s pericenter within the earth’s orbit Lucchesi:2010zzb . Based on the analysis of tracking data spanning a period of 13 years, the PPN level factor $\epsilon_{\omega}$ is estimated to be $\epsilon_{\omega}=1+\left(0.28\pm 2.14\right)\times 10^{-3}$. When $\epsilon_{\omega}=1$, the situation reverts back to the case of GR. When compared to the excessive coefficient in Eq. (63), it results in the constraints on $r_{0}$ and $\lambda$ as follows:

$\displaystyle 0<r_{0}<1.44\times 10^{-2},\quad 0<\lambda<8.55\times 10^{-2}.$

(65)

Additionally, the observations of star S2 orbit around $\mathrm{Sgr\ A^{*}}$, the closest massive black hole candidate at the centre of the Milky Way, provide an alternative means of testing the Schwarzchild precession (SP) GRAVITY:2020gka . GR predicts a precession advance angle of $\Delta{\phi}_{\rm{GR}}=12.1^{\prime}$ per orbital period. The data analysis by the GRAVITY collaboration yields a PPN-like parameter $f_{\rm{SP}}$. In the context of GR, it is anticipated that this parameter would have a value of $1$. Nevertheless, a fiducial value with uncertainty $f_{\rm{SP}}=1.10\pm 0.19$ was determined. Comparing these results with Eq. (63), the corresponding upper bounds on $r_{0}$ and $\lambda$ are as follows:

$\displaystyle 0<r_{0}<0.93,\quad 0<\lambda<2.59.$

(66)

This is in agree with the result given in Balali:2023ccr , where $\lambda\in\left[0,2.65\right]$ is obtained by using the $\mathrm{Sgr\ A^{*}}$ shadow’s angular diameter.

In this subsection, our attention turns to another test of GR known as geodetic precession. This test serves the purpose of probing the spacetime geometry and place constraints on the LQG parameter.

We commence by studying the motion of the spin of a point test particle in free fall Schiff:1960gi ; hartle2003gravity . We assume that the test particle moves along a timelike geodesic, whose four-velocity vectors $u^{\mu}=\dot{x}^{\mu}={\mathrm{d}x^{\mu}}/{\mathrm{d}\tau}$, governed by the geodesic equation:

$\displaystyle\frac{\mathrm{d}u^{\mu}}{\mathrm{d}\tau}+\Gamma_{\ \nu\alpha}^{\mu}u^{\nu}u^{\alpha}$

$\displaystyle=$

$\displaystyle 0,$

(67)

where $\Gamma_{\enspace\nu\alpha}^{\mu}$ represents the four-dimensional Christoffel symbol. The evolution of its spin four-vector $S^{\mu}$ along the geodesic is described as follows:

$\displaystyle\frac{\mathrm{d}S^{\mu}}{\mathrm{d}\tau}+\Gamma_{\ \nu\alpha}^{\mu}S^{\nu}u^{\alpha}=0.$

(68)

This equation is commonly referred to as the gyroscope equation or parallel transport equation walker1935note ; synge1960relativity ; Misner:1973prb . In addition, we will use the orthogonality and normalization conditions, which are expressed as:

	$\displaystyle u^{\mu}S_{\mu}$	$\displaystyle=$	$\displaystyle 0,$		(69)
	$\displaystyle S^{\mu}S_{\mu}$	$\displaystyle=$	$\displaystyle 1.$		(70)

To simplify the calculation, we assume that the trajectory of the test particle follows a circular orbit and is confined to the equatorial plane, where $\theta={\pi}/{2}$. Here, we introduce the effective potential $V_{\mathrm{ep}}$ to analyze the stability of the test particle’s orbit. From Eq. (29), we can derive:

$\displaystyle\dot{r}^{2}+V_{\mathrm{ep}}=E^{2},$

(71)

where the effective potential is

$\displaystyle V_{\mathrm{ep}}=E^{2}-\left(-1+\frac{E^{2}}{f\left(r\right)}-\frac{l^{2}}{r^{2}}\right)f\left(r\right)g\left(r\right).$

(72)

To have a stable circular orbit in the equatorial plane, both the radial velocity and radial acceleration need to be zero simultaneously, which means $V_{\mathrm{ep}}=E^{2}$ and $\mathrm{d}V_{\mathrm{ep}}/{\mathrm{d}r}=0$. These conditions yield:

$\displaystyle{E}=\left(\frac{2f^{2}\left(r\right)}{2f\left(r\right)-f^{\prime}\left(r\right)r}\right)^{\frac{1}{2}},\quad{l}=\left(\frac{r^{3}f^{\prime}\left(r\right)}{2f\left(r\right)-f^{\prime}\left(r\right)r}\right)^{\frac{1}{2}}.$

(73)

Hence, the related four-velocity vectors can be recast as:

$\displaystyle u^{t}=\left(\frac{2}{2f\left(r\right)-f^{\prime}\left(r\right)r}\right)^{\frac{1}{2}},\quad u^{\phi}=\left(\frac{f^{\prime}\left(r\right)}{2rf\left(r\right)-f^{\prime}\left(r\right)r^{2}}\right)^{\frac{1}{2}}.$

(74)

By definition, we find that:

$\displaystyle\Omega\equiv\frac{\mathrm{d}\phi}{\mathrm{d}t}=\frac{u^{\phi}}{u^{t}}=\left(\frac{f^{\prime}\left(r\right)}{2r}\right)^{\frac{1}{2}},$

(75)

where $\Omega$ is the orbital angular velocity of the test particle.

Based on these results, we can form the parallel transport equations as follows:

	$\displaystyle\frac{\mathrm{d}S^{t}}{\mathrm{d}\tau}+\frac{1}{2}\frac{f^{\prime}\left(r\right)}{f\left(r\right)}S^{r}u^{t}=0,$		(76)
	$\displaystyle\frac{\mathrm{d}S^{r}}{\mathrm{d}\tau}+\frac{1}{2}{f\left(r\right)g\left(r\right)}{f^{\prime}\left(r\right)}S^{t}u^{t}-{rf\left(r\right)g\left(r\right)}S^{\phi}u^{\phi}=0,$		(77)
	$\displaystyle\frac{\mathrm{d}S^{\theta}}{\mathrm{d}\tau}=0,$		(78)
	$\displaystyle\frac{\mathrm{d}S^{\phi}}{\mathrm{d}\tau}+\frac{1}{r}S^{r}u^{\phi}=0.$		(79)

For convenience, we substitute the derivatives with respect to coordinate time $t$ for derivatives with respect to proper time $\tau$. It is worth noting that $\mathrm{d}\tau=\mathrm{d}t/u^{t}$, and this leads to the following expressions for the spin vectors:

	$\displaystyle S^{t}\left(t\right)=-\frac{f^{\prime}\left(r\right)}{2\omega}\sqrt{\frac{g\left(r\right)}{f\left(r\right)}}\sin\left(\omega t\right),$		(80)
	$\displaystyle S^{r}\left(t\right)=\sqrt{f\left(r\right)g\left(r\right)}\cos\left(\omega t\right),$		(81)
	$\displaystyle S^{\theta}\left(t\right)=0,$		(82)
	$\displaystyle S^{\phi}\left(t\right)=-\frac{\Omega}{r\omega}\sqrt{f\left(r\right)g\left(r\right)}\sin\left(\omega t\right).$		(83)

where

$\displaystyle\omega=\Omega\left(f\left(r\right)g\left(r\right)-\frac{r}{2}g\left(r\right)f^{\prime}\left(r\right)\right)^{\frac{1}{2}},$

(84)

is the angular velocity of the spin vector. Here we have assumed that the spin vector is radial directed at $t=0$, i.e. $S^{t}\left(0\right)=S^{\theta}\left(0\right)=S^{\phi}\left(0\right)=0$. The coefficients can be gained by the conditions (69) and (70). Hence, we can utilize the discrepancy between $\Omega$ and $\omega$ to detect the geodetic effect.

Upon completing one orbit, where $\phi$ changes from 0 to $2\pi$, the corresponding coordinate time is $\delta{t}=2\pi/\Omega$. Consequently, the geodetic precession angle per revolution can be expressed as:

$\displaystyle\Delta{\Phi}_{\mathrm{geo}}=2\pi-\omega\delta{t}=2\pi\left(1-\frac{\omega}{\Omega}\right).$

(85)

To obtain experimental constraints, substituting Eq. (84) into Eq. (85) and expanding it as power series in terms of $M/r$ up to first order gives:

$\displaystyle\Delta{\Phi}_{\mathrm{geo}}\approx\frac{3\pi M}{r}\left(1+\frac{2r_{0}}{6-3r_{0}}\right)=\Delta{\Phi}_{\mathrm{GR}}\left(1+\frac{2r_{0}}{6-3r_{0}}\right).$

(86)

We can conclude that the geodetic precession angle increases with the LQG parameter when $r_{0}>0$.

Detecting such phenomena can be quite challenging. Fortunately, the Gravity Probe B (GP-B) mission, equipped with four nearly perfect spherical gyroscopes and a star-tracking telescope, operates on a polar orbit around earth at an altitude of 642 km Everitt:2011hp . GP-B measures the geodetic drift rate in the north-south direction, with the GR prediction being $-6066.1$ milliarcseconds per year. The analysis of data from the four gyroscopes reveals a geodetic drift rate of $\Delta{\Phi}=\left(-6066.8\pm 18.3\right)$ milliarcseconds per year. This measurement provides bounds for $r_{0}$ and $\lambda$:

$\displaystyle 0<r_{0}<6.34\times 10^{-3},\quad 0<\lambda<5.65\times 10^{-2}.$

(87)

Additionally, lunar laser ranging (LLR) has proven to be one of the most powerful tools for rigorously testing GR theory with high level of precision muller2019lunar . The Earth-Moon system can be deemed as a gyroscope moving around the sun, the geodetic precession is manifested through the change of the lunar orbit which has already reach a level that can be observed by laser ranging. Thus, by measuring the lunar orbit within the Earth-Moon system’s dynamic in the weak field of the sun, LLR provides a relative deviation of geodetic precession from GR value, yielding $K_{\mathrm{gp}}=-0.0019\pm 0.0064$. Based on this result, we therefore obtain upper limits for the parameter $r_{0}$ and $\lambda$ as:

$\displaystyle 0<r_{0}<1.34\times 10^{-2},\quad 0<\lambda<8.24\times 10^{-2}.$

(88)