Quantum Fluctuations in the Interior of Black Holes and Backreactions

Hassan<EMAIL_ADDRESS>
School of Astronomy, Institute for Research in Fundamental Sciences (IPM)
P. O. Box 19395-5531, Tehran, Iran

###### Abstract

We study the propagation of the quantum field perturbations in the interior of
the Schwarzschild black hole. The interior of the black hole is like an
anisotropic cosmological background which expands in one extended direction
while contracting in azimuthal directions. Solving the quantum mode functions
approximately, we calculate the expectation values of physical quantities such
as the energy density $\langle\rho\rangle$ and pressure $\langle P\rangle$ as
measured by an observer in the interior of the black hole. We employ a
combination of the $\zeta$ function regularization and dimensional
regularization schemes to regularize the UV divergences in the energy momentum
tensor associated to these quantum perturbations. By solving the Einstein
field equations, we calculate the quantum backreactions induced in the
background geometry. We speculate that the effects of quantum fluctuations in
the interior of the black hole can be measured by the exterior observer.

##  1 Introduction

The event horizon surrounding a black hole (BH) separates it into two causally
disconnected regions. To a remote observer outside the BH, it represents a
point-like singularity which may be observed via their gravitational effects
such as in merging and generation of gravitational wave as observed in
LIGO/Virgo/Advanced Virgo collaborations [1, 2]. Most of the studies
concerning BHs are naturally devoted to the exterior region as viewed by the
exterior observers. This is mainly motivated from the simple fact that the
events occurring inside the BH are inaccessible to the exterior observer.
Having said this, the interior of the BH reveals interesting physical
phenomena which deserve more attention. Specifically, the interior of BH looks
like a cosmological background in which the singularity is spacelike in the
form of big crunch singularity. In simple setups such as in Schwarzschild
solution, the interior of BH looks like a homogenous but anisotropic
cosmological background. The spatial part has the topology of $R\times S^{2}$
known as the Kantowski-Sachs spacetime [3]. This anisotropic cosmological
background contains two scale factors, and correspondingly, two “Hubble”
expansion rates. The scale factor along the extended spatial direction expands
while the scale factor associated to the azimuthal directions contracts.
Eventually, near the singularity, the scale factor of the azimuthal direction
reaches zero, signalling the onset of the big crunch singularity.

In the past there have been works to treat the interior of a BH as a
cosmological background. One interesting idea is that the interior of BH is
replaced by a non-singular dS spacetime, as studied in [4, 5, 6, 7, 8, 9, 10,
11, 12, 13, 14, 15] or the proposal of gravastars [16, 17, 18]. The main
motivation in these studies was to remove the singularity of BH. On the
physical grounds one expects that the singularity of BH is a shortcoming of
the classical general relativity. It is expected that on very small scales,
say the Planck scale, the quantum gravity effects play important roles to
resolve the classical singularity. In particular, the proposal of maximum
curvature of spacetime [19, 6, 7] is an interesting idea in this direction.

The question of quantum field perturbations in the exterior of the BH is
vastly studied. Among the non-trivial issues in these studies is the fact that
there is no unique vacuum in a curved background, specially in a BH background
which has an event horizon. This has lead to the celebrated Hawking radiation
effect [20, 21], for some reviews of Hawking radiation see [22, 23, 24, 25,
26]. While QFTs effects are vastly studied in the exterior region, they are
less investigated for the interior region. In this work we study the
propagation of quantum field perturbations in the interior of BH. Intuitively
speaking, the perturbations enter the horizon after crossing the event horizon
and propagate towards the singularity. Since the interior of the BH is like an
anisotropic cosmological background, the propagation of the quantum
perturbations in the interior region is similar to the evolution of
perturbations in FLRW cosmology. While this similarity is very helpful for the
understanding of the nature of quantum perturbations in the interior of BH,
but there are a number of important differences as well. First, compared to
isotropic FLRW cosmology, we have less symmetry, the former enjoys an $O(3)$
symmetry while the latter has the smaller $O(2)$ symmetry. The other important
difference is that in FLRW cosmology the quantum perturbations are generated
near the time of big bang, as in inflationary cosmology. Here, however, the
perturbations are viewed to be generated near the event horizon and terminated
at the moment of big crunch.

In this study we consider a massless scalar field in the interior of a
Schwarzschild background. After reviewing the background cosmology associated
to the interior of the BH (section 2), we consider the quantum perturbations.
Upon decomposing the perturbations into appropriate bases of the creation and
annihilation operators, we solve the mode function approximately (section 3).
Our main interest is to calculate the expectation values of physical
quantities such as the energy density $\langle\rho\rangle$ and pressure
$\langle P\rangle$. As in many QFT analysis, we inevitably encounter UV
divergences which should be regularized appropriately (section 4). Our goal is
to calculate the backreactions of the regularized energy momentum tensor
associated to these quantum perturbations on the background geometry. Upon
solving the Einstein field equations, we obtain the corrections in the metric
to leading order in perturbations and look for their physical implications
(section 5).

##  2 Background Geometry

In this section we briefly review basic properties of the Schwarzschild
solution which will be employed for our treatment of the quantum fields in
this background.

###  2.1 Black Hole Preliminaries

We consider the Schwarzschild background with the following metric,

$\displaystyle\differential
s^{2}=-\big{(}1-\frac{2GM}{r}\big{)}\differential{\bar{t}\,}^{2}+\frac{\differential
r^{2}}{\big{(}1-\frac{2GM}{r}\big{)}}+r^{2}\differential\Omega^{2}\,,$ (2.1)

in which $G$ is the Newton constant, $M$ is the mass of the BH as measured by
a distant observer and $\differential\Omega^{2}$ represents the metric on a
unit two-sphere. The coordinate system $(r,\bar{t})$ covers only the exterior
of the whole manifold which becomes singular on the BH event horizon at
$r=R_{S}\equiv 2GM$. Note that we use the notation $\bar{t}$ for the time
coordinate employed by the exterior observer while we reserve $t$ for a
timelike coordinate employed for the interior of the black hole. Naively
speaking, for the interior of the BH the roles of $\bar{t}$ and $r$
coordinates are switched in which $\bar{t}$ becomes spacelike while $r$
becomes timelike. This also suggests that the interior of BH is dynamical, as
in cosmological backgrounds.

Since the coordinate $(r,\bar{t})$ is singular at the horizon, to cover the
entire manifold we can use the Kruskal-Szekeres coordinate $(T,R)$ which are
related to the original coordinate $(r,\bar{t})$ by

$\displaystyle T^{2}-R^{2}=e^{r/2GM}\big{(}1-\frac{r}{2GM}\big{)}\,,$ (2.2)

and

$\displaystyle\frac{T}{R}=\tanh(\frac{\bar{t}}{4GM})\,.$ (2.3)

Furthermore, it is very convenient to use the Kruskal-Szekeres coordinate in
its lightcone base $(U,V)$ defined via

$\displaystyle U\equiv GM(T-R),\quad\quad V\equiv GM(T+R)\,.$ (2.4)

We use the convention such that the coordinates $(U,V)$ carry the dimension of
length (or time).

A conformal diagram of the entire Schwarzschild manifold is depicted in Fig.
1. The exterior of the BH is the region $V>0,U<0$ while the interior of the
BH, which we are interested in here, is the region $U,V>0$ bounded by the
future singularity $r=0$. The white hole (WH) region is given by $U,V<0$
bounded by the past singularity $r=0$. Both the interiors of the BH and the WH
share the common property that their backgrounds are dynamical, corresponding
to anisotropic cosmological setups. However, the crucial difference is that
the singularity of the BH is in future, i.e. it represents a big crunch
singularity while the singularity of the WH is in the past, i.e. it represents
a big bang singularity. Mathematically speaking, the WH solution is realized
from the time reversal symmetry of the BH solution in classical general
relativity.

The generation and propagation of the quantum perturbations in WH region were
studied in [27] as a realization of Hawking radiation measured by an observer
far outside the BH. More specifically, it was argued in [27] that the quantum
fluctuations which are generated at past singularity $r=0$ (big bang) will
propagate towards the past horizon $V=0,U<0$ and eventually towards ${\cal
I}^{+}$ which is measured by an observer at future infinity in the exterior
part of the BH. It was shown that while the observer deep inside the WH and
the observer far outside the BH both share the same vacuum, but their vacuum
is different than the one defined for a freely falling observer near the
horizon. In particular, the WH observer and the observer far outside the BH
detect a spectrum of particle with the Planck distribution having the Hawking
temperature $T_{H}=\frac{1}{8\pi GM}$. Intuitively speaking, the conclusion of
[27] implies that if a BH is not entirely black due to Hawking radiation, then
the WH is not entirely white either!

In this works, however, we consider the quantum fluctuations which enters
inside the BH from the future event horizon $U=0,V>0$ and propagate towards
the future singularity $r=0$ (big crunch) as shown schematically by the
arrowed lines in Fig. 1.

Figure 1:  The Kruskal-Szekeres diagram. Region I (II) denotes the exterior
(interior) of BH while the region III denotes the WH part. The future event
horizon is located at $U=0,V>0$. The arrowed lines represent the propagation
of the quantum perturbations towards the interior region, heading to future
singularity (big crunch) at $r=0$.

###  2.2 Cosmology of the Interior of the Black Hole

As mentioned above, the interior of the BH represents an anisotropic
cosmological setup known as the Kantowski-Sachs background [3, 28] with the
spatial part having the structure of $R\times S^{2}$. To see this more
clearly, let us start with the usual tortoise coordinate

$\displaystyle\differential r_{*}=(1-\frac{2GM}{r})^{-1}\differential r\,,$
(2.5)

which for the interior of BH ($r<2GM$) yields the following relation,

$\displaystyle
r_{*}=r+2GM\ln\left(1-\frac{r}{2GM}\right),\qquad\quad({r<2GM})\,.$ (2.6)

We note that $-\infty<r_{*}\leq 0$ in which $r_{*}=-\infty$ corresponds to the
future event horizon $r=2GM,U=0$ with $V>0$ while $r_{*}=0$ corresponds to the
future singularity at $r=0$.

Let us define the future directed time coordinate
$\differential\tau=\differential r_{*}$ so

$\displaystyle\tau=r+2GM\ln\left(1-\frac{r}{2GM}\right)\qquad\quad({r<2GM})\,.$
(2.7)

Correspondingly, the BH singularity is at $\tau=0$ while the future event
horizon $V>0,U=0$ is mapped to $\tau=-\infty$. Proceeding further, let us
define the spacelike coordinate $\differential x\equiv\differential\bar{t}$
with $-\infty<x<+\infty$. In the new coordinate $(\tau,x)$, the metric for the
interior of BH takes the following cosmological form

$\displaystyle\differential
s^{2}=a(\tau)^{2}(-\differential\tau^{2}+\differential
x^{2})+r(\tau)^{2}\differential\Omega^{2}\,,$ (2.8)

in which the scale factor $a(\tau)$ is defined via

$\displaystyle
a(\tau)\equiv\Big{(}\frac{2GM}{r(\tau)}-1\Big{)}^{\frac{1}{2}}\,.$ (2.9)

Note that from Eq. (2.7), one can solve for $r$ as a function of $\tau$.

The metric (2.8) represents an anisotropic cosmological background with two
scale factors $a(\tau)$ and $r(\tau)$. Near the horizon, $a(\tau)\rightarrow
0$ while $r(\tau)\rightarrow R_{S}$. As time goes by, the space along the $x$
direction expands while the scale factor along the two-sphere shrinks. At the
moment of big crunch $\tau=0$, $a(\tau)\rightarrow\infty$ while
$r(\tau)\rightarrow 0$.

It is insightful to look at the metric for two limiting times,
$\tau\rightarrow-\infty$ and $\tau\rightarrow 0$ more closely. Using Eq.
(2.7), at later time when $r\rightarrow 0$ we have

$\displaystyle r^{2}\simeq-4GM\tau\,,\quad\quad
a(\tau)\simeq\big{(}\frac{GM}{-\tau}\big{)}^{\frac{1}{4}}\quad\quad\quad(r\rightarrow
0)\,,$ (2.10)

so the metric takes the following form

$\displaystyle\differential
s^{2}\simeq\sqrt{\frac{GM}{-\tau}}(-\differential\tau^{2}+\differential
x^{2})-4GM\tau~{}\differential\Omega^{2}\quad\quad(\tau\rightarrow 0^{-})\,.$
(2.11)

The contraction along the $x$ direction and the expansion along the azimuthal
directions is evident from the above line element.

On the other hand, near the past horizon $r=2GM,U=0$, we have,

$\displaystyle r(\tau)\simeq
2GM\Big{[}1-\exp\big(\frac{\tau}{2GM}\big{missing})\Big{]}\quad\quad(\tau\rightarrow-\infty)\,,$
(2.12)

so the metric takes the following form,

$\displaystyle\differential
s^{2}\simeq\exp\big({\frac{\tau}{2GM}}\big{missing})(-\differential\tau^{2}+\differential
x^{2})+4G^{2}M^{2}\differential\Omega^{2}\quad\quad(\tau\rightarrow-\infty)\,.$
(2.13)

Note that near the horizon, the two-sphere reaches the fixed radius
$r\rightarrow R_{S}$ while the scale factor in the $x$-direction falls off
exponentially. Although the scale factor $a(\tau)$ falls off, but the
spacetime is regular near the horizon. This is because the original coordinate
$(r,\bar{t})$ was singular near the horizon. More specifically, calculating
the Riemann tensor associated to the metric (2.13) we find that
$R_{\theta\phi\theta\phi}=4G^{2}M^{2}\sin(\theta)^{2}$ while all other
components are zero so the background is regular near the horizon as expected.

