# Rainbow Black Hole From Quantum Gravitational Collapse

Aliasghar Parvizi<EMAIL_ADDRESS>Department of Physics, University of
Tehran, 14395-547 Tehran, Iran  School of Particles and Accelerators,
Institute for Research in Fundamental Sciences (IPM)
P.O. Box 19395-5531, Tehran, Iran  Tomasz Pawłowski
<EMAIL_ADDRESS>Institute for Theoretical Physics, Faculty of
Physics and Astronomy, University of Wrocław, pl. M. Borna 9, 50-204 Wrocław,
Poland  Yaser Tavakoli<EMAIL_ADDRESS>Department of Physics,
University of Guilan, 41335-1914 Rasht, Iran  Jerzy Lewandowski
<EMAIL_ADDRESS>Faculty of Physics, University of Warsaw,
Pasteura 5, 02-093 Warsaw, Poland

###### Abstract

Quantum evolution of a scalar field’s modes propagating on quantum spacetime
of a collapsing homogeneous dust ball is written effectively, as an evolution
of the same quantum modes on a (semiclassical) dressed geometry. When the
backreaction of the field is discarded, the classical spacetime singularity is
resolved due to quantum gravity effects and is replaced by a quantum bounce on
the dressed collapse background. In the presence of backreaction, the emergent
(interior) dressed geometry becomes mode dependent and the energy density
associated with the backreaction of each mode scales as a radiation fluid.
Semiclassical dynamics of this so-called rainbow dressed background is
analyzed. It turns out that the backreaction effects speed up the occurrence
of the bounce in comparison to the case where only a dust fluid is present. By
matching the interior and exterior regions at the boundary of dust, a mode-
dependent black hole geometry emerges as the exterior spacetime. Properties of
such a rainbow black hole are discussed. That mode dependence causes, in
particular, a chromatic aberration in the gravitational lensing process of
which maximal magnitude is estimated via calculation of the so-called Einstein
angle.

###### pacs:

04.60.-m, 04.60.Pp, 98.80.Qc, 04.60.Bc

##  I Introduction

In contrast to classical general relativity (GR), quantum gravity models based
on discrete spacetime structure often predict that the local Lorentz
invariance may be modified or broken at sufficiently high energy [1, 2]. This
leads, in particular, to the deformation of dispersion relations for the
propagation of particles [1, 2, 3, 4, 5]. Consequently, it has some
phenomenological implications that can provide an empirical ground to test
quantum gravity theories [6, 4] (for a review on such an issue, see [7]). The
effects of possible Lorentz invariance violation are expected, in particular,
to be present within frameworks implementing the quantum theory of fields
propagating on a quantized background spacetime.

An example of such a framework, where significant progress has been achieved
recently, is a quantum field theory (QFT) on a spherically symmetric quantum
spacetime described by loop quantum gravity (LQG) [8, 9, 10]. It was shown
that the evolution of quantum fields in a quantum spacetime leads to emerging
an effective (semiclassical) dressed background metric. The components of this
dressed metric depend on the fluctuations of the background quantum geometry.
In the presence of the backreaction of the fields, the emergent dressed
metric’s components depend further on the energy of the field modes [11, 10],
which is called “rainbow metric” in the literature [5, 3]. Propagation of
electromagnetic signals (or massless scalar perturbations) on this rainbow
background is superluminal, which violates the local Lorentz symmetry [10].
[See other scenarios for violation of the Lorentz symmetry, e.g., due to an
emerging rainbow metric from a massive quantum field on a loop quantum
cosmology (LQC) spacetime [12], or due to polymer quantization of the field
[13].]

The highest observable energies in the Universe are provided by cosmic gamma
rays and cosmic rays. Considering the fact that all long-duration gamma ray
bursts (GRBs) are physically connected with the core-collapse supernovae (SNe)
[6, 14, 15, 16], it is not unreasonable to expect that when observing GRBs, we
directly observe the gravitational collapse of a massive and compact star
core. Therefore, it is of pertinence to search for Lorentz violation signals
in the gravitational collapse of a massive star. Gravitational collapse of a
fluid with a variety of matter fields has been well studied within the
framework of LQG [17, 18, 19, 20]. It turns out that, when considering a
background quantum spacetime for gravitational collapse of a spherically
symmetric star, a rainbow interior region would emerge which effects can be
carried out to the exterior region through appropriate junction conditions at
the boundary of the star. This provides a fruitful scenario for the formation
of a rainbow black hole as the final state of gravitational collapse in
quantum gravity. If such black holes exist in nature, high-energy
astrophysical observations from them can raise the possibility of tests of the
Lorentz symmetry and quantum gravity theories.

The present work is concerned with the phenomenological issues of quantum
gravity in the context of gravitational collapse. It is organized as follows.
In Sec. II we present the dynamics of the gravitational collapse of a
spherically symmetric dust cloud. In Sec. III and IV we study the quantum
theory of a massless scalar field on the background spacetime of a collapsing
dust ball. Then, we quantize the background due to LQG and show that the
theory of quantum field on this quantum background corresponds to a quantum
theory of the same field on an effective, dressed background. By considering
the backreaction of the field, we will show that the components of the dressed
background metric depend on the energy of the field. Next, we expand the
dynamics of the emerging dressed background spacetime by means of the higher-
order quantum corrections provided by fluctuations due to moments of the
quantum spacetime state through a semiclassical regime. In Sec. V, we match
the interior spacetime to a convenient exterior geometry. We will show that
the quantum gravity effects in the interior region are carried out to the
collapse of radiation exterior spacetime by matching and a rainbow black hole
will emerge. We will discuss some optical properties of the emerging rainbow
black hole in Sec. VI. Finally, in the Sec. VII, we will present the
conclusion of our work.

##  II Gravitational collapse of a dust field

Our purpose in this section is to construct a classical model of gravitational
collapse with an interior region filled with an irrotational dust field $T$,
such that in the late time stages of the collapse (cf. next sections), when
the fluid enters the Planck regime, the quantum gravity effects could alter
the nature of singularity or/and development of trapped surfaces in spacetime.
Thus, we consider a homogeneous, isotropic barotropic fluid for the matter
content of the interior collapse background which can have a Friedmann-
Lemaître-Robertson-Walker (FLRW) metric, equipped with the coordinates
$(x_{0},\mathbf{x})$ and a scale factor $a(x_{0})$. We assume that
$x_{0}\in\mathbb{R}$ is a generic time coordinate and
$\mathbf{x}\in\mathbb{T}^{3}$ is the spatial coordinates [$\mathbb{T}^{3}$ is
the three-torus with coordinates $x^{j}\in(0,\ell)$].

For the background matter source, being the irrotational dust field $T$, the
Lagrangian density is given by [21]

$\mathcal{L}_{T}\,=\,-\tfrac{1}{2}\sqrt{-g}\rho_{T}\big{(}g^{\mu\nu}\partial_{\mu}T\partial_{\nu}T+1\big{)},$
(1)

where $\rho_{T}$ is a multiplier enforcing the gradient of the dust field to
be timelike (e.g., see also [22, 23]). We further consider a (inhomogeneous)
massless scalar field perturbation $\phi(t,\mathbf{x})$, with the Lagrangian

$\mathcal{L}_{\phi}\,=\,-\tfrac{1}{2}\sqrt{-g}\,g^{\mu\nu}\partial_{\mu}\phi\partial_{\nu}\phi,$
(2)

which propagates on the background spacetime of the collapsing cloud. The
corresponding action for the background geometry, coupled to dust and the
scalar perturbation, reads

$\displaystyle S=\int d^{4}x\sqrt{-g}\Big{[}\frac{\mathcal{R}}{16\pi
G}+\mathcal{L}_{T}\Big{]}+S_{\phi}\,.$ (3)

Then, on the full phase space, total Hamiltonian density is written as

$\displaystyle{\cal H}$ $\displaystyle=$ $\displaystyle{\cal H}_{\rm
grav}+{\cal H}_{T}+{\cal H}_{\phi},$ (4)

where, ${\cal H}_{\rm grav}$, ${\cal H}_{T}$ and ${\cal H}_{\phi}$ are
respectively, the Hamiltonian densities of the gravitational sector, dust and
the scalar field.

Using a general background metric in Arnowitt-Deser-Misner decomposition,

${\rm d}s^{2}=-N^{2}{\rm d}x_{0}^{2}+q_{ab}(N^{a}{\rm d}x_{0}+{\rm
d}x^{a})(N^{b}{\rm d}x_{0}+{\rm d}x^{b}),$ (5)

the Hamiltonian densities of the scalar perturbation, ${\cal H}_{\phi}$, and
the background dust field, ${\cal H}_{T}$, are written respectively, as

$\displaystyle{\cal
H}_{\phi}=\frac{N}{2}\left[\frac{P_{\phi}^{2}}{\sqrt{q}}+\sqrt{q}q^{ab}(\nabla_{a}\phi)(\nabla_{b}\phi)\right],\quad\quad$
(6)

and

$\displaystyle{\cal
H}_{T}=\frac{1}{2}\left[\frac{p_{T}^{2}}{\rho_{T}\sqrt{q}}+\frac{\rho_{T}\sqrt{q}}{p_{T}^{2}}\left(p_{T}^{2}+q^{ab}C^{D}_{a}C^{D}_{b}\right)\right].\quad\quad$
(7)

In the above equations, $N$ and $N^{a}$ denote, respectively, the lapse
function and the shift vector, and $q_{ab}$ is the (spatial) three-metric with
the conjugate momentum $\pi^{ab}$. Moreover, $C^{D}_{a}=-p_{T}\partial_{a}T$,
where $p_{T}$ is the momentum conjugate to $T$, given by

$\displaystyle p_{T}=(\sqrt{q}\rho_{T}/N)(\partial_{0}T+N^{a}\partial_{a}T).$
(8)

The stress-energy tensor for the dust field can then be obtained as

$\displaystyle{\rm T}_{\mu\nu}^{\rm(D)}=\frac{2}{\sqrt{-g}}\frac{\delta
S_{T}}{\delta g^{\mu\nu}}=\rho_{T}u_{\mu}u_{\nu}\,,$ (9)

where $u_{\mu}=\partial_{\mu}T$ is the four-velocity field in the dust
coordinate frame which satisfies the condition $g^{\mu\nu}u_{\mu}u_{\nu}=-1$.
From the equation of motion for $\rho_{T}$, we obtain [21]

$\displaystyle\rho_{T}=\frac{1}{\sqrt{q}}\frac{p_{T}^{2}}{\sqrt{p_{T}^{2}+q^{ab}C^{D}_{a}C^{D}_{b}}}\,.$
(10)

Equation (9) represents the stress-energy tensor of a perfect fluid with the
energy density $\rho_{T}$ and a vanishing pressure, so $\rho_{T}$ is regarded
as the dust energy density. Now, by substituting Eq. (10) in ${\cal H}_{T}$
and fixing time gauge $x_{0}=T$, the total Hamiltonian density (4) becomes

$\displaystyle{\cal H}={\cal H}_{\rm grav}+p_{T}+{\cal H}_{\phi}.$ (11)

The dynamics of the gravity-matter system is then generated by the regulated
integral

$\displaystyle NH=\int_{\cal{V}}N\left({\cal H}_{\rm grav}+p_{T}+{\cal
H}_{\phi}\right),$ (12)

where, for simplicity, we have considered $\cal{V}$ to be a cell of unit
volume.

We will consider the internal spacetime model in the classical regime given by
a marginally bound ($k=0$) case,

$ds^{2}=-dx_{0}^{2}+a^{2}(x_{0})\,d\mathbf{x}^{2},$ (13)

so that the dust fluid begins to collapse from a very large physical radius.
Then, the physical trajectories lie on the surface of the Hamiltonian
constraint

$\displaystyle-p_{T}=H_{\rm grav}+H_{\phi},$ (14)

where $H_{\rm grav}$ is given by

$\displaystyle H_{\rm grav}=-\frac{3\pi G}{2\alpha_{o}}b^{2}|v|.$ (15)

In LQC, the gravitational part of the phase space is conveniently
coordinatized by a canonically conjugate pair $\\{b,v\\}=2$, where
$v=a^{3}/\alpha_{o}$ is the oriented volume and $b$ is the Hubble parameter
$b=\gamma(\dot{a}/a)$. As usual, a “dot” refers to a derivative with respect
to the proper time $x_{0}$. Moreover,
$\alpha_{o}=2\pi\gamma\sqrt{\Delta}\;\ell_{\rm Pl}^{2}$, in which
$\Delta\equiv 4\sqrt{3}\pi\gamma\ell_{\rm Pl}^{2}$ is the LQC area gap [24].

By solving the constraint equation (14), we find the evolution equation for
the collapse as $(b/\gamma)^{2}=(\dot{a}/a)^{2}=(8\pi G/3)\rho$, being the
standard Friedmann equation, where $\rho$ is the total density of the system
including the energy densities of dust and scalar perturbation. Moreover,
notice that for the collapsing process here, $\dot{a}<0$. For the homogeneous
spatial slices here we get $C^{D}_{a}=0$, so
$p_{T}=a^{3}\rho_{T}(\partial_{0}T)$. Therefore, the energy density $\rho_{T}$
in Eq. (10) reduces to $\rho_{T}=p_{T}/a^{3}$. From the Hamilton equation of
motion $\dot{p}_{T}=\\{p_{T},H\\}=0$, it turns out that $p_{T}$ is a constant
of motion, so as expected, the energy density of the dust becomes
$\rho_{T}=p_{T}\,a^{-3}$.

In order to fix sufficient initial conditions for the collapse, we assume that
$\rho_{0}$ and $a_{0}$ are, respectively, the total energy density and the
scale factor of the collapsing cloud at the initial time, $x_{0}=0$. On large
scales in the classical region, the energy density of the scalar perturbation
is negligible so that, at the initial configuration of the collapse, the
energy density of the dust field dominates; we thus assume that
$\rho\approx\rho_{T}$. Nevertheless, as the collapse enters the quantum
regime, the backreaction of the scalar field, $\phi$, on the background
quantum geometry will become important. Our aim in the next sections,
therefore, will be to investigate the effects of this backreaction on the
evolution of trapped surfaces and the emergence of mode-dependent exterior
spacetime. In fact, the scenario that we will consider is to discuss, by
taking a full quantum system, the effects of the backreaction of scalar
perturbation on the evolution of the background quantum geometry, and to
explore a suitable (semiclassical) geometry for the exterior region. In
particular, we shall assume that the homogenous classical sector of the scalar
field, $\hat{\phi}$, has no effect on the background; by definition
$\hat{\phi}=\langle\hat{\phi}\rangle\openone+\delta\hat{\phi}$, this means
that we assume $\langle\hat{\phi}\rangle=0$ for the vacuum expectation value
of $\hat{\phi}$, while $\delta\hat{\phi}$ describes the inhomogeneous part of
$\hat{\phi}$. So, there is no backreaction on the background geometry caused
by the homogenous part and we will consider only the excitations and particle
states of the massless field. At the semiclassical level, these excitations
give only an effective description (sum over all modes) for spacetime.

In the classical region, we write the energy density of the collapsing dust
cloud in the form

$\displaystyle\rho={\rho_{T}}_{0}(a_{0}/a)^{3},$ (16)

where $p_{T}={\rho_{T}}_{0}a_{0}^{3}$. As the collapse evolves, the energy
density of the dust grows and ultimately diverges at $a=0$. Therefore, a
singularity will form at the end state of the collapse. This singularity will
be covered by a Schwarzschild horizon during the dynamical evolution of the
collapse which is extracted through a suitable matching at the boundary of the
dust cloud. At a given time $x_{0}$ and for a fixed shell with the radius $r$,
the mass of the dust cloud reads $M=(4\pi/3)\rho R^{3}$, where $R=ra$ is the
physical radius of the collapsing shell. It turns out that, for the dust
energy density (16), the mass $M$ is constant and equals the initial mass,
$M_{0}=(4\pi/3)\rho_{0}r^{3}a^{3}_{0}$, at $x_{0}=0$. This is the mass of the
exterior Schwarzschild black hole with the horizon radius $R_{\rm S}=2GM$.

##  III QFT in quantized interior background

In this section, we will first show that the quantum theory of the scalar
field perturbation, $\phi$, propagating on the interior quantum background of
the dust cloud, corresponds to an emerging quantum theory for the field on an
effective, dressed background spacetime.

###  III.1 Perturbed background

In canonical quantum gravity coupled with dust field [21], in the case when
the effect of a scalar perturbation (denoted by $\hat{\mathcal{H}}_{\phi}$) in
Eq. (4) is negligible, the total Hamiltonian operator, $\hat{H}_{\rm
geo}=\hat{H}_{\rm grav}+\hat{H}_{T}$, of the system is well defined on
$\mathscr{H}_{\rm kin}^{o}=\mathscr{H}_{\rm grav}\otimes\mathscr{H}_{T}$,
where $\mathscr{H}_{\rm grav}$ is a suitable Hilbert space for the gravity
sector and $\mathscr{H}_{T}$ is the dust sector of the kinematical Hilbert
space, which is quantized according to the Schrödinger picture with the
Hilbert space $L^{2}(\mathbb{R},dT)$. The physical states
$\Psi_{o}(v,T)\in\mathscr{H}_{\rm kin}^{o}$ are those lying on the kernel of
$\hat{H}_{\rm geo}$. Therefore, the states $\Psi_{o}(v,T)$ are solutions to
the self-adjoint Hamiltonian constraint $\hat{H}_{\rm geo}\Psi_{o}=0$, so that

$\displaystyle i\hbar\partial_{T}\Psi_{o}(v,T)=\hat{H}_{\rm
grav}\Psi_{o}(v,T).$ (17)

Here, $\hat{H}_{\rm grav}$ is a well-defined, self-adjoint operator acting on
$\mathscr{H}_{\rm grav}$. In quantum theory, the polymer representation of the
Poisson algebra of $v$ and $b$ is characterized by the Hilbert space
$\mathscr{H}_{\rm grav}=L^{2}(\bar{\mathbb{R}},{\rm
d}\mu_{\operatorname*{Bohr}})$, where $\bar{\mathbb{R}}$ is the Bohr
compactification of the real line and $d\mu_{\rm Bohr}$ is the Haar measure on
it [25]. Thereby, the gravitational Hamiltonian operator is expressed as [26]

$\hat{H}_{\rm grav}=\frac{3\pi
G}{8\alpha_{o}}\sqrt{\hat{v}}\big{(}\hat{N}^{2}-\hat{N}^{-2}\big{)}^{2}\sqrt{\hat{v}}\,,$
(18)

where, $\hat{v}|v\rangle=v|v\rangle$, and the operator
$\hat{N}\equiv\widehat{\exp(ib/2)}$ acts on the basis $\\{|v\rangle\\}$, i.e.
the eigenstates of $\hat{v}$, as $\hat{N}|v\rangle=|v+1\rangle$, so that,
$[\hat{b},\hat{v}]=2i\hbar$.

When the contribution of the scalar field $\phi$ in the Hamiltonian constraint
(4) is significant, the kinematical Hilbert space for the full, quantized
gravity-matter system (dust plus scalar perturbation) becomes
$\mathscr{H}_{\rm kin}=\mathscr{H}_{\rm
grav}\otimes\mathscr{H}_{T}\otimes\mathscr{H}_{\phi}$, where the perturbation
sector is quantized due to the Schrödinger picture with
$\mathscr{H}_{\phi}=L^{2}(\mathbb{R},d\phi)$. Now, the total states
$\Psi\in\mathscr{H}_{\rm kin}$ of the system are different from the pure
geometrical states $\Psi_{o}$ and are solutions to a new evolution equation

$\displaystyle i\hbar\partial_{T}\Psi(v,\phi,T)=\big{(}\hat{H}_{\rm
grav}+\hat{H}_{\phi}\big{)}\Psi(v,\phi,T).$ (19)

On the quantized background here, the gravitational sectors of the quantized
Hamiltonian (6) of the perturbation turn out to be operators on
$\mathscr{H}_{\rm kin}$; thus, the quantum Hamiltonian of the massless scalar
perturbation becomes