We can define the “cosmic time” $\hat{t}$ related to the conformal time via
$\differential\hat{t}\equiv a(\tau)~{}\differential\tau$ so

$\displaystyle\differential\hat{t}=\frac{-\differential
r}{\sqrt{\frac{2GM}{r}-1}}\,.$ (2.14)

Choosing the onset of big crunch singularity to be at $\hat{t}=0$, we obtain

$\displaystyle\hat{t}=\sqrt{r(2GM-r)}+MG\Big{[}-\frac{\pi}{2}+\mathrm{arc}\sin\Big(1-\frac{r}{GM}\Big{missing})\Big{]}\,,$
(2.15)

so $-{\pi}GM<\hat{t}<0$ inside the BH. From the inverse of the above
expression, one can obtain $r=r(\hat{t})$ so the metric in term of cosmic time
is given by

$\displaystyle\differential
s^{2}=-\differential\hat{t}^{2}+\Big{(}\frac{2GM}{r(\hat{t}\,)}-1\Big{)}\differential
x^{2}+r(\hat{t}\,)^{2}\differential\Omega^{2}.$ (2.16)

It is instructive to look at the geometry near the big crunch singularity when
expressed in the form of metric (2.16). Near $r=0$, we have

$\displaystyle
r\simeq(GM)^{\frac{1}{3}}\bigg{(}\frac{-3\hat{t}}{\sqrt{2}}\bigg{)}^{\frac{2}{3}}\quad\quad(\hat{t}\rightarrow
0^{-})\,,$ (2.17)

so the metric (2.16) takes the following form,

$\displaystyle\differential
s^{2}\simeq-\differential\hat{t}^{2}+(GM)^{\frac{2}{3}}\Big{(}c_{1}{\hat{t}}^{-\frac{2}{3}}\differential
x^{2}+c_{2}{\hat{t}}^{\frac{4}{3}}\differential\Omega^{2}\Big{)}\,,$ (2.18)

in which $c_{1}$ and $c_{2}$ are two numerical factors. The above geometry is
in the form of the Kasner metric. More specifically, the Kasner metric is the
solution of the Einstein field equation in the vacuum with the following form,

$\displaystyle\differential
s^{2}=-{\differential\hat{t}\,}^{2}+\sum_{i=1}^{3}{\hat{t}\,}^{2p_{i}}\differential
x_{i}^{2}\,,$ (2.19)

where the exponents $p_{i}$ satisfy the following constraints:

$\displaystyle\sum_{i=1}^{3}p_{i}=\sum_{i=1}^{3}p_{i}^{2}=1\,.$ (2.20)

From the above constraints one concludes that at least one of the exponents
$p_{i}$ should be negative. In our example of BH cosmology near the big crunch
singularity with the metric (2.18), one has $p_{1}=-\frac{1}{3}$ while
$p_{2}=p_{3}=\frac{2}{3}$. It can be checked that these exponents satisfy the
constraints (2.20) with the contraction being along the $x$ direction.

We can also define the “Hubble rate” associated to metric (2.16). Since there
are two scale factors $a(\hat{t})$ and $r(\hat{t})$, correspondingly we should
define two Hubble rates. Defining the Hubble rate associated to the scale
factor $a(\hat{t})$ by $H_{1}$, we obtain

$\displaystyle H_{1}\equiv\frac{1}{a}\frac{\differential
a}{\differential\hat{t}}=\frac{GM}{{r(\hat{t})}^{2}}\Big{(}\frac{2GM}{r(\hat{t}\,)}-1\Big{)}^{-\frac{1}{2}}\,.$
(2.21)

Since the scale factor $a$ is expanding, the Hubble rate $H_{1}$ is positive.
Furthermore, $H_{1}$ diverges both near the BH big crunch singularity and near
the horizon. The former divergence is genuine as the spacetime is singular at
$r=0$ so physical quantities, such as the energy density or Hubble rate,
diverge. On the other hand, the divergence of $H_{1}$ near the event horizon
is an artifact of coordinate singularity as suggested from the relation
between $\differential\hat{t}$ and $\differential r$ in Eq. (2.14).

On the other hand, denoting the Hubble rate associated to the scale factor
$r(\hat{t})$ by $H_{2}$ we have,

$\displaystyle H_{2}\equiv\frac{1}{r(\hat{t})}\frac{\differential
r(\hat{t})}{\differential\hat{t}}=\frac{-1}{r(\hat{t})}\Big{(}\frac{2GM}{r(\hat{t}\,)}-1\Big{)}^{\frac{1}{2}}\,.$
(2.22)

Since the scale factor $r(\hat{t})$ is contracting, $H_{2}<0$. In addition,
$H_{2}$ diverges near the singularity as expected but it vanishes near the
event horizon. This is expected since near the horizon
($\tau\rightarrow-\infty$) the radius of the two-sphere approaches a constant
size so $H_{2}$ vanishes in this limit. From the expressions for $H_{1}$ and
$H_{2}$, one can check that the following relation holds,

$\displaystyle H_{1}H_{2}=-\frac{GM}{r(\hat{t}\,)^{3}}\,,$ (2.23)

so the two Hubble rates have opposite signs as required from the structure of
Kasner geometry (i.e. expansion in one direction and contractions in the
remaining directions).

##  3 Quantum Fields Inside the Black Hole

Here we study the quantum field perturbations in the interior of BHs. A
similar analysis were performed for the quantum fluctuations in the WH
backgrounds in [27], see also [29] for a recent study of the quantum field in
the interior of the BH.

We consider a real massless scalar field $\Phi$ which propagates in the
interior of the BH. In this section, we work in the test field limit where the
backreaction of the field on the background is negligible. However, in the
next section we take into account the effects of the backreactions on the
background geometry. We use the metric (2.8) with the conformal time
$\tau\equiv r_{*}$ while $x\equiv\bar{t}$ represents the extended spatial
direction (as mentioned before, in the interior of the BH, $\bar{t}$ is
spacelike while $r$ is timelike).

To perform the QFT regularizations, we employ the dimensional regularization
scheme [30, 31, 32, 33, 34, 35] in which the extended directions ($R^{2}$
parts of the whole $R^{2}\times S^{2}$ manifold) is generalized to $d$
dimensions so the whole manifold now is $R^{d}\times S^{2}$. The $R^{d}$ part
has one timelike coordinate along $\tau$ as before while the $d-1$ spatial
directions are denoted by $\bf{x}$. To perform the regularization, we take
$d=2-\epsilon$ and take $\epsilon\rightarrow 0$ at the end. The divergent
terms with the poles like $\epsilon^{-1}$ in physical quantities such as the
energy density or pressure are discarded as these UV divergences are taken
care of by appropriate counter terms. The regularized contributions are the
finite $\epsilon^{0}$ terms.

With these discussions in mind, the action of the scalar field in $d+2$
dimension is given by

$\displaystyle S$ $\displaystyle=$
$\displaystyle-\frac{1}{2}\int\differential\tau\differential
x^{d-1}\differential\Omega\,a(\tau)^{d}r(\tau)^{2}\partial^{\mu}\Phi\partial_{\mu}\Phi\,$
(3.1) $\displaystyle=$
$\displaystyle\frac{1}{2}\int\differential\tau\differential
x^{d-1}\differential\Omega\,a(\tau)^{d-2}r(\tau)^{2}\Big{[}\big{(}\partial_{\tau}\Phi\big{)}^{2}-\big{(}\partial_{\bf
x}\Phi\big{)}^{2}-\frac{a(\tau)^{2}}{r(\tau)^{2}}\gamma^{ab}\partial_{a}\Phi\partial_{b}\Phi\Big{]}\,,$

in which $\gamma_{ab}=\mathrm{diag}(1,\sin^{2}\theta)$ is the metric along the
azimuthal directions $\\{a,b\\}=\\{\theta,\phi\\}$. Note that the factor
$a(\tau)^{d}$ in the first line appears because the part of the background
metric (2.8) endowed with the scale factor $a(\tau)$ is $d$-dimensional.

Since the background enjoys the rotation symmetry in the azimuthal directions,
we expand the quantum field $\Phi$ in terms of the spherical harmonics
$Y_{\ell m}(\theta,\phi)$ as follows,

$\displaystyle\Phi=f(\tau)\sum_{\ell=0}^{\infty}\sum_{m=-\ell}^{\ell}Z_{\ell
m}(\tau,{\bf x})Y_{\ell m}(\theta,\phi)\,,\quad\quad
f(\tau)\equiv\frac{a(\tau)^{\frac{2-d}{2}}}{r(\tau)},$ (3.2)

in which $Z_{\ell m}(\tau,{\bf x})$ should be viewed as the quantum mode
operator. We have pulled out a factor $f(\tau)=a^{(2-d)/2}/r$ so $Z_{\ell m}$
will be the canonically normalized field in the perturbation analysis. Note
that when $d=2$, we recover the usual convention that $f=1/r$ as employed in
BH quasi-normal mode analysis. Finally, the reality of $\Phi$ requires that
$Z_{\ell m}^{*}=(-1)^{m}Z_{\ell\,-m}$.

Plugging the decomposition (3.2) for the quantum field in the action (3.1) we
obtain the following canonical form for the action,

$\displaystyle S=\frac{1}{2}\int\differential\tau\differential
x^{d-1}\sum_{\ell,m}\Big{[}\big{|}\frac{\partial Z_{\ell
m}}{\partial\tau}\big{|}^{2}-\big{|}\frac{\partial Z_{\ell m}}{\partial{\bf
x}}\big{|}^{2}+\Big{(}\frac{2f^{\prime
2}}{f^{2}}-\frac{f^{\prime\prime}}{f}+\frac{\ell(\ell+1)}{r^{2}}\big{(}1-\frac{2GM}{r}\big{)}\Big{)}|Z_{\ell
m}|^{2}\Big{]},$ (3.3)

in which a prime denotes a derivative with respect to conformal time
$\tau=r_{*}$. To obtain the above result, the following relations concerning
the spherical harmonics have been used,

$\displaystyle\int\differential\Omega\,Y_{\ell
m}(\theta,\phi)Y_{\ell^{\prime}m^{\prime}}^{*}(\theta,\phi)=\delta_{\ell\ell^{\prime}}\delta_{mm^{\prime}}\,,$
(3.4)

and,

$\displaystyle
Y_{lm}(\theta,\phi)^{*}=(-1)^{m}Y_{l\,-m}(\theta,\phi)\,,\quad\quad\sum_{m=-\ell}^{\ell}|Y_{\ell
m}(\theta,\phi)|^{2}=\frac{2\ell+1}{4\pi}\,.$ (3.5)

Since the spatial coordinate $-\infty<{\bf x}<\infty$ has the translation
invariance we can Fourier expand $Z_{\ell m}\propto e^{-i{\bf k}\cdot{\bf
x}}$. Correspondingly, the quantum mode $Z_{\ell m}(\tau,x)$ is expanded in
$(d-1)$-dimensional Fourier space as follows,

$\displaystyle Z_{\ell m}(\tau,{\bf
x})=\int_{-\infty}^{+\infty}\frac{\differential^{d-1}{\bf
k}}{(2\pi)^{(d-1)/2}}\left[e^{i{\bf k}\cdot{\bf x}}g_{k}^{\ell}(\tau)b_{{\bf
k}}^{\ell m}+(-1)^{m}e^{-i{\bf k}\cdot{\bf x}}{g_{k}^{\ell}}^{*}(\tau){b_{{\bf
k}}^{\ell\,-m}}^{\dagger}\right]\,.$ (3.6)

Here $b_{\bf k}^{\ell m}$ and ${b_{{\bf k}}^{\ell m}}^{\dagger}$ are the
annihilation and the creation operators defined for the observer inside the BH
which satisfy the usual commutation relations,

$\displaystyle[b_{\bf k}^{\ell m},{b_{{\bf
k}^{\prime}}^{\ell^{\prime}m^{\prime}}}^{\dagger}]=\delta_{\ell\ell^{\prime}}\delta_{mm^{\prime}}\delta({\bf
k}-{\bf k}^{\prime})\,,\quad[b_{\bf k}^{\ell m},b_{{\bf
k}^{\prime}}^{\ell^{\prime}m^{\prime}}]=[{b_{{\bf k}}^{\ell
m}}^{\dagger},{b_{{\bf
k}^{\prime}}^{\ell^{\prime}m^{\prime}}}^{\dagger}]=0\,.$ (3.7)

Note that the factor $(-1)^{m}$ appears in the second term of the mode
expansion to enforce the reality of $\Phi$.

From the action (3.3) the equation of the mode function in Fourier space is
obtained as follows,

$\displaystyle\frac{\differential^{2}g^{\ell}_{k}(\tau)}{\differential\tau^{2}}+\Big{(}k^{2}-V_{\mathrm{eff}}(\tau)\Big{)}g^{\ell}_{k}(\tau)=0\,,$
(3.8)

in which the effective potential is given by,

$\displaystyle V_{\mathrm{eff}}(\tau)=\frac{2f^{\prime
2}}{f^{2}}-\frac{f^{\prime\prime}}{f}+\frac{\ell(\ell+1)}{r^{2}}\big{(}1-\frac{2GM}{r}\big{)}\,.$
(3.9)

Eq. (3.8) is similar to the standard Regge-Wheeler perturbation equation for
the field in the exterior region of the BH [36]. However, in our setup where
the interior of BH is treated as a cosmological background, the above equation
can be interpreted as the extension of the Sasaki-Mukhanov equation to the
anisotropic cosmological background. Technically speaking, the mode function
in our analysis is time-dependent (i.e. the differential equation is with
respect to time $\tau$) while in the standard Regge-Wheeler equation the
relevant coordinate to solve the field equation is spacelike, the coordinate
$r_{*}$ in the exterior of BH.

With some efforts, and using the relations (2.5), one can show that

$\displaystyle
V_{\mathrm{eff}}(\tau)=\frac{2GM}{r(\tau)^{3}}\Big{(}1-\frac{(6-d)(2+d)GM}{8r(\tau)}\Big{)}+\frac{\ell(\ell+1)}{r(\tau)^{2}}\big{(}1-\frac{2GM}{r(\tau)}\big{)}\,.$
(3.10)

One can check that for $d=2$, this yields the standard form of
$V_{\mathrm{eff}}$ as in BH quasi-normal mode analysis [37, 38].

Note that $V_{\mathrm{eff}}(\tau)<0$. It goes to zero near the horizon
corresponding to $\tau\rightarrow-\infty$ and diverges near the singularity at
$\tau\rightarrow 0^{-}$. We try to solve Eq. (3.8) approximately as it can not
be solved analytically for the entire region inside BH for $-\infty<\tau<0$.
For the regions near the horizon the effective potential falls off rapidly and
the dominant term in the bracket in Eq. (3.8) is simply $k^{2}$. The mode
functions then are the plane wave solutions $g_{k}(\tau)\propto e^{\mp
i\omega\tau}$ in which $\omega\equiv|k|\geq 0$. Since we assume that the
perturbations are generated near the horizon and move forward towards the big
crunch singularity (future singularity), we take the positive frequency
solution $g_{k}(\tau)\propto e^{-i\omega\tau}$ as the incoming solution. On
the other hand, near the singularity the contribution of the term containing
$\ell(\ell+1)$ can be neglected compared to the first term in
$V_{\mathrm{eff}}(\tau)$ which is more divergent. As a result, as argued in
[27], we can neglect the contribution of the term $\ell(\ell+1)$ and the term
containing $r^{-3}$ in the effective potential and consider the most divergent
term $r^{-4}$.

Near the singularity $r^{2}\simeq-4GM\tau$ so within our approximation, Eq.
(3.8) simplifies to

$\displaystyle
g_{k}^{\prime\prime}(\tau)+\left(k^{2}+\frac{(6-d)(2+d)}{64\tau^{2}}\right)g_{k}(\tau)=0\,\quad\quad(\tau\rightarrow
0)\,,$ (3.11)

in which a prime denotes the derivative with respect to $\tau$. To simplify
the notation, we have dropped the superscript $\ell$ in $g_{k}^{\ell}(\tau)$
so effectively we work in the $s$-wave limit with $\ell=m=0$.

As mentioned before, Eq. (3.11) is similar to the Sasaki-Mukhanov equation in
cosmological perturbation theory. Therefore, interesting insights can be
obtained when compared to perturbations in FLRW background. We have two scale
factors $r(\tau)$ and $a(\tau)$, in which the latter is expanding as time goes
by (i.e. as we move forward towards the big crunch singularity). Since we have
decomposed the perturbations into spherical harmonics, Eq. (3.11) represents
the dynamics of perturbations in a 1+1 expanding universe with the scale
factor $a(\tau)$ given in Eq. (2.9). Similar to the studies of perturbations
in FLRW cosmology, we can divide the perturbations into long “superhorizon”
and short “subhorizon” perturbations. Perturbations which satisfy $k|\tau|\gg
1$ are initially inside the horizon. As we move towards big crunch and
$\tau\rightarrow 0$, $k|\tau|\ll 1$ and the mode becomes superhorizon.
Eventually, at the moment of big crunch, all modes become superhorizon.

Fortunately Eq. (3.11) can be solved with the solution,

$\displaystyle
g_{k}(\tau)=\sqrt{-\tau}\left[C_{1}(\omega)H_{\nu}^{(1)}(-\omega\tau)+C_{2}(\omega)H_{\nu}^{(2)}(-\omega\tau)\right]\,,\quad\quad\nu\equiv\frac{d-2}{8}\,,$
(3.12)

in which $H_{\nu}^{(1,2)}$ are the Hankel functions and
$C_{1}(\omega),C_{2}(\omega)$ are constants of integrations. As discussed
before, the solution (3.12) is an approximate solution since we have neglected
various subleading terms in the effective potential. However, if we manage to
connect the incoming plane wave solution to the near singularity solution
(3.12) then we can be hopeful that the resultant solution mimics the true
solution to reasonable accuracy. This is possible because the Hankel functions
asymptotically approach to the plane wave solutions $e^{\mp i\omega\tau}$.
This suggests that the solution (3.12) is a good approximation for both
regions $\tau\rightarrow 0$ and $\tau\rightarrow-\infty$ with some
modifications in the intermediate region. With this approximation, we take the
analytic solution (3.12) to be qualitatively valid for the entire region
$-\infty<\tau<0$.

Our job is now to fix the constants $C_{1}(\omega),C_{2}(\omega)$. This can be
achieved by imposing the equal time quantum commutation relation between the
field $Z$ and its conjugate momentum $\partial_{\tau}Z$,

$\displaystyle[Z(\tau,x),\partial_{\tau}Z(\tau,x^{\prime})]=i\delta(x-x^{\prime})\,.$
(3.13)

Imposing the above quantization condition in our mode function (3.12) yields
the following relation,

$\displaystyle\big{|}C_{1}(\omega)\big{|}^{2}-\big{|}C_{2}(\omega)\big{|}^{2}=\frac{\pi}{4}\,.$
(3.14)

The above constraint does not uniquely fix $C_{1}(\omega)$ and
$C_{2}(\omega)$. However, the above equation suggests the natural possibility
that $C_{1}=\sqrt{\pi}/2$ and $C_{2}=0$. Indeed one can fix $C_{1}$ and
$C_{2}$ definitely only considering the initial conditions near the horizon
where the perturbations are generated. The initial condition was fixed by the
positive frequency mode with $g_{k}(\tau)\propto e^{-i\omega\tau}$. Using the
asymptotic form of the Hankel function,

$\displaystyle
H_{\nu}^{(1)}(-\omega\tau)\rightarrow\sqrt{\frac{-2}{\pi{\omega}\tau}}e^{-i\big{(}\omega\tau+\frac{\pi}{4}{(1+2\nu)}\big{)}}\quad\quad(\tau\rightarrow-\infty)\,,$
(3.15)

one can check that only the term containing $C_{1}$ in Eq. (3.12) can match to
the positive frequency mode while the term containing $C_{2}$ match to the
negative frequency mode $e^{i\omega\tau}$. This justifies our conclusion that
$C_{2}=0$. This conclusion is similar to the standard “Bunch-Davies” vacuum in
a dS background which carries the lowest energy among all possible vacua. Of
course, if we consider a mixed initial condition with both positive and
negative frequencies, then one can consider the general case in which both
$C_{1}(\omega)$ and $C_{2}(\omega)$ are present, similar to the “non Bunch-
Davies” vacuum in inflation. With $C_{2}=0$, we finally obtain the mode
function as follows

$\displaystyle
g_{k}(\tau)=\sqrt{\frac{-\pi\tau}{4}}H_{\nu}^{(1)}(-\omega\tau)\,.$ (3.16)

In summary, putting all together, the quantum perturbation $\Phi(x^{\mu})$
from Eq. (3.2) is obtained to be

$\displaystyle\Phi(x^{\mu})=\frac{a(\tau)^{\frac{2-d}{2}}}{r(\tau)}\sum_{\ell,m}\int\frac{\differential^{d-1}{\bf
k}}{(2\pi)^{(d-1)/2}}\left[e^{i{\bf k}\cdot{\bf x}}g_{k}(\tau)b_{{\bf
k}}^{\ell m}+(-1)^{m}e^{-i{\bf k}\cdot\bf{x}}g^{*}_{k}(\tau){b_{{\bf
k}}^{\ell\,-m}}^{\dagger}\right]Y_{\ell m}(\theta,\phi)\,,$ (3.17)

with the mode function $g_{k}(\tau)$ given in Eq. (3.16).

##  4 The Regularization

Armed with the mode function given in Eq. (3.17), we can calculate the energy
momentum tensor associated to quantum perturbations. Technically speaking,
these quantum corrections correspond to one-loop correction from vacuum bubble
diagrams. Similar studies were performed for dS spacetime in [39, 40].

The energy momentum tensor for the perturbations $\Phi(x^{\mu})$ is given by

$\displaystyle
T^{\mu}_{\nu}=\partial^{\mu}\Phi\partial_{\nu}\Phi-\frac{1}{2}\delta^{\mu}_{\nu}\,\partial^{\alpha}\Phi\partial_{\alpha}\Phi\,.$
(4.1)

In this section as our emphasis is on the quantum properties and the
regularizations of the UV divergences, we consider only the energy density
$\rho=-T^{0}_{0}$. The remaining components of the anisotropic energy density
will be considered in details in next section where we study the backreactions
induced by these quantum perturbations.

Constructing the energy density, it is given by three contributions,

$\displaystyle\rho=\rho_{1}+\rho_{2}+\rho_{3}\,,$ (4.2)

in which

$\displaystyle\rho_{1}\equiv\frac{1}{2a(\tau)^{2}}\big{(}\partial_{\tau}\Phi(x^{\mu})\big{)}^{2}\,,$
(4.3) $\displaystyle\rho_{2}\equiv\frac{1}{2a(\tau)^{2}}\big{(}\partial_{\bf
x}\Phi(x^{\mu})\big{)}^{2}\,,$ (4.4)

and

$\displaystyle\rho_{3}\equiv\frac{1}{2r(\tau)^{2}}\gamma^{ab}\partial_{a}\Phi(x^{\mu})\partial_{b}\Phi(x^{\mu})\,.$
(4.5)

Before going into technicalities of calculating each of the above components,
it is instructive to compare the magnitude of the above components near the
singularity region, $r\rightarrow 0$. Since both $\rho_{1}$ and $\rho_{2}$
deal with derivatives along the extended direction then on physical grounds we
expect that $\rho_{1}\sim\rho_{2}$. On the other hand, $\rho_{3}$ involves
derivative along the compact directions so intuitively speaking
$\rho_{3}\sim(a^{2}/r^{2})\rho_{2}$. Using the relations (2.10) we conclude
that $\rho_{3}\sim a(\tau)^{2}r(\tau)^{-2}\rho_{2}\sim(-\tau/GM)^{1/2}$.
Correspondingly, we expect that near the singularity $\rho_{3}$ to be
suppressed compared to $\rho_{1}$ or $\rho_{2}$ by a factor
$(-\tau/GM)^{1/2}$.

We are interested in quantum expectation values such as $\langle\rho\rangle$.
As is well known, there is no unique notion of the vacuum in a curved
background such as in BH background where the spacetime contains an event
horizon. In our analysis, we restrict ourselves to the interior of the BH with
the mode function given in Eq. (3.17). The vacuum for the observer inside the
BH is defined by $b_{\bf k}|0\rangle=0$. With our choice $C_{2}=0$, this
vacuum is like the Bunch-Davies vacuum in inflationary setup. Of course, since
the coordinate $(r,\bar{t})$ or equivalently $(\tau,x)$ do not cover the whole
manifold, this vacuum is not well-defined globally. For example, one may
instead use the vacuum associated to the Kruskal-Szekeres coordinate, i.e. the
Hartle-Hawking vacuum which is well-defined on both sides of the horizon. The
difference in vacua employed by various observers is the prime reason for the
existence of Hawking radiation in BH background. In this work, we exclusively
concentrate to the interior region with the vacuum $|0\rangle$ described above
which is well-defined for all interior observers.

Let us start with $\langle\rho_{2}\rangle$ which is easier to calculate. Also,
to see the structure of divergences, let us for the moment take $d=2$ (i.e. no
dimensional regularization yet). Using the mode function (3.17) and performing
the summation over $m$ as given in Eq. (3.5), we obtain

$\displaystyle\langle\rho_{2}(\tau)\rangle=\frac{1}{2\pi
r(\tau)^{2}a(\tau)^{2}}\sum_{\ell=0}^{\infty}(2\ell+1)\int_{0}^{\infty}\frac{d\omega}{2\pi}\omega^{2}|g_{\omega}(-\omega\tau)|^{2}\quad\quad(d=2)\,.$
(4.6)

As usual in QFT analysis, the above expectation value suffers from the UV
divergences. Since we decomposed the spatial directions into the form $R\times
S^{2}$, we have two different types of UV divergences. The first divergence is
the infinite sum $\sum_{\ell=0}^{\infty}(2\ell+1)$ associated to the azimuthal
directions while the other divergence is due to the integral of the high
frequency mode $\omega$. Fortunately, we are able to regularize both
divergences as follows.

The divergence associated to the infinite sum
$\sum_{\ell=0}^{\infty}(2\ell+1)$ can be dealt with using the zeta function
regularization [41, 42] scheme. More specifically, we have

$\displaystyle\sum_{\ell=0}^{\infty}(2\ell+1)=-\zeta(-1)=\frac{1}{12}\,.$
(4.7)

The UV divergence over the high frequency mode can be regulated via
dimensional regularization as we advertised before. For this, we replace