$\displaystyle\hat{H}_{\phi}=\tfrac{1}{2}\Big{[}\hat{V}^{-1}\otimes\hat{P}_{\phi}^{2}+\hat{V}^{1/3}\otimes(\nabla_{i}\hat{\phi})^{2}\Big{]},$
(20)

where $V$, defined as $V=\ell^{3}a^{3}=\alpha_{o}v$, is the physical volume of
the collapsing cloud. For convenience, we set $\ell=1$ throughout this section
and will bring it back again into our formalism in Sec. V.

By using the Fourier expansion, we can rewrite $\hat{H}_{\phi}$ in Eq. (20) as
an assembly of the Hamiltonians of decoupled harmonic oscillators, each
represented by a pair of canonically conjugate variables
$(Q_{\mathbf{k}},P_{\mathbf{k}})$ [8]. It reads

$\hat{H}_{\phi}:=\sum_{\mathbf{k}\in{\cal
L}}\hat{H}_{\mathbf{k}}=\frac{1}{2}\sum_{\mathbf{k}\in{\cal
L}}\Big{[}\hat{V}^{-1}\otimes\hat{P}_{\mathbf{k}}^{2}+k^{2}\,{V}^{1/3}\otimes\hat{Q}_{\mathbf{k}}^{2}\Big{]},$
(21)

where, the wave vectors $\mathbf{k}~{}(\in 2\pi\mathbb{Z})$ span a three-
dimensional lattice ${\cal L}$. We will focus on a linear response theory and
will study the quantum theory of a single mode $\mathbf{k}$ of the scalar
field on the background quantum spacetime. Thereby, Eqs. (19) and (21) for a
given mode $\mathbf{k}$ yield

$i\hbar\partial_{T}\Psi_{\mathbf{k}}(v,Q_{\mathbf{k}},T)=\big{(}\hat{H}_{\rm
grav}+\hat{H}_{\mathbf{k}}\big{)}\Psi_{\mathbf{k}}(v,Q_{\mathbf{k}},T).\quad\quad$
(22)

In order to find the quantum theory of the test field $Q_{\mathbf{k}}$ on an
emergent background spacetime, we should simplify the Schrödinger equation
(22) and represent it as an effective equation for the state
$\psi_{\mathbf{k}}(Q_{\mathbf{k}})\in\mathscr{H}_{\mathbf{k}}$ only
($\mathscr{H}_{\mathbf{k}}$ represents the Hilbert space of each mode
$\mathbf{k}$). To do so, we employ the following algorithm:

  1. i)

To decompose the heavy (gravity) and light (scalar perturbation) degrees of
freedom in the wave function, as
$\Psi_{\mathbf{k}}(v,Q_{\mathbf{k}},T)=\Psi(v,T)\otimes\psi_{\mathbf{k}}(Q_{\mathbf{k}},T)$,
we will employ the Born-Oppenheimer (BO) approximation. This approximation
enables us to take into account the backreaction between the field and the
geometry.

  2. ii)

To make our resulting evolution comparable to that of a quantum field on a
classical dynamical background, instead of working in the “interaction
picture” used in Ref. [8], we will trace out the heavy and light degrees of
freedom in Eq. (22) to drive the evolution equation for the scalar
perturbation and the background states, respectively, as

$\displaystyle\qquad i\hbar\partial_{T}\psi_{\mathbf{k}}(Q_{\mathbf{k}},T)$
$\displaystyle=\hat{H}_{\mathbf{k}}\,\psi_{\mathbf{k}}(Q_{\mathbf{k}},T),$
(23) $\displaystyle\quad\quad i\hbar\partial_{T}\Psi(v,T)$
$\displaystyle=\big{(}\hat{H}_{\rm
grav}+\langle\hat{H}_{\mathbf{k}}(\hat{v})\rangle\big{)}\Psi(v,T),$ (24)

where,
$\langle\hat{H}_{\mathbf{k}}(\hat{v})\rangle=\langle\psi_{\mathbf{k}}(T)|\hat{H}_{\mathbf{k}}(\hat{v})|\psi_{\mathbf{k}}(T)\rangle$.
In this pattern, quantum geometry and field are described using the
Schrödinger picture, in which expectation values evolve over the time
parameter and give results of Ref. [8] as the mean (test) field approximation.

  3. iii)

We plan to solve Eqs. (23) and (24) step by step perturbatively. At the first
step, we will consider the mean field limit, where
$\langle\hat{H}_{\mathbf{k}}(\hat{v})\rangle=0$, and find the unperturbed
geometry state $\Psi_{o}(v,T)$. This state is then used to evaluate the
expectation values of the geometry component operators present as components
of $\hat{H}_{k}$, which allows us to evaluate the solutions to Eq. (23) to
find the scalar field’s eigenfunctions,
$\chi^{n}_{\mathbf{k}}(Q_{\mathbf{k}};g^{(1)},g^{(2)})$, parametrized by said
geometry expectation values, $g^{(1)}$ and $g^{(2)}$, given by

$\displaystyle
g^{(1)}\equiv(1/2)\langle\hat{V}^{-1}\rangle_{o}\quad\text{and}\quad
g^{(2)}\equiv(1/2)\langle\hat{V}^{1/3}\rangle_{o}.$ (25)

A procedure analogous to this step has already been performed in Ref. [8]. In
the next step, in order to construct the backreacted state $\Psi_{1}(v,T)$, we
will use the obtained eigenfunctions
$\chi^{n}_{\mathbf{k}}(Q_{\mathbf{k}};g^{(1)},g^{(2)})$ to evaluate
$\langle\hat{H}_{\mathbf{k}}(\hat{v})\rangle$. The result is then put into Eq.
(24); however, here the variables $g^{(1)}$ and $g^{(2)}$ are again promoted
to geometry operators (powers of $\hat{V}$) according to the form specified in
Eq. (25). This, in turn, allows us to find the modified eigenfunctions
$\xi^{\mu}_{\mathbf{k}}(v)$ of the geometry. This step is the second-order
modification to the so-called test field approximation, presented in Ref. [8],
wherein the backreaction effects were discarded.

In BO approximation, the total wave function consists of the products of two
sets of eigenstates: The first one is the (discrete) field mode’s eigenstate
$\chi^{n}_{\mathbf{k}}$, being the solution to the stationary state equation

$\displaystyle\hat{\tilde{H}}_{\mathbf{k}}\,\chi^{n}_{\mathbf{k}}\big{(}Q_{\mathbf{k}};g^{(1)},g^{(2)}\big{)}=\epsilon^{n}_{\mathbf{k}}\big{(}g^{(1)},g^{(2)}\big{)}\,\chi^{n}_{\mathbf{k}}\big{(}Q_{\mathbf{k}};g^{(1)},g^{(2)}\big{)},$
(26)

where $\hat{\tilde{H}}_{\mathbf{k}}\in\mathscr{H}_{\mathbf{k}}$ is the
Hamiltonian operator of the scalar field propagating on the specified (fixed)
quantum geometry. The second one is a chosen (usually semiclassical) state of
the background wave function $\Psi(v,T)$, being a solution to Eq. (24). More
precisely, Eq. (26) is constructed by a partial tracing over the geometry
degrees of freedom, Eq. (22) [defined on the full Hilbert space with the
product state $\Psi(v,T)\otimes\psi_{\mathbf{k}}(Q_{\mathbf{k}},T)$]. In this
approximation, $\hat{\epsilon}^{n}_{\mathbf{k}}$, the energy eigenvalue of the
test field $Q_{\mathbf{k}}$, is still an operator on the gravitational Hilbert
space, $\mathscr{H}_{\rm grav}$.

Partial tracing of Eq. (22) over the geometry degrees of freedom (d.o.f)
yields

$\displaystyle\hat{\tilde{H}}_{\mathbf{k}}$
$\displaystyle=g^{(1)}\,\hat{P}_{\mathbf{k}}^{2}+k^{2}g^{(2)}\,\hat{Q}_{\mathbf{k}}^{2}\,,$
(27a) $\displaystyle\epsilon^{n}_{\mathbf{k}}(g^{(1)},g^{(2)})$
$\displaystyle=\langle\Psi_{o}(T)|\,\epsilon^{n}_{\mathbf{k}}\big{(}\hat{g}^{(1)},\hat{g}^{(2)}\big{)}|\Psi_{o}(T)\,\rangle,$
(27b)

where $\hat{g}^{(1)}$ and $\hat{g}^{(2)}$ are composite operators expressed
via $\hat{v}$. Thus, the scalar perturbation behaves as that of a (geometry-
dependent) quantum harmonic oscillator. The resulting eigenvalue problem,
thus, takes the form

$\displaystyle
g^{(1)}\,\frac{d^{2}|\chi^{n}\rangle}{dQ_{\mathbf{k}}^{2}}+k^{2}g^{(2)}\,Q_{\mathbf{k}}^{2}\,|\chi^{n}\rangle$
$\displaystyle=$ $\displaystyle\epsilon^{n}_{\mathbf{k}}\,|\chi^{n}\rangle.$
(28)

The solutions to differential equation (28) are well known, given by

$\displaystyle\epsilon^{n}_{\mathbf{k}}$
$\displaystyle=\sqrt{g^{(1)}g^{(2)}}\left(n+\tfrac{1}{2}\right)\hbar k,$ (29a)
$\displaystyle|\chi^{n}\rangle$
$\displaystyle=a_{n}k^{1/4}B^{1/8}\exp\left(-x^{2}/2\right)H_{n}(x),$ (29b)

where, we have defined

$\displaystyle a_{n}$ $\displaystyle\equiv$ $\displaystyle\left(\frac{1}{\pi
2^{2n}(n!)^{2}}\right)^{\frac{1}{4}},$ $\displaystyle B$ $\displaystyle\equiv$
$\displaystyle\frac{g^{(2)}}{g^{(1)}\hbar^{2}}\quad{\rm and}\quad x\equiv
k^{1/2}B^{1/4}Q_{\mathbf{k}}.$ (30)

In further treatment, we would like to follow the procedure introduced in Ref.
[27]. There, an essential step was treating the emergent description of the
matter field (in our case, the scalar field) as parametrized by a single
geometric variable, the volume
$V^{\prime}=\langle\Psi_{o}(T)|\hat{V}|\Psi_{o}(T)\rangle$. In subsequent
steps of the considered procedure that parameter was promoted back to quantum
operator. Here, however, the scalar field description involves two expectation
values: $\langle\hat{V}^{-1}\rangle$ and $\langle\hat{V}^{1/3}\rangle$. In
order to introduce the parametrization analogous to that in Ref. [27], we note
that the functions $g^{(i)}$ can be expanded in terms of the central Hamburger
moments corresponding to the volume [28], namely,

$\langle\hat{V}^{\alpha}\rangle=\langle\hat{V}\rangle^{\alpha}+\sum_{i=1}^{\infty}\left({\begin{array}[]{c}\alpha\\\
i\\\ \end{array}}\right)\langle\hat{V}\rangle^{\alpha-i}G^{i00}\,,$ (31)

where $G^{i00}=\langle(\delta\hat{V})^{i}\rangle$. Since the background state
$\Psi_{o}$ is _chosen_ , both $V^{\prime}$ and $G^{i00}$ are determined as
functions of $T$. Had $V^{\prime}(T)$ been invertible (which happens, for
example, in geometrodynamics if we restrict ourselves to the post-big-bang
epoch), we would be able to define

$G^{i00}(V^{\prime})=G^{i00}(T(V^{\prime})).$ (32)

In LQC, however, there are two reasons preventing us from doing so:

  1. 1.

The dynamics of the background state features a bounce (see, for example,
[26]), thus, the function $V^{\prime}(T)$ is not globally invertible. One
could, in principle, choose a state symmetric with respect to the bounce, that
is, such that

$\displaystyle V^{\prime}(T_{B}+\delta T)$
$\displaystyle=V^{\prime}(T_{B}-\delta T),$ (33a) $\displaystyle
G^{i00}(T_{B}+\delta T)$ $\displaystyle=G^{i00}(T_{B}-\delta T),$ (33b)

where, $T_{B}$ is the time of the bounce, however, this choice is a fine-
tuning and there is no physical reason distinguishing it.

  2. 2.

The expectation value of $\hat{V}$ never drops below
$V^{\prime}_{B}\propto-\langle\hat{H}_{\rm grav}\rangle$; thus $T(V^{\prime})$
is not defined for $V^{\prime}<V^{\prime}_{B}$.

As a consequence, in order to introduce the desired parametrization, we have
to neglect at this step all the second- and higher-order quantum corrections
(encoded in $G^{i00}$) of the background state, leaving only the quantum
imprint on the trajectory. Then, $g^{(i)}$’s in Eq. (25) reduce to

$g^{(1)}=(1/2)\langle\hat{V}\rangle^{-1}_{o},\quad
g^{(2)}=(1/2)\langle\hat{V}\rangle^{1/3}_{o},$ (34)

and in consequence

$\displaystyle\hat{\epsilon}^{n}_{\mathbf{k}}$
$\displaystyle=\left(b_{\mathbf{k}}^{\dagger
n}b_{\mathbf{k}}^{n}+\frac{1}{2}\right)\langle\hat{V}\rangle^{-1/3}_{o}\,\hbar
k\,,$ (35a) $\displaystyle|\chi^{n}\rangle$
$\displaystyle=a_{n}\left(\frac{k}{\hbar}\right)^{\frac{1}{4}}\langle\hat{V}\rangle^{1/12}_{o}\exp\left(-\frac{k\langle\hat{V}\rangle^{1/3}_{o}}{2\hbar}Q_{\mathbf{k}}^{2}\right)\,H_{n}\left(\sqrt{\frac{k}{\hbar}}\langle\hat{V}\rangle^{1/6}_{o}Q_{\mathbf{k}}\right).$
(35b)

Having defined $|\chi^{n}_{\mathbf{k}}\rangle$ as the eigenfunctions of
$\hat{\tilde{{H}}}_{\mathbf{k}}$ with eigenvalues $\epsilon^{n}_{\mathbf{k}}$,
we can now turn back to the geometry eigenfunction components
$\xi^{\mu}_{\mathbf{k}}(g^{(1)},g^{(2)})$ of the perturbed quantum geometry
state (24). After tracing out the scalar field degrees of freedom, the
evolution equation for eigenfunctions of geometry is obtained as

$\big{[}\hat{H}_{\rm
gr}+\epsilon^{n}_{\mathbf{k}}(\hat{v})\big{]}\xi^{\mu}_{\mathbf{k}}(v)=E^{\mu}_{\mathbf{k}}\;\xi^{\mu}_{\mathbf{k}}(v),$
(36)

where

$\displaystyle\epsilon^{n}_{\mathbf{k}}(\hat{v})$
$\displaystyle=\left(n+\dfrac{1}{2}\right)\hbar k\hat{V}^{-1/3}$
$\displaystyle=:N_{k}\hbar k\,\hat{V}^{-1/3},$ (37)

represents the energy of the mode of test field state, in which $\hat{N}_{k}$,
in the adiabatic regime, reduces to the number operator of a harmonic
oscillator [29].

The eigenvalue problem (36) differs from the one of the background state
(studied in Ref. [26]) only by the bounded potential quickly decaying to zero
as $v$ increases. As a consequence, the operator on its left-hand side will
share the spectral property of the background one: Its spectrum is
nondegenerate and continuous and consists of the entire real line. The
eigenvectors $\xi^{\mu}_{\mathbf{k}}(v)$ can be found by numerical means via
methods used in Refs. [30, 31]. The properties of these eigenvectors are quite
similar to those of the background Hamiltonian, $\hat{H}_{\rm grav}$. Each has
a form of the reflected wave further featuring the region of exponential
suppression around $v$ of the size depending on $\mu$ and $\mathbf{k}$. In the
next subsection we will find eigenfunctions $\xi^{\mu}_{\mathbf{k}}(v)$ of Eq.
(36) and show that, for large $v$, they feature the following asymptotic
behavior:

$\displaystyle\xi^{\mu}_{\mathbf{k}}(v)$
$\displaystyle=\frac{C}{v^{1/4}}(1+f(\mu,\mathbf{k}))\cos\left(\mu
v^{1/2}+\varphi(\mu,\mathbf{k})\right)+\mathcal{O}(v^{-9/4}),$ (38)

where $\varphi(\mu,\mathbf{k})$ is a phase shift and $C$ is a normalization
factor. This asymptotic actually provides for us a precise definition of the
label $\mu$ in the choice of which we have a freedom due to the continuity of
the spectrum of the studied operator. The form of the asymptotic implies that
$\xi^{\mu}_{\mathbf{k}}$ are Dirac-delta normalizable. We can thus form out of
them an orthonormal basis (for each value of ${\mathbf{k}}$ independently),
setting

$(\xi^{\mu}_{\mathbf{k}}|\xi^{\mu^{\prime}}_{\mathbf{k}})=\delta(\mu-\mu^{\prime}).$
(39)

###  III.2 Rate of convergence of bases

To solve Eq. (36) numerically, we need to explicitly show that the convergence
rate (38) exists. This specific rate has applications in numerical
calculations of LQC [31], results of which will be presented in a separate
paper [32]. To study the rate of convergence, we compare
$\xi^{\mu}_{\mathbf{k}}(v)$ with eigenfunctions
$\underline{e}^{\mu}_{\mathbf{k}}(v)$ of the Wheeler-DeWitt (WDW) analog of
Eq. (19) at the asymptotic region.

The quantum Hamiltonian constraint in WDW theory can be expressed as a
differential analog of LQC evolution operator $\Theta$, where an action of the
operator $\Theta$ equals

$\begin{split}[\Theta\psi](v)&=f_{-}(v)\Psi(v-4)-f_{o}(v)\Psi(v)+f_{+}(v)\Psi(v+4),\end{split}$
(40)

with

$\displaystyle f_{\pm}(v)$ $\displaystyle=(3\pi G/8\alpha_{o})(v\pm
4)^{1/2}v^{1/2},$ $\displaystyle f_{o}(v)$ $\displaystyle=(3\pi
G/4\alpha_{o})v-N_{k}\hbar k\alpha_{o}^{-1/3}v^{-1/3}.$ (41)

To arrive to the WDW equation, we select the factor ordering consistent with
the one of Eq. (18), so we get

$\displaystyle i\partial_{T}\underline{\Psi}(v,\phi)$
$\displaystyle=\underline{\Theta}\,\underline{\Psi}(v,\phi)$
$\displaystyle:=\frac{6\pi
G}{\alpha_{o}}\,|v|^{1/2}\partial_{v}\partial_{v}|v|^{1/2}\,\underline{\Psi}(v,\phi).$
(42)

One can build a basis in $\underline{\mathscr{H}}_{\operatorname*{grav}}$
(Hilbert space of WDW theory) out of the eigenfunctions
$\underline{e}_{\mu}(v)$ corresponding to non-negative eigenvalues

$[\underline{\Theta}\;\underline{e}_{\mu}](v)=-\omega^{2}(\mu)\underline{e}_{\mu}(v),$
(43)

where $\omega(\mu)=\sqrt{3\pi G/2\alpha_{o}}\mu$ with $\mu>0$. A general
solution to Eq. (43) is

$\underline{e}_{\mu}(v)=c_{1}J_{1}\left(2\sqrt{v}\omega\right)+2c_{2}Y_{1}\left(2\sqrt{v}\omega\right),$
(44)

where $J_{1}\left(2\sqrt{v}\omega\right)$ and
$Y_{1}\left(2\sqrt{v}\omega\right)$ are Bessel functions of the first and
second kind, respectively. These solutions asymptotically tend to the
orthonormal basis:

$\underline{e}^{\pm}_{\mu}(v)\approx\dfrac{|v|^{-1/4}}{\sqrt{4\pi}}\;e^{\pm
i\mu|v|^{1/2}}.$ (45)

To verify asymptotes of $\xi^{\mu}_{\mathbf{k}}(v)$, we start with rewriting
Eq. (40), the second-order difference equation, in a first-order form,
introducing the vector notation

$\vec{\xi}^{\mu}_{\mathbf{k}}(v):=\left(\begin{array}[]{l}\xi^{\mu}_{\mathbf{k}}(v)\\\
\xi^{\mu}_{\mathbf{k}}(v-4)\end{array}\right).$ (46)

Using it, Eq. (40) turns to

$\vec{\xi}^{\mu}_{\mathbf{k}}(v+4)=\mathbf{A}(v)\,\vec{\xi}^{\mu}_{\mathbf{k}}(v),$
(47)

where the matrix $\mathbf{A}$ can be expressed as

$\mathbf{A}(v)=\left(\begin{array}[]{cc}\frac{f_{o}(v)-\omega^{2}(\mu)}{f_{+}(v)}&-\frac{f_{-}(v)}{f_{+}(v)}\\\
1&0\end{array}\right).$ (48)

To relate $\xi^{\mu}_{\mathbf{k}}(v)$ with $\underline{e}^{\pm}_{\mu}$, we
note that the value of $\xi^{\mu}_{\mathbf{k}}(v)$ at each pair of consecutive
points $v$ and $v+4$ can be encoded as a linear combination of the WDW
components of $\underline{e}^{\pm}_{\mu}$, that is,

$\vec{\xi}^{\mu}_{\mathbf{k}}(v+4)=\mathbf{B}^{\mu}_{\mathbf{k}}(v)\vec{\chi}^{\mu}_{\mathbf{k}}(v+4),$
(49)

where the transformation matrix $\mathbf{B}^{\mu}_{\mathbf{k}}$ is defined as
follows:

$\mathbf{B}^{\mu}_{\mathbf{k}}(v):=\left(\begin{array}[]{ll}\underline{e}^{+}_{\mu}(v+4)&\underline{e}^{-}_{\mu}(v+4)\\\
\underline{e}^{+}_{\mu}(v)&\underline{e}^{-}_{\mu}(v)\end{array}\right).$ (50)

Using the objects defined above, we can rewrite Eq. (47) as the iterative
equation for the vectors of coefficients $\vec{\chi}^{\mu}_{\mathbf{k}}$:

$\displaystyle\vec{\chi}^{\mu}_{\mathbf{k}}(v+4)$
$\displaystyle=\mathbf{B}^{\mu^{-1}}_{\mathbf{k}}(v)\mathbf{A}(v)\mathbf{B}^{\mu}_{\mathbf{k}}(v-4)\vec{\chi}^{\mu}_{\mathbf{k}}(v)$
$\displaystyle=:\mathbf{M}^{\mu}_{\mathbf{k}}(v)\vec{\chi}^{\mu}_{\mathbf{k}}(v).$
(51)

The exact elements of the matrix $\mathbf{M}^{\mu}_{\mathbf{k}}(v)$ can be
calculated explicitly for the coefficients of the evolution operator (40). By
straightforward calculations, one can find that it has the following
asymptotic behavior:

$\mathbf{M}^{\mu}_{\mathbf{k}}(v)=\openone+\mathcal{O}(v^{-1/3}).$ (52)

This result does not grant the rate of convergence needed for Eq. (38). We can
improve the level of convergence by replacing the components
$\underline{e}^{\pm}_{\mu}(v)$ in Eq. (45) with functions,

$\displaystyle\underline{e}^{\pm}_{\mu,\mathbf{k}}(v)$
$\displaystyle=\frac{|v|^{-1/4}}{\sqrt{4\pi}}\left[1+\sum_{n=1}^{5}a_{n}\;|v^{-n/3}|\right]\exp\left[\pm
i\mu|v|^{1/2}\left(1+\sum_{n=1}^{7}b_{n}|v^{-n/3}|\right)\right],$ (53)

where coefficients $a_{n}$ and $b_{n}$ are presented in Appendix B. Direct
inspection of the asymptotics of $\mathbf{M}^{\mu}_{\mathbf{k}}(v)$ shows that

$\mathbf{M}^{\mu}_{\mathbf{k}}(v)=\openone+\mathcal{O}(v^{-3}),$ (54)

which now admits the level of convergence needed for Eq. (38). Modified
eigenfunctions (53) will be used in a subsequent paper for the normalization
procedure in LQC numerical calculations [32].

###  III.3 Emerging mode-dependent dressed cosmological background

At this point, we have at our disposal the matter field basis eigenfunctions
$|\chi^{n}_{\mathbf{k}}\rangle$, corresponding to the discrete eigenvalues
$\epsilon^{n}_{\mathbf{k}}$, parametrized by $n$ and a family of bases of the
gravitational Hilbert space formed of eigenfunctions
$|\xi^{\mu}_{\mathbf{k}})$ which can be determined numerically and normalized
using Eq. (53). To construct the complete wave function $|\Psi_{1}\rangle$, we
need to determine the spectral profiles $c_{\mathbf{k}}(\mu)$ which may differ
from the profile of the original background state. In order to determine
$\Psi_{1}(v,T)$ and describe the eigenfunctions $|\xi^{\mu}_{\mathbf{k}})$ in
a convenient manner, which is more suitable for our perturbation treatment, we
expand $|\xi^{\mu}_{\mathbf{k}})$ as follows, distinguishing the hierarchy of
corrections by

$\displaystyle|\xi^{\mu}_{\mathbf{k}})\,=:\,\underline{N}\Big{[}|\xi^{\mu}_{o})+|\delta\xi^{\mu}_{\mathbf{k}})\Big{]},$
(55)

where $\underline{N}$ is the overall normalization factor determined from
orthonormality of the bases. A first-order solution to Eq. (24) can be
constructed using profile $c_{\mathbf{k}}(\mu)$ and eigenfunctions
$\xi^{\mu}_{\mathbf{k}}(v)$ as [31, 24]

$\displaystyle\Psi_{1}(v,T)$ $\displaystyle=$
$\displaystyle\int_{\mu\in\mathbb{R}}d\mu\,c_{\mathbf{k}}(\mu)\,\xi^{\mu}_{\mathbf{k}}(v)\,e^{i\omega(k)T}.\qquad$
(56)

Here, we are interested in considering only backreaction effects of the field
on geometry and ignoring any correlation effects between these two. Within
this approximation, the perturbed wave function can be expressed [by
substituting Eq. (55) in the wave function (56)] as

$\displaystyle\Psi_{1}(v,T)\,$
$\displaystyle=\,\int_{\mu\in\mathbb{R}}d\mu\,c(\mu)\,\xi^{\mu}_{o}(v)\,e^{i\omega(k)T}+\int_{\mu\in\mathbb{R}}d\mu\,c(\mu)\,\delta\xi^{\mu}_{\mathbf{k}}(v)\,e^{i\omega(k)T}$
$\displaystyle=:\,\Psi_{o}(v,T)+\delta\Psi_{\mathbf{k}}(v,T),$ (57)

where, we used the same profile for the perturbed and unperturbed states. The
first term on the right-hand side above denotes the unperturbed wave function,
while the second term represents the corrections in the geometry quantum state
induced by backreaction of each mode of the field on the geometry. To build
first the order total wave function, following BO approximation, we focus on
the situation where the geometry and field components of the above
backreaction term are uncorrelated (separable), that is,

$\Psi^{1}_{\mathbf{k}}(v,Q_{\mathbf{k}},T)\,=\,\Psi_{1}(v,T)\otimes\psi_{\mathbf{k}}(Q_{\mathbf{k}},T),\quad$
(58)

where, $\psi_{\mathbf{k}}(Q_{\mathbf{k}},T)$ is constructed using
eigenfunctions $\chi^{n}_{\mathbf{k}}(Q_{\mathbf{k}};v)$. This indicates that,
the wave function $\Psi_{1}(v,T)$ of the background, perturbed by field’s
backreaction, depends on the energy of the mode $\mathbf{k}$. In other words,
due to backreaction, each mode of the field induces different changes in the
geometry state and probes a specific background geometry which depends on
$k=|\mathbf{k}|$.

In zeroth order, in mean field limit
$\langle\hat{H}_{\mathbf{k}}(\hat{v})\rangle=0$, Eq. (24) leads us to an
eigenvalue equation for the geometry part, which is

$\hat{H}_{\rm grav}|\xi^{\mu}_{o})\,=\,E^{\mu}_{o}|\xi^{\mu}_{o}).$ (59)

Tracing out the geometrical d.o.f. in Eq. (23), using the state
$\Psi_{o}(v,T)$ which is constructed from eigenfunctions (59), yields the
unperturbed Schrodinger-like equation

$\displaystyle
i\hbar\partial_{T}\psi_{\mathbf{k}}=\frac{1}{2}\Big{[}\langle\hat{V}^{-1}\rangle_{o}\,\hat{P}_{\mathbf{k}}^{2}+k^{2}\,\langle\hat{V}^{\frac{1}{3}}\rangle_{o}\,\hat{Q}_{\mathbf{k}}^{2}\Big{]}\psi_{\mathbf{k}}.\qquad\quad$
(60)

Having found the perturbed eigenfunctions $|\xi^{\mu}_{\mathbf{k}})$ of the
geometry eigenvalue equation (36) and constructing the perturbed state
$\Psi_{1}(v,T)$, we get the following Schrodinger-like equation for each mode
of the scalar field:

$\displaystyle
i\hbar\partial_{T}\psi_{\mathbf{k}}=\tfrac{1}{2}\Big{[}\langle\hat{V}^{-1}\rangle\,\hat{P}_{\mathbf{k}}^{2}+k^{2}\,\langle\hat{V}^{\frac{1}{3}}\rangle\,\hat{Q}_{\mathbf{k}}^{2}\Big{]}\psi_{\mathbf{k}},\qquad\quad$
(61)

in which we have defined the expectation values $\langle\cdot\rangle$ with
respect to the total perturbed state $\Psi_{1}(v,T)$, i.e.,

$\langle\hat{V}^{-1}\rangle:=\langle\Psi_{1}(v,T)|\hat{V}^{-1}|\Psi_{1}(v,T)\rangle$
(62a) and
$\langle\hat{V}^{1/3}\rangle:=\langle\Psi_{1}(v,T)|\hat{V}^{1/3}|\Psi_{1}(v,T)\rangle.$
(62b)

(Henceforth, we will drop the subscript index for the expectation value of any
operator with respect to the perturbed background quantum state.) Now, by
substituting the decomposition (57) into Eq. (61), we obtain the following
equation for the effects of the perturbed geometry state on the evolution
equation of the field:

$\displaystyle i\hbar\partial_{T}\psi_{\mathbf{k}}$
$\displaystyle=\tfrac{1}{2}\Big{[}\Big{(}\langle\hat{V}^{-1}\rangle_{o}+\langle\hat{V}^{-1}\rangle_{\delta}\Big{)}\hat{P}_{\mathbf{k}}^{2}+k^{2}\Big{(}\langle\hat{V}^{1/3}\rangle_{o}+\langle\hat{V}^{1/3}\rangle_{\delta}\Big{)}\hat{Q}_{\mathbf{k}}^{2}\Big{]}\psi_{\mathbf{k}}\,,\quad$
(63)

in which, we have defined

$\displaystyle\langle\hat{V}^{-1}\rangle_{\delta}$
$\displaystyle:=\langle\Psi_{o}|\hat{V}^{-1}|\delta\Psi_{\mathbf{k}}\rangle+\langle\delta\Psi_{\mathbf{k}}|\hat{V}^{-1}|\Psi_{o}\rangle+\langle\delta\Psi_{\mathbf{k}}|\hat{V}^{-1}|\delta\Psi_{\mathbf{k}}\rangle,$
(64) $\displaystyle\langle\hat{V}^{1/3}\rangle_{\delta}$
$\displaystyle:=\langle\Psi_{o}|\hat{V}^{1/3}|\delta\Psi_{\mathbf{k}}\rangle+\langle\delta\Psi_{\mathbf{k}}|\hat{V}^{1/3}|\Psi_{o}\rangle+\langle\delta\Psi_{\mathbf{k}}|\hat{V}^{1/3}|\delta\Psi_{\mathbf{k}}\rangle.$
(65)

Here, $\langle\hat{V}^{-1}\rangle_{\delta}$ and
$\langle\hat{V}^{1/3}\rangle_{\delta}$ are modifications to the background
dressed metric, being probed in a BO approximation due to backreaction
effects. The state $\delta\Psi_{\mathbf{k}}$ is expanded in terms of the
eigenfunctions $\delta\xi^{\mu}_{\mathbf{k}}(v)$ on the right-hand side of Eq.
(56). A numerical analysis111A thorough numerical investigation in the LQG
context will be presented in an upcoming paper [32]. for computing this state
was performed in the Appendix B. Therein, for a normalization procedure, we
have used the perturbative eigenfunctions (53) of the WDW theory. It turns out
that the leading-order terms in these eigenfunctions depend on $k$. Therefore,
the correction terms $\langle\hat{V}^{-1}\rangle_{\delta}$ and
$\langle\hat{V}^{1/3}\rangle_{\delta}$ depend also on $k$: hence, they are
mode dependent.

The effective equation (61) corresponds to an evolution equation for the
scalar perturbation’s state, $\psi_{\mathbf{k}}$, on a dressed background
metric

$\displaystyle\tilde{g}_{ab}dx^{a}dx^{b}=-\tilde{N}(T)dT^{2}+\tilde{a}^{2}(T)d\mathbf{x}^{2}.$
(66)

By comparison, we find the following relations between the components of the
emerging dressed metric and the expectation values of quantum operators of the
original spacetime metric:

$\displaystyle\tilde{N}\tilde{a}^{-3}$ $\displaystyle=$
$\displaystyle\langle\hat{V}^{-1}\rangle_{o}\left(1+\delta_{1}\right),$ (67)
$\displaystyle\tilde{N}\tilde{a}$ $\displaystyle=$
$\displaystyle\langle\hat{V}^{1/3}\rangle_{o}\left(1+\delta_{2}\right),$ (68)

where

$\displaystyle\delta_{1}(k,T)=\frac{\langle\hat{V}^{-1}\rangle_{\delta}}{\langle\hat{V}^{-1}\rangle_{o}},\quad\quad\delta_{2}(k,T)=\frac{\langle\hat{V}^{1/3}\rangle_{\delta}}{\langle\hat{V}^{1/3}\rangle_{o}}.\quad\qquad$
(69)

By solving Eqs. (67) and (68), we obtain

$\displaystyle\tilde{N}(k,T)$ $\displaystyle=\bar{N}(T)\,f(k,T)\,,$ (70)
$\displaystyle\tilde{a}(k,T)$ $\displaystyle=\bar{a}(T)\,q(k,T)\,,\quad\quad$
(71)

where

$\displaystyle f(k,T)$ $\displaystyle=$
$\displaystyle\big{(}1+\delta_{1}\big{)}^{1/4}\big{(}1+\delta_{2}\big{)}^{3/4},$
(72) $\displaystyle q(k,T)$ $\displaystyle=$
$\displaystyle\left(\frac{1+\delta_{2}}{1+\delta_{1}}\right)^{1/4},\quad\quad$
(73)

are mode-dependent functions representing the backreaction effects in the
emerged dressed metric $\tilde{g}$. In the absence of backreaction, $f(k,T)$
and $q(k,T)$ tend to unity. Moreover, $\bar{N}_{T}$ and $\bar{a}$ are
components of the dressed metric given in a test field approximation (where no
backreaction is taken into account):

$\displaystyle\bar{N}_{T}(T)$ $\displaystyle=$
$\displaystyle\left[\big{\langle}\hat{V}^{-1}\big{\rangle}_{o}~{}\big{\langle}\hat{V}^{1/3}\big{\rangle}_{o}^{3}\right]^{\frac{1}{4}},$
(74) $\displaystyle\bar{a}(T)$ $\displaystyle=$
$\displaystyle\left[\big{\langle}\hat{V}^{1/3}\big{\rangle}_{o}\,\big{\langle}\hat{V}^{-1}\big{\rangle}_{o}^{-1}\right]^{\frac{1}{4}}.\quad$
(75)

Equations. (70) and (71) present the mode-dependent components of the dressed
metric $\tilde{g}$ emerged in the herein quantum gravity regime. This implies
that a “rainbow” metric emerges in the interior region of the collapse
background spacetime, due to backreaction effects.

##  IV Effective dynamics of the dressed metric

From the point of view of a semiclassical observer, we intend to find a
corresponding evolution equation for the dressed metric component $\tilde{a}$
with respect to a new time coordinate $\tau$ with $d\tau=\tilde{N}(T,k)dT$. In
particular, we compute the Friedmann equation corresponding to the dressed
scale factor $\tilde{a}$ as
$\tilde{H}\equiv\partial_{\tau}\tilde{a}/\tilde{a}$:

$\displaystyle\tilde{H}=\frac{1}{4}\left(\frac{\partial_{\tau}\langle\hat{v}^{1/3}\rangle}{\langle\hat{v}^{1/3}\rangle}-\frac{\partial_{\tau}\langle\hat{v}^{-1}\rangle}{\langle\hat{v}^{-1}\rangle}\right),\quad$
(76)

where $\langle\hat{v}^{1/3}\rangle=\alpha_{o}^{-1/3}\langle\hat{a}\rangle$ and
$\langle\hat{v}^{-1}\rangle=\alpha_{o}\langle\hat{a}^{-3}\rangle$ and
expectation values are computed with respect to the backreacted background
states. To the first-order quantum corrections, one can show that the
backreacted dressed Hubble rate (76) reduces to222The details of the
calculations are presented in Appendix A.

$\displaystyle\frac{\partial_{\tau}\tilde{a}}{\tilde{a}}\,\approx\,\langle\hat{H}\rangle\,=\,\frac{1}{3}\langle\widehat{:\dot{v}/v:}\rangle.\quad$
(77)

The right-hand side of the equation above is given by [26]

$\displaystyle\langle\hat{H}\rangle$ $\displaystyle=$
$\displaystyle\frac{i}{3\hbar}\big{\langle}\hat{v}^{-1/2}\big{[}\hat{{H}}_{\rm
grav},\hat{v}\big{]}\hat{v}^{-1/2}\big{\rangle}\,=\,\frac{\pi
G}{\alpha_{o}}\langle\hat{h}\rangle,\quad\quad$ (78)

where the operator $\hat{h}$ is defined in Eq. (151).

By using the central moments $C^{abc}$, defined in Eq. (154), we get

$\langle\hat{H}\rangle^{2}=\langle\hat{H}^{2}\rangle-(\pi
G/\alpha_{o})^{2}G^{002},$ (79)

where $G^{002}\equiv\langle(\delta\hat{h})^{2}\rangle$. Taking the expectation
value of the total Hamiltonian constraint (14) with respect to the background
quantum state yields

$\displaystyle\big{\langle}\widehat{:{V}^{-1}{H}_{\rm
grav}:}\big{\rangle}+\big{\langle}\widehat{:{V}^{-1}{H}_{T}:}\big{\rangle}+\big{\langle}\widehat{:{V}^{-1}{H}_{\mathbf{k}}:}\big{\rangle}=0.\quad\quad$
(80)

Also, an energy density related to backreaction can be defined as

$\rho_{k}:=\big{\langle}\widehat{:{V}^{-1}{H}_{\mathbf{k}}:}\big{\rangle}=\big{\langle}\widehat{:{V}^{-1}\epsilon^{n}_{\mathbf{k}}:}\big{\rangle}=N_{k}\hbar
k\big{\langle}\widehat{V^{-\frac{4}{3}}}\big{\rangle}.$ (81)

Then, from Eqs. (79)-(81) we get

$\displaystyle\langle\hat{\rho}_{T}\rangle+N_{k}\hbar
k\big{\langle}\widehat{V^{-\frac{4}{3}}}\big{\rangle}=\frac{3\pi
G}{2\alpha_{o}^{2}}\,\langle\hat{r}\rangle=:\langle\hat{\rho}\rangle,$ (82)

where $\hat{r}$ is defined in Eq. (151) and $\langle\hat{\rho}\rangle$ is the
expectation value of the total energy density operator, which turns out to be
the sum of the energy densities of dust field and backreaction.

From Eq. (78) we have that $\langle\hat{H}^{2}\rangle=(\pi
G/\alpha_{o})^{2}\langle\hat{h}^{2}\rangle$. Now, by setting this into Eq.
(79), we get