$\displaystyle\differential^{d-1}{\bf
k}\rightarrow\omega^{d-2}\mu^{2-d}\differential\omega\,\differential^{d-2}\Omega\,,$
(4.8)

in which $\differential^{d-2}\Omega$ represents the volume element of a unit
$(d-2)$-dimensional sphere with the volume

$\displaystyle
V_{d-2}\equiv\int\differential^{d-2}\Omega=\frac{2\pi^{\frac{d-1}{2}}}{\Gamma(\frac{d-1}{2})}\,,$
(4.9)

and $d=2-\epsilon$ in which $\epsilon\rightarrow 0$ is the small parameter in
dimensional regularization expansion. The mass scale $\mu$ is introduced to
keep track of the correct dimension of the physical quantities. Putting all
together, we obtain

$\displaystyle\langle\rho_{2}(\tau)\rangle=\frac{-\zeta(-1)\mu^{2-d}}{2\pi
r(\tau)^{2}a(\tau)^{d}}V_{d-2}\int_{0}^{\infty}\frac{d\omega}{(2\pi)^{d-1}}\omega^{d}|g_{\omega}(-\omega\tau)|^{2}\,.$
(4.10)

The above integral can be performed analytically, yielding†††We use the Maple
computational software to calculate the integrals.

$\displaystyle\langle\rho_{2}(\tau)\rangle=\frac{-\zeta(-1)\,(d-1)\cos(\pi\nu)}{2^{d+2}\pi^{2+\frac{d}{2}}a(\tau)^{d}r(\tau)^{2}(-\tau)^{d}}\mu^{2-d}\Gamma\big{(}\frac{-d}{2}\big{)}\Gamma\big{(}\frac{5d+2}{8}\big{)}\Gamma\big{(}\frac{3d+6}{8}\big{)}\,.$
(4.11)

The divergence of the above expression appears in
$\Gamma\big{(}\frac{-d}{2}\big{)}$ when $d\rightarrow 2$. Now, expanding to
leading order in $d=2-\epsilon$ we obtain

$\displaystyle\langle\rho_{2}(\tau)\rangle=\frac{\zeta(-1)}{32a(\tau)^{2}r(\tau)^{2}\pi^{2}\tau^{2}}\Big{(}\frac{1}{\epsilon}+\ln\big(-\mu\tau
a(\tau)\big{missing})\Big{)}+{\cal O}(\epsilon)\,,$ (4.12)

in which we have absorbed some numerical constants into the mass scale
parameter $\mu$ in the final result. As expected, there is a pole at
$\epsilon=0$ which highlights the UV divergence which should be regularized by
appropriate counter terms as in standard QFT analysis. In addition, the
appearance of the logarithmic term $\ln\big(-\mu\tau\big{missing})$ is
expected which is the hallmark of the dimensional regularization scheme.

Discarding the divergence term via regularization, the regularized value of
$\langle\rho_{2}(\tau)\rangle$ is obtained to be

$\displaystyle\langle\rho_{2}(\tau)\rangle_{\mathrm{reg}}=\frac{\zeta(-1)}{32a(\tau)^{2}r(\tau)^{2}\pi\tau^{2}}\ln(-\mu\tau
a(\tau))\,.$ (4.13)

The above expression is valid for all region inside the BH. We are mostly
interested in the value of the energy momentum tensor near the singularity
$\tau\rightarrow 0$ where the expectation values like $\langle\rho\rangle$
grow with negative power of $\tau$. Using near singularity relations (2.10),
we obtain

$\displaystyle\langle\rho_{2}(\tau)\rangle_{\mathrm{reg}}\equiv
c_{2}\ln(-\mu\tau)(-\tau)^{-\frac{5}{2}}\,,$ (4.14)

in which $c_{2}\propto(GM)^{-\frac{3}{2}}$ is a constant involving $\zeta(-1)$
and other numerical factors. In addition, we have absorbed a factor
$(GM)^{1/4}$ into $\mu$ which is an arbitrary scale at this level. We observe
that near the singularity $\langle\rho_{2}(\tau)\rangle_{\mathrm{reg}}$ grows
like $(-\tau)^{-\frac{5}{2}}\sim r^{-5}$ with logarithmic running. The factor
$\ln(-\mu\tau)$ is the hallmark of the quantum corrections.

Now we calculate $\langle\rho_{1}(\tau)\rangle$ which from Eq. (4.3) is given
by

$\displaystyle\langle\rho_{1}(\tau)\rangle=\frac{-\zeta(-1)\mu^{2-d}}{2\pi
a(\tau)^{2}}V_{d-2}\int_{0}^{\infty}\frac{d\omega\,\omega^{d-2}}{(2\pi)^{d-1}}\Bigg{|}\frac{\differential}{\differential\tau}\Big{(}\frac{a(\tau)^{\frac{2-d}{2}}}{r(\tau)}g_{\omega}(-\omega\tau)\Big{)}\Bigg{|}^{2}\,,$
(4.15)

in which we have used Eq. (4.7) to perform the sum over $\ell$ as before.
Performing the integral using the dimensional regularization approach with the
mode function given in Eq. (3.16), and considering near the singularity
region, we obtain

$\displaystyle\langle\rho_{1}(\tau)\rangle=\frac{3(2+d)}{6-5d}\langle\rho_{2}(\tau)\rangle\,.$
(4.16)

Regularizing $\langle\rho_{1}(\tau)\rangle$ by discarding the divergent
$\epsilon^{-1}$ term we obtain

$\displaystyle\langle\rho_{1}(\tau)\rangle_{\mathrm{reg}}=c_{1}\ln(-\mu\tau)(-\tau)^{-\frac{5}{2}}\,,$
(4.17)

in which $c_{1}$ is a numerical constant with the dimension $c_{1}\sim
c_{2}\sim(GM)^{-\frac{3}{2}}$. Naively, Eq. (4.16) may suggest that
$c_{1}=-3c_{2}$ when $d=2$. However, this is misleading since we should
actually set $d=2-\epsilon$ which brings an additional contribution. More
specifically, we have

$\displaystyle\langle\rho_{1}(\tau)\rangle_{\mathrm{reg}}+3\langle\rho_{2}(\tau)\rangle_{\mathrm{reg}}=\frac{3\zeta(-1)}{128\pi^{2}(GM)^{\frac{3}{2}}(-\tau)^{\frac{5}{2}}}\,.$
(4.18)

Now we calculate $\langle\rho_{3}\rangle$. For this purpose, first note that
$\gamma^{ab}\nabla_{a}\nabla_{b}\Phi=-\ell(\ell+1)\Phi$ in which $\gamma_{ab}$
is the metric on a unit two-sphere and $\\{a,b\\}=\\{\theta,\phi\\}$. This
yields

$\displaystyle\langle\rho_{3}\rangle$ $\displaystyle=$
$\displaystyle\frac{1}{2r(\tau)^{2}}\gamma^{ab}\langle\partial_{a}\Phi\partial_{b}\Phi\rangle$
(4.19) $\displaystyle=$ $\displaystyle\frac{1}{2\pi
r(\tau)^{4}a(\tau)^{d-2}}\sum_{\ell=0}^{\infty}(2\ell+1)\ell(\ell+1)\int_{0}^{\infty}\frac{d\omega\omega^{d-2}}{(2\pi)^{d-1}}|g_{\omega}(-\omega\tau)|^{2}\,.$

There are a number of important differences compared to the case of
$\langle\rho_{2}\rangle$. First, we have a non-trivial sum over $\ell$ which
requires more sophisticated form of $\zeta$ function regularization scheme
[43, 42]. However, the important difference is that we have a weaker
dependence on $a(\tau)$ while the dependence on $r(\tau)$ is stronger. This is
because $a(\tau)$ is the scale factor along the extended spatial direction
while $r(\tau)$ represents the scale factor along the compact dimensions. Now,
using the scalings of $a(\tau)$ and $r(\tau)$ near the singularity given in
Eq. (2.10), we find that $\langle\rho_{3}\rangle\propto(GM)^{-2}\tau^{-2}$.
Upon performing the zeta regularization over the sum over $\ell$ and
dimensional regularization over the frequency,
$\langle\rho_{3}(\tau)\rangle_{\mathrm{reg}}$ is obtained to have the
following form,

$\displaystyle\langle\rho_{3}(\tau)\rangle_{\mathrm{reg}}=c_{3}\ln(-\mu\tau)\,\tau^{-2}\,,$
(4.20)

in which $c_{3}\propto(GM)^{-2}$.

Comparing the scaling of $\langle\rho_{3}\rangle$ to $\langle\rho_{1}\rangle$
or $\langle\rho_{2}\rangle$ we observe that the $\tau$-dependence in
$\langle\rho_{3}\rangle$ is milder compared to latter by a factor
$(-\tau/GM)^{1/2}$ as concluded previously based on qualitative arguments.
Therefore, near the singularity region we can neglect the contributions of
$\langle\rho_{3}\rangle$ compared to those of $\langle\rho_{1}\rangle$ and
$\langle\rho_{2}\rangle$.

Armed with the results for $\langle\rho_{i}\rangle$ obtained above, we can
calculate their backreactions on the background metric.

##  5 Anisotropic $T^{\mu}_{\nu}$ and Backreactions

In the previous section we have calculated
$\langle\rho\rangle_{\mathrm{reg}}$. To complete the study, we need to
calculate other components of the energy momentum tensor. Since we are done
with the regularization, we set $d=2$.

As the background is homogenous but anisotropic, the energy momentum tensor
induced by quantum fluctuations takes the form of an imperfect fluid, given by
[44]

$\displaystyle
T_{\mu\nu}=(\rho+P)u_{\mu}u_{\nu}+Pg_{\mu\nu}+q_{\mu}u_{\nu}+q_{\nu}u_{\mu}+\pi_{\mu\nu}\,,$
(5.1)

with the additional conditions

$\displaystyle u^{\mu}u_{\mu}=-1,\quad
q_{\mu}u^{\mu}=0,\quad\pi^{\mu}_{\mu}=0,\quad\pi_{\mu\nu}=\pi_{\nu\mu},\quad\pi_{\mu\nu}u^{\mu}=0\,.$
(5.2)

In the above decomposition, $\rho$ is the energy density, $P$ is the isotropic
pressure, $u^{\mu}$ is the fluid’s comoving four velocity, $\pi_{\mu\nu}$ is
the anisotropic pressure (stress) and $q^{\mu}$ is the heat conduction which
also measures the energy flux relative to $u^{\mu}$. In our setup we expect no
heat flow for the fluid so $q_{\mu}=0$.

From the symmetry of the cosmological background, we obtain

$\displaystyle u^{\mu}=(a(\tau),0,0,0)\,.$ (5.3)

In the absence of the heat flow, the energy momentum tensor takes the diagonal
form

$\displaystyle
T^{\mu}_{\nu}=\mathrm{diag}\left(-\rho,P+\pi^{1}_{1},P+\pi^{2}_{2},P+\pi^{3}_{3}\right)\,.$
(5.4)

As the compact dimensions are spherically symmetric, we impose
$\pi^{2}_{2}=\pi^{3}_{3}$. Combined with the traceless condition,
$\pi^{\mu}_{\mu}=0$, we obtain

$\displaystyle
T^{\mu}_{\nu}=\mathrm{diag}\left(-\rho,P-2\Pi,P+\Pi,P+\Pi\right)\,,$ (5.5)

where the anisotropic pressure is now controlled by a single scalar quantity,
$\Pi\equiv\pi^{2}_{2}$.

We are interested in vacuum expectation values of $\langle
T^{\mu}_{\nu}\rangle$. From the above form of $T^{\mu}_{\nu}$ we need to
calculate $\langle\rho\rangle$, $\langle P\rangle$ and $\langle\Pi\rangle$.
The first quantity was already calculated in previous section as

$\displaystyle\langle\rho\rangle=\langle\rho_{1}\rangle+\langle\rho_{2}\rangle+\langle\rho_{3}\rangle\,.$
(5.6)

As we shall see, $\langle P\rangle$ and $\langle\Pi\rangle$ can be expressed
in terms of $\langle\rho_{i}\rangle$ as well.

Let us consider the azimuthal components of $\langle T^{\mu}_{\nu}\rangle$.
From the spherical symmetry of the azimuthal directions we have
$\langle\nabla_{a}\Phi\nabla_{b}\Phi\rangle=\kappa\gamma_{ab}$ in which
$\kappa$ is a constant [45]. To find the value of $\kappa$, contract this
relation by the inverse metric $\gamma^{ab}$, yielding

$\displaystyle\langle\nabla_{a}\Phi\nabla_{b}\Phi\rangle=\gamma_{ab}\langle\rho_{3}\rangle\,.$
(5.7)

From this relation, we obtain

$\displaystyle\langle
T^{a}_{b}\rangle=\delta^{a}_{b}\Big{(}\langle\rho_{1}\rangle-\langle\rho_{2}\rangle\Big{)}\,,\quad\quad\\{a,b\\}=\\{\theta,\phi\\}\,.$
(5.8)

On the other hand, for $\langle T^{x}_{x}\rangle$ component we have

$\displaystyle\langle
T^{x}_{x}\rangle=\langle\rho_{1}\rangle+\langle\rho_{2}\rangle-\langle\rho_{3}\rangle\,.$
(5.9)

Comparing Eqs. (5.9) and (5.8) with the general form of anisotropic
$T^{\mu}_{\nu}$ given in Eq. (5.5), we obtain

$\displaystyle\langle
P\rangle=\langle\rho_{1}\rangle-\frac{1}{3}\langle\rho_{2}\rangle-\frac{1}{3}\langle\rho_{3}\rangle\,,$
(5.10)