$\displaystyle\langle\hat{H}\rangle^{2}$ $\displaystyle=$
$\displaystyle\frac{8\pi
G}{3}\langle\hat{\rho}\rangle\left(1-\frac{\langle\hat{\rho}\rangle}{\rho_{\rm
cr}}\right)-\frac{2\pi G}{3}\rho_{\rm cr}\big{[}G^{002}+4G^{020}\big{]},$ (83)

where $\rho_{\rm cr}\equiv 3\pi G/(2\alpha_{o}^{2})$ [24] and $G^{abc}$’s are
defined in terms of the expectation values of the central moments given in Eq.
(154). Then, by substituting $\langle\hat{\rho}\rangle$ from Eq. (82) into Eq.
(83), the dressed Hubble rate Eq. (77) to the leading order terms becomes

$\displaystyle\tilde{H}^{2}$ $\displaystyle\approx\frac{8\pi
G}{3}\left(\langle\hat{\rho}_{T}\rangle+\rho_{k}\right)\left[1-\frac{\langle\hat{\rho}_{T}\rangle+\rho_{k}}{\rho_{\rm
cr}}\right].$ (84)

The modified Friedmann equation (84) for the dressed metric
$(\tilde{N},\tilde{a})$ will be sufficient for our purpose in the rest of this
paper. Interestingly, the energy density of the backreaction behaves as a
radiation fluid. Clearly, a quantum bounce still occurs at the collapse final
state. However, this bounce may occur much earlier, because, in the presence
of a radiationlike backreaction, the total energy density of the collapse
grows faster and the critical energy density will be reached earlier than the
case where the backreaction is absent. It should be noticed that, the modified
Friedmann equation above is the evolution equation for the perturbed
background which is explored by the $\mathbf{k}$th mode of the scalar
perturbation.

##  V Exterior geometry: emergence of a rainbow black hole

To explore the whole spacetime structure of the herein model of gravitational
collapse, we need to find a suitable exterior geometry to be matched with the
interior dressed spacetime at the boundary of dust. If the pressures at the
boundary of the cloud vanish (as for a pure dust model), then it is always
possible to match the interior collapsing spacetime with an empty
Schwarzschild exterior. However, in the present model, an effective nonzero
pressure emerges due to the radiationlike behavior, induced by backreaction,
and the LQG effects (included in terms proportional to $1/\rho_{\rm cr}$). In
the following, we will therefore match the interior region with a generic
nonstatic exterior presented by a generalized Vaidya geometry which
establishes a radiating generalization of the Schwarzschild geometry [33, 34].

###  V.1 Matching conditions

Let us rewrite the interior metric as

$\displaystyle
d\tilde{s}_{-}^{2}=-d\tau^{2}+\tilde{a}^{2}(\tau)dr^{2}+r^{2}\tilde{a}^{2}(\tau)d\Omega^{2}.$
(85)

Now, we define the matching surface $\Sigma$ as the (interior) boundary shell,
$\partial M^{-}$, of the collapsing cloud with the radius $r=r_{\rm b}=\rm
const$. The unit normal vector to the matching surface, $\Sigma$, is
$n^{a}_{-}=\tilde{a}^{-1}(\partial/\partial r)^{a}$, and the induced metric on
$\Sigma$ is given by

$\displaystyle h_{ab}^{-}dx^{a}dx^{b}$ $\displaystyle=$
$\displaystyle-d\tau^{2}+r^{2}\tilde{a}^{2}d\Omega^{2}.$

The extrinsic curvature $K_{ab}=\frac{1}{2}{\cal L}_{n}h_{ab}$ for $\Sigma$ at
the interior boundary is obtained by

$\displaystyle K_{ab}^{-}$ $\displaystyle=$
$\displaystyle\tfrac{1}{2}\left(n_{-}^{c}\partial_{c}h_{ab}^{-}+h_{cb}^{-}\partial_{a}n_{-}^{c}+h_{ac}^{-}\partial_{b}n_{-}^{c}\right).$

Now, let the exterior geometry have a Vaidya form333Here the mass $\tilde{M}$,
unlike the mass given by the Schwarzschild metric, is not necessarily a
constant. Then, our choice of the metric may constitute the simplest non-
static generalization of the non-radiative Schwarzschild solution to an
effective Einstein’s field equation. Therefore, we consider a generalized
Vaidya metric in which the mass parameter is extended from a constant to a
function of the corresponding null coordinates $u,X$, as $\tilde{M}(u,X)$ [35,
36, 37]., which in Eddington-Finkelstein coordinates its line element reads

$ds^{2}_{+}=-F(u,X)du^{2}+2dudX+X^{2}d\Omega^{2},\quad$ (86)

where $F(u,X)$ is the boundary function given by

$F(u,X)=1-\frac{2G\tilde{M}(u,X)}{X}\,$ (87)

and $\tilde{M}(u,X)$ is the generalized Vaidya mass. We take only a region
with $X>X(u)$ as an exterior region of the collapsing cloud, to be matched to
the interior dressed FLRW geometry. Once the relations for $u(\tau)$ and
$X(\tau)$ at $r_{\rm b}$ are derived, the form of matching surface in the
exterior region, $\mathcal{M}^{+}$, is determined. So, we consider a general
matching surface for $\mathcal{M}^{+}$ being parametrized by
$(u(\tau),X(\tau))$ in terms of $\tau$, which has to be identified with the
interior proper time later. The metric on this surface reads

$h_{ab}^{+}dx^{a}dx^{b}=-\Big{[}F(\partial_{\tau}u)^{2}-2(\partial_{\tau}u)(\partial_{\tau}X)\Big{]}d\tau^{2}+X^{2}d\Omega^{2}.$
(88)

The unit normal vector components on this surface are given by

$\displaystyle n_{+}^{u}$ $\displaystyle=$
$\displaystyle\frac{1}{\sqrt{1-2G\tilde{M}/X-2(\partial_{\tau}X)/(\partial_{\tau}u)}}\
,$ $\displaystyle n_{+}^{X}$ $\displaystyle=$
$\displaystyle\frac{1-2G\tilde{M}/X-(\partial_{\tau}X)/(\partial_{\tau}u)}{\sqrt{1-2G\tilde{M}/X-2(\partial_{\tau}X)/(\partial_{\tau}u)}}\
.$ (89)

Similarly to the interior surface, using these vectors, we can find the
explicit components of the extrinsic curvature for the exterior surface.

Following the relations above, we are prepared now to formulate the junction
conditions at $\Sigma$:

  1. i)

From the condition $h_{\theta\theta}^{-}=h_{\theta\theta}^{+}$, we obtain

$\displaystyle X(\tau)|_{\Sigma}=r_{\rm
b}\tilde{a}(\tau)\equiv\tilde{R}(\tau).$ (90)

That is, for a given interior scale factor, $\tilde{a}(\tau)$, we can
determine $X$ on the exterior matching surface for the given shell $r_{\rm
b}$.

  2. ii)

From the condition $h_{\tau\tau}^{-}=h_{\tau\tau}^{+}$, a differential
equation is found for $X(u)$ as

$\displaystyle(dX/du)^{2}=(\partial_{\tau}X)^{2}\big{(}1-2G\tilde{M}/X-2dX/du\big{)},$
(91)

which determines the relations between the exterior coordinates (for the given
shell $r_{\rm b}$).

  3. iii)

The condition $K_{\tau\tau}^{-}=K_{\tau\tau}^{+}$ leads to a differential
equation for $X(u)$, as

$\quad\quad K_{uu}^{+}+2K_{uX}^{+}(dX/du)+K_{XX}^{+}(dX/du)^{2}=0,$ (92)

which turns out to be automatically satisfied when given the other junction
conditions.

  4. iv)

From matching $K_{\theta\theta}^{-}=K_{\theta\theta}^{+}$, by setting
$X=r_{\rm b}\tilde{a}$, an equation is obtained as

$\quad\quad\quad\Big{(}1-\frac{2G\tilde{M}}{X}-\frac{dX}{du}\Big{)}^{2}=1-\frac{2G\tilde{M}}{X}-2\frac{dX}{du},$
(93)

by expanding of which, and using Eq. (91), we obtain

$\displaystyle\frac{2G\tilde{M}}{X}\Big{|}_{\Sigma}=(\partial_{\tau}X)^{2}\big{|}_{\Sigma}=r_{\rm
b}^{2}(\partial_{\tau}\tilde{a})^{2}.$ (94)

Thus, we obtain a relation between the generalized Vaidya mass function
$\tilde{M}(u,X)$ and the $\tau$-time derivative of the interior dressed scale
factor $\tilde{a}$ at $\Sigma$.

The right-hand side of Eq. (94) is already determined by the modified
(dressed) Hubble parameter in Eq. (84). However, it would be more convenient
to rewrite the exterior mass function, $\tilde{M}$, in terms of a quantum mass
$\hat{M}$ induced by the interior quantum-gravity-inspired operators. We will
evaluate such a quantum mass in the what follows.

From Eqs. (90) and (71), the physical radius of the collapse reads

$\tilde{R}\,=\,r_{\rm b}\tilde{a}\,\approx\,r_{\rm
b}\langle\hat{V}\rangle^{1/3}.$ (95)

Note that, the time dependence in $\tilde{a}$ is encoded in the expectation
value $\langle\hat{V}\rangle^{1/3}$ with respect to the perturbed state
$\Psi_{1}$. Using this, we can write the physical volume of the spherical
cloud as

$\tilde{V}=(4\pi/3)\tilde{R}^{3}=(4\pi/3)r_{\rm
b}^{3}\,\langle\hat{V}\rangle.$ (96)

We note that, in classical theory, we have $\tilde{V}=V=\ell^{3}a^{3}$, where
$q(k,T)=1$ and $\tilde{a}=a$, so that $(4\pi/3)r_{\text{b}}^{3}=\ell^{3}$.
Recall that we set $\ell=1$ throughout previous sections, which yields
$r_{\text{b}}=(4\pi/3)^{-1/3}$.

Let $M$ be the classical mass of the collapsing cloud given by $M=\rho_{T}V$.
A quantum mass $\hat{M}$ can then be introduced for the collapsing cloud,
generated now by the background quantum (dust) matter source plus quantum
corrections induced by backreaction. Since both the matter and geometry are
quantized in quantum gravity, the quantized mass is given by
$\hat{M}=\widehat{\rho V}=-\hat{{H}}_{\rm grav}$. Now, from Eq. (184), we
write the expectation value of the quantum mass as

$\langle\hat{M}\rangle=\sum_{n=0}^{\infty}\beta_{n}\langle\hat{v}\rangle^{-n}\Big{[}\langle\hat{\rho}\rangle\langle\hat{{V}}\rangle\,G^{n00}+\rho_{\rm
cr}\langle\hat{V}\rangle G^{n10}\Big{]}.$ (97)

It turns out that, up to zeroth order (i.e., $n=0$), the expectation value of
the total quantum mass, $\langle\hat{M}\rangle_{0}$, is obtained from
expectation value of the quantum mass associated to the quantized dust cloud,
and that of the quantum backreaction of each scalar mode, $\mathbf{k}$. It
reads

$\displaystyle\langle\hat{M}\rangle_{0}\,:=\,\langle\hat{\rho}\rangle\langle\hat{V}\rangle=\langle\hat{M}_{T}\rangle+M_{k}\,,$
(98)

where
$\langle\hat{M}_{T}\rangle=\langle\hat{\rho}_{T}\rangle\langle\hat{{V}}\rangle$
is the mass of the dust field and the second term on the right-hand side is
the mass generated by backreaction of the $\mathbf{k}$th mode,

$M_{k}=\rho_{k}\tilde{V}=\frac{r_{\rm b}N_{k}\hbar k}{\tilde{R}}.$ (99)

In leading-order terms, by ignoring the moments $G^{abc}$, the time evolution
of the backreacted, dressed scale factor $\tilde{a}$ reads
$(\partial_{\tau}\tilde{a})^{2}=\tilde{a}^{2}\tilde{H}^{2}\approx\langle\hat{V}\rangle^{2/3}\tilde{H}^{2}$,
where $\tilde{H}$ is given by Eq. (84). Having that $r_{\rm
b}^{2}(\partial_{\tau}\tilde{a})^{2}=(\partial_{\tau}X)^{2}|_{\Sigma}$ [cf.
Eq. (94)], we get

$\displaystyle(\partial_{\tau}X)^{2}|_{\Sigma}\,\approx\,\frac{2G}{\tilde{R}}\left(\langle\hat{M}_{T}\rangle+M_{k}\right)\left[1-\frac{3}{4\pi\rho_{\rm
cr}}\frac{\langle\hat{M}_{T}\rangle+M_{k}}{\tilde{R}^{3}}\right].$ (100)

For convenience, in what follows we will set $\tilde{R}\approx r_{\rm
b}\langle\hat{V}\rangle^{1/3}\equiv q(k,T)\bar{R}$ to distinguish the area
radius $\tilde{R}$ in the presence of the backreaction from $\bar{R}\equiv
r_{\rm b}\langle\hat{V}\rangle_{o}^{1/3}$, in which the backreaction is absent
(i.e., $q(k,T)=1$). So, ”tilde” and ”bar” refer to high-energy and low-energy
observers, respectively. Now, from Eq. (94) we obtain an effective (dressed)
mass $\tilde{M}_{\text{b}}$:

$\displaystyle\tilde{M}_{\text{b}}$
$\displaystyle=\langle\hat{M}\rangle_{0}\left[1-\frac{3}{4\pi\rho_{\rm
cr}}\frac{\langle\hat{M}\rangle_{0}}{\tilde{R}^{3}}\right]$
$\displaystyle=\left(\langle\hat{M}_{T}\rangle+\frac{r_{\rm b}N_{k}\hbar
k}{\tilde{R}}\right)\left[1-\frac{3}{4\pi\rho_{\rm
cr}}\left(\frac{\langle\hat{M}_{T}\rangle}{\tilde{R}^{3}}+\frac{r_{\rm
b}N_{k}\hbar k}{\tilde{R}^{4}}\right)\right].$ (101)

From Eq. (101) it is clear that, unlike the relativistic collapse of a dust
matter, the effective mass $\tilde{M}_{\text{b}}$ here is not a constant and
depends on the (dressed backreacted) physical radius $\tilde{R}$ (at the
boundary $\Sigma$) and the mode $k$ of the scalar perturbation. This is a
consequence of the fact that the energy density growth of the (semiclassical
interior) matter cloud is accompanied by a negative pressure which does not
vanish on the boundary surface.

It should be noticed that $\tilde{M}_{\text{b}}$ given in Eq. (101) is not the
total mass in the Vaidya region (i.e., $X\geq\tilde{R}$) but the mass on the
boundary surface $\Sigma$. That is why we have denoted it by a subscript ”b”.
Therefore, to find the total Vaidya mass $\tilde{M}(u,X)$, one needs to
determine the exact form of the energy-momentum tensor, satisfying the
(effective) Einstein field equations and the energy conditions in the exterior
region. Indeed, the reason for calling this exterior geometry a generalized
Vaidya is that, in a region sufficiently close to the boundary, we expect that
quantum gravitational effects could be effectively described by an effective
energy-momentum tensor in GR.

###  V.2 Rainbow horizons at the boundary surface

Let us investigate a qualitative behavior of horizons in the vicinity of the
matter shells due to the quantum gravity effects. The exterior region can be a
generalized Vaidya spacetime, supported by finite density and pressures, which
vanishes rapidly at large distances. Such a spacetime can be consistently
matched with a quantum-gravity-corrected interior region using the discussed
boundary conditions. These conditions are effectively satisfied in our model
for the dust collapse in the quantum regime. The question of effective
dynamics of the exterior spacetime will be answered in the subsequent
subsection. Here we address the formation of dynamical horizon, when we know
only the dynamics of matching surface $\Sigma$.

When the condition $F(u,\tilde{R})=0$ holds in Eq. (87), which is equivalent
to $2G\tilde{M}_{\rm b}=X$ for the known mass function (101), a dynamical
horizon forms and intersects the boundary surface $\Sigma$ [38, 39]. From the
point of view of the interior spacetime parameters, the right-hand side of Eq.
(94) is already determined by the modified (dressed) Hubble parameter in Eq.
(84). Thus, when $\partial_{\tau}\tilde{a}$ reaches the value
$|\partial_{\tau}\tilde{a}|=r_{\rm b}^{-1}$, a dynamical horizon forms. In
classical GR, $\partial_{\tau}\tilde{a}$ is unbounded and diverges at the
singularity where $a\to 0$ (blue and orange curve in Fig. 1). However, in the
presence of quantum gravity effects, $\partial_{\tau}\tilde{a}$ changes from a
finite initial condition and vanishes at some point where a quantum bounce
occurs (red and green curves in Fig. 1).

Figure 1: Qualitative behavior of Eq. (84) is plotted. Horizontal lines
correspond to different values of $r_{\rm b}$. The colored curves represent
different effective collapse scenarios; a situation where neither backreaction
nor loop correction are applied (orange curve); the case where only loop
correction is present (green curve); the case where only backreaction effect
is considered (blue curve); and the case where both backreaction and loop
corrections (i.e., both quantum gravity effects) are applied (red curve).

In this case, there is a turning point in $|\dot{a}|$. Based on the initial
values for $r_{\rm b}$ (i.e. $r_{\rm b1}$, $r_{\rm b2}$, etc.), a horizon may
or may not form. In the presence of backreaction (red curve), the threshold
$a_{\rm th}$ for the horizon formation, given from $|\dot{a}|_{\rm max}(a_{\rm
th})=r_{\rm b}^{-1}$, is changed and would happen earlier than the case where
the backreaction is ignored (green curve). This is depicted in Fig. 1,
$a_{2}>a_{1}$, where $a_{2}$ is the threshold in the presence of both
backreaction and loop corrections (i.e., terms proportional to $1/\rho_{\rm
cr}$), whereas $a_{1}$ is the threshold for the case of a pure loop
correction.

###  V.3 Exterior rainbow geometry

So far we have found the boundary mass $\tilde{M}_{\text{b}}(\tilde{R})$ in
the Vaidya region due to matching with the interior spacetime [cf. Eq. (101)].
However, in order to specify the exterior geometry, the total mass
$\tilde{M}(u,X)$ should be exactly determined in the exterior Vaidya region
(i.e., $X\geq\tilde{R}$). This requires the knowledge about the modified
Einstein field equations in the exterior region.

Let us assume that the (effective) energy-momentum tensor in the exterior
Vaidya region can be written as [40]

$\displaystyle T_{\mu\nu}=\sigma
N_{\mu}N_{\nu}+(\varrho+p)(N_{\mu}L_{\nu}+N_{\nu}L_{\mu})+pg_{\mu\nu},$ (102)

with the help of two null vectors,

$\displaystyle N_{\mu}=\delta_{\mu}^{\,0}\quad\text{and}\quad
L_{\nu}=\tfrac{1}{2}F(u,X)\delta_{\nu}^{\,0}+\delta_{\nu}^{\,1},\quad\quad$
(103)

such that, $N_{\lambda}N^{\lambda}=L_{\lambda}L^{\lambda}=0$ and
$N_{\lambda}L^{\lambda}=-1$. Here, $\varrho$ and $p$ are, respectively, the
energy density and pressure associated to $T_{\mu\nu}$, the effective energy-
momentum tensor in the exterior region that should satisfy the (effective)
Einstein field equations. Using these parameters, the nonvanishing components
of Einstein field equations are

$\displaystyle\sigma(X,u)$ $\displaystyle=$ $\displaystyle\frac{2}{\kappa
X^{2}}\frac{d\tilde{M}}{du},$ (104) $\displaystyle\varrho(X,u)$
$\displaystyle=$ $\displaystyle\frac{2}{\kappa X^{2}}\frac{d\tilde{M}}{dX},$
(105) $\displaystyle p(X,u)$ $\displaystyle=$ $\displaystyle-\frac{1}{\kappa
X}\frac{d^{2}\tilde{M}}{d^{2}X}.\quad\quad$ (106)

If we project $T_{\mu\nu}$ to the following orthonormal basis:

$\displaystyle\textbf{e}_{(a)}^{~{}~{}\mu}=\begin{pmatrix}\frac{1}{\sqrt{2}}\left(\frac{3}{2}-\frac{\tilde{M}}{X}\right)&-\frac{1}{\sqrt{2}}&0&0\\\
\frac{1}{\sqrt{2}}\left(\frac{3}{2}+\frac{\tilde{M}}{X}\right)&\frac{1}{\sqrt{2}}&0&0\\\
0&0&\frac{1}{X}&0\\\ 0&0&0&\frac{1}{X\sin\theta}\end{pmatrix},\quad\quad$
(107)

where
$T_{(a)(b)}=T_{\mu\nu}\;\textbf{e}^{~{}~{}\mu}_{(a)}\;\textbf{e}^{~{}~{}\nu}_{(b)}$,
we find,

$\displaystyle
T^{(a)(b)}=\begin{pmatrix}\frac{\sigma}{2}+\varrho&\frac{\sigma}{2}&0&0\\\
\frac{\sigma}{2}&\frac{\sigma}{2}-\varrho&0&0\\\ 0&0&p&0\\\
0&0&0&p\end{pmatrix}.\quad\quad$ (108)

Assuming an equation of state (EOS) $p=w\varrho$ with $w=\rm const.$, for the
fluid and replacing it into Eqs. (105) and (106), we find the following
solution for the mass function:

$\displaystyle\tilde{M}\left(u,X\right)=\alpha F_{1}(u)X^{\beta}+F_{2}(u),$
(109)

where $\beta\neq 1/2$ and $F_{1}(u),\,F_{2}(u)$ are two arbitrary functions.
Constants $\beta$ and $\alpha$ are related to each other by

$\displaystyle\beta(w)=1-2w=:\frac{1}{\alpha(w)}\,.\quad\ $ (110)

Using the mass function (109) along with Eqs. (104)-(106), the effective
profiles for the energy-momentum tensor take the form

$\displaystyle\sigma(u,X)$ $\displaystyle=$ $\displaystyle\frac{2}{\kappa
X^{2}}\left(\frac{{dF}_{2}(u)}{du}+\alpha
X^{\beta}\frac{d{F}_{1}(u)}{du}\right),$ (111) $\displaystyle\varrho(u,X)$
$\displaystyle=$ $\displaystyle\frac{2}{\kappa}\,F_{1}(u)X^{\beta-3},\quad$
(112) $\displaystyle p(u,X)$ $\displaystyle=$
$\displaystyle\frac{1-\beta}{\kappa}\,F_{1}(u)X^{\beta-3}.\quad\quad\quad\quad$
(113)

Arbitrary functions $F_{1}(u)$ and $F_{2}(u)$ can be obtained through matching
conditions and should be chosen such that the energy conditions are satisfied.
Once they are found, a physical energy-momentum tensor is established for the
Vaidya region.

In the absence of quantum effects (i.e., $\rho_{\rm crit}\to 0$) and
backreaction (i.e., $k=0$), the exterior solution should reduce to the
standard Schwarzschild geometry. This corresponds to a vacuum region with a
vanishing energy-momentum tensor, requiring that $\sigma=\varrho=p=0$, thus,

$\displaystyle\alpha
X^{\beta}\frac{d{F}_{1}(u)}{du}+\frac{d{F}_{2}(u)}{du}=0,$ (114a)
$\displaystyle F_{1}(u)X^{\beta-3}=0,$ (114b) $\displaystyle
F_{1}(u)(1-\beta)X^{\beta-3}=0.$ (114c)

The above equations have solution if $F_{1}(u)$ and $d{F}_{2}(u)/du$ vanish in
the relativistic limit.

On the other hand, as mentioned above, the energy conditions will establish
additional constraints on $F_{1}$ and $F_{2}$ [41, 40]. In particular, the
weak energy condition ($\sigma\geq 0,\,\varrho\geq 0,\,\varrho+p\geq 0$)
requires that

$F_{1}(u)\geq 0\quad\text{and}\quad\beta\leq 2.$ (115)

The weak and strong energy conditions together ($\sigma\geq 0,\;\varrho\geq
0,\;p\geq 0$) lead to

$F_{1}(u)\geq 0\quad\text{and}\quad\beta\leq 1.$ (116)

Finally, the dominant energy condition ($\sigma\geq 0,\;\varrho\geq
0,\;\varrho\geq|p|$) yields

$F_{1}(u)\geq 0\quad\text{and}\quad 0\leq\beta\leq 2.$ (117)

In summary, the energy conditions acquire the physical ranges $F_{1}(u)\geq 0$
and $0\leq\beta\leq 1$ for the arbitrary parameters $F_{1}$ and $\beta$ of the
model.

On the boundary surface $\Sigma$, where $X=\tilde{R}=r_{\text{b}}\tilde{a}$,
we have

$\displaystyle\tilde{M}(u,X)|_{\Sigma}$
$\displaystyle=\tilde{M}_{\text{b}}(\tilde{R}).$ (118)

This yields

$\displaystyle\alpha F_{1}(u)\tilde{R}^{\beta}+F_{2}(u)$
$\displaystyle=\langle\hat{M}\rangle_{0}\left(1-\frac{3}{4\pi\rho_{\rm
cr}}\frac{\langle\hat{M}\rangle_{0}}{\tilde{R}^{3}}\right),$ (119)

where, $\langle\hat{M}\rangle_{0}$ is given by Eq. (98):

$\langle\hat{M}\rangle_{0}=\langle\hat{M}_{T}\rangle+\frac{r_{\rm b}N_{k}\hbar
k}{\tilde{R}}\,.$ (120)

Furthermore, in the absence of backreaction (when $k=0$) and quantum gravity
effects (where $\rho_{\rm cr}\to\infty$), the standard Schwarzschild solution
should be retrieved. This, together with Eqs. (119) and (120) gives the the
unknown functions as

$\displaystyle F_{1}(u)$
$\displaystyle=\frac{\beta}{\tilde{R}^{\beta}}\left(\tilde{M}_{\rm
b}-\langle\hat{M}_{T}\rangle\right)$
$\displaystyle=\frac{\beta\langle\hat{M}\rangle_{0}}{\tilde{R}^{\beta}}\left(1-\frac{3}{4\pi\rho_{\rm
cr}}\frac{\langle\hat{M}\rangle_{0}}{\tilde{R}^{3}}\right)-\frac{\beta\langle\hat{M}_{T}\rangle}{\tilde{R}^{\beta}}\,,$
(121) $\displaystyle F_{2}(u)$ $\displaystyle=\langle\hat{M}_{T}\rangle.$
(122)

Once the free functions $F_{1}$ and $F_{2}$ are found, we can extract the
fluid profiles $\sigma,\varrho,p$ and the mass function $\tilde{M}(u,X)$ in
the Vaidya region.

For the $\mathbf{k}$th mode of the interior scalar field, the mass
$\langle\hat{M}\rangle_{0}$ [cf. Eq. (120)] clearly depends on $k$. Likewise,
the boundary radius $\tilde{R}$ depends implicitly on $k$. Therefore, $F_{1}$,
which is a function of $\langle\hat{M}\rangle_{0}$ and $\tilde{R}$, will also
depend on the mode $k$. In other words, for a fixed mode $k$ on the interior
region of the collapsing ball, $F_{1}$ is unique and depends on the value of
$k$, so different $k$’s acquire distinct functions $F_{1}$. Thereby, the fluid
profiles (111)-(113), defined as functions of $F_{1}$ and $F_{2}$, are
uniquely obtained for a fixed value of $k$, as

$\displaystyle\sigma_{k}$ $\displaystyle=$ $\displaystyle\frac{2}{\kappa
X^{2}}\frac{d\langle\hat{M}_{T}\rangle}{du}+\frac{2\alpha\beta}{\kappa
X^{2-\beta}}\frac{d}{du}\left(\frac{\tilde{M}_{\rm
b}(k)-\langle\hat{M}_{T}\rangle}{\tilde{R}^{\beta}}\right),$ (123)
$\displaystyle\varrho_{k}$ $\displaystyle=$
$\displaystyle\frac{2\beta}{\kappa\tilde{R}^{\beta}}\,\,\frac{\tilde{M}_{\rm
b}(k)-\langle\hat{M}_{T}\rangle}{X^{3-\beta}},\quad$ (124) $\displaystyle
p_{k}$ $\displaystyle=$
$\displaystyle\frac{\beta(1-\beta)}{\kappa\tilde{R}^{\beta}}\,\frac{\tilde{M}_{\rm
b}(k)-\langle\hat{M}_{T}\rangle}{X^{3-\beta}}.\quad\quad\quad\quad$ (125)

These induced profiles establish a unique energy-momentum tensor (102) in the
exterior Vaidya region as

$\displaystyle
T_{\mu\nu}^{(k)}=\sigma_{k}N_{\mu}N_{\nu}+(\varrho_{k}+p_{k})(N_{\mu}L_{\nu}+N_{\nu}L_{\mu})+p_{k}g_{\mu\nu}^{(k)}.$
(126)

This implies that, for matching between the inner FLRW region and the outer
Vaidya region to be unique at the boundary, associated with each interior mode
$k$ there should exist a unique energy-momentum tensor $T_{\mu\nu}^{(k)}$ in
the exterior region. Interestingly, for the specified $k$, the exterior metric
$g_{\mu\nu}^{(k)}$, which appeared in the last term of (126), is a unique
solution to the Einstein field equation outside the dust ball. We will find
this unique exterior metric in what follows.

Using the derived $F_{1}(u)$ and $F_{2}(u)$ in Eq. (109), the mass function
$\tilde{M}(u,X)$ in the Vaidya region is achieved.444There can be other
choices for the fluid, such as one with EoS $p=w\varrho^{\gamma}$, which carry
the interior quantum gravity corrections to the exterior Vaidya region. In
that case mass function $\tilde{M}\left(u,X\right)$ will be a polynomial of
$X$. In the present model, we are interested in considering gravitational
collapse of a radiating fluid for the exterior spacetime, so the choice
$\gamma=1$ suffices for our purpose. It is clear that the loop correction,
encoded in terms proportional to $1/\rho_{\rm cr}$, appears in the first term
of $F_{1}(u)$, whereas the backreaction effect is included in
$\langle\hat{M}\rangle_{0}$ and $\tilde{R}$. Having found $\tilde{M}(u,X)$,
the exterior function (87) takes the form

$\displaystyle
F_{k}(u,X)=1-\frac{2G}{X}\langle\hat{M}_{T}\rangle-\frac{2G}{X^{1-\beta}}\tilde{R}^{-\beta}\left[\langle\hat{M}\rangle_{0}\left(1-\frac{3}{4\pi\rho_{\rm
cr}}\frac{\langle\hat{M}\rangle_{0}}{\tilde{R}^{3}}\right)-\langle\hat{M}_{T}\rangle\right].$
(127)

As expected, the boundary function $F_{k}(u,X)$, induced by the fluid profiles
in the Vaidya region, is uniquely defined for a given value of $k$ and the
fluid EOS, $\beta$. Note that, in order to show this $k$ dependency clearly,
we have set the subscript $k$ in $F_{k}(u,X)$.

By changing the time coordinate through $du=dt+dX/F_{k}(u,X)$ in Eq. (86), the
exterior generalized Vaidya metric becomes

$\displaystyle
ds^{2}_{+}=-F_{k}(t,X)dt^{2}+F_{k}^{-1}(t,X)dX^{2}+X^{2}d\Omega^{2}.\quad$
(128)

Here, the physical radius $X$ belongs to the interval
$r_{\text{b}}\left(\alpha_{o}v_{\rm m}(k)\right)^{1/3}\leq X<+\infty$, where
$v_{\rm m}(k)$ is the minimum volume of the collapse at the quantum bounce. A
physical interpretation of the emergent exterior geometry (128) is as follows.
The interior quantum gravity effects (provided by the loop correction and the
backreaction of field modes) induce a fluid with the energy-momentum
$T_{\mu\nu}$, in the exterior region of the collapsing ball. In fact, each
mode $\mathbf{k}$ of the field in the inner region induces an energy-momentum
$T_{\mu\nu}^{(k)}$ associated with the external fluid, so that there exists a
one-to-one correspondence between the discrete interior energy-momentum tensor
(of each mode $\mathbf{k}$) and the discretized induced exterior one,
$T_{\mu\nu}^{(k)}$. The full $T_{\mu\nu}$ of the exterior fluid would consist
of the full set of the interior modes $\mathbf{k}$ on the same lattice
$\mathcal{L}$. The $\mathbf{k}$th mode of the exterior fluid constructs the
$k$-dependent spacetime geometry (128) in the exterior region so that,
different modes of the external fluid experience different geometries; a
rainbow exterior spacetime background emerges.

By setting $F_{k}(t,X)=0$ in Eq. (128), one finds a set of solutions for
horizons. It turns out that, if such solutions exist for a given value of $k$
associated with the exterior fluid profile, with a specified EOS parameter
$\beta$, a particular set of horizons emerges. Then, different modes of the
fluid would experience different horizons, so that a refraction of black hole
horizons can occur in the exterior region. Therefore, a rainbow black hole can
emerge in the framework of the distant observer.

From Eq. (127) it is clear that, even at large distances far from the Planck
region ($X\gg\ell_{\rm Pl}$), where the loop effects [i.e., terms proportional
to $1/\rho_{\rm cr}$ in Eq. (127)] are negligible, a quantum gravity effect
will still exist due to backreaction of the $\mathbf{k}$th mode:

$\displaystyle
F_{k}(t,X)=1-\frac{2G}{X}\langle\hat{M}_{T}\rangle-\frac{2Gr_{\rm
b}}{\tilde{R}^{1+\beta}}\,\frac{N_{k}\hbar k}{X^{1-\beta}}\,.$ (129)

We are interested in ranges far from the bounce, where the backreaction
effects are still significant while loop corrections are negligible.
Therefore, in the rest of the paper, we will explore the physical consequences
of the solution (127) with no loop correction included.

##  VI Gravitational lensing effect

So far, we have seen that the quantum gravity effects of the interior dust
ball can be carried out to the exterior region through suitable junction
conditions. For each mode $\mathbf{k}$ of the interior matter field, a
generalized Vaidya geometry then emerges in the exterior region through an
induced fluid with energy-momentum tensor (126). We assume that, the external
matter is a photonic fluid with discrete energy-momentum $T_{\mu\nu}^{(k)}$
and an EoS $w=1/3$ (i.e., radiation). Consequently, the exterior spacetime
will be modified by the backreaction of the induced photons, thus, photons
with different energies experience different spacetime geometries555The
backreaction of each mode of photons on the exterior background changes the
spacetime metric, so different modes will induce different geometries which in
turn, can be probed by the same modes. This is similar to the refraction of
light wavelengths in a medium with refractive indices which depend on those
wavelengths. Indeed, this is a consequence of the backreaction of light modes
on the medium so that different modes of light feel different refractive
indices and hence pass through different trajectories. This leads to a rainbow
feature, as emerging from a prism.. In this section, we will investigate the
phenomenological implications of the propagation of such photons on the
exterior induced spacetime, by studying their gravitational lensing and
calculate, perturbatively, the effects of backreaction on the Einstein angle.

We are interested only in effects of backreaction on the exterior spacetime.
Therefore, we consider large distances from the Planck scale where loop
effects are negligible, whereas backreaction effects through semiclassical
gravity are still present. In this case, the exterior boundary function is
given from Eq. (129) as

$F(t,X)=1-\frac{X_{\rm S}}{X}-\frac{R_{k}^{2/3}(t)}{X^{2/3}},$ (130)

where, we have defined a Schwarzschild radius as $X_{\rm S}\equiv
2G\langle\hat{M}_{T}\rangle$ and
$R_{k}^{2/3}(t)\equiv{2r_{\text{b}}N_{k}k\ell_{\rm
Pl}^{2}/\tilde{R}^{4/3}(t)}$. Location of the apparent horizon, $X_{\rm AH}$,
can be obtained by solving $F_{k}(t,X_{\rm AH})=0$:

$\displaystyle(X-X_{\rm S})^{3}-R_{k}^{2}X=0,$ (131)

which has the following solution:

$X_{\rm AH}(k)=X_{\rm
S}+\frac{(\frac{2}{3})^{1/3}R_{k}^{2}}{\mathcal{W}(X_{\rm
S};R_{k})}+\frac{\mathcal{W}(X_{\rm S};R_{k})}{18^{1/3}}.$ (132)

where $\mathcal{W}(X_{\rm S};R_{k})\equiv\left(9X_{\rm
S}R_{k}^{2}+\sqrt{81X_{\rm S}^{2}R_{k}^{4}-12R_{k}^{6}}\right)^{1/3}$.

Now, we intend to compute the deflection angle with respect to the exterior
metric (128) with the boundary function (130). The generic geodesic equation
can be written as

$\frac{dv^{i}}{d\vartheta}+\Gamma^{i}_{jk}v^{j}v^{k}=0,$ (133)

where $v^{i}\equiv dx^{i}/d\vartheta$ is the tangent vector to the geodesic
curve and $\vartheta$ is an affine parameter. For null curves, with
$g_{ij}v^{i}v^{j}=0$, on the exterior background (128) with the boundary
function (133), the full set of geodesic equations becomes

$\displaystyle
t^{\prime\prime}+\frac{1}{F_{k}}\frac{dF_{k}}{dX}\,t^{\prime}X^{\prime}-\frac{1}{2F_{k}^{3}}\frac{dF_{k}}{dt}\,{X^{\prime}}^{2}+\frac{1}{2F_{k}}\frac{dF_{k}}{dt}\,{t^{\prime}}^{2}=0,\quad$
(134a)
$\displaystyle\varphi^{\prime\prime}+\frac{2}{R}\,R^{\prime}\,\varphi^{\prime}+2\cot\theta\,\varphi^{\prime}\,\theta^{\prime}=0,$
(134b)
$\displaystyle\theta^{\prime\prime}+\frac{2}{R}\,\theta^{\prime}\,R^{\prime}-\sin\theta\cot\theta\left(\varphi^{\prime}\right)^{2}=0,$
(134c) $\displaystyle
X^{\prime\prime}-\frac{1}{2F_{k}}\frac{dF_{k}}{dX}\left(X^{\prime}\right)^{2}+\frac{1}{2}F_{k}\frac{dF_{k}}{dX}\left(t^{\prime}\right)^{2}$
$\displaystyle
XF_{k}\left[\left(\theta^{\prime}\right)^{2}+\sin^{2}\theta\left(\varphi^{\prime}\right)^{2}\right]+\frac{1}{F_{k}}\frac{dF_{k}}{dt}\,t^{\prime}X^{\prime}=0,$
(134d)

where a prime stands for a derivative with respect to $\vartheta$.

To study lensing in the exterior spacetime, we need to make a stationarity
assumption. We assume that, at an early stage of the collapse, the crossing
time of photons, $t_{p}$, is much smaller than the timescale of variation of
the lens, $t_{l}$ (i.e., $t_{p}\ll t_{l}$). Then we will derive the lensing
observables through perturbative methods. In the weak field limit, we will
expand quantities in terms of the expansion parameters $\epsilon_{m}\equiv
X_{\rm S}/X$ and $\epsilon_{k}\equiv(R_{k}/X)^{2/3}$ so that the time
derivative terms in Eqs. (134) will become of next order; for example,
$dF_{k}/dt\propto\epsilon_{k}\,d\tilde{R}/dt\propto\epsilon_{k}\,d\tilde{a}/dt\propto\epsilon_{k}\,\epsilon_{m}^{1/2}$.
Being interested in first-order effects, i.e., terms proportional to
$\epsilon_{m}$ or $\epsilon_{k}$, we ignore terms with time derivatives in Eq.
(134). In other words, we will consider regimes in which stationary features
of spacetime are more important than nonstationary ones. Without loss of
generality, we will work on the equatorial plane $\theta=\pi/2$. Then, Eqs.
(134) reduce to [42],

$\displaystyle t^{\prime}\simeq C/F_{k},$ (135a) $\displaystyle
X^{2}\,\varphi^{\prime}=b,$ (135b) $\displaystyle
F^{-1}_{k}\left(X^{\prime}\right)^{2}+b^{2}/X^{2}-1/F_{k}\simeq-\lambda,$
(135c)