and

$\displaystyle\langle\Pi\rangle=-\frac{2}{3}\langle\rho_{2}\rangle+\frac{1}{3}\langle\rho_{3}\rangle\,.$
(5.11)

Our results for the values of $\langle T^{\mu}_{\nu}\rangle$ were general so
far. Now, we consider our limit of interests where $\tau\rightarrow 0$ with
$\langle\rho_{i}\rangle_{\mathrm{reg}}$ as obtained in the previous section.
In this limit, the regularized components of the $\langle
T^{\mu}_{\nu}\rangle$ are given by

$\displaystyle\langle\rho\rangle_{\mathrm{reg}}=\Big{[}(c_{1}+c_{2})(-\tau)^{\frac{-5}{2}}+c_{3}\tau^{-2}\Big{]}\ln(-\mu\tau)\,,$
(5.12) $\displaystyle\langle
P\rangle_{\mathrm{reg}}=\Big{[}(c_{1}-\frac{c_{2}}{3})(-\tau)^{\frac{-5}{2}}-\frac{c_{3}}{3}\tau^{-2}\Big{]}\ln(-\mu\tau)\,,$
(5.13)

and

$\displaystyle\langle\Pi\rangle_{\mathrm{reg}}=\Big{[}-\frac{2c_{2}}{3}(-\tau)^{\frac{-5}{2}}+\frac{c_{3}}{3}\tau^{-2}\Big{]}\ln(-\mu\tau)\,.$
(5.14)

###  5.1 Backreactions

Having calculated the components of the energy momentum tensor generated from
the quantum perturbations, we can calculate the backreactions of this energy
momentum tensor on the geometry.

The Einstein field equations in the presence of the quantum backreactions are
given as usual by [46]

$\displaystyle R^{\mu}_{\nu}=8\pi G\langle
S^{\mu}_{\nu}\rangle_{\mathrm{reg}}\,,\quad\quad S^{\mu}_{\nu}\equiv
T^{\mu}_{\nu}-\frac{1}{2}\delta^{\mu}_{\nu}T\,.$ (5.15)

In the absence of backreaction, $T^{\mu}_{\nu}=0$ so one deals with the vacuum
Einstein field equations which of course yield the background Schwarzschild
solution (2.1). In the presence of backreactions, we expect this solution to
be modified. As we still assume the spherical symmetry to hold in the presence
of backreactions, it is convenient to use a metric ansatz which takes the
spherical symmetry into account. Motivated by the metric (2.1) for the
exterior of the BH, we propose the following time-dependent spherical
symmetric metric for the interior of the BH (see also [28, 47] for a similar
ansatz),

$\displaystyle ds^{2}=-\frac{dt^{2}}{A(t)}+B(t)dx^{2}+t^{2}d\Omega^{2}\,.$
(5.16)

In the above ansatz $t$ plays the roles of $r$ in the exterior region. In
particular, in the absence of backreactions, the solution for the functions
$A(t)$ and $B(t)$ are

$\displaystyle A(t)=B(t)=-1+\frac{2GM}{t}\,,\quad\quad(T^{\mu}_{\nu}=0)\,.$
(5.17)

Of course, instead of the metric ansatz (5.16), we could have chosen a
different cosmological ansatz, such as

$\displaystyle
ds^{2}=a_{1}(\tau)^{2}(-d\tau^{2}+dx^{2})+a_{2}(\tau)^{2}d\Omega^{2}\,.$
(5.18)

This metric is the extension of the metric (2.8) to the case where the
backreactions are included with two unknown scale factors $a_{1}(\tau)$ and
$a_{2}(\tau)$. However, it turns out that the field equations are easier with
the ansatz (5.16) which we consider in the following analysis.

From the combinations $2R^{\theta}_{\theta}-R^{x}_{x}-R^{t}_{t}$ and
$R^{x}_{x}-R^{t}_{t}$ of the Einstein fields equation associated to the metric
(5.16) we obtain the following equations,

$\displaystyle\frac{1}{t^{2}}\Big{(}t\dot{A}+A+1\Big{)}$ $\displaystyle=$
$\displaystyle 8\pi G\Big{(}2\langle\Pi\rangle_{\mathrm{reg}}-\langle
P\rangle_{\mathrm{reg}}\Big{)},$ (5.19)
$\displaystyle\frac{A}{t}\Big{(}\frac{\dot{B}}{B}-\frac{\dot{A}}{A}\Big{)}$
$\displaystyle=$ $\displaystyle 8\pi
G\Big{(}\langle\rho\rangle_{\mathrm{reg}}+\langle
P\rangle_{\mathrm{reg}}-2\langle\Pi\rangle_{\mathrm{reg}}\Big{)}\,.$ (5.20)

One can easily check that in the absence of the backreaction the solutions of
the above equations are given by Eq. (5.17) as expected. We can solve the
above equations perturbatively as follows:

$\displaystyle A(t)=\bar{A}{(t)}+\delta A(t)\,,\quad\quad
B(t)=\bar{B}{(t)}+\delta B(t)\,,$ (5.21)

in which $\bar{A}=\bar{B}$ are the background solutions given in Eq. (5.17).
Starting with Eq. (5.19), we obtain $\delta A$ as follows,

$\displaystyle\delta A=\frac{8\pi G}{t}\int^{t}dt^{\prime}t^{\prime
2}\Big{[}2\langle\Pi(t^{\prime})\rangle_{\mathrm{reg}}-\langle
P(t^{\prime})\rangle_{\mathrm{reg}}\Big{]}\,.$ (5.22)

Now using this solution into Eq. (5.20) we obtain

$\displaystyle\delta B=\delta A+8\pi
G\bar{A}(t)\int^{t}dt^{\prime}\frac{t^{\prime}}{\bar{A}(t^{\prime})}\Big{[}\langle\rho(t^{\prime})\rangle_{\mathrm{reg}}+\langle
P(t^{\prime})\rangle_{\mathrm{reg}}-2\langle\Pi(t^{\prime})\rangle_{\mathrm{reg}}\Big{]}\,.$
(5.23)

The above solutions for $\delta A$ and $\delta B$ were general with no
assumption on the source terms in the right hand sides of Eqs. (5.22) and
(5.23). Now, we consider our case of interests with
$\langle\rho\rangle_{\mathrm{reg}},\langle P\rangle_{\mathrm{reg}}$ and
$\langle\Pi\rangle_{\mathrm{reg}}$ given in Eqs. (5.12), (5.13) and (5.14). To
use these expressions for the source terms, we need to express $\tau$ in terms
of $t$. From Eq. (2.10) we note that $\tau\simeq-t^{2}/4GM$ (note that in our
treatment of the interior of the BH, $t$ plays the same roles as the
coordinate $r$ outside the BH). Correspondingly, the source terms are obtained
to be

$\displaystyle 2\langle\Pi\rangle_{\mathrm{reg}}-\langle
P\rangle_{\mathrm{reg}}=\Big{[}-(c_{1}+c_{2})(-\tau)^{-5/2}+c_{3}\tau^{-2}\Big{]}\ln(-\mu\,\tau)\simeq\alpha\ln(\mu\,t)\,t^{-5}\,,$
(5.24)

in which $\alpha\propto(c_{1}+c_{2})\sim GM$ is a constant and we have
rescaled the undetermined scale $\mu$ as well. We have neglected the
subleading term containing $c_{3}$ as long as $c_{1}\neq-c_{2}$ which is the
case in our analysis. Similarly, for the other source term, we obtain

$\displaystyle\langle\rho\rangle_{\mathrm{reg}}+\langle
P\rangle_{\mathrm{reg}}-2\langle\Pi\rangle_{\mathrm{reg}}=2(c_{1}+c_{2})\ln(-\mu\,\tau)(-\tau)^{-5/2}=-2\alpha\ln(\mu\,t)\,t^{-5}\,.$
(5.25)

With the above source terms at hand, we can calculate $\delta A$ and $\delta
B$. Considering the near singularity limit $t\rightarrow 0$, from Eq. (5.22)
we obtain,

$\displaystyle\delta A\simeq\frac{8\pi
G\alpha}{t}\int^{t}dt^{\prime}\,\ln(\mu\,t^{\prime})t^{\prime-3}\simeq-\frac{4\pi
G\alpha}{t^{3}}\ln(\mu t)\,,$ (5.26)

in which we have absorbed some numerical constant in the scale $\mu$ as well.
Using this solution in equation (5.23), we obtain

$\displaystyle\delta B\simeq\delta A-\frac{16\pi
G\alpha}{t}\int^{t}dt^{\prime}t^{\prime-3}\simeq\frac{4\pi
G\alpha}{t^{3}}\ln(\mu t)\,.$ (5.27)

In conclusion, the metric functions to leading order in backreaction effects
are given by

$\displaystyle A(t)$ $\displaystyle\simeq$
$\displaystyle-1+\frac{2GM}{t}-\frac{4\pi G\alpha}{t^{3}}\ln(\mu t)\,,$ (5.28)
$\displaystyle B(t)$ $\displaystyle\simeq$
$\displaystyle-1+\frac{2GM}{t}+\frac{4\pi G\alpha}{t^{3}}\ln(\mu t)\,.$ (5.29)

with $\alpha\propto GM$. As mentioned before, the appearance of the
logarithmic correction $\ln(\mu t)$ is reassuring, highlighting the nature of
the regularized quantum corrections.

###  5.2 Physical Implications

Having obtained the backreaction in the metric, one may ask for some physical
implications.

One question to ask is what is the threshold where the backreaction
contributions become comparable to the background solution, $\delta
A\sim\bar{A}$? For this purpose, let us consider the metric along the
azimuthal direction with the line element $t^{2}d\Omega^{2}$. This may define
the shortest distance (radius of two-sphere near the singularity) that we can
trust our approximate solution before the backreactions overcome the
background solution. Otherwise, this can define the scale in which new physics
may enter to resolve the singularity. Defining this radius (actually the time
scale towards singularity) by $t_{c}$, it is given by

$\displaystyle\frac{2GM}{t_{c}}\sim\frac{4\pi G\alpha}{t_{c}^{3}}\ln(\mu
t_{c})\quad\rightarrow\quad
t_{c}\sim\big{(}\frac{\alpha}{M}\big{)}^{\frac{1}{2}}\,,$ (5.30)

in which in obtaining the final estimation, we have neglected the subleading
logarithmic correction involving the scale $\mu$. With $\alpha\sim GM$, this
timescale is obtained to be $t_{c}\sim\sqrt{G}\sim 1/M_{P}$, i.e. the Planck
time. This is somewhat expected as $M_{P}$ defines the scale of quantum
gravity. Of course the above estimation is rough as we have considered only
near the singularity limit of our approximate solution. Treating the analysis
more carefully, one may obtain a somewhat larger scale for $t_{c}$ prior to
the singularity when the backreaction becomes important. This requires a
careful investigation which is beyond the scope of the current analysis.

The other interesting question is what the exterior observer interpret the
corrections in the metric? Based on our starting discussion, for an exterior
observer the roles of the timelike and spacelike coordinates are switched.
This suggests the speculation that for the exterior observer in which
$t\rightarrow r$ and $x\rightarrow\bar{t}$, the line element with the quantum
corrections is given by

$\displaystyle ds^{2}=-{\cal
B}(r)\differential{\bar{t}\,}^{2}+\frac{\differential r^{2}}{{\cal
A}(r)}+r^{2}d\Omega^{2}\,\quad\quad\quad(r>2GM)\,,$ (5.31)

with the functions ${\cal A}(r)$ and ${\cal B}(r)$ are given as follows

$\displaystyle{\cal A}(r)=-A(r)$ $\displaystyle\simeq$ $\displaystyle
1-\frac{2GM}{r}+\frac{4\pi G\alpha}{r^{3}}\ln(\mu r)\,,$ (5.32)
$\displaystyle{\cal B}(r)=-B(r)$ $\displaystyle\simeq$ $\displaystyle
1-\frac{2GM}{r}-\frac{4\pi G\alpha}{r^{3}}\ln(\mu r)\,.$ (5.33)

Interestingly, the above speculative results are in agreement with the results
for the quantum backreactions for the exterior region [48, 49, 50, 51] who
obtained $1/r^{3}$ profile with possible logarithmic running for the metric
components. This similarity is encouraging, indicating that the effects of
backreactions from the quantum fields in the interior of the BH are genuine.

An interesting corollary from the above result is that the quantum corrections
in the interior of BH can be probed by the exterior observer. Indeed, the
exterior observer equipped with the quantum modified metric (5.31) can
observationally distinguish this metric from the classical Schwarzschild
metric. This is intriguing, implying that the effects occurring inside the
event horizon is accessible to the exterior observer. Of course, like in
Hawking radiation, this is a pure quantum effect. Having said this, the
quantum corrections in Eqs. (5.32) and (5.33) are so small that they seem
impossible to be probed observationally.

With the quantum corrections in background metric included, it is insightful
to look at various of its physical effects. First, the quantum corrections to
the metric will modify the position of the event horizon. The new position of
the event horizon $R_{S}^{\prime}$ from the solution of ${\cal B}(r)=0$ is
given by

$\displaystyle R_{S}^{\prime}\simeq R_{S}-\frac{\pi\alpha}{GM^{2}}\ln(2\mu
GM)\,.$ (5.34)

With $\alpha\sim GM$ the fractional corrections in the position of the event
horizon is at the order $(MR_{S})^{-1}\sim(M_{P}/M)^{2}$ which is extremely
small and practically unobservable for astrophysical BHs. Since the fractional
correction induced by parameter $\alpha$ is tiny, the quantum corrections can
not modify the number of event horizons so there is only one event horizon.

The surface gravity associated to the metric (5.31) is given by [52, 53]

$\displaystyle\kappa=\frac{{\cal B}^{\prime}}{2}\sqrt{\frac{{\cal A}}{{\cal
B}}}\,\Big{|}_{R_{S}}\,,$ (5.35)

which to leading order in $\alpha$, is obtained to be

$\displaystyle\kappa\simeq\frac{1}{4GM}+\frac{\pi\alpha}{8G^{3}M^{4}}\big{(}3\ln(2\mu
GM)-1\big{)}\,.$ (5.36)

Correspondingly, the fractional corrections in surface gravity is at the order
$\alpha/G^{2}M^{3}\sim(M_{P}/M)^{2}$ which is extremely small for
astrophysical BHs. Related to this issue, note that since the fractional
corrections in event horizon and surface gravity is very small, we can not
hope to have a extremal black hole with a vanishing surface gravity.
Physically, this is because the parameter $\alpha$ is not independent of $M$
as $\alpha\sim GM$. Therefore, we can not expect to trade $M$ with $\alpha$ as
an independent parameter to obtain an extremal black hole.

We can calculate the ADM mass associated to the metric (5.31) as well. Since
the quantum correction falls off like $\ln(\mu r)/r^{3}$, it is expected that
its contribution in surface integrals at infinity in calculating the ADM mass
to vanish. Indeed, performing the analysis [54], we confirmed that the ADM
mass is simply $M$.

Another question is whether the quantum corrections can resolve the
singularity at the center of BH. Since we have treated the quantum effects as
perturbative backreactions, we do not expect that they can resolve the
singularity at this stage. Indeed, calculating the Ricci scalar associated to
the interior metric (5.16) with the metric functions (5.28) and (5.29) for
$t\rightarrow 0$ we obtain

$\displaystyle R\simeq 4\pi G\alpha\big{(}3-2\ln(\mu t)\big{)}t^{-5}\,.$
(5.37)

Similarly, calculating the Kretschmann scalar near the center, we obtain

$\displaystyle R_{\mu\nu\kappa\rho}R^{\mu\nu\kappa\rho}\simeq
16\pi^{2}G^{2}\alpha^{2}\big{(}53-192\ln(\mu t)+184\ln(\mu
t)^{2}\big{)}t^{-10}\,.$ (5.38)

Both of the above scalar invariants diverge for $t\rightarrow 0$ so the
perturbative quantum corrections do not resolve the singularity at the center
of the BH.

Before closing this section, we comments that the metric (5.31) is obtained in
the semi-classical approximation in which the matter field is treated quantum
mechanically while the geometry is still at the classical level. In this view,
this solution can not capture the whole quantum gravity effects. However, this
result may be useful towards the understanding of a full quantum gravity
solution. It is speculated that the final remnant of the BH after the
completion of Hawking radiation is a quantum primordial BH with the mass
around $M_{P}$. This can have interesting observational implications, for
example as the source of dark matter [55, 56]. In addition, as we discussed
previously, the backreactions by the quantum field modify the position of the
event horizon and the BH surface area. This in turn induces corrections in BH
entropy and its thermodynamics properties [57, 58, 59, 60].However, as we have
seen, the fractional correction is physical quantities such as the surface
gravity and the position of the event horizon scales like
$\alpha/G^{2}M^{3}\sim(M_{P}/M)^{2}$ which is extremely small for
astrophysical BHs. Therefore, we expect that the corrections in BH
thermodynamics are small and unobservable for astrophysical BHs.

##  6 Summary and Discussions

In this work we have studied the propagation of quantum field perturbations in
the interior of a Schwarzschild black hole. As reviewed in section 2.2, the
interior of the BH looks like an anisotropic cosmological background with the
structure of $R^{2}\times S^{2}$. As time moves forward, the extended
direction with the scale factor $a(\tau)$ is expanding while the compact
direction, with the scale factor $r(\tau)$, is shrinking. Eventually, we
encounter the big crunch singularity when
$r(\tau)\sim(-GM\tau)^{1/2}\rightarrow 0$.

In our treatment, the quantum perturbations are initially generated near the
horizon $\tau\rightarrow-\infty$ with the positive frequency solution
$e^{-i\omega\tau}$. In the interior region, this plane wave encounters the
effective potential $V_{\mathrm{eff}}$ which has a non-trivial time-dependent
profile. We are able to find the analytic solution for the high frequency mode
where the effects of the term $\ell(\ell+1)$ can be neglected. In this limit
the solution is given in terms of the Hankel functions $H_{\nu}^{(1)}$ and
$H^{(2)}_{\nu}$ of which only the former can match to the incoming positive
frequency mode near the horizon. With the profile of the mode function at
hand, we have calculated the vacuum expectation values of the components of
the anisotropic energy momentum tensor such as $\langle\rho\rangle$, $\langle
P\rangle$ and $\langle\Pi\rangle$. To regularize the UV divergences, we have
employed the $\zeta$ function regularization for the infinite sum over $\ell$
and the dimensional regularization scheme for the integral over ${\bf k}$.
Upon performing the regularizations, we have found that near the big crunch
singularity $\langle\rho\rangle_{\mathrm{reg}}\sim\langle
P\rangle_{\mathrm{reg}}\sim\langle\Pi\rangle_{\mathrm{reg}}\sim(GM)^{-3/2}(-\tau)^{\frac{-5}{2}}\ln(-\mu\tau)$.
The appearance of the logarithmic term $\ln(-\mu\tau)$ is reassuring,
indicating the quantum nature of the corrections.

Having obtained $\langle T^{\mu}_{\nu}\rangle_{\mathrm{reg}}$, we have
calculated the backreactions of the quantum perturbations on the background
metric. We have found the interesting result that $\delta A=-\delta B\sim
G^{2}M\ln(\mu t)/t^{3}$. To trust our perturbative treatment, we require that
the backreactions on the background solution to be small, requiring that
$\delta A\ll\bar{A}$. This in turn yields the time scale (or the radius of the
two-sphere) $t_{c}\sim\sqrt{G}$ beyond which the backreactions become strong.
It is possible that the backreactions become important sooner (i.e. at larger
radius) than what we obtained here since we have imposed various
simplification assumptions to solve the system analytically. Alternatively,
one may argue that the quantum corrections, like what are calculated here, may
play important roles to resolve the singularity. Whether or not these quantum
corrections can resolve the classical singularity require more sophisticated
QFT analysis with higher curvature terms included in the gravitational action
as well.

One intriguing implication of our result is that the corrections induced on
the metric by the quantum effects in the interior of the BH are, in principle,
accessible to the exterior observer. This is somewhat counterintuitive as it
is usually believed that events occurring in the interior of the BHs are
hidden to the exterior observer. However, this picture seems to be valid only
classically and in the presence of quantum effects, as in the case of Hawking
radiation, one may have to revise the classical interpretation about the roles
of the event horizon. Having said this, the fractional correction in the
metric and in the position of the horizon is at the order $(M_{P}/M)^{2}$
which is extremely small for astrophysical BHs.

There are a number of directions in which the current analysis can be
extended. In the current analysis we have considered the quantum perturbations
of a real scalar field in the interior of the BH. However, this analysis can
be extended to other fields such as the $U(1)$ gauge (vector) fields or tensor
perturbations as well. As in classical quasi-normal mode analysis, the
effective potential $V_{\mathrm{eff}}$ will now depend on the spin of the
field $s$. We would like to come back to the study of quantum perturbations
for the vector ( $s=1$) and tensor perturbations ($s=2$) in future. Finally,
we have calculated the vacuum expectation values with respect to the vacuum
$|0\rangle$ defined for the interior observer. On the other hand, the vacuum
associated to the Kruskal observer $|0\rangle_{K}$ (or the Hartle-Hawking
vacuum) has the advantage that it is well-defined for both interior and
exterior observers. These two vacua are related to each other via the
Bogolyubov transformations. It is therefore an interesting question to examine
what the Kruskal observer, equipped with the Hartle-Hawking vacuum, measures
for the expectation values.

Acknowledgments: We thank Mohammad Ali Gorji, Raihaneh Moti, Haidar
Sheikhahmadi and Alireza Talebian for useful discussions and correspondences.
This work is supported by the INSF of Iran under the grant number 4022911.

## References

  * [1] LIGO Scientific, Virgo collaboration, B. P. Abbott et al.__;_Phys. Rev. Lett._ 116 (2016) 061102, [1602.03837].
  * [2] LIGO Scientific, Virgo collaboration, B. P. Abbott et al.__;_Phys. Rev. Lett._ 116 (2016) 241103, [1606.04855].
  * [3] R. Kantowski and R. K. Sachs __;_J. Math. Phys._ 7 (1966) 443.
  * [4] A. D. Sakharov __;_Sov. Phys. JETP_ 22 (1966) 241.
  * [5] E. Poisson and W. Israel __;_Class. Quant. Grav._ 5 (1988) L201–L205.
  * [6] V. P. Frolov, M. A. Markov and V. F. Mukhanov __;_Phys. Lett. B_ 216 (1989) 272–276.
  * [7] V. P. Frolov, M. A. Markov and V. F. Mukhanov __;_Phys. Rev. D_ 41 (1990) 383.
  * [8] D. A. Easson and R. H. Brandenberger __;_JHEP_ 06 (2001) 024, [hep-th/0103019].
  * [9] H. Firouzjahi __;1610.03767.
  * [10] S. Wenda and S. Zhu __;_Physics Letters A_ 126 (1988) 229–232.
  * [11] D. Yoshida and R. H. Brandenberger __;_JCAP_ 07 (2018) 022, [1801.05070].
  * [12] R. Brandenberger, L. Heisenberg and J. Robnik __;_JHEP_ 05 (2021) 090, [2103.02842].
  * [13] M. de Cesare, S. S. Seahra and E. Wilson-Ewing __;_JCAP_ 07 (2020) 018, [2002.11658].
  * [14] E. Gaztanaga __;_Hal Open Science_ 7 (2021) 03783110.
  * [15] E. Gaztanaga __;_Universe_ 8 (2022) 257, [2204.11608].
  * [16] P. O. Mazur and E. Mottola __;_Universe_ 9 (2023) 88, [gr-qc/0109035].
  * [17] M. Visser and D. L. Wiltshire __;_Class. Quant. Grav._ 21 (2004) 1135–1152, [gr-qc/0310107].
  * [18] C. Cattoen, T. Faber and M. Visser __;_Class. Quant. Grav._ 22 (2005) 4189–4202, [gr-qc/0505137].
  * [19] M. Markov __;_JETP Lett.(Engl. Transl.);(United States)_ 36 (1982) .
  * [20] S. W. Hawking __;_Commun. Math. Phys._ 43 (1975) 199–220.
  * [21] W. G. Unruh __;_Phys. Rev. D_ 14 (1976) 870.
  * [22] V. Mukhanov and S. Winitzki, _Introduction to quantum effects in gravity_.  Cambridge University Press, 6, 2007.
  * [23] P. K. Townsend __;gr-qc/9707012.
  * [24] A. Fabbri and J. Navarro-Salas, _Modeling black hole evaporation_.  2005\.
  * [25] M. Visser __;_Int. J. Mod. Phys. D_ 12 (2003) 649–661, [hep-th/0106111].
  * [26] T. Jacobson __; in _School on Quantum Gravity_ , pp. 39–89, 8, 2003, gr-qc/0308048, DOI.
  * [27] H. Firouzjahi and A. Talebian __;_Phys. Rev. D_ 106 (2022) 123505, [2210.15186].
  * [28] R. Doran, F. S. N. Lobo and P. Crawford __;_Found. Phys._ 38 (2008) 160–187, [gr-qc/0609042].
  * [29] P. Gusin, A. Radosz, A. T. Augousti, J. Polonyi, O. B. Zaslavskii and R. J. Sciborski __;_Universe_ 9 (2023) 299, [2306.15549].
  * [30] G. ’t Hooft and M. J. G. Veltman __;_Nucl. Phys. B_ 44 (1972) 189–213.
  * [31] G. ’t Hooft and M. J. G. Veltman __;_Ann. Inst. H. Poincare Phys. Theor. A_ 20 (1974) 69–94.
  * [32] C. G. Bollini and J. J. Giambiagi __;_Nuovo Cim. B_ 12 (1972) 20–26.
  * [33] S. Deser and P. van Nieuwenhuizen __;_Phys. Rev. D_ 10 (1974) 401.
  * [34] J. S. Dowker and R. Critchley __;_Phys. Rev. D_ 13 (1976) 3224.
  * [35] A. O. Barvinsky and G. A. Vilkovisky __;_Phys. Rept._ 119 (1985) 1–74.
  * [36] T. Regge and J. A. Wheeler __;_Phys. Rev._ 108 (1957) 1063–1069.
  * [37] H.-P. Nollert __;_Class. Quant. Grav._ 16 (1999) R159–R216.
  * [38] K. D. Kokkotas and B. G. Schmidt __;_Living Rev. Rel._ 2 (1999) 2, [gr-qc/9909058].
  * [39] H. Firouzjahi and H. Sheikhahmadi __;_Phys. Rev. D_ 108 (2023) 065002, [2307.00977].
  * [40] H. Firouzjahi and H. Sheikhahmadi __;2403.10463.
  * [41] S. W. Hawking __;_Commun. Math. Phys._ 55 (1977) 133.
  * [42] E. Elizalde, S. D. Odintsov, A. Romeo, A. A. Bytsenko and S. Zerbini, _Zeta regularization techniques with applications_.  World Scientific Publishing, Singapore, 1994, 10.1142/2065.
  * [43] E. Elizalde __;_J. Math. Phys._ 35 (1994) 3308–3321, [hep-th/9308048].
  * [44] G. F. R. Ellis and H. van Elst __;_NATO Sci. Ser. C_ 541 (1999) 1–116, [gr-qc/9812046].
  * [45] H. Firouzjahi __;_Phys. Rev. D_ 106 (2022) 045015, [2205.06561].
  * [46] S. Weinberg, _Gravitation and Cosmology: Principles and Applications of the General Theory of Relativity_.  John Wiley and Sons, New York, 1972.
  * [47] M. A. Gorji, A. Allahyari, M. Khodadi and H. Firouzjahi __;_Phys. Rev. D_ 101 (2020) 124060, [1912.04636].
  * [48] M. J. Duff __;_Phys. Rev. D_ 9 (1974) 1837–1839.
  * [49] N. E. J. Bjerrum-Bohr, J. F. Donoghue and B. R. Holstein __;_Phys. Rev. D_ 68 (2003) 084005, [hep-th/0211071].
  * [50] J. F. Donoghue, B. R. Holstein, B. Garbrecht and T. Konstandin __;_Phys. Lett. B_ 529 (2002) 132–142, [hep-th/0112237].
  * [51] X. Calmet and B. K. El-Menoufi __;_Eur. Phys. J. C_ 77 (2017) 243, [1704.00261].
  * [52] S. M. Carroll, _Spacetime and Geometry: An Introduction to General Relativity_.  Cambridge University Press, 7, 2019, 10.1017/9781108770385.
  * [53] M. Visser __;_Phys. Rev. D_ 46 (1992) 2445–2451, [hep-th/9203057].
  * [54] E. Poisson, _A Relativist’s Toolkit: The Mathematics of Black-Hole Mechanics_.  Cambridge University Press, 12, 2009, 10.1017/CBO9780511606601.
  * [55] J. H. MacGibbon __;_Nature_ 329 (1987) 308–309.
  * [56] G. Domènech and M. Sasaki __;_Class. Quant. Grav._ 40 (2023) 177001, [2303.07661].
  * [57] A. Bonanno and M. Reuter __;_Phys. Rev. D_ 62 (2000) 043008, [hep-th/0002196].
  * [58] Z. W. Feng, H. L. Li, X. T. Zu and S. Z. Yang __;_Eur. Phys. J. C_ 76 (2016) 212, [1604.04702].
  * [59] Z.-W. Feng, S.-Z. Yang, H.-L. Li and X.-T. Zu __;_Adv. High Energy Phys._ 2016 (2016) 2341879, [1607.04114].
  * [60] Y. Xiao and Y. Tian __;_Phys. Rev. D_ 105 (2022) 044013, [2104.14902].