where $C$, $b$ and $\lambda$ are constants of integration. For photons, we set
$\lambda=0$, and for simplicity we set $C=1$. Putting everything together in
Eqs. (135), we get the following geodesic for photons:

$\displaystyle\varphi(X_{\rm so})-\varphi(X_{\rm ob})$
$\displaystyle=\int\,\frac{dX}{X^{2}}\left[\frac{1}{b^{2}}-\frac{1}{X^{2}}+\frac{X_{\rm
S}}{X^{3}}+\frac{R_{k}^{2/3}}{X^{8/3}}\right]^{-1/2},$ (136)

where, $X_{\rm so}$ and $X_{\rm ob}$ are the locations of source and observer,
respectively. We consider a collection of photons (as part of the fluid) which
start emitting from the source in $X_{\rm so}$, moving toward the turning
point on $X_{0}$ (the closest distance to the lens, i.e. the dust ball here),
where $dX/d\varphi=0$, and keep propagating until they reach the observer on
$X_{\rm ob}$. The deflection angle of the trajectory with respect to a
straight line can then be obtained as

$\Delta\varphi(X_{0})=|\varphi(X_{\rm so})-\varphi(X_{\rm ob})|-\pi.$ (137)

In turning point, we have

$\frac{F_{k}(X_{0})}{X_{0}^{2}}=\frac{1}{b^{2}},$ (138)

whence Eq. (136), in terms of the expansion parameters $\epsilon_{m}$ and
$\epsilon_{k}$ and a new variable $x=X_{0}/X$, becomes

$\displaystyle\varphi(X_{\rm so})-\varphi(X_{\rm ob})$
$\displaystyle=\left(\int_{x=\frac{X_{0}}{X_{\rm
so}}}^{1}dx+\int_{x=\frac{X_{0}}{X_{\rm
ob}}}^{1}dx\right)(1-x^{2})^{-\frac{1}{2}}$
$\displaystyle\quad\quad\times\left[1-\frac{\epsilon_{m}(1-x^{3})}{1-x^{2}}-\frac{\epsilon_{k}(1-x^{8/3})}{1-x^{2}}\right]^{-\frac{1}{2}}.$
(139)

Integration in the right-hand side of Eq. (139) can be done perturbatively.
Then, up to first order in $X_{0}/X_{\rm so}$, $X_{0}/X_{\rm ob}$,
$\epsilon_{m}$ and $\epsilon_{k}$, we get

$\displaystyle\Delta\varphi(X_{0})$ $\displaystyle=\varphi(X_{\rm
so})-\varphi(X_{\rm ob})-\pi$ $\displaystyle\simeq
2\,\epsilon_{m}+\epsilon_{k}\frac{\sqrt{\pi}\;\Gamma\left(11/6\right)}{\Gamma\left(4/3\right)}-(1+\epsilon_{m})\left(\frac{X_{0}}{X_{\rm
so}}+\frac{X_{0}}{X_{\rm
ob}}\right)-\frac{\epsilon_{k}}{2}\left(\frac{X_{0}}{X_{\rm
so}}+\frac{X_{0}}{X_{\rm ob}}\right)+\mathcal{O}(x^{2}).$ (140)

Deflection angle (140) contains local and nonlocal terms. Since the metric
function (130) is not asymptotically flat [43], we cannot isolate the whole
gravitational system, so it would be better to take into account the local
terms. To rewrite Eq. (140) in terms of the constant of motion $b$, we solve
Eq. (138) for $X_{0}$. Taking leading-order terms for $X_{\rm S}/b$ and
$R_{k}/b$, we get [44],

$X_{0}\simeq b\left(1-\frac{X_{\rm
S}}{2b}-\frac{R_{k}^{2/3}}{2b^{2/3}}\right).$ (141)

Then, the deflection angle becomes

$\displaystyle\Delta\varphi(k)$ $\displaystyle\simeq 2\frac{X_{\rm
S}}{b}+\left(\frac{R_{k}}{b}\right)^{2/3}\frac{\sqrt{\pi}\;\Gamma\left(11/6\right)}{\Gamma\left(4/3\right)}-b\left(\frac{1}{X_{\rm
so}}+\frac{1}{X_{\rm ob}}\right)+\mathcal{O}(x^{2}).$ (142)

Since the backreaction effects will be important in the high-energy regimes,
using the perturbative approach to solve the lens equation is not suitable.
However, we apply a perturbative treatment to get a sense of quantum gravity
effects on light rays propagating on the herein mode-dependent emergent
spacetime666For strong gravity and high energy regimes and including time
variation of spacetime into consideration, an independent study is required by
using the numerical techniques; this will be studied in a separated work
[45]..

ObserverLensA${\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\rm
Image}$$\varphi$$\theta$$X_{0}$$b$$D_{\rm ol}$$D_{\rm ls}$$D_{\rm os}$ Figure
2: Lensing configuration: deflection of different modes by a point mass, where
lens, observer and source are highly aligned. Different modes probe different
curvatures and create images at different angular positions, this chromatic
gravitational aberration can be called “ Quantum Gravitational Prism”.

It is well known that when the observer, source and lens (gravitational field)
are aligned, it gives rise to the formation of an Einstein ring. The third
term in deflection angle (142) can be canceled according to initial values for
$\varphi(X_{\rm so})$ and $\varphi(X_{\rm ob})$ at the flat spacetime, when
there are no curvature terms. In this case, we make use of lens equation
$\theta_{E}D_{\rm os}=\Delta\phi D_{\rm ls}$, in which $\theta_{E}$, $D_{\rm
os}$ and $D_{\rm ls}$ are the Einstein angle of the image, observer-source and
lens-source distances, respectively. In the configuration depicted in Fig. 2,
we look for the first-order effects, i.e., sub-leading corrections in weak
field approximation. We consider terms up to first orders in $\epsilon_{m}$
and $\epsilon_{k}$ in the deflection angle (142) and then apply a perturbative
approach to find a solution for the lens equation. The aim of this
investigation is to get a sense of the order of corrections arising from
backreactions on the Einstein angle.

In our case, the black hole spacetime can be described by the metric function
(130). Thus, angular diameter distances are different from Schwarzschild at
spatial infinity and are different from radial coordinates. To calculate the
Einstein angle, we use asymptotically flat spacetime relations $b\simeq D_{\rm
ol}\,\theta+\mathcal{O}(\epsilon_{k}^{n})$, $D_{\rm ls}\simeq X_{\rm
so}+\mathcal{O}(\epsilon_{k}^{n})$ and $D_{\rm ol}\simeq X_{\rm
ob}+\mathcal{O}(\epsilon_{k}^{n})$ which bring errors of the order
$\mathcal{O}(\epsilon_{k}^{n+1})$ in the Einstein angle. As stated before, we
are not looking for exact values, but instead, we are interested in finding
the order of corrections using a solution like
$\theta_{E}\approx\theta_{0}+\lambda_{\rm pert}\theta_{1}$ for the lens
equation [46]; $\lambda_{\rm pert}$ is the perturbation parameter that can be
read from Eq. (142) as,

$\displaystyle\Delta\varphi(k)$
$\displaystyle\simeq\frac{1}{Z}+\frac{\lambda_{\rm pert}}{Z^{2/3}}\,,$ (143)
$\displaystyle\lambda_{\rm pert}$ $\displaystyle=\left(\frac{R_{k}}{2X_{\rm
S}}\right)^{2/3}\frac{\sqrt{\pi}\;\Gamma\left(11/6\right)}{\Gamma\left(4/3\right)},$
(144)

where we introduced $Z=b/(2X_{\rm S})$. Using the lens equation, up to the
first order in perturbation parameter $\lambda_{\rm pert}$ we find

$\displaystyle\theta_{0}$ $\displaystyle=$ $\displaystyle\left(2X_{\rm
S}\frac{D_{\rm ls}}{D_{\rm os}D_{\rm ol}}\right)^{1/2},$ (145)
$\displaystyle\theta_{1}$ $\displaystyle=$
$\displaystyle\frac{1}{2}\left(\frac{D_{\rm ls}}{D_{\rm
os}}\right)^{2/3}\left(\frac{2X_{\rm S}}{D_{\rm ol}}\right)^{1/3},$ (146)

where $D_{\rm ol}$ is the observer-lens distance. Equation (146) indicates
that different modes create images at different angular positions which leads
to chromatic aberration for the gravitational lens. Considering more terms in
expansion of $b$ and $\theta$ adds next-order corrections to $\theta_{0}$ and
$\theta_{1}$. For better precision, it is more convenient to evaluate
$\Delta\varphi(k)$ numerically and use the lens equation to find the Einstein
angle; these investigations will be presented in a consequent paper [45].

###  VI.1 An astronomical example

Here, we give an example for gravitational chromatic aberration in GRBs
propagating from near a rainbow black hole spacetime. This example will show
how the free parameters in the backreaction term of metric function (130) can
be practically interpreted. Each degree of freedom of the electromagnetic
field can be treated as a massless scalar field and has an identical
backreaction term [10]. In particular, here we consider astronomical data of
Cyg X-1 for a stellar black hole candidate in our Milky Way galaxy, to study
the gravitational chromatic aberration effects within our herein model.

GRBs can be created by gravitational collapse of a star releasing energy
(typically $\sim 10^{44-47}$J) as electromagnetic waves [47]. Table 1 depicts
two sets of data, “optimistic” and “realistic” sets of expected values for
aberration effects of GRBs when they are passing by Cyg X-1. In the optimistic
part, we have presented numerical values made from observations at short
distances to Cyg X-1, a few billions of kilometers, whereas in the realistic
part, we considered observations from long distances to Cyg X-1, a few
kiloparsecs. Although we made many simplifying assumptions in our perturbative
analysis, we can still trust the order of corrections for the weak field
regime.

Optimistic set, $D_{\rm ol}\simeq 10^{-3}$pc, $\mathcal{D}:=D_{\rm ls}/D_{\rm
os}=0.5$
$E_{k}$ $\lambda_{\rm pert}$ $\theta_{E}$ (arcsec) 0 0 9.02453456 $\rm 1\,keV$
$3.29064131\times 10^{-10}$ 9.02453460 $\rm 1\,MeV$ $3.29064131\times 10^{-7}$
9.02456801
Realistic set, $D_{\rm ol}\simeq$2 kpc, $\mathcal{D}:=D_{\rm ls}/D_{\rm
os}=0.0005$
$E_{k}$ $\lambda_{\rm pert}$ $\theta_{E}$ ($\mu$ arcsec) 0 0 201.79472758 $\rm
1\,keV$ $3.29064131\times 10^{-10}$ 201.79473023 $\rm 1\,MeV$
$3.29064131\times 10^{-7}$ 201.79738211

Table 1: Optimistic and realistic values for the image positions of the
Einstein ring with the source position $\beta=0$, due to lensing by a stellar
black hole with quantum backreaction effects (a rainbow black hole with
chromatic aberration effects). The lens is the Milky Way black hole candidate
Cyg X-1 with mass $M\simeq 20\times M_{\odot}$ [48]. We take
$r_{\text{b}}/R_{\rm sch}=1$ with stationarity assumption
$\tilde{R}(t)\sim\tilde{R}\sim r_{\rm b}$ at an early stage of collapse. Here,
$E_{k}$ represents the energy of the massless particle in electron volts with
particle number $N_{k}\sim 10^{55}$; $E=0$ shows results for the classical
case where $R_{k}=0$.

The modification term $R_{k}^{2/3}\propto N_{k}\ell_{\rm Pl}^{2}$ suffers from
a huge suppression of the order of $\sim\ell_{\rm Pl}^{2}$, whereas the free
parameter $N_{k}$ can compensate the squared Planck length, and play the role
of amplification parameter that enhances the backreaction effects. In the
numerical calculations, we have taken the radius of the boundary shell,
$r_{\text{b}}$, to equal the Schwarzschild radius (note that, this is the
least value for the shell radius). Moreover, $N_{k}$ can be interpreted as the
number of photons in an adiabatic regime. Based on the energy released by
GRBs, we take the optimistic value $N_{k}\sim 10^{55}$ as an approximate value
for the number of photons in one energy bound. In Table 1, we have considered
two energy bands in keV and MeV for bursts probing a rainbow black hole. In
the optimistic case where rainbow effects can be observed from a nearby
source, modification to Einstein ring can be of the order of
$\delta\theta_{E}\sim 100\;\mu$arcsec, while for the realistic case, in which
observation is made from Earth at large distances from the source, i.e.,
distances of the order of $\sim 2$kps, rainbow effects in $\theta_{E}$ are
minuscule and at best can be of the order of $\delta\theta_{E}\sim
10^{-3}\mu$arcsec.

##  VII Conclusion and discussion

In this paper, we considered the gravitational collapse of a (homogeneous)
spherically symmetric dust cloud plus a (inhomogeneous) massless scalar
perturbation. Classically, when discarding the effects of the homogeneous
sector of the scalar field on the background, this model leads to the
formation of a Schwarzschild black hole in the exterior region. In quantum
theory in the interior region, the dust field $T$ plays the role of internal
time which represents the evolution of the physical Hamiltonian of the
gravitational system coupled to the scalar field $\phi$. For each mode of the
scalar perturbation propagating on this quantized background, the evolution
equation corresponds to evolution for the same field’s mode on an effective
dressed background. The components of this dressed metric depend on
fluctuations of the background quantum geometry. When the backreaction of the
quantum modes is taken into account, the emergent dressed background turns out
to be mode dependent; a rainbow metric (with components depending on the
energy of the field modes) emerges.

The semiclassical behavior of the interior dressed geometry was presented by
employing the quantum backreactions and higher-order corrections due to
quantum fluctuations of the spacetime geometry. The nonclassical features of
the interior spacetime were carried out to the exterior region due to
convenient matching conditions at the boundary of the dust cloud: An exterior
(nonstatic) black hole geometry [cf. Eq. (128)] could emerge whose components
depend on the mode of the induced fluid in the outer region. Properties of the
interior and the induced exterior geometries are summarized as follows.

  1. i)

At the late stage of collapse inside the dust ball, LQG effects in the
interior region are large; i.e., the terms $\propto 1/\rho_{\rm cr}$ in Eq.
(84) are dominant. This indicates that the classical singularity is removed
and is replaced by a quantum bounce at the final stage of the collapse. There
is an additional correction, $\propto k\langle V\rangle^{-4/3}$, from the
backreaction of each mode $\mathbf{k}$ of the scalar perturbation on the
interior quantum spacetime. This term represents the energy density of a
radiation fluid that appears to be dominant in very short distances, so that,
as the collapse proceeds, the total energy density of the interior region
grows faster, compared with the case where a pure dust field is considered. A
thorough numerical analysis would confirm our herein results 777These analyzes
will be presented in an upcoming paper [32]. which imply that a bounce still
occurs in our model, but the backreaction effects speed up its occurrence.

  2. ii)

This radiationlike effect leads to a unique evolution associated with each
mode inside the dust ball. Thus, for each scalar field mode in the interior
region of the collapse, a unique dressed, classical-like metric emerges. In
other words, different modes explore different backgrounds and thus, a rainbow
geometry emerges in the interior region.

  3. iii)

Outside the dust ball, a generalized Vaidya spacetime, with a suitable choice
of (external) fluids, can be matched consistently to each interior mode-
dependent spacetime. Therefore, corresponding to each mode $\mathbf{k}$ in the
interior region, there exists a unique Vaidya geometry, $g_{\mu\nu}^{+}(k)$
(labeled by the number $k$), provided by an external fluid with an energy-
momentum tensor $T^{(k)}_{\mu\nu}$ satisfying the Einstein field equations
outside the dust ball. This leads to a one-to-one correspondence between the
interior field mode $\mathbf{k}$ and an exterior fluid with profiles labeled
by $k$. As a consequence, different modes and labels, $k^{\prime}$, of the
external fluid explore different Vaidya backgrounds. This is equivalent to
saying that the components of the exterior Vaidya metric depend on the mode or
label of the external fluid, being a rainbow geometry. Such a geometry may
feature “rainbow horizons” with optical properties different from the
classical black holes.

At distances much larger than the scale of the bounce, where the loop effect
(i.e., terms proportional to $1/\rho_{\rm cr}$) is negligible, the quantum
gravity effects are still significant due to the backreaction effects, through
a term $\propto 1/X^{2/3}$ in the exterior metric (129). Assuming that the
external matter is a radiation fluid, the spacetime corresponding to each mode
of this fluid can be probed as a source of gravitational lensing; each mode
can provide its own particular Einstein’s ring, leading to a chromatic
aberration in the gravitational lensing process. Therefore, different modes
see different rings, so that a rainbowlike collection of rings can be detected
from astrophysical observation of such spacetimes (cf. Fig. 2).

Finally, we should emphasize that the results we achieved within this paper
are not limited to the LQG approach only. It can be shown that similar
conclusions might arise from other approaches to quantum gravity such as the
geometrodynamics approach.

## Acknowledgments

The work of A.P. was supported in part by the Ministry of Science, Research
and Technology of Iran. He is grateful for the support and kind hospitality of
University of Warsaw and University of Wrocław, where part of this work was
completed and he wishes to thank Professor Mohammad Nouri-Zonoz for comments
and helpful discussions concerning this project. T.P. acknowledges the support
by the Polish Narodowe Centrum Nauki (NCN) grants 2012/05/E/ST2/03308 and
2020/37/B/ST2/03604. The work of Y.T. was supported by the Research deputy of
University of Guilan. He also acknowledges the financial support from the
Polish Narodowe Centrum Nauki (NCN) through the grant 2012/05/E/ST2/03308.
J.L. was supported by the Polish Narodowe Centrum Nauki, Grant No.
2011/02/A/ST2/00300. This paper is based upon work from European Cooperation
in Science and Technology (COST) action CA18108 – Quantum gravity
phenomenology in the multi-messenger approach, supported by COST.

##  Appendix A Derivation of dressed Hubble rate

In order to compute $\partial_{\tau}\tilde{a}/\tilde{a}$ in Eq. (76), we
should look for $\partial_{\tau}\langle\hat{v}^{\alpha}\rangle$,
$\alpha=-1,1/3$. This is given by

$\partial_{\tau}\langle\hat{v}^{\alpha}\rangle=\frac{\big{\langle}\big{[}\hat{H}_{\rm
grav},\hat{v}^{\alpha}\big{]}\big{\rangle}}{-i\hbar}\,,$ (147)

in which the scalar field Hamiltonian cancels out in the right-hand side of
the equation, because it depends only on the volume operator. The main tool
available at the dynamical level is that of effective equations which describe
the evolution of expectation values for a dynamical state. Thus, quantum
fluctuations and higher moments act on the evolution of expectation values,
described in effective equations by coupling classical and quantum degrees of
freedom.

In LQC coupled with dust, the Hamiltonian $\hat{H}_{\rm grav}$ is given by

$\displaystyle\hat{H}_{\rm grav}\ =\ \frac{3\pi
G}{8\alpha_{o}}\sqrt{|\hat{v}|}\left(\hat{\mathsf{N}}^{2}-\hat{\mathsf{N}}^{-2}\right)^{2}\sqrt{|\hat{v}|}\
.\quad$ (148)

Volume operator $\hat{v}$ is defined as $\hat{v}|v\rangle=v|v\rangle$ and
$\hat{\mathsf{N}}|v\rangle=e^{i\hat{b}/2}|v\rangle=|v+1\rangle$ with
$[\hat{b},\hat{v}]=2i$. The operators $\hat{\mathsf{N}}=\exp(i\hat{b}/2)$,
$\hat{v}$ and $\hat{b}$ satisfy the relations

$\displaystyle\big{[}\hat{\mathsf{N}}^{n},\hat{v}\big{]}=-n\hat{\mathsf{N}}^{n}\
,\quad\quad\big{[}\sin^{2}(\hat{b}),\hat{v}\big{]}=2i\sin(2\hat{b})\,.\quad\quad$
(149)

By substituting $\hat{\mathsf{N}}^{2}-\hat{\mathsf{N}}^{-2}=2i\sin(\hat{b})$
in Eq. (148) we have