(original image) uses a spatial correlation to detect sources. The top row
DAWIS fractions are therefore better (less contaminated, see
$e^{\rm{cont}}_{\rm{ICL}}$ in Figure 8) than the other methods, showcasing its
capacity to detect faint and extended features.

##  6 Example on real astronomical data

Figure 11: From left to right: original astronomical image, complete restored
image, residuals and restored ICL profile. The images are $1.6\times 1.6$
arcmin2. The green circle has a radius of 50 kpc, and the cyan circle has a
radius of 100 kpc.

In this section we demonstrate the capabilities of DAWIS on a galaxy cluster
image at optical wavelengths. The targeted galaxy cluster is LCDCS 0541, which
is part of both the Las Campanas Distant Cluster Survey (LCDCS; Gonzalez et
al., 2001) and the Dark energy American French Team / French American DArk
energy Team (DAFT/FADA; Guennou et al., 2010) surveys. The redshift of this
cluster is 0.542, and a deep 2x2 arcmin2 HST Advanced Camera for Survey (ACS)
image in the F814W filter is available. The photometric depth of the image was
estimated (Guennou et al., 2012, hereafter G12) to be around
$\mu_{\rm{F814W}}\sim$ 24 mag.arcsec-2 for extended sources with SExtractor
($\sim$ 1.25 arcsec, $1.8\sigma$ detection limit). Details about the data and
the reduction procedure are available in G12.

We chose LCDCS 0541 and these data because this target was the case study of
G12, where a multiscale analysis was performed to detect ICL in the images of
ten galaxy clusters using OV_WAV, better suited for this task than SExtractor.
By running DAWIS on the same exact cluster image, we can compare the results
of DAWIS to those obtained by this earlier multiscale analysis method (for
more details on the differences between both approaches, see Section 2 for
traditional multiscale analysis methods, and Section 3 for an in-depth
description of DAWIS).

In G12, OV_WAV was run two consecutive times on the LCDCS 0541 image. The
first run allowed the modeling of sources with characteristic sizes up to
$j=10$ (1024 pixels), and the second run up to $j=7$ (128 pixels). The
detected and restored sources were removed from the original scene, and the
ICL was detected by applying a 2.5$\sigma$ detection limit on the residuals.
The background level for the detection limit was computed in a wide annulus in
the outskirts of the residuals. The total flux in the detected ICL source was
measured, and the associated absolute magnitude was computed assuming a
$k$-correction of an early-type galaxy: this gave an absolute magnitude of
$-20.4$. Note that with this procedure, the multiscale analysis made by the
authors is analogous to a PF detection method, as the ICL sources are detected
by applying a detection limit in the direct space residuals after fitting and
removing all bright and small sources (MW stars and galaxies).

We run DAWIS on the same LCDCS 0541 F814W image, and compare both the absolute
magnitude value and morphological properties of the resulting ICL source. In
our case, we resize the image to a FoV of $2048\times 2048$ pixel size,
corresponding to $\sim 1.6\times 1.6$ arcmin2. The CLEAN factor is set to
$\delta=0.5$, and the relative threshold parameter is set to $\tau=0.8$. As we
are only interested here in detecting the diffuse ICL source in the central
part of the image, we do not apply any particular criterion for the
classification step, except for a separation based on the wavelet scale of the
interscale maximum $j_{\rm{sep}}$, and for the fact that the restored sources
composing the atoms associated to ICL must be centered on the BCG position. We
set $j_{\rm{set}}=8$ here, corresponding to a characteristic size of 256
pixels ($\sim$ 80 kpc at $z=0.542$).

The resulting restored astronomical field, residual image and restored ICL
distribution are displayed in Figure 11. The restored ICL distribution is
anisotropic, with a boxy shape of $\sim 200\times 100$ kpc2 and has an average
apparent surface brightness of $\sim$ 26.3 mag.arcsec-2 (corresponding to a
mean surface brightness of $\sim$ 24.5 mag.arcsec-2 when correcting for the
cosmological dimming). The morphology of our restored ICL distribution agrees
with the apparent morphology of the detection made by G12 (see their Figure
4), but ours is continuous and seems slightly more extended. We compute the
associated absolute magnitude using the same $k$-correction as G12, and find a
value of $-22.4$. This value is two orders of magnitude brighter than the
source detected in G12. We adopt the same cosmology than G12 to compute these
numbers. This detection shows the capabilities of DAWIS to detect these faint
and extended sources where other multiscale methods could not before, and this
in an automated way.

We now discuss the implications of this comparison. We detect with DAWIS more
flux than G12, resulting in a brighter absolute magnitude. This can be
attributed to two things: DAWIS detecting more flux due to the detection of
the ICL component directly in the wavelet space, and differences in
classification. First, it is expected that DAWIS, due to the synthesis-by-
analysis and the meticulous removal of all bright sources before detecting
faint ones, detects more flux in faint and extended sources than what is
detected with a simple detection limit in the direct space residuals. Second,
such a careful approach allows to clearly assess and classify each restored
source. In the case of G12, all the sources in an OV_WAV pass are restored at
once, meaning that no classification is done and that signal from ICL may have
been removed by this approach (especially in the first pass, that allows the
removal of objects with characteristic sizes up to 1024 pixels). This also
explains such a discrepancy in flux and absolute magnitude.

Of course, the question of the classification criteria, and more specifically
of the separation wavelet scale $j_{\rm{sep}}$ and its value must be brought
up again. We chose $j_{\rm{sep}}=8$ because a lower value would imply a
significant increase in the flux associated to ICL and of its average apparent
surface brightness (jumping to $\mu\sim 25.7$ mag.arcsec-2). The source we
detected has already a bright average surface brightness, compared to results
from literature. More importantly, the bright inner profile of the BCG of
LCDCS 0541 has a spatial extent of $\sim$ 40 to 50 kpc. Associating to ICL
sources with characteristic sizes of at least twice this spatial extent seems
to be an accurate choice. On the other hand, a value of $j_{\rm{sep}}=7$
corresponding to a typical size of $\sim$ 40 kpc sets a precarious separation.
In any case, the question of whether this extended LSB source should actually
be called ICL or is rather the outskirts of the BCG intensity profile is still
pending, and cannot really be answered with certainty with only monochromatic
data.

##  7 Conclusions

We presented in Section 3 our novel detection algorithm DAWIS, based on a
synthesis-by-analysis approach with an operating mode based on an isotropic
wavelet dictionary and on interscale connectivity analysis. The algorithm is
developed keeping in mind the limitations and defects of previous wavelet
algorithm packages, and therefore integrates a semi-greedy structure with
parallelized modules. For an input image, DAWIS computes in this iterative
procedure the synthesised dictionary composed of atoms (e.g. the sum of
restored source profiles) and the associated denoised restored image of the
original field. We showed that there is a great range of possibilities to
separate objects and compose different synthesised dictionaries based on
astrophysical priors and on the goals of the performed analysis.

To estimate the efficiency of the algorithm, we simulated in Section 4 mock
images of galaxy clusters with noise using the modeling package Galsim.
Galaxies are represented by Sérsic profiles, and the ICL by a large and faint
exponential profile, the properties of which are controlled by few parameters,
notably an input ICL fraction. The images are only emulating the MegaCam
instrument properties, but the same study with simulations of other
instruments should be led in the future. We take advantage of this work to
test in Section 5 the efficiency of other detection methods that we roughly
classified in two categories: surface brightness threshold and fitting
algorithms. We apply these two methods in the same way as DAWIS to the
simulations, and compare the resulting measured fluxes and ICL fractions.

We showed in Section 5 the limitation of the three detection methods,
especially the SBT method which applies a constant threshold in surface
brightness to the image to separate ICL from galaxies. As already shown in
Rudick et al. (2011), galaxy light profiles and ICL are overlapping and
superimposed, and the use of a single threshold to separate these two
components leads to strong contamination of the detected ICL (see Section
5.5). While we managed to list at least part of the contamination and
degradation effects (masks, contamination by galaxy outskirt profiles,
contamination by faint undetected galaxies, loss of ICL flux under the
$3\sigma$ detection limit), these effects are co-added, and future studies
should be led to isolate and constrain them in a more refined way.

The PF method, acting as an upper limit to the fitting method quality, allowed
us to showcase in Section 5.4 the amount of flux lost under the detection
limit when applying a standard detection procedure in the mock images. The
differences on the measures of ICL fractions for faint and extended ICL
profiles between this method and DAWIS come mainly from the fact that DAWIS
detects sources in the wavelet scales where spatial correlation information is
used. This gives food for thoughts on the importance of the methodology used
to estimate the limiting depth and detection level in images, a topic already
addressed by Mihos (2019) and not specific to ICL, but rather touching all low
surface brightness observations. We advocate here that in order to better
constrain ICL detection, future studies should concentrate (in parallel with
methods to separate galaxies from ICL) on the effects of sky background noise
estimation in images, specially in more complex and realistic cases than the
flat background of our images.

We showed that in most cases, DAWIS does a good job to measure ICL fractions
as little contaminated as possible (see Section 5.5). However, a limitation of
the priors we use to separate galaxies from ICL (based on the scale of the
interscale maxima as detailed in Section 5.3) becomes apparent when DAWIS
tries to process images with bright and concentrated ICL halos; part of the
synthesis atoms that should be associated to ICL are in this case wrongly
classified and associated to galaxies. The classification operator used in
this work is simplistic in nature and based mainly on the characteristic size
and spatial information of the processed atoms. This fact is accentuated by
our results on LCDCS 0541, for which we discussed the relevance of such a
prior to separate ICL from the BCG intensity distribution in real astronomical
data. The large discrepancy in absolute magnitude between our detection and
the detection made by G12 could be in part explained by the different
classifications made in both studies. In future studies, we plan to implement
a more complete and robust classification for the ICL, using additional
criteria such as morphology or granularity for example.

As all methods show flaws in the detection of ICL, we want to emphasize here
the notion of bias compensation. Some ICL fractions measured with the PF or
the SBT method appear accurate (see Section 5.5) not because the measured ICL
and galaxy fluxes are correct, but because different contamination and
degradation effects are compensating each other: because ICL fractions are
fractions of integrated values, an effect of bias compensation is occurring.
This effect of bias compensation cannot really be seen or estimated in real
astronomical images, which is problematic. This would definitely have
worrisome effects in multi-band studies (ICL spectral energy distribution for
example), where evolution of ICL fractions with wavelength could be
misinterpreted. Our conclusion here is that the ICL fraction taken
individually, is not a good metric to characterise ICL. We advocate for the
use of physical flux values in future studies for both galaxies and ICL, in
order to better put in evidence the results of different approaches to detect
ICL.