$\displaystyle\hat{\mathcal{H}}_{\rm grav}\ =\ -\frac{3\pi
G}{2\alpha_{o}}\sqrt{|\hat{v}|}\sin^{2}(\hat{b})\sqrt{|\hat{v}|}\,.\quad$
(150)

It is convenient to define operators

$\displaystyle\hat{r}\ :=\ \sin^{2}(\hat{b})\quad\quad{\rm
and}\quad\quad\hat{h}\ :=\ \sin(2\hat{b})\ ,$ (151)

with the following commutation relations:

$\displaystyle[\hat{r},\hat{v}]=2i\hbar\hat{h}\quad\quad{\rm
and}\quad\quad[\hat{h},\hat{v}]=4i(1-2\hat{r})\ .$ (152)

To determine the evolution equation for observables
$\langle\widehat{v^{\alpha}}\rangle$ under Hamiltonian $\langle\hat{H}_{\rm
grav}\rangle$, we make use of a background-dependent expansion method. A
combination of operators $\hat{D}\big{(}\hat{v},\hat{r},\hat{h}\big{)}$ can be
expanded as (by considering a convenient choice of symmetric ordering)

$\displaystyle\hat{D}$ $\displaystyle=$ $\displaystyle
D\big{(}\langle\hat{v}\rangle,\langle\hat{r}\rangle,\langle\hat{h}\rangle\big{)}+\sum_{a,b,c=0}^{\infty}\frac{1}{a!b!c!}\frac{\partial^{a+b+c}~{}D^{klm}}{\partial\langle\hat{v}\rangle^{a}~{}\partial\langle\hat{r}\rangle^{b}~{}\partial\langle\hat{h}\rangle^{c}}C^{abc}.\quad\quad\quad$
(153)

Note $\hat{D}\big{(}\hat{v},\hat{r},\hat{h}\big{)}$ is different than the
multiplication of expectation values of operators $\hat{v},\hat{r}$ and
$\hat{h}$; to describe the system completely, infinite central fluctuation
operators $C^{abc}$ are needed, where $C^{abc}$ are the (symmetric ordered)
central moments defined by

$\displaystyle C^{abc}$ $\displaystyle:=$
$\displaystyle(\delta\hat{v})^{a}(\delta\hat{r})^{b}(\delta\hat{h})^{c}$ (154)

and $\delta\hat{v}=\hat{v}-\langle\hat{v}\rangle\mathbb{I}$,
$\delta\hat{r}=\hat{r}-\langle\hat{r}\rangle\mathbb{I}$, and
$\delta\hat{h}=\hat{h}-\langle\hat{h}\rangle\mathbb{I}$ are fluctuations of
the operators $\hat{v}$, $\hat{r}$ and $\hat{h}$ around the background state
with expectation values $\langle\hat{v}\rangle$, $\langle\hat{r}\rangle$ and
$\langle\hat{h}\rangle$, respectively.

Now, following the definitions above, we can find the expansion of
$\langle\hat{v}^{\alpha}\rangle$ and $\langle\hat{H}_{\rm grav}\rangle$ in
terms of operators (154). To do so, we can expand any operator
$\hat{v}^{\alpha}$ as

$\displaystyle\hat{v}^{\alpha}$ $\displaystyle=$
$\displaystyle\big{(}\langle\hat{v}\rangle\mathbb{I}+\delta\hat{v}\big{)}^{\alpha}=\sum_{n=0}^{\infty}\left({\begin{array}[]{c}\alpha\\\
n\\\ \end{array}}\right)\langle\hat{v}\rangle^{\alpha-n}C^{n00},\quad\quad$
(157)

then, by taking its expectation value, we find

$\displaystyle\langle\hat{v}^{\alpha}\rangle$ $\displaystyle=$
$\displaystyle\sum_{k=0}^{\infty}\left({\begin{array}[]{c}\alpha\\\ k\\\
\end{array}}\right)\langle\hat{v}\rangle^{\alpha-k}G^{k00}\ ,$ (160)

where we have defined the moments $G^{k00}$ as

$\displaystyle G^{k00}$ $\displaystyle=$ $\displaystyle\langle
C^{k00}\rangle,\quad C^{k00}=(\delta\hat{v})^{k}.$ (161)

Moreover, in terms of moments (161), by using Eq. (160), we can now expand the
dressed scale factor as

$\displaystyle\tilde{a}$ $\displaystyle=$
$\displaystyle\alpha_{o}^{-1/3}\left(\frac{\langle\hat{v}^{1/3}\rangle}{\langle\hat{v}^{-1}\rangle}\right)^{1/4}$
(166) $\displaystyle=$
$\displaystyle\alpha_{o}^{-1/3}\langle\hat{v}\rangle^{1/3}\left(\sum_{k=0}^{\infty}\left({\begin{array}[]{c}-1\\\
k\\\
\end{array}}\right)\langle\hat{v}\rangle^{-k}G^{k00}\right)^{-1/4}\left(\sum_{m=0}^{\infty}\left({\begin{array}[]{c}\frac{1}{3}\\\
m\\\ \end{array}}\right)\langle\hat{v}\rangle^{-m}G^{m00}\right)^{1/4}.\quad$

To the zeroth order in quantum fluctuations, $k,m=0$, the dressed scale factor
$\tilde{a}$ reduces to $\alpha_{o}^{-1/3}\langle\hat{v}\rangle^{1/3}$ (it is
worth noting that expectation values are taken with respect to the backreacted
states, solutions to the full quantum Hamiltonian constraint (14) ).

In order to compute the multimoment expansion for the gravitational
Hamiltonian:

$\displaystyle\hat{H}_{\rm grav}$ $\displaystyle=$ $\displaystyle-\frac{3\pi
G}{2\alpha_{o}}\big{(}\hat{v}^{1/2}\hat{r}\hat{v}^{1/2}\big{)}\ ,$ (167)

one needs to expand the term $\hat{v}^{1/2}\hat{r}\hat{v}^{1/2}$. It should be
noted that the status of the symmetry here is similar to that given by the
symmetric operator $(\hat{v}\hat{r}+\hat{r}\hat{v})/2$ with integer powers of
$\hat{v}$ and $\hat{r}$. Therefore, we expect that the expansion of the herein
Hamiltonian operator (of LQC coupled with dust) will be symmetric
automatically and no reordering procedure is required. Despite this analogy,
in the latter case, the binomial expansion of the symmetric operator
$(\hat{v}\hat{r}+\hat{r}\hat{v})/2$ will lead to the finite terms around the
background state with expectation values $\langle\hat{v}\rangle$ and
$\langle\hat{r}\rangle$, and will be truncated to a certain order of quantum
corrections provided by
$\langle\delta\hat{v}\delta\hat{r}+\delta\hat{r}\delta\hat{v}\rangle/2$ once
their expectation values are taken. However, in our case, because of one-half
power of the volume operator, expansion of $\hat{v}^{1/2}\hat{r}\hat{v}^{1/2}$
involves infinitely many terms of binomial series. So, the Hamiltonian
operator can be expanded as

$\displaystyle\hat{H}_{\rm grav}$ $\displaystyle=-\frac{3\pi
G}{2\alpha_{o}}\sum_{a=0}^{\infty}\sum_{b=0}^{\infty}\left({\begin{array}[]{c}\frac{1}{2}\\\
a\\\ \end{array}}\right)\left({\begin{array}[]{c}\frac{1}{2}\\\ b\\\
\end{array}}\right)\langle\hat{v}\rangle^{1-a-b}\Big{[}\langle\hat{r}\rangle(\delta\hat{v})^{a+b}+(\delta\hat{v})^{a}\delta\hat{r}(\delta\hat{v})^{b}\Big{]}.$
(172)

The rhs of Eq. (172) indicates that, the Hamiltonian operator $\hat{H}_{\rm
grav}$ on the full Hilbert space is totally symmetric; that is, it constitutes
all possible reorderings of the operator $\delta\hat{v}$ on both sides of the
operator $\delta\hat{r}$. The expectation value of $\hat{H}_{\rm grav}$ can be
written now as

$\displaystyle\langle\hat{H}_{\rm grav}\rangle$ $\displaystyle=$
$\displaystyle-\frac{3\pi
G}{2\alpha_{o}}\sum_{n=0}^{\infty}\sum_{m=0}^{n}\left({\begin{array}[]{c}\frac{1}{2}\\\
m\\\ \end{array}}\right)\left({\begin{array}[]{c}\frac{1}{2}\\\ n-m\\\
\end{array}}\right)\langle\hat{v}\rangle^{1-n}\Big{[}\langle\hat{r}\rangle\big{\langle}(\delta\hat{v})^{n}\big{\rangle}+\big{\langle}(\delta\hat{v})^{m}\delta\hat{r}(\delta\hat{v})^{n-m}\big{\rangle}\Big{]}.\quad\quad$
(177)

By defining the central moments $G^{n10}$ as

$\displaystyle G^{n10}$ $\displaystyle:=$
$\displaystyle\frac{1}{\beta_{n}}\sum_{m=0}^{n}\left({\begin{array}[]{c}\frac{1}{2}\\\
m\\\ \end{array}}\right)\left({\begin{array}[]{c}\frac{1}{2}\\\ n-m\\\
\end{array}}\right)\big{\langle}C^{n10}\big{\rangle},$ (182) $\displaystyle
C^{n10}$ $\displaystyle:=$
$\displaystyle(\delta\hat{v})^{m}\delta\hat{r}(\delta\hat{v})^{n-m},$ (183)

together with $G^{n00}$ [defined in Eq. (161)], we can describe the
expectation value of the Hamiltonian of the quantum system completely by

$\displaystyle\langle\hat{H}_{\rm grav}\rangle=-\frac{3\pi
G}{2\alpha_{o}}\sum_{n=0}^{\infty}\beta_{n}\langle\hat{v}\rangle^{1-n}\Big{[}\langle\hat{r}\rangle
G^{n00}+G^{n10}\Big{]},\quad\quad$ (184)

where $\beta_{n}$ is a normalization constant defined by

$\displaystyle\beta_{n}$ $\displaystyle:=$
$\displaystyle\sum_{m=0}^{n}\tilde{\beta}_{nm}\,;\
\quad\tilde{\beta}_{nm}=\left({\begin{array}[]{c}\frac{1}{2}\\\ m\\\
\end{array}}\right)\left({\begin{array}[]{c}\frac{1}{2}\\\ n-m\\\
\end{array}}\right).\quad\quad$ (189)

Now, following Eq. (147), in order to obtain the time evolution of
$\langle\hat{v}^{\alpha}\rangle$ we should compute commutators between central
fluctuation operators $C^{j00}$ and $C^{n10}$. In particular, we have

$\displaystyle\langle\big{[}\hat{H}_{\rm
grav}\;,\hat{v}^{\alpha}\big{]}\rangle\,=\,-\frac{3\pi
G}{2\alpha_{o}}\sum_{n=0}^{\infty}\sum_{j=1}^{\infty}\sum_{m=0}^{n}\left({\begin{array}[]{c}\alpha\\\
j\\\
\end{array}}\right)\langle\hat{v}\rangle^{\alpha+1-n-j}\tilde{\beta}_{nm}\Big{\langle}\langle\hat{r}\rangle\big{[}C^{n00},C^{j00}\big{]}+\big{[}C^{n10},C^{j00}\big{]}\Big{\rangle}.\quad\quad\quad$
(192)

The first bracket on the rhs of the equation above is zero, so our task will
be computing only the second bracket. We get

$\displaystyle\sum_{m=0}^{n}\tilde{\beta}_{nm}\big{\langle}\big{[}C^{n10},\,C^{j00}\big{]}\big{\rangle}$
$\displaystyle=$
$\displaystyle\sum_{m=0}^{n}\tilde{\beta}_{nm}\big{\langle}\big{[}(\delta\hat{v})^{m}\delta\hat{r}(\delta\hat{v})^{n-m},\,(\delta\hat{v})^{j}\big{]}\big{\rangle}$
$\displaystyle=$
$\displaystyle\frac{1}{2}\sum_{m=0}^{n}\sum_{a=0}^{m}\sum_{b=0}^{n-m}\sum_{c=1}^{j}\tilde{\beta}_{nm}\left({\begin{array}[]{c}m\\\
a\\\ \end{array}}\right)\left({\begin{array}[]{c}n-m\\\ b\\\
\end{array}}\right)\left({\begin{array}[]{c}j\\\ c\\\ \end{array}}\right)$
$\displaystyle\quad\quad\times(-1)^{n+j-a-b-c}\langle\hat{v}\rangle^{n+j-a-
b-c}\big{\langle}\big{[}\hat{v}^{a}\hat{r}\hat{v}^{b},\,\hat{v}^{c}\big{]}\big{\rangle}$
$\displaystyle+\frac{1}{2}\sum_{m=0}^{n}\sum_{a=0}^{n-m}\sum_{b=0}^{m}\sum_{c=1}^{j}\tilde{\beta}_{nm}\left({\begin{array}[]{c}n-m\\\
a\\\ \end{array}}\right)\left({\begin{array}[]{c}m\\\ b\\\
\end{array}}\right)\left({\begin{array}[]{c}j\\\ c\\\ \end{array}}\right)$
$\displaystyle\quad\quad\times(-1)^{n+j-a-b-c}\langle\hat{v}\rangle^{n+j-a-
b-c}\big{\langle}\big{[}\hat{v}^{a}\hat{r}\hat{v}^{b},\,\hat{v}^{c}\big{]}\big{\rangle}.$

In deriving the equation above, we have replaced the expectation value of an
(nonsymmetric) operator

$\displaystyle\hat{A}_{n,m}$ $\displaystyle:=$
$\displaystyle\sum_{a=0}^{m}\sum_{b=0}^{n-m}\left({\begin{array}[]{c}m\\\ a\\\
\end{array}}\right)\left({\begin{array}[]{c}n-m\\\ b\\\
\end{array}}\right)(-1)^{n-a-b}\langle\hat{v}\rangle^{n-a-b}\,\hat{v}^{a}\hat{r}\hat{v}^{b},\quad\quad$
(211)

by its symmetric counterpart as

$\displaystyle\langle\hat{B}_{n,m}\rangle$ $\displaystyle=$
$\displaystyle\tfrac{1}{2}\langle\hat{A}_{n,m}+\hat{A}_{n,n-m}\rangle\,.$
(212)

Using linearity and the Leibniz rule, we get

$\displaystyle\sum_{m=0}^{n}\tilde{\beta}_{nm}\big{\langle}\big{[}C^{n10},\,C^{j00}\big{]}\big{\rangle}$
$\displaystyle=\frac{1}{2}\sum_{m=0}^{n}\sum_{a=0}^{m}\sum_{b=0}^{n-m}\sum_{c=1}^{j}\tilde{\beta}_{nm}\left({\begin{array}[]{c}m\\\
a\\\ \end{array}}\right)\left({\begin{array}[]{c}n-m\\\ b\\\
\end{array}}\right)\left({\begin{array}[]{c}j\\\ c\\\ \end{array}}\right)$
(219) $\displaystyle\quad\quad\times(-1)^{n+j-a-
b-c}\langle\hat{v}\rangle^{n+j-a-
b-c}\big{\langle}\hat{v}^{a}\big{[}\hat{r},\,\hat{v}^{c}\big{]}\hat{v}^{b}\big{\rangle}$
$\displaystyle\quad+\frac{1}{2}\sum_{m=0}^{n}\sum_{a=0}^{n-m}\sum_{b=0}^{m}\sum_{c=1}^{j}\tilde{\beta}_{nm}\left({\begin{array}[]{c}n-m\\\
a\\\ \end{array}}\right)\left({\begin{array}[]{c}m\\\ b\\\
\end{array}}\right)\left({\begin{array}[]{c}j\\\ c\\\ \end{array}}\right)$
(226) $\displaystyle\quad\quad\times(-1)^{n+j-a-
b-c}\langle\hat{v}\rangle^{n+j-a-
b-c}\big{\langle}\hat{v}^{a}\big{[}\hat{r},\,\hat{v}^{c}\big{]}\hat{v}^{b}\big{\rangle}$
$\displaystyle=\sum_{m=0}^{n}\sum_{c=1}^{j}\tilde{\beta}_{nm}\left({\begin{array}[]{c}j\\\
c\\\
\end{array}}\right)(-1)^{j-c}\langle\hat{v}\rangle^{j-c}\Big{(}\big{\langle}(\delta\hat{v})^{m}\big{[}\hat{r},\,\hat{v}^{c}\big{]}(\delta\hat{v})^{n-m}\big{\rangle}\Big{)}\quad\quad$
(229)
$\displaystyle=\sum_{m=0}^{n}\sum_{c=1}^{j}\tilde{\beta}_{nm}\left({\begin{array}[]{c}j\\\
c\\\
\end{array}}\right)(-1)^{j-c}\langle\hat{v}\rangle^{j-c}\Big{(}\sum_{l=0}^{c-1}\left({\begin{array}[]{c}c-1\\\
l\\\
\end{array}}\right)\big{\langle}(\delta\hat{v})^{m}\hat{v}^{l}\big{[}\hat{r},\,\hat{v}\big{]}\hat{v}^{c-1-l}(\delta\hat{v})^{n-m}\big{\rangle}\Big{)}$
(234)
$\displaystyle=2i\hbar\sum_{m=0}^{n}\sum_{c=1}^{j}\tilde{\beta}_{nm}\left({\begin{array}[]{c}j\\\
c\\\
\end{array}}\right)(-1)^{j-c}\langle\hat{v}\rangle^{j-1}\Big{(}\langle\hat{h}\rangle\;\sum_{l=0}^{c-1}\left({\begin{array}[]{c}c-1\\\
l\\\ \end{array}}\right)\sum_{k=0}^{c-1}\left({\begin{array}[]{c}c-1\\\ k\\\
\end{array}}\right)\langle\hat{v}\rangle^{-k}G^{(n+k)00}$ (241)
$\displaystyle\quad+\sum_{l=0}^{c-1}\sum_{k=0}^{l}\sum_{d=0}^{c-1-l}\left({\begin{array}[]{c}c-1\\\
l\\\ \end{array}}\right)\left({\begin{array}[]{c}l\\\ k\\\
\end{array}}\right)\left({\begin{array}[]{c}c-1-l\\\ d\\\
\end{array}}\right)\langle\hat{v}\rangle^{-k-d}\big{\langle}(\delta\hat{v})^{m+k}\delta\hat{h}(\delta\hat{v})^{n-m+d}\big{\rangle}\Big{)}.\quad$
(248)

Now we can rewrite Eq. (147) as

$\displaystyle\frac{\big{\langle}\big{[}\hat{\mathcal{H}}_{\rm
grav},\hat{v}^{\alpha}\big{]}\big{\rangle}}{-i\hbar}$ $\displaystyle=:$
$\displaystyle\frac{3\pi G}{\alpha_{o}}\left({\cal
A}(\alpha)\langle\hat{h}\rangle+{\cal B}(\alpha)\right),\quad$ (249)

where ${\cal A}(\alpha)$ and ${\cal B}(\alpha)$ are given, respectively, by

$\displaystyle{\cal A}(\alpha)$ $\displaystyle=$
$\displaystyle\sum_{n=0}^{\infty}\sum_{j=1}^{\infty}\sum_{c=1}^{j}\sum_{k=0}^{c-1}\sum_{l=0}^{c-1}\left({\begin{array}[]{c}c-1\\\
l\\\ \end{array}}\right)\left({\begin{array}[]{c}\alpha\\\ j\\\
\end{array}}\right)\left({\begin{array}[]{c}j\\\ c\\\
\end{array}}\right)\left({\begin{array}[]{c}c-1\\\ k\\\
\end{array}}\right)\beta_{n}(-1)^{j-c}\langle\hat{v}\rangle^{\alpha-
n-k}G^{(n+k)00},\quad\quad$ (258) $\displaystyle{\cal B}(\alpha)$
$\displaystyle=$
$\displaystyle\sum_{n=0}^{\infty}\sum_{j=1}^{\infty}\sum_{m=0}^{n}\sum_{c=1}^{j}\sum_{l=0}^{c-1}\sum_{k=0}^{l}\sum_{d=0}^{c-1-l}\tilde{\beta}_{nm}\left({\begin{array}[]{c}\alpha\\\
j\\\ \end{array}}\right)\left({\begin{array}[]{c}j\\\ c\\\
\end{array}}\right)\left({\begin{array}[]{c}c-1\\\ l\\\
\end{array}}\right)\left({\begin{array}[]{c}l\\\ k\\\
\end{array}}\right)\left({\begin{array}[]{c}c-1-l\\\ d\\\ \end{array}}\right)$
(270) $\displaystyle\times(-1)^{j-c}\langle\hat{v}\rangle^{\alpha-n-
k-d}\big{\langle}(\delta\hat{v})^{m+k}\delta\hat{h}(\delta\hat{v})^{n-m+d}\big{\rangle}.$