Our simulations in this work do not integrate several other components or
contamination effects in the detection of LSB sources: large-scale background
spatial variation, bright foreground stars, scattered light halos, presence of
tidal streams just to cite a few. While the isotropic wavelet dictionary is a
good tool to detect and restore quasi-circular sources such as galaxies, these
other features are not particularly well expressed in the associated
representation space. However, the core operating mode of DAWIS is versatile
enough to allow future upgrades for these components, not only on the atom
classification step to limit contamination effects, but also in the analysis
dictionaries used for the detection or the restoration of sources. Including
analyses based on other representation spaces and function bases adequate for
these signal shapes is also a good direction to take when trying to capture
the low surface brightness content of astronomical images.

###### Acknowledgements.

We acknowledge long-term support from CNES. N. Martinet acknowledges support
from a CNES fellowship. This work has made use of the CANDIDE Cluster at the
Institut d’Astrophysique de Paris and was made possible by grants from the
PNCG and the DIM-ACAV.

## References

  * Abraham & van Dokkum (2014) Abraham, R. G. & van Dokkum, P. G. 2014, PASP, 126, 55
  * Adami et al. (2013) Adami, C., Durret, F., Guennou, L., & Da Rocha, C. 2013, A&A, 551, A20
  * Adami et al. (2005) Adami, C., Slezak, E., Durret, F., et al. 2005, A&A, 429, 39
  * Anscombe (1948) Anscombe, F. J. 1948, Biometrika, 35, 246
  * Atkinson et al. (2013) Atkinson, A. M., Abraham, R. G., & Ferguson, A. M. N. 2013, ApJ, 765, 28
  * Beckwith et al. (2006) Beckwith, S. V. W., Stiavelli, M., Koekemoer, A. M., et al. 2006, AJ, 132, 1729
  * Bertin & Arnouts (1996) Bertin, E. & Arnouts, S. 1996, A&AS, 117, 393
  * Bijaoui & Rué (1995) Bijaoui, A. & Rué, F. 1995, Signal Processing, 46, 345
  * Burke et al. (2012) Burke, C., Collins, C. A., Stott, J. P., & Hilton, M. 2012, MNRAS, 425, 2058
  * Burke et al. (2015) Burke, C., Hilton, M., & Collins, C. 2015, Monthly Notices of the Royal Astronomical Society, 449, 2353
  * Candes et al. (2006) Candes, E. J., Romberg, J., & Tao, T. 2006, IEEE Transactions on Information Theory, 52, 489
  * Chen et al. (2001) Chen, S. S., Donoho, D. L., & Saunders, M. A. 2001, SIAM Rev., 43, 129
  * Da Rocha & Mendes de Oliveira (2005) Da Rocha, C. & Mendes de Oliveira, C. 2005, MNRAS, 364, 1069
  * Da Rocha et al. (2008) Da Rocha, C., Ziegler, B. L., & Mendes de Oliveira, C. 2008, MNRAS, 388, 1433
  * Dabbech et al. (2015) Dabbech, A., Ferrari, C., Mary, D., et al. 2015, A&A, 576, A7
  * Daubechies (1988) Daubechies, I. 1988, IEEE Transactions on Information Theory, 34, 605
  * Daubechies (1990) Daubechies, I. 1990, IEEE Transactions on Information Theory, 36, 961
  * Daubechies (1992) Daubechies, I. 1992, in CBMS-NSF regional conference series in applied mathematics
  * DeMaio et al. (2018) DeMaio, T., Gonzalez, A. H., Zabludoff, A., et al. 2018, MNRAS, 474, 3009
  * Donoho (2006) Donoho, D. L. 2006, IEEE Transactions on Information Theory, 52, 1289
  * Donzelli et al. (2011) Donzelli, C. J., Muriel, H., & Madrid, J. P. 2011, ApJS, 195, 15
  * Duc et al. (2011) Duc, P.-A., Cuillandre, J.-C., Serra, P., et al. 2011, MNRAS, 417, 863
  * Durret et al. (2019) Durret, F., Tarricq, Y., Márquez, I., Ashkar, H., & Adami, C. 2019, A&A, 622, A78
  * Elad et al. (2007) Elad, M., Milanfar, P., & Rubinstein, R. 2007, Inverse Problems, 23, 947
  * Ellien et al. (2019) Ellien, A., Durret, F., Adami, C., et al. 2019, A&A, 628, A34
  * Euclid Collaboration et al. (2019) Euclid Collaboration, Martinet, N., Schrabback, T., et al. 2019, A&A, 627, A59
  * Feldmeier et al. (2004) Feldmeier, J. J., Mihos, J. C., Morrison, H. L., et al. 2004, ApJ, 609, 617
  * Feldmeier et al. (2002) Feldmeier, J. J., Mihos, J. C., Morrison, H. L., Rodney, S. A., & Harding, P. 2002, ApJ, 575, 779
  * Ferrarese et al. (2012) Ferrarese, L., Côté, P., Cuilland re, J.-C., et al. 2012, ApJS, 200, 4
  * Fletcher & Reeves (1964) Fletcher, R. & Reeves, C. M. 1964, The Computer Journal, 7, 149
  * Gonzalez et al. (2005) Gonzalez, A. H., Zabludoff, A. I., & Zaritsky, D. 2005, ApJ, 618, 195
  * Gonzalez et al. (2001) Gonzalez, A. H., Zaritsky, D., Dalcanton, J. J., & Nelson, A. 2001, ApJS, 137, 117
  * Guennou et al. (2012) Guennou, L., Adami, C., Da Rocha, C., et al. 2012, A&A, 537, A64
  * Guennou et al. (2010) Guennou, L., Adami, C., Ulmer, M. P., et al. 2010, A&A, 523, A21
  * Gwyn (2012) Gwyn, S. D. J. 2012, AJ, 143, 38
  * Holmberg (1958) Holmberg, E. 1958, Meddelanden fran Lunds Astronomiska Observatorium Serie II, 136, 1
  * Holschneider et al. (1989) Holschneider, M., Kronland-Martinet, R., Morlet, J., & Tchamitchian, P. 1989, in Wavelets. Time-Frequency Methods and Phase Space, ed. J.-M. Combes, A. Grossmann, & P. Tchamitchian, 286
  * Ibata et al. (2017) Ibata, R. A., McConnachie, A., Cuilland re, J.-C., et al. 2017, ApJ, 848, 128
  * Ivezić et al. (2019) Ivezić, Ž., Kahn, S. M., Tyson, J. A., et al. 2019, ApJ, 873, 111
  * Jiménez-Teja & Benítez (2012) Jiménez-Teja, Y. & Benítez, N. 2012, ApJ, 745, 150
  * Jiménez-Teja & Dupke (2016) Jiménez-Teja, Y. & Dupke, R. 2016, ApJ, 820, 49
  * Jiménez-Teja et al. (2018) Jiménez-Teja, Y., Dupke, R., Benítez, N., et al. 2018, ApJ, 857, 79
  * Kluge et al. (2020) Kluge, M., Neureiter, B., Riffeser, A., et al. 2020, The Astrophysical Journal Supplement Series, 247, 43
  * Klypin et al. (2016) Klypin, A., Yepes, G., Gottlöber, S., Prada, F., & Heß, S. 2016, MNRAS, 457, 4340
  * Ko & Jee (2018) Ko, J. & Jee, M. J. 2018, ApJ, 862, 95
  * Krick & Bernstein (2007) Krick, J. E. & Bernstein, R. A. 2007, AJ, 134, 466
  * Krick et al. (2006) Krick, J. E., Bernstein, R. A., & Pimbblet, K. A. 2006, AJ, 131, 168
  * Laigle et al. (2016) Laigle, C., McCracken, H. J., Ilbert, O., et al. 2016, ApJS, 224, 24
  * Lotz et al. (2017) Lotz, J. M., Koekemoer, A., Coe, D., et al. 2017, ApJ, 837, 97
  * Mallat (1989) Mallat, S. 1989, IEEE Trans. Pattern Anal. Mach. Intell., 11, 674
  * Mallat (2008) Mallat, S. 2008, Academic, Burlington, Mass
  * Mallat & Zhang (1993) Mallat, S. G. & Zhang. 1993, IEEE Transactions on Signal Processing, 41, 3397
  * Mihos (2019) Mihos, J. C. 2019, arXiv e-prints, arXiv:1909.09456
  * Mihos et al. (2017) Mihos, J. C., Harding, P., Feldmeier, J. J., et al. 2017, ApJ, 834, 16
  * Montes (2019) Montes, M. 2019, arXiv e-prints, arXiv:1912.01616
  * Montes & Trujillo (2018) Montes, M. & Trujillo, I. 2018, MNRAS, 474, 917
  * Montes & Trujillo (2019) Montes, M. & Trujillo, I. 2019, MNRAS, 482, 2838
  * Morishita et al. (2017) Morishita, T., Abramson, L. E., Treu, T., et al. 2017, ApJ, 846, 139
  * Navarro et al. (1997) Navarro, J. F., Frenk, C. S., & White, S. D. M. 1997, ApJ, 490, 493
  * Peng et al. (2002) Peng, C. Y., Ho, L. C., Impey, C. D., & Rix, H.-W. 2002, AJ, 124, 266
  * Puchwein et al. (2010) Puchwein, E., Springel, V., Sijacki, D., & Dolag, K. 2010, MNRAS, 406, 936
  * Refregier (2003) Refregier, A. 2003, MNRAS, 338, 35
  * Remus et al. (2017) Remus, R.-S., Dolag, K., Naab, T., et al. 2017, MNRAS, 464, 3742
  * Rowe et al. (2015) Rowe, B. T. P., Jarvis, M., Mandelbaum, R., et al. 2015, Astronomy and Computing, 10, 121
  * Rudick et al. (2009) Rudick, C. S., Mihos, J. C., Frey, L. H., & McBride, C. K. 2009, ApJ, 699, 1518
  * Rudick et al. (2011) Rudick, C. S., Mihos, J. C., & McBride, C. K. 2011, ApJ, 732, 48
  * Scheick & Kuhn (1994) Scheick, X. & Kuhn, J. R. 1994, ApJ, 423, 566
  * Seigar et al. (2007) Seigar, M. S., Graham, A. W., & Jerjen, H. 2007, MNRAS, 378, 1575
  * Shensa (1992) Shensa, M. J. 1992, IEEE Trans. Signal Processing, 40, 2464
  * Slezak et al. (1990) Slezak, E., Bijaoui, A., & Mars, G. 1990, A&A, 227, 301
  * Slezak et al. (1994) Slezak, E., Durret, F., & Gerbal, D. 1994, AJ, 108, 1996
  * Starck et al. (2007) Starck, J.-L., M. Fadili, J., & Murtagh, F. 2007, IEEE transactions on image processing : a publication of the IEEE Signal Processing Society, 16, 297
  * Starck et al. (1998) Starck, J.-L., Murtagh, F. D., & Bijaoui, A. 1998, Image Processing and Data Analysis, 297
  * Steinhardt et al. (2020) Steinhardt, C. L., Jauzac, M., Acebron, A., et al. 2020, arXiv e-prints, arXiv:2001.09999
  * Tang et al. (2018) Tang, L., Lin, W., Cui, W., et al. 2018, ApJ, 859, 85
  * Uson et al. (1991) Uson, J. M., Boughn, S. P., & Kuhn, J. R. 1991, ApJ, 369, 46
  * Valls-Gabaud & MESSIER Collaboration (2017) Valls-Gabaud, D. & MESSIER Collaboration. 2017, in IAU Symposium, Vol. 321, Formation and Evolution of Galaxy Outskirts, ed. A. Gil de Paz, J. H. Knapen, & J. C. Lee, 199–201
  * Zibetti et al. (2005) Zibetti, S., White, S. D. M., Schneider, D. P., & Brinkmann, J. 2005, MNRAS, 358, 949

##  Appendix A Varying SBT

Different values have been tested for the SBT method with
$25\leq\mu_{\rm{SBT}}\leq 28$. We applied exactly the same methodology as in
Section 5.1, and measured the ICL fractions in the GAL+ICL+NOISE images. The
results are shown in Figure 12, with the same format as Figure 10. As expected
and consistent with the literature (Rudick et al. 2009; Tang et al. 2018), the
choice of SBT influences greatly the measured ICL fraction, with values
ranging from $\sim$ 0% to $\sim$ 80% in some cases for the same input ICL
fraction. Most of the times however, thresholds $\mu_{\rm{SBT}}=26.5$
mag.arcsec-2 or $\mu_{\rm{SBT}}=27$ mag.arcsec-2, which are the generic SBTs
considered in the literature, give results that are more consistent with the
input ICL fraction. However, this trend needs to be contrasted with our in
depth analysis of biases in the measure of $F_{\rm{cluster}}$ and
$F_{\rm{ICL}}$ (see Section 5.4 and Section 5.5).

Figure 12: ICL fractions displayed for the different ICL profile parameters
(half-light radius $r_{e}$, input ICL fraction $f^{\rm{input}}_{\rm{ICL}}$ and
redshift $z$) as measured with different SBT. The circles with continuous
lines are the ICL fractions $f_{\rm{ICL}}$computed from each of the 270
GAL+ICL+NOISE MegaCam-type images, and the displayed values correspond to the
average of the ten images with the same input parameters. The error bars
correspond to standard deviations across ten different clusters. The input ICL
fraction is also displayed as a black horizontal dashed line for better
visualization.