By dividing Eq. (249) to $\langle\hat{v}^{\alpha}\rangle$ [given in Eq. (160)]
we obtain

$\displaystyle\frac{1}{3}\frac{\partial_{\tau}\langle\hat{v}^{\alpha}\rangle}{\langle\hat{v}^{\alpha}\rangle}$
$\displaystyle=$ $\displaystyle\frac{\pi G}{\alpha_{o}}\Big{(}{\cal
F}(\alpha)\langle\hat{h}\rangle+{\cal G}(\alpha)\Big{)}\,,\quad$ (271)

where

$\displaystyle{\cal F}(\alpha)$ $\displaystyle:=$
$\displaystyle\left(\sum_{k=0}^{\infty}\left({\begin{array}[]{c}\alpha\\\ k\\\
\end{array}}\right)\langle\hat{v}\rangle^{\alpha-k}G^{k00}\right)^{-1}{\cal
A}(\alpha)\,,$ (274) $\displaystyle{\cal G}(\alpha)$ $\displaystyle:=$
$\displaystyle\left(\sum_{k=0}^{\infty}\left({\begin{array}[]{c}\alpha\\\ k\\\
\end{array}}\right)\langle\hat{v}\rangle^{\alpha-k}G^{k00}\right)^{-1}{\cal
B}(\alpha)\,.\quad\quad\quad$ (277)

Inserting the definition (78) into Eq. (271) we obtain

$\displaystyle\frac{1}{3}\frac{\partial_{\tau}\langle\hat{v}^{\alpha}\rangle}{\langle\hat{v}^{\alpha}\rangle}$
$\displaystyle=$ $\displaystyle{\cal F}(\alpha)\langle\hat{H}\rangle+\frac{\pi
G}{\alpha_{o}}{\cal G}(\alpha)\,.\quad$ (278)

By computing (278) for values $\alpha=-1$ and $\alpha=1/3$ and replacing them
in Eq. (76) we obtain the dressed Friedmann equation $\tilde{H}$. Let us
define

$\displaystyle{\cal F}_{1}$ $\displaystyle:=$ $\displaystyle{\cal
F}({\alpha=-1})\,,\quad\quad{\cal G}_{1}\ :=\ {\cal
G}({\alpha=-1})\,,\quad\quad$ $\displaystyle{\cal F}_{2}$ $\displaystyle:=$
$\displaystyle{\cal F}({\alpha=1/3})\,,\quad\quad{\cal G}_{2}\ :=\ {\cal
G}({\alpha=1/3})\,,\quad\quad$ (279)

and

$\displaystyle\delta{\cal F}$ $\displaystyle=$ $\displaystyle{\cal
F}_{2}-{\cal F}_{1}\,,\quad\quad\delta{\cal G}\ =\ {\cal G}_{2}-{\cal
G}_{1}\,.$ (280)

Using these definitions we derive the dressed Friedmann equation (76) as

$\displaystyle\tilde{H}$ $\displaystyle=$
$\displaystyle\frac{3}{4}\left(\delta{\cal
F}\cdot\langle\hat{H}\rangle+\frac{\pi G}{\alpha_{o}}\delta{\cal G}\right).$
(281)

To the leading order, we have

$\displaystyle{\cal F}_{1}$ $\displaystyle\approx$
$\displaystyle-1+\mathcal{O}(G^{100}/\langle\hat{v}\rangle),\quad{\cal
G}_{1}\approx 0+\mathcal{O}(G^{101}/\langle\hat{v}\rangle)$
$\displaystyle{\cal F}_{2}$ $\displaystyle\approx$
$\displaystyle\frac{1}{3}+\mathcal{O}(G^{100}/\langle\hat{v}\rangle),\quad{\cal
G}_{2}\approx 0+\mathcal{O}(G^{101}/\langle\hat{v}\rangle).$

In Eq. (281), the quantum fluctuations included in the function ${\cal G}$ are
given by moments of the order of $G^{(n+k+d)01}$ which are very small for
large volumes. Therefore, the second term is negligible, and the squared
dressed Hubble rate can be approximated as
$\tilde{H}^{2}\approx(9/16)(\delta{\cal F})^{2}\langle\hat{H}\rangle^{2}$.
Considering only the leading-order correction terms, where
$\delta\mathcal{F}\approx 4/3$, the dressed Hubble rate reduces to
$\tilde{H}^{2}\approx\langle\hat{H}\rangle^{2}$. In this approximation, the
dressed Hubble rate has a similar form as the one provided by the effective
dynamics of LQC [21]; however, in the present case, the backreaction should be
included. This leading-order term of the modified Friedmann equation for the
dressed metric $(\tilde{N},\tilde{a})$ will be sufficient for our purpose in
this paper.

##  Appendix B Sub-leading terms in eigenfunctions
${\underline{e}}^{\pm}_{\mu,\mathbf{k}}(v)$

The rate of convergence for eigenfunctions is a few orders weaker than what is
needed for numerical calculations. We found the following corrections
improving the convergence rate of functions (53). $b_{n}$ and $a_{n}$ are
corrections to the phase and amplitude of the eigenfunctions, respectively:

$\displaystyle b_{1}$ $\displaystyle=\frac{\alpha_{o}l}{\pi
G\mu^{2}},\quad\quad\quad
b_{2}=\frac{\alpha_{o}^{2}l^{2}}{6\pi^{2}G^{2}\mu^{4}},\quad\quad\quad
b_{3}=\frac{-36\pi^{3}G^{3}\mu^{8}+81\pi^{3}G^{3}\mu^{4}-4\alpha_{o}^{3}l^{3}}{216\pi^{3}G^{3}\mu^{6}},$
$\displaystyle b_{4}$
$\displaystyle=\frac{-108\pi^{3}G^{3}\mu^{8}\alpha_{o}l-21\pi^{3}G^{3}\mu^{4}\alpha_{o}l+5\alpha_{o}^{4}l^{4}}{1080\pi^{4}G^{4}\mu^{8}},\quad\quad
b_{5}=\frac{-108\pi^{3}G^{3}\mu^{8}\alpha_{o}^{2}l^{2}-237\pi^{3}G^{3}\mu^{4}\alpha_{o}^{2}l^{2}-14\alpha_{o}^{5}l^{5}}{9072\pi^{5}G^{5}\mu^{10}},$
$\displaystyle b_{6}$
$\displaystyle=\frac{-11664\pi^{6}G^{6}\mu^{16}+29160\pi^{6}G^{6}\mu^{12}-76545\pi^{6}G^{6}\mu^{8}+480\pi^{3}G^{3}\mu^{8}\alpha_{o}^{3}l^{3}+14600\pi^{3}G^{3}\mu^{4}\alpha_{o}^{3}l^{3}+280\alpha_{o}^{6}l^{6}}{466560\pi^{6}G^{6}\mu^{12}},$
$\displaystyle b_{7}$
$\displaystyle=\frac{-11664\pi^{6}G^{6}\mu^{16}\alpha_{o}l+12456\pi^{6}G^{6}\mu^{12}\alpha_{o}l+53967\pi^{6}G^{6}\mu^{8}\alpha_{o}l-72\pi^{3}G^{3}\mu^{8}\alpha_{o}^{4}l^{4}-9310\pi^{3}G^{3}\mu^{4}\alpha_{o}^{4}l^{4}-88\alpha_{o}^{7}l^{7}}{342144\pi^{7}G^{7}\mu^{14}},$
(282)

and for $a_{n}$ we find

$\displaystyle a_{1}$ $\displaystyle=-\frac{\alpha_{o}l}{6\pi
G\mu^{2}},\quad\quad\quad
a_{2}=\frac{5\alpha_{o}^{2}l^{2}}{72\pi^{2}G^{2}\mu^{4}},\quad\quad\quad
a_{3}=\frac{36\pi^{3}G^{3}\mu^{8}+27\pi^{3}G^{3}\mu^{4}-5\alpha_{o}^{3}l^{3}}{144\pi^{3}G^{3}\mu^{6}},$
$\displaystyle a_{4}$
$\displaystyle=\frac{432\pi^{3}G^{3}\mu^{8}\alpha_{o}l-380\pi^{3}G^{3}\mu^{4}\alpha_{o}l+65\alpha_{o}^{4}l^{4}}{3456\pi^{4}G^{4}\mu^{8}},$
$\displaystyle a_{5}$
$\displaystyle=\frac{-216\pi^{3}G^{3}\mu^{8}\alpha_{o}^{2}l^{2}+670\pi^{3}G^{3}\mu^{4}\alpha_{o}^{2}l^{2}-221\alpha_{o}^{5}l^{5}}{20736\pi^{5}G^{5}\mu^{10}},$
(283)

where $l=\alpha_{o}^{-1/3}N_{k}\hbar k$.

## References

  * Gambini and Pullin [1999] R. Gambini and J. Pullin, Nonstandard optics from quantum space-time, Phys. Rev. D59, 124021 (1999), arXiv:gr-qc/9809038 [gr-qc] .
  * Alfaro  _et al._ [2000] J. Alfaro, H. A. Morales-Tecotl, and L. F. Urrutia, Quantum gravity corrections to neutrino propagation, Phys. Rev. Lett. 84, 2318 (2000), arXiv:gr-qc/9909079 [gr-qc] .
  * Magueijo and Smolin [2004] J. Magueijo and L. Smolin, Gravity’s rainbow, Class. Quant. Grav. 21, 1725 (2004), arXiv:gr-qc/0305055 [gr-qc] .
  * Mattingly [2005] D. Mattingly, Modern tests of Lorentz invariance, Living Rev. Rel. 8, 5 (2005), arXiv:gr-qc/0502097 [gr-qc] .
  * Lafrance and Myers [1995] R. Lafrance and R. C. Myers, Gravity’s rainbow, Phys. Rev. D51, 2584 (1995), arXiv:hep-th/9411018 [hep-th] .
  * Amelino-Camelia  _et al._ [1998] G. Amelino-Camelia, J. R. Ellis, N. E. Mavromatos, D. V. Nanopoulos, and S. Sarkar, Tests of quantum gravity from observations of gamma-ray bursts, Nature 393, 763 (1998), arXiv:astro-ph/9712103 [astro-ph] .
  * Addazi  _et al._ [2021] A. Addazi _et al._ , Quantum gravity phenomenology at the dawn of the multi-messenger era – A review, (2021), arXiv:2111.05659 [hep-ph] .
  * Ashtekar  _et al._ [2009] A. Ashtekar, W. Kaminski, and J. Lewandowski, Quantum field theory on a cosmological, quantum space-time, Phys. Rev. D79, 064030 (2009), arXiv:0901.0933 [gr-qc] .
  * Dapor  _et al._ [2013] A. Dapor, J. Lewandowski, and J. Puchta, QFT on quantum spacetime: a compatible classical framework, Phys. Rev. D87, 104038 (2013), [Erratum: Phys. Rev.D87,no.12,129904(2013)], arXiv:1302.3038 [gr-qc] .
  * Lewandowski  _et al._ [2017] J. Lewandowski, M. Nouri-Zonoz, A. Parvizi, and Y. Tavakoli, Quantum theory of electromagnetic fields in a cosmological quantum spacetime, Phys. Rev. D96, 106007 (2017), arXiv:1709.04730 [gr-qc] .
  * Dapor  _et al._ [2012] A. Dapor, J. Lewandowski, and Y. Tavakoli, Lorentz Symmetry in QFT on Quantum Bianchi I Space-Time, Phys. Rev. D86, 064013 (2012), arXiv:1207.0671 [gr-qc] .
  * Assanioussi  _et al._ [2015] M. Assanioussi, A. Dapor, and J. Lewandowski, Rainbow metric from quantum gravity, Phys. Lett. B751, 302 (2015), arXiv:1412.6000 [gr-qc] .
  * Garcia-Chung  _et al._ [2021] A. Garcia-Chung, J. B. Mertens, S. Rastgoo, Y. Tavakoli, and P. Vargas Moniz, Propagation of quantum gravity-modified gravitational waves on a classical FLRW spacetime, Phys. Rev. D 103, 084053 (2021), arXiv:2012.09366 [gr-qc] .
  * Bernardini  _et al._ [2015] M. G. Bernardini _et al._ , Comparing the spectral lag of short and long gamma-ray bursts and its relation with the luminosity, Mon. Not. Roy. Astron. Soc. 446, 1129 (2015), arXiv:1410.5216 [astro-ph.HE] .
  * Yi  _et al._ [2006] T.-F. Yi, E.-W. Liang, Y.-P. Qin, and R.-J. Lu, On the spectral lags of the short gamma-ray bursts, Mon. Not. Roy. Astron. Soc. 367, 1751 (2006), arXiv:astro-ph/0512270 .
  * Meszaros [2006] P. Meszaros, Gamma-Ray Bursts, Rept. Prog. Phys. 69, 2259 (2006), arXiv:astro-ph/0605208 .
  * Tavakoli  _et al._ [2013] Y. Tavakoli, J. Marto, A. H. Ziaie, and P. Vargas Moniz, Semiclassical collapse with tachyon field and barotropic fluid, Phys. Rev. D87, 024042 (2013).
  * Tavakoli  _et al._ [2014] Y. Tavakoli, J. Marto, and A. Dapor, Semiclassical dynamics of horizons in spherically symmetric collapse, Int. J. Mod. Phys. D23, 1450061 (2014), arXiv:1303.6157 [gr-qc] .
  * Goswami  _et al._ [2006] R. Goswami, P. S. Joshi, and P. Singh, Quantum evaporation of a naked singularity, Phys. Rev. Lett. 96, 031302 (2006), arXiv:gr-qc/0506129 [gr-qc] .
  * Bojowald  _et al._ [2005] M. Bojowald, R. Goswami, R. Maartens, and P. Singh, A Black hole mass threshold from non-singular quantum gravitational collapse, Phys. Rev. Lett. 95, 091302 (2005), arXiv:gr-qc/0503041 [gr-qc] .
  * Husain and Pawlowski [2012] V. Husain and T. Pawlowski, Time and a physical Hamiltonian for quantum gravity, Phys. Rev. Lett. 108, 141301 (2012), arXiv:1108.1145 [gr-qc] .
  * Brown and Kuchar [1995] J. D. Brown and K. V. Kuchar, Dust as a standard of space and time in canonical quantum gravity, Phys. Rev. D51, 5600 (1995), arXiv:gr-qc/9409001 [gr-qc] .
  * Bojowald [2010] M. Bojowald,  _Canonical Gravity and Applications: Cosmology, Black Holes, and Quantum Gravity_ (Cambridge University Press, 2010).
  * Ashtekar  _et al._ [2006a] A. Ashtekar, T. Pawlowski, and P. Singh, Quantum nature of the big bang, Phys. Rev. Lett. 96, 141301 (2006a), arXiv:gr-qc/0602086 [gr-qc] .
  * Ashtekar  _et al._ [2003] A. Ashtekar, M. Bojowald, and J. Lewandowski, Mathematical structure of loop quantum cosmology, Adv. Theor. Math. Phys. 7, 233 (2003), arXiv:gr-qc/0304074 [gr-qc] .
  * Husain and Pawlowski [2011] V. Husain and T. Pawlowski, Dust reference frame in quantum cosmology, Class. Quant. Grav. 28, 225014 (2011), arXiv:1108.1147 [gr-qc] .
  * Giesel  _et al._ [2009] K. Giesel, J. Tambornino, and T. Thiemann, Born-Oppenheimer decomposition for quantum fields on quantum spacetimes, (2009), arXiv:0911.5331 [gr-qc] .
  * Bojowald and Skirzewski [2006] M. Bojowald and A. Skirzewski, Effective equations of motion for quantum systems, Rev. Math. Phys. 18, 713 (2006), arXiv:math-ph/0511043 .
  * Berger [1975] B. K. Berger, Scalar Particle Creation in an Anisotropic Universe, Phys. Rev. D12, 368 (1975).
  * Ashtekar  _et al._ [2006b] A. Ashtekar, T. Pawlowski, and P. Singh, Quantum Nature of the Big Bang: Improved dynamics, Phys. Rev. D 74, 084003 (2006b), arXiv:gr-qc/0607039 .
  * Mena Marugan  _et al._ [2011] G. A. Mena Marugan, J. Olmedo, and T. Pawlowski, Prescriptions in Loop Quantum Cosmology: A comparative analysis, Phys. Rev. D 84, 064012 (2011), arXiv:1108.0829 [gr-qc] .
  * Parvizi  _et al._ [shed] A. Parvizi, T. Pawlowski, and Y. Tavakoli, (To be published).
  * Vaidya [1951] P. C. Vaidya, Nonstatic Solutions of Einstein’s Field Equations for Spheres of Fluids Radiating Energy, Phys. Rev. 83, 10 (1951).
  * Joshi [1993] P. S. Joshi,  _Global aspects in gravitation and cosmology_ (Clarendon Press, Oxford, 1993).
  * Vaidya [1999a] P. C. Vaidya, The External Field of a Radiating Star in Relativity, Gen. Rel. Grav. 31, 119 (1999a).
  * Vaidya [1999b] P. C. Vaidya, The Gravitational Field of a Radiating Star, Gen. Rel. Grav. 31, 121 (1999b).
  * Vaidya [1953] P. C. Vaidya, Newtonian Time in General Relativity, Nature 171, 260 (1953).
  * Hayward [1993] S. A. Hayward, Marginal surfaces and apparent horizons, (1993), arXiv:gr-qc/9303006 .
  * Ashtekar and Krishnan [2002] A. Ashtekar and B. Krishnan, Dynamical horizons: Energy, angular momentum, fluxes and balance laws, Phys. Rev. Lett. 89, 261101 (2002), arXiv:gr-qc/0207080 .
  * Ziaie and Tavakoli [2020] A. H. Ziaie and Y. Tavakoli, Null Fluid Collapse in Rastall Theory of Gravity, Annalen Phys. 532, 2000064 (2020), arXiv:1912.08890 [gr-qc] .
  * Hawking and Ellis [2011] S. W. Hawking and G. F. R. Ellis,  _The Large Scale Structure of Space-Time_, Cambridge Monographs on Mathematical Physics (Cambridge University Press, 2011).
  * Weinberg [1972] S. Weinberg,  _Gravitation and Cosmology: Principles and Applications of the General Theory of Relativity_ (John Wiley and Sons, New York, 1972).
  * Husain [1996] V. Husain, Exact solutions for null fluid collapse, Phys. Rev. D 53, R1759 (1996).
  * Keeton and Petters [2005] C. R. Keeton and A. O. Petters, Formalism for testing theories of gravity using lensing by compact objects. I. Static, spherically symmetric case, Phys. Rev. D 72, 104006 (2005), arXiv:gr-qc/0511019 .
  * Nouri-Zonoz  _et al._ [shed] M. Nouri-Zonoz, A. Parvizi, and Y. Tavakoli, , (To be published).
  * Sereno [2004] M. Sereno, Weak field limit of Reissner-Nordstrom black hole lensing, Phys. Rev. D 69, 023002 (2004), arXiv:gr-qc/0310063 .
  * Piron [2016] F. Piron, Gamma-Ray Bursts at high and very high energies, Comptes Rendus Physique 17, 617 (2016), arXiv:1512.04241 [astro-ph.HE] .
  * Miller-Jones  _et al._ [2021] J. C. A. Miller-Jones _et al._ , Cygnus X-1 contains a 21–solar mass black hole—Implications for massive star winds, Science 371, 1046 (2021), arXiv:2102.09091 [astro-ph.HE] .