# Statistics of boundary encounters by a particle diffusing outside a compact
planar domain

Denis S. Grebenkov<EMAIL_ADDRESS>Laboratoire de Physique
de la Matière Condensée (UMR 7643),
CNRS – Ecole Polytechnique, IP Paris, 91128 Palaiseau, France

###### Abstract

We consider a particle diffusing outside a compact planar set and investigate
its boundary local time $\ell_{t}$, i.e., the rescaled number of encounters
between the particle and the boundary up to time $t$. In the case of a disk,
this is also the (rescaled) number of encounters of two diffusing circular
particles in the plane. For that case, we derive explicit integral
representations for the probability density of the boundary local time
$\ell_{t}$ and for the probability density of the first-crossing time of a
given threshold by $\ell_{t}$. The latter density is shown to exhibit a very
slow long-time decay due to extremely long diffusive excursions between
encounters. We briefly discuss some practical consequences of this behavior
for applications in chemical physics and biology.

###### pacs:

02.50.-r, 05.40.-a, 02.70.Rr, 05.10.Gg

††: J. Phys. A: Math. Gen.

Keywords: Boundary local time; Reflected Brownian motion; Diffusion-influenced
reactions; Surface reactivity; Robin boundary condition; Heterogeneous
catalysis

##  1 Introduction

When a Brownian particle diffuses in a geometric confinement, its encounters
with the reflecting boundary can be characterized by the boundary local time
$\ell_{t}$, which plays the central role in the theory of stochastic processes
[1, 2, 3]. In the basic setting of ordinary diffusion, one considers reflected
Brownian motion $\bm{X}_{t}$, released at time $t=0$ from a fixed point
$\bm{x}_{0}$ and diffusing inside an Euclidean domain
$\Omega\subset{\mathbb{R}}^{d}$ with diffusion coefficient $D$ and normal
reflections on the smooth boundary $\partial\Omega$. The boundary local time
$\ell_{t}$ of this process on a subset of the boundary,
$\Gamma\subset\partial\Omega$, is defined as

$\ell_{t}=\lim\limits_{a\to
0}\frac{D}{a}\underbrace{\int\limits_{0}^{t}dt^{\prime}\,\Theta(a-|\bm{X}_{t^{\prime}}-\Gamma|)}_{\textrm{residence
time}},$ (1)

where $\Theta(z)$ is the Heaviside step function, and $|\bm{x}-\Gamma|$
denotes the distance between a point $\bm{x}$ and the boundary region
$\Gamma$. In this definition, the integral is the residence time of reflected
Brownian motion $\bm{X}_{t}$ up to time $t$ inside a thin boundary layer of
width $a$ around $\Gamma$:
$\Gamma_{a}=\\{\bm{x}\in\Omega~{}:~{}|\bm{x}-\Gamma|<a\\}$ (see [4, 5, 6, 7,
8, 9, 10, 11, 12, 13, 14, 15, 16, 17] and references therein). In the limit
$a\to 0$, the residence time vanishes but its rescaling by $a$ yields a well-
defined nontrivial limit $\ell_{t}$. According to Eq. (1), the boundary local
time $\ell_{t}$ is a non-decreasing process that remains constant when
$\bm{X}_{t}$ is the bulk, and increases only when $\bm{X}_{t}$ hits the
boundary. Alternatively, the boundary local time $\ell_{t}$ can also be
written as

$\ell_{t}=\lim\limits_{a\to 0}a\,\mathcal{N}_{t}^{a},$ (2)

where $\mathcal{N}_{t}^{a}$ is the number of downcrossings of the thin
boundary layer $\Gamma_{a}$ up to time $t$ (Fig. 1), i.e., a regularized
version of the number of encounters of the process with the region $\Gamma$
(see [18] for further discussion). For a fixed time $t$, $\ell_{t}$ is a
random variable, which can be characterized by the probability density
$\rho(\ell,t|\bm{x}_{0})$. Note that $\ell_{t}$ has units of length, whereas
$\ell_{t}/D$ has units of time per length, reflecting the rescaling by $a$ in
Eq. (1).

Figure 1:  A simulated trajectory of reflected Brownian motion in the exterior
of a disk of radius $R$:
$\Omega=\\{\bm{x}\in{\mathbb{R}}^{2}~{}:~{}|\bm{x}|>R\\}$. Pink region denotes
a thin boundary layer $\Gamma_{a}$ of width $a$ near the boundary:
$\Gamma=\partial\Omega$. The particle is released in a close vicinity of the
boundary (blue diamond) and diffuses until the first encounter with that
boundary (filled circle enumerated by $1$). From that point, the particle is
released a distance $a$ above the circle and resumes its motion until the next
encounter, and so on. Such a regularization with small $a>0$ allows one to
define the number of encounters $\mathcal{N}_{t}^{a}$ and to split the
trajectory into excursions between encounters (such three excursions are drawn
by different colors).

While general stochastic properties of the boundary local time were thoroughly
investigated in the past (see [20, 21, 22] and references therein), its
geometry-specific properties are less known. For instance, how does the
distribution of the number of encounters depend on the shape of the confining
domain and evolve with time? When does this number exceed a prescribed
threshold? To answer these questions, general spectral expansions for
$\rho(\ell,t|\bm{x}_{0})$ and for the probability density
$U(\ell,t|\bm{x}_{0})$ of the first-crossing time ${\mathcal{T}}_{\ell}$ of a
prescribed threshold $\ell$ by the process $\ell_{t}$, were established in
[23, 18]. These spectral expansions rely on the eigenmodes of the Dirichlet-
to-Neumann operator, which are known explicitly only for some simple domains
[24]. In particular, closed analytical formulas were derived for diffusion in
a half-line,

$\displaystyle\rho(\ell,t|x_{0})$ $\displaystyle=$
$\displaystyle\mathrm{erf}\left(\frac{x_{0}}{\sqrt{4Dt}}\right)\delta(\ell)+\frac{\exp\bigl{(}-\frac{(x_{0}+\ell)^{2}}{4Dt}\bigr{)}}{\sqrt{\pi
Dt}}\,,$ (3) $\displaystyle U(\ell,t|x_{0})$ $\displaystyle=$
$\displaystyle(\ell+x_{0})\frac{e^{-(\ell+x_{0})^{2}/(4Dt)}}{\sqrt{4\pi
Dt^{3}}},$ (4)

and in the exterior of a ball of radius $R$:

$\displaystyle\rho(\ell,t|\bm{x}_{0})$ $\displaystyle=$
$\displaystyle\left(1-\frac{R}{|\bm{x}_{0}|}\mathrm{erfc}\left(\frac{|\bm{x}_{0}|-R}{\sqrt{4Dt}}\right)\right)\delta(\ell)$
$\displaystyle+$
$\displaystyle\frac{e^{-\ell/R}}{|\bm{x}_{0}|}\left(\mathrm{erf}\left(\frac{|\bm{x}_{0}|-R+\ell}{\sqrt{4Dt}}\right)+\frac{R\,e^{-(|\bm{x}_{0}|-R+\ell)^{2}/(4Dt)}}{\sqrt{\pi
Dt}}\right)\,,$ $\displaystyle U(\ell,t|\bm{x}_{0})$ $\displaystyle=$
$\displaystyle\frac{R\,e^{-\ell/R}}{|\bm{x}_{0}|}\,\frac{|\bm{x}_{0}|-R+\ell}{\sqrt{4\pi
Dt^{3}}}e^{-(|\bm{x}_{0}|-R+\ell)^{2}/(4Dt)}\,,$ (6)

where $\delta(\ell)$ is the Dirac distribution (see [25, 24] for details).

In turn, the analysis of diffusion in the exterior of a compact planar domain
turns out to be more subtle. As planar diffusion is recurrent, the diffusing
particle never escapes to infinity and repeatedly returns to the boundary so
that the boundary local time $\ell_{t}$ grows to infinity as $t\to\infty$ and
thus crosses any threshold with probability $1$. However, the probability of
not hitting the boundary up to time $t$ is known to decay logarithmically
slowly with $t$ [26], so that each return to the boundary may take abnormally
long time. In this paper, we focus on diffusion outside a disk and derive
explicit integral representations for both probability densities
$\rho(\ell,t|\bm{x}_{0})$ and $U(\ell,t|\bm{x}_{0})$. We use then these
presentations to analyze their asymptotic behavior. We also characterize the
relative contributions of random trajectories with different exploration
sizes. Finally, we discuss extensions to more general planar domains, describe
some applications in physics and chemistry, and outline several open
questions.

##  2 Formal solution

In this Section, we recall the formal general solution for
$\rho(\ell,t|\bm{x}_{0})$ and $U(\ell,t|\bm{x}_{0})$ from [23, 18]. The
central relation is the expression for the moment-generating function of the
boundary local time $\ell_{t}$:

${\mathbb{E}}_{\bm{x}_{0}}\\{e^{-q\ell_{t}}\\}=S_{q}(t|\bm{x}_{0}),$ (7)

where $S_{q}(t|\bm{x}_{0})$ is the survival probability, which satisfies the
(backward) diffusion equation

$\partial_{t}S_{q}(t|\bm{x}_{0})=D\Delta
S_{q}(t|\bm{x}_{0})\qquad(\bm{x}_{0}\in\Omega),$ (8)

with the initial condition $S_{q}(0|\bm{x}_{0})=1$ and the mixed Robin-Neumann
boundary condition:

$\displaystyle\left.(\partial_{n}+q)S_{q}(t|\bm{x}_{0})\right|_{\Gamma}$
$\displaystyle=$ $\displaystyle 0,$ (9)
$\displaystyle\left.\partial_{n}S_{q}(t|\bm{x}_{0})\right|_{\partial\Omega\backslash\Gamma}$
$\displaystyle=$ $\displaystyle 0$ (10)

(for unbounded domains, the regularity condition $S_{q}(t|\bm{x}_{0})\to 1$ as
$|\bm{x}_{0}|\to\infty$ is also imposed). Here $\Delta$ is the Laplace
operator acting on $\bm{x}_{0}$, $\partial_{n}$ is the normal derivative at
the boundary oriented outward the domain, and $q\geq 0$ is a parameter in Eq.
(7), which can be related to the surface reactivity of the subset $\Gamma$
[18]. On the other hand, the moment-generating function is defined via the
probability density of $\ell_{t}$:

${\mathbb{E}}_{\bm{x}_{0}}\\{e^{-q\ell_{t}}\\}=\int\limits_{0}^{\infty}d\ell\,e^{-q\ell}\,\rho(\ell,t|\bm{x}_{0}).$
(11)

Inverting this Laplace transform, one formally gets

$\rho(\ell,t|\bm{x}_{0})={\mathcal{L}}_{q,\ell}^{-1}\\{S_{q}(t|\bm{x}_{0})\\}.$
(12)

The first-crossing time ${\mathcal{T}}_{\ell}$ of a given threshold $\ell\geq
0$ by the boundary local time $\ell_{t}$ is defined as

${\mathcal{T}}_{\ell}=\inf\\{t>0~{}:~{}\ell_{t}>\ell\\}.$ (13)

As the boundary local time is a non-decreasing process, the cumulative
distribution function of the first-crossing time is determined as

$\displaystyle Q(\ell,t|\bm{x}_{0})$ $\displaystyle=$
$\displaystyle{\mathbb{P}}_{\bm{x}_{0}}\\{{\mathcal{T}}_{\ell}<t\\}={\mathbb{P}}_{\bm{x}_{0}}\\{\ell_{t}>\ell\\}=1-{\mathcal{L}}_{q,\ell}^{-1}\bigl{\\{}S_{q}(t|\bm{x}_{0})/q\bigr{\\}}\,,$
(14)

where we used Eq. (12). The probability density of the first-crossing time
follows:

$U(\ell,t|\bm{x}_{0})=\partial_{t}Q(\ell,t|\bm{x}_{0})={\mathcal{L}}_{q,\ell}^{-1}\left\\{\frac{-\partial_{t}S_{q}(t|\bm{x}_{0})}{q}\right\\}.$
(15)

Recalling that $H_{q}(t|\bm{x}_{0})=-\partial_{t}S_{q}(t|\bm{x}_{0})$ is the
probability density of the first-passage time to a partially reactive target
$\Gamma$ (see, e.g., [26]), one can rewrite Eq. (15) as

$H_{q}(t|\bm{x}_{0})=\int\limits_{0}^{\infty}d\ell\,q\,e^{-q\ell}\,U(\ell,t|\bm{x}_{0}).$
(16)

In other words, the probability density $U(\ell,t|\bm{x}_{0})$ determines the
whole family of the first-passage times described by $H_{q}(t|\bm{x}_{0})$
[18].

Even though the probability densities $\rho(\ell,t|\bm{x}_{0})$ and
$U(\ell,t|\bm{x}_{0})$ are formally determined via Eqs. (12, 15), their
analysis requires the Laplace transform inversion,
${\mathcal{L}}_{q,\ell}^{-1}$, which may be numerically unstable [27].
Moreover, as the parameter $q$ enters through the Robin boundary condition
(9), the dependence of the survival probability $S_{q}(t|\bm{x}_{0})$ on $q$
is implicit that renders its analysis even more difficult. The aim of the
present paper is to derive fully explicit representations for both densities
in the case of diffusion outside a disk and then to use them for the
asymptotic analysis.

##  3 Diffusion outside a disk

In this section, we consider reflected Brownian motion in the exterior of the
disk of radius $R$, $\Omega=\\{\bm{x}\in{\mathbb{R}}^{2}~{}:~{}|\bm{x}|>R\\}$,
and focus on the boundary local time $\ell_{t}$ on the boundary of that disk:
$\Gamma=\partial\Omega$ (Fig. 1). The rotational invariance of the domain
implies that the probability densities and related quantities do not depend on
the angular coordinate. In the following, we replace the starting point
$\bm{x}_{0}$ by its radial coordinate $r_{0}=|\bm{x}_{0}|$.

We start by recalling the integral representation for the survival probability
$S_{q}(t|r_{0})$, which can be derived from the classical solution of a
similar heat conduction problem [28, 29]

$\displaystyle\qquad S_{q}(t|r_{0})$ $\displaystyle=$
$\displaystyle\frac{2qR}{\pi}\int\limits_{0}^{\infty}\frac{dz}{z}\,e^{-z^{2}Dt/R^{2}}$
(17) $\displaystyle\times$
$\displaystyle\frac{Y_{0}(zr_{0}/R)(zJ_{1}(z)+qRJ_{0}(z))-J_{0}(zr_{0}/R)(zY_{1}(z)+qRY_{0}(z))}{(zJ_{1}(z)+qRJ_{0}(z))^{2}+(zY_{1}(z)+qRY_{0}(z))^{2}}\,,$

where $J_{\nu}(z)$ and $Y_{\nu}(z)$ are the Bessel functions of the first and
second kind, respectively. To render the dependence of this expression on $q$
even more explicit, we represent it as

$S_{q}(t|r_{0})=\frac{4}{\pi}\int\limits_{0}^{\infty}\frac{dz}{z}e^{-z^{2}Dt/R^{2}}{\rm
Re}\left(\frac{qA(z,r_{0}/R)}{B(z)/R+q}\right),$ (18)

where

$\displaystyle A(z,r)$ $\displaystyle=$
$\displaystyle\frac{i}{2}\,\frac{(J_{0}(z)+iY_{0}(z))(J_{0}(zr)-iY_{0}(zr))}{J_{0}^{2}(z)+Y_{0}^{2}(z)}\,,$
(19) $\displaystyle B(z)$ $\displaystyle=$
$\displaystyle\frac{z(J_{0}(z)J_{1}(z)+Y_{0}(z)Y_{1}(z))+i\frac{2}{\pi}}{J_{0}^{2}(z)+Y_{0}^{2}(z)}\,.$
(20)

Substituting Eq. (18) into Eq. (12) and evaluating the inverse Laplace
transform with respect to $q$, we get

$\displaystyle\rho(\ell,t|r_{0})$ $\displaystyle=$ $\displaystyle
S_{\infty}(t|r_{0})\delta(\ell)$ $\displaystyle-$ $\displaystyle\frac{4}{\pi
R}\int\limits_{0}^{\infty}\frac{dz}{z}e^{-z^{2}Dt/R^{2}}\,{\rm
Re}\biggl{(}A(z,r_{0}/R)B(z)e^{-B(z)\ell/R}\biggr{)},$

where

$S_{\infty}(t|r_{0})=\frac{4}{\pi}\int\limits_{0}^{\infty}\frac{dz}{z}e^{-z^{2}Dt/R^{2}}\,{\rm
Re}\bigl{(}A(z,r_{0}/R)\bigr{)}$ (22)

is obtained as the limit of Eq. (18) when $q\to\infty$. As in Eqs. (3, 1), the
first term in Eq. (3) accounts for trajectories that did not hit the boundary
up to time $t$, for which the boundary local time $\ell_{t}$ remained $0$. The
positive-order moments of the boundary local time are analyzed in A.
Similarly, Eq. (14) implies that the inverse Laplace transform of
$S_{q}(t|r_{0})/q$ with respect to $q$ yields

$Q(\ell,t|r_{0})=1-\frac{4}{\pi}\int\limits_{0}^{\infty}\frac{dz}{z}e^{-z^{2}Dt/R^{2}}\,{\rm
Re}\left(A(z,r_{0}/R)e^{-B(z)\ell/R}\right).$ (23)

The time derivative of this expression gives the probability density
$U(\ell,t|r_{0})$ of the first-crossing time:

$U(\ell,t|r_{0})=\frac{4D}{\pi
R^{2}}\int\limits_{0}^{\infty}dz\,z\,e^{-z^{2}Dt/R^{2}}\,{\rm
Re}\left(A(z,r_{0}/R)e^{-B(z)\ell/R}\right).$ (24)

Note also that setting $\ell=0$ in Eq. (24) yields

$\displaystyle U(0,t|r_{0})$ $\displaystyle=$ $\displaystyle\frac{4D}{\pi
R^{2}}\int\limits_{0}^{\infty}dz\,z\,e^{-z^{2}Dt/R^{2}}\,{\rm
Re}\bigl{(}A(z,r_{0}/R)\bigr{)}=H_{\infty}(t|r_{0}),$ (25)

where we used Eq. (22) and thus retrieved the probability density of the
first-passage time to the disk.

The explicit integral representations (3, 23, 24) are the main analytical
results of the paper. Thanks to the exponential factor $e^{-z^{2}Dt/R^{2}}$,
these integrals rapidly converge for large $z$. In turn, the integrals in Eqs.
(3, 23) exhibit logarithmically slow convergence at small $z$. A practical
solution of this issue is discussed in B.

The boundary local time $\ell_{t}$ remains zero until the first encounter with
the boundary. As a consequence, the first-crossing time ${\mathcal{T}}_{\ell}$
can be decomposed into two contributions: the first-passage time from
$\bm{x}_{0}$ to the circle,
${\mathcal{T}}_{0,r_{0}}=\inf\\{t>0~{}:~{}\ell_{t}>0~{}|~{}|\bm{X}_{0}|=r_{0}\\}$,
and the first-crossing time of the level $\ell$ after starting from the
circle,
${\mathcal{T}}_{\ell,R}=\inf\\{t>0~{}:~{}\ell_{t}>\ell~{}|~{}|\bm{X}_{0}|=R\\}$:

${\mathcal{T}}_{\ell}={\mathcal{T}}_{0,r_{0}}+{\mathcal{T}}_{\ell,R}.$ (26)

The strong Markovian character of reflected Brownian motion and of the
boundary local time, as well as the rotational invariance of the problem imply
that these two contributions are independent variables. As a consequence, the
probability density of ${\mathcal{T}}_{\ell}$ can be obtained by convolving
the densities of ${\mathcal{T}}_{0,r_{0}}$ and ${\mathcal{T}}_{\ell,R}$. The
probability density of the first-passage time to the circle,
$H_{\infty}(t|r_{0})$, has been studied long ago (see [26, 30, 29] and
references therein). As a consequence, one can focus on the probability
density $U(\ell,t|R)$ of the second contribution ${\mathcal{T}}_{\ell,R}$, for
which $A(z,1)=i/2$, and Eq. (24) becomes

$U(\ell,t|R)=-\frac{2D}{\pi
R^{2}}\int\limits_{0}^{\infty}dz\,z\,e^{-z^{2}Dt/R^{2}}\,{\rm
Im}\left(e^{-B(z)\ell/R}\right).$ (27)

In the next section, we discuss the asymptotic behavior of the density
$U(\ell,t|r_{0})$.

##  4 Asymptotic analysis

To investigate the asymptotic behavior of the probability density
$U(\ell,t|r_{0})$, it is convenient to change the integration variable in Eq.
(24) as

$\displaystyle U(\ell,t|r_{0})$ $\displaystyle=$ $\displaystyle\frac{4}{\pi
t}\int\limits_{0}^{\infty}dz\,z\,e^{-z^{2}}\,{\rm
Re}\biggl{(}A\bigl{(}zR/\sqrt{Dt},r_{0}/R\bigr{)}e^{-B(zR/\sqrt{Dt})\,\ell/R}\biggr{)}.$
(28)

###  4.1 Short-time asymptotic behavior

The limit $t\to 0$ corresponds to the large-$z$ expansions of $A(z,r)$ and
$B(z)$:

$\displaystyle A(z,r)$ $\displaystyle\simeq$
$\displaystyle\frac{i}{2\,r^{\frac{1}{2}}}e^{-iz(r-1)}\qquad(z\to\infty),$
(29) $\displaystyle B(z)$ $\displaystyle\simeq$
$\displaystyle\frac{1}{2}+iz+O(1/z)\qquad(z\to\infty).$ (30)

Substituting these approximations into Eq. (28), we get

$U(\ell,t|r_{0})\simeq\frac{e^{-\ell/(2R)}}{(r_{0}/R)^{\frac{1}{2}}}\,\frac{(r_{0}-R+\ell)\,e^{-(r_{0}-R+\ell)^{2}/(4Dt)}}{\sqrt{4\pi
Dt^{3}}}\qquad(t\to 0).$ (31)

Apart from the factor $e^{-\ell/(2R)}\,(r_{0}/R)^{-\frac{1}{2}}$, this
expression is identical to Eq. (4) for diffusion in the half-line, if
$x_{0}=r_{0}-R$ denotes the distance to the boundary. This is not surprising
because the circle looks locally flat at short times, and the behavior should
be close to that of the half-plane. In turn, the supplementary factor
$e^{-\ell/(2R)}\,(r_{0}/R)^{-\frac{1}{2}}$ accounts for the curvature of the
boundary. Note that in the case $r_{0}=R$, the same result could alternatively
be derived by analyzing Eq. (49) in the short-time limit (i.e., $p\to\infty$).

In the limit $t\to 0$, the particle does not spend enough time near the
boundary to allow for the boundary local time $\ell_{t}$ to cross a given
threshold $\ell$, and the probability density of too short first-crossing
times, $t\ll\ell^{2}/(4D)$, is extremely small, even if the particle starts on
the boundary. At first thought, this statement may sound to contradict the
well-known property that (reflected) Brownian motion that crossed a smooth
boundary returns infinitely many times to that boundary during an infinitely
short period [31]. This seeming contradiction reflects the crucial difference
between the boundary local time, $\ell_{t}$, and the number of encounters
$N_{t}^{a}$ with a thin boundary layer of width $a$ (i.e., the number of
returns up to $t$). The latter is defined for any $a>0$ but diverges in the
limit $a\to 0$. This divergence, $N_{t}^{a}\to N_{t}^{0}=+\infty$, is evoked
in the above statement about infinitely many returns. In turn, the boundary
local time $\ell_{t}$, which is obtained by rescaling $N_{t}^{a}$ by $a$ in
Eq. (2), remains a meaningful characteristics of the encounters with the
boundary in the limit $a\to 0$.

###  4.2 Long-time asymptotic behavior

The large-$t$ limit corresponds to the small-$z$ behavior of $A(z,r)$ and
$B(z)$:

$\displaystyle A(z,r)$ $\displaystyle\simeq$
$\displaystyle\frac{i}{2}\bigl{(}1-B(z)\ln r\bigr{)}\qquad(z\to 0),$ (32)
$\displaystyle B(z)$ $\displaystyle\simeq$ $\displaystyle\frac{1}{\ln(1/z)+\ln
2-\gamma-\frac{\pi i}{2}}\qquad(z\to 0)\,,$ (33)

where $\gamma\approx 0.5772\ldots$ is the Euler constant. We substitute these
expressions in Eq. (28) and then neglect $\ln z$ as a slowly varying function
to get

$\displaystyle\qquad U(\ell,t|r_{0})$ $\displaystyle\simeq$
$\displaystyle\frac{1}{\pi
t}\exp\left(-\frac{(\ell/R)b_{t}}{b_{t}^{2}+(\pi/2)^{2}}\right)\left[\sin\left(\frac{(\ell/R)\pi/2}{b_{t}^{2}+(\pi/2)^{2}}\right)\left(1-\frac{b_{t}\ln(r_{0}/R)}{b_{t}^{2}+(\pi/2)^{2}}\right)\right.$
$\displaystyle+$
$\displaystyle\left.\cos\left(\frac{(\ell/R)\pi/2}{b_{t}^{2}+(\pi/2)^{2}}\right)\frac{(\pi/2)\ln(r_{0}/R)}{b_{t}^{2}+(\pi/2)^{2}}\right],$
(34)

where the remaining integral over $z$ gave the factor $1/2$, and we set
$b_{t}=\ln(\sqrt{4Dt}/R)-\gamma$. When $b_{t}\gg\pi/2$ and
$b_{t}\gg\sqrt{\ell/R}$, we finally get

$U(\ell,t|r_{0})\simeq
2\frac{\ell/R+\ln(r_{0}/R)}{t\,[\ln(Dt/R^{2})]^{2}}\qquad(t\to\infty).$ (35)

For comparison, we also present the long-time asymptotic behavior of the
probability density $H_{q}(t|r_{0})$ of the first-passage time to a partially
reactive target, which was derived in [29]:

$H_{q}(t|r_{0})\simeq\frac{2\biggl{(}\frac{1}{qR}+\ln(r_{0}/R)\biggr{)}}{t\biggl{(}\pi^{2}+\left[\ln(R^{2}/(4Dt))+2\gamma-\frac{2}{qR}\right]^{2}\biggr{)}}\,.$
(36)

In the very long time limit, one gets

$H_{q}(t|r_{0})\simeq
2\frac{\frac{1}{qR}+\ln(r_{0}/R)}{t\bigl{[}\ln(Dt/R^{2})\bigr{]}^{2}}\,.$ (37)

In the particular case $q=\infty$, this behavior coincides with Eq. (35) at
$\ell=0$, in agreement with Eq. (25). More generally, both densities
$U(\ell,t|r_{0})$ and $H_{q}(t|r_{0})$ that are related by Eq. (16), exhibit
the same dependence on time in the limit $t\to 0$. Curiously, the relation
(16) is satisfied even for the asymptotic forms in Eqs. (35, 37).

One can see that the probability density $U(\ell,t|r_{0})$ exhibits a very
heavy tail at large $t$. In particular, any positive moment of the first-
crossing time, ${\mathbb{E}}_{\bm{x}_{0}}\\{{\mathcal{T}}_{\ell}^{k}\\}$, with
$k\in{\mathbb{R}}_{+}$, is infinite. This divergence is caused by the
contribution of very long trajectories that explore the plane and can move far
away from the target but unavoidably return to it due to the recurrent
character of planar diffusion. The same mechanism leads to the infinite mean
first-crossing time for diffusion in the half-line. However, the $t^{-3/2}$
decay of $U(\ell,t|x_{0})$ in Eq. (4) ensures that at least the moments
${\mathbb{E}}_{x_{0}}\\{{\mathcal{T}}_{\ell}^{k}\\}$ with $k<1/2$ exist in the
one-dimensional setting.

Figure 2 illustrates the behavior of the probability density $U(\ell,t|R)$ for
several values of the threshold $\ell$. For cross-validation, this density was
computed by two independent methods: (i) the numerical evaluation of the
integral in Eq. (27), and (ii) the numerical computation of the inverse
Laplace transform in Eq. (52) by the Talbot algorithm. Both methods show an
excellent agreement (circles versus crosses). At $\ell/R=0.1$ and $\ell/R=1$,
the short-time asymptotic relation (31) accurately captures the behavior of
$U(\ell,t|R)$ at short and even moderate times, $t\lesssim R^{2}/D$, but fails
at long times, as expected. In particular, one can use the simple form of Eq.
(31) to estimate the most probable first-crossing time: $t_{\rm
mp}\simeq(r_{0}-R+\ell)^{2}/(6D)$. As the mean value is infinite, the most
probable time plays the role of a natural time scale. At long times, $t\gg
R^{2}/D$, Eq. (35) correctly captures the long-time behavior but approaches it
slowly, particularly for large $\ell$. In contrast, the approximate relation
(34) turns out to be surprisingly accurate for large and even intermediate
times.

The quality of both asymptotic relations (31, 35) is much lower in the case
$\ell/R=10$ (Fig. 2(c)). On one hand, Eq. (31) does not reproduce correctly
the maximum of $U(\ell,t|R)$ and captures only the steep decay at small $t$.
In other words, in the range of validity of this relation (i.e., $t\lesssim
R^{2}/D$), the density $U(\ell,t|R)$ is already very small and thus not much
useful in practice. On the other hand, the approach to Eq. (35) is very slow
so that this asymptotic relation becomes accurate only at extremely large
times. In contrast, one can still rely on the approximate Eq. (34).

Figure 2:  The probability density $U(\ell,t|R)$ as a function of $t$ for
diffusion outside the disk of radius $R=1$, with $D=1$, $\ell=0.1$ (a),
$\ell=1$ (b), and $\ell=10$ (c). Empty circles show the exact solution
obtained by numerical computation of the integral in Eq. (27); dashed line
presents the short-time asymptotic behavior (31); solid and dash-dotted lines
show the long-time asymptotic relations (34) and (35), respectively; crosses
present the numerical computation of the inverse Laplace transform in Eq. (52)
by the Talbot algorithm.

###  4.3 Exploration size of trajectories

In order to quantify the relative contributions of trajectories with different
spatial extents into the heavy tail of the probability density $U(\ell,t|R)$,
we introduce an absorbing circle of radius $L>R$ and compute the probability
density of the first-crossing time ${\mathcal{T}}_{\ell}$ for the trajectories
that remained inside that circle up to ${\mathcal{T}}_{\ell}$. In order words,
while computing the statistics of the first-crossing times, we discard all
trajectories that hit the circle of radius $L$ up to time
${\mathcal{T}}_{\ell}$. In D, we describe an extension of the spectral
approach from [23, 18, 24, 25] to compute the density $U_{L}(\ell,t|r_{0})$ in
the presence of the absorbing circle. Expectedly, the density
$U_{L}(\ell,t|r_{0})$ is not normalized to $1$, and
$1-\int\nolimits_{0}^{\infty}dt\,U_{L}(\ell,t|r_{0})$ is the probability that
the trajectory $\bm{X}_{t}$ has crossed the circle of radius $L$ before the
associated boundary local time could reach the level $\ell$. In this setting,
we could not derive a fully explicit representation for $U_{L}(\ell,t|r_{0})$
and thus performed the Laplace transform inversion in Eq. (56) numerically by
the Talbot algorithm.

Figure 3 shows the probability densities $U_{L}(\ell,t|R)$ for several values
of $L$ (with $L=\infty$ corresponding to the former case without absorbing
circle). As the “restricted” trajectories explore a bounded domain between two
concentric circles (of radii $R$ and $L$), the particle returns to the inner
circle (target) more often, and the boundary local time increases faster. As a
consequence, the probability density $U_{L}(\ell,t|r_{0})$ of the associated
first-crossing time is expected to decay much faster as $t\to\infty$, as
confirmed by Fig. 3. This faster decay is also related to the fact that
keeping long trajectory within a bounded region becomes more and more unlikely
as time goes on. This figure explicitly illustrates that the heavy tail of the
density $U(\ell,t|R)$ in Eq. (35) is caused by very long trajectories that
explore the unbounded planar space and thus create long stalling periods,
during which the boundary local time does not change.

Figure 3:  The probability density $U(\ell,t|R)$ (solid line) as a function of
$t$ for diffusion outside the disk of radius $R=1$, with $D=1$, $\ell=1$ (a)
and $\ell=10$ (b). This density is compared to the probability densities
$U_{L}(\ell,t|R)$ in the presence of the absorbing circle of radius $L$, with
three values of $L$ as indicated in the plot. All these densities were
obtained via the numerical computation of the inverse Laplace transform in Eq.
(56) by the Talbot algorithm. Note that two last curves in the panel (a) were
truncated at large times to avoid numerical instabilities of the Laplace
transform inversion.

##  5 Discussion and conclusion

In this paper, we considered diffusion in the exterior of a disk and derived
the explicit integral representations (3) and (24) for the probability
densities $\rho(\ell,t|r_{0})$ and $U(\ell,t|r_{0})$ of the boundary local
time $\ell_{t}$ and of the first-crossing time ${\mathcal{T}}_{\ell}$,
respectively. As the former density $\rho(\ell,t|r_{0})$ was investigated in
[23], we mainly focused on the latter one. With the aid of the integral
representation (24), we deduced Eqs. (31, 35) describing the short-time and
the long-time asymptotic behavior of $U(\ell,t|r_{0})$, respectively.
Moreover, Eq. (34) was checked to be an accurate approximation of
$U(\ell,t|r_{0})$ for moderate and long times. We illustrated that the
probability density $U(\ell,t|r_{0})$ vanishes at short times, reaches the
maximum at intermediate times $\sim\frac{1}{D}(r_{0}-R+\ell)^{2}$, and
exhibits a heavy long-time tail such that all positive-order moments are
infinite. To quantify the effect of long excursions between encounters, we
also analyzed the probability density $U_{L}(\ell,t|r_{0})$ in the presence of
an absorbing circle of radius $L$.

As discussed in Secs. 1 and 2, the boundary local time $\ell_{t}$ is a proxy
of the number of encounters, $\mathcal{N}_{t}^{a}$, of a diffusing particle
with the boundary layer $\Gamma_{a}$ of width $a$ up to time $t$ (here,
$\Gamma=\partial\Omega$). Consequently, the first-crossing time
${\mathcal{T}}_{\ell}$ describes the moment when the number of encounters
exceeds a prescribe threshold. While the random variables $\ell_{t}$ and
${\mathcal{T}}_{\ell}$ characterize the reflecting boundary, they form a
natural ground for incorporating surface reactivity and to describe various
diffusion-mediated surface phenomena such as diffusion-influenced chemical
reactions, permeation across biological membranes, surface relaxation in
nuclear magnetic resonance, etc. [19, 18]. For instance,
$U(0,t|r_{0})=H_{\infty}(t|r_{0})$ describes the first-passage time to the
perfectly reactive boundary [32], whereas $H_{q}(t|r_{0})$, expressed through
Eq. (16) in terms of $U(\ell,t|r_{0})$, determines the reaction time on a
partially reactive boundary [33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44,
45, 46, 47, 48, 49]. More sophisticated surface reaction mechanisms can also
be implemented with the help of $U(\ell,t|r_{0})$, as discussed in [18].

Interestingly, the probability densities $\rho(\ell,t|r_{0})$ and
$U(\ell,t|r_{0})$ also describe the first-encounter statistics of two
independent particles diffusing in the plane. In fact, associating the
coordinate frame with one of the particles, the diffusive dynamics of two
circular particles of radii $R_{1}$ and $R_{2}$ with diffusivities $D_{1}$ and
$D_{2}$ can be mapped onto planar diffusion of a single point-like particle
with diffusivity $D=D_{1}+D_{2}$ toward a static circular target of radius
$R=R_{1}+R_{2}$. As a consequence, $\ell_{t}/a$ is the number of encounters of
such two particles, where $a$ is the separation between two particles at each
encounter. In other words, a single encounter is counted when two particles
hit each other and then moved away and become separated by a distance
exceeding $a$. In turn, ${\mathcal{T}}_{aN}$ is the first moment when the
number of such encounters exceeds $N$. As an encounter of two species is the
necessary step for most bimolecular reactions, reproduction in animals, and
virus spreading in humans, both quantities can be used in describing various
encounters in chemistry, biology and ecology (see [50, 51, 52, 53, 54, 55, 56,
57, 58] and references therein).

Even though the explicit representations (3) and (24) were derived for the
particular case of a disk, the asymptotic behavior of the probability density
$U(\ell,t|\bm{x}_{0})$ is expected to remain valid for diffusion in the
exterior of a general compact planar domain with smooth boundary. Indeed, at
short times, $t\ll R_{c}^{2}/D$, the boundary looks locally flat with respect
to the radius of curvature $R_{c}$ of the boundary in a vicinity of the
starting point. As a consequence, the time dependence of the density
$U(\ell,t|\bm{x}_{0})$ is expected to be given by Eq. (4) for the one-
dimensional setting, where $x_{0}=|\bm{x}_{0}-\Gamma|$ is the distance between
the starting point $\bm{x}_{0}$ and the target $\Gamma$. In turn, the
correction factor $e^{-\ell/(2R)}\,(r_{0}/R)^{-\frac{1}{2}}$ from Eq. (31) is
specific to the disk. Finding the shape-dependent correction factor for other
domains presents an interesting perspective. Similarly, at large times, long
excursions between encounters are expected to lead to a heavy tail of the
probability density $U(\ell,t|\bm{x}_{0})$ as $t\to\infty$. The characteristic
$1/(t\ln^{2}(t))$ decay may remain universal but the dependence on $\ell$ and
$\bm{x}_{0}$ could be shape-specific. Validation and specification of these
conjectures present another exciting direction for future research. This
problem is tightly related to the asymptotic behavior of the spectrum of the
Dirichlet-to-Neumann operator associated with the modified Helmholtz equation
[23] (see also D).

## Acknowledgments

The author acknowledges a partial financial support from the Alexander von
Humboldt Foundation through a Bessel Research Award.

##  Appendix A Long-time behavior of the mean boundary local time

According to Eq. (3), the probability density $\rho(\ell,t|r_{0})$ exhibits a
rapid decay for large $\ell$, ensuring the existence of all positive-order
moments of the boundary local time $\ell_{t}$:

${\mathbb{E}}_{\bm{x}_{0}}\\{\ell_{t}^{k}\\}=\int\limits_{0}^{\infty}d\ell\,\ell^{k}\,\rho(\ell,t|\bm{x}_{0})\qquad(k>0).$
(38)

The explicit representation (3) can be used to compute these moments. However,
the direct exchange of the order of integrals over $\ell$ in Eq. (38) and over
$z$ in Eq. (3) is not possible as the resulting integral over $z$ would be
divergent. To overcome this limitation, one can first express the time
derivative of the moment ${\mathbb{E}}_{\bm{x}_{0}}\\{\ell_{t}^{k}\\}$ as

$\partial_{t}{\mathbb{E}}_{\bm{x}_{0}}\\{\ell_{t}^{k}\\}=\frac{4DR^{k-2}\Gamma(k+1)}{\pi}\int\limits_{0}^{\infty}dz\,z\,e^{-z^{2}Dt/R^{2}}\,{\rm
Re}\left(\frac{A(z,r_{0}/R)}{[B(z)]^{k}}\right).$ (39)

Rescaling the integral variable yields

$\partial_{t}{\mathbb{E}}_{\bm{x}_{0}}\\{\ell_{t}^{k}\\}=\frac{4R^{k}\Gamma(k+1)}{\pi
t}\int\limits_{0}^{\infty}dz\,z\,e^{-z^{2}}\,{\rm
Re}\left(\frac{A(zR/\sqrt{Dt},r_{0}/R)}{\bigl{[}B(zR/\sqrt{Dt})\bigr{]}^{k}}\right).$
(40)

In the long-time limit, one uses the small-$z$ expansions (33, 32) to get

$\displaystyle\partial_{t}{\mathbb{E}}_{\bm{x}_{0}}\\{\ell_{t}^{k}\\}$
$\displaystyle\simeq$ $\displaystyle\frac{2R^{k}\Gamma(k+1)}{\pi
t}\int\limits_{0}^{\infty}dz\,z\,e^{-z^{2}}\,{\rm
Re}\left(i\biggl{[}(b_{t}-\ln z-\pi i/2)^{k}\right.$ (41) $\displaystyle-$
$\displaystyle\left.\ln(r_{0}/R)(b_{t}-\ln z-\pi i/2)^{k-1}\biggr{]}\right),$

where $b_{t}=\ln(\sqrt{4Dt}/R)-\gamma$. Integrating this expression over $t$,
we get

$\displaystyle{\mathbb{E}}_{\bm{x}_{0}}\\{\ell_{t}^{k}\\}$
$\displaystyle\simeq$
$\displaystyle\frac{2R^{k}\Gamma(k+1)}{\pi}\int\limits_{0}^{\infty}dz\,e^{-z}\,{\rm
Im}\left(-\frac{(b_{t}-\frac{1}{2}\ln z-\frac{1}{2}\pi i)^{k+1}}{k+1}\right.$
(42) $\displaystyle+$
$\displaystyle\left.\ln(r_{0}/R)\frac{(b_{t}-\frac{1}{2}\ln z-\frac{1}{2}\pi
i)^{k}}{k}\right),$

where we also changed the integration variable $z^{2}\to z$. Using the
binomial expansion, one can evaluate these integrals term by term. In
particular, we get the mean and the variance in the long-time limit:

$\displaystyle{\mathbb{E}}_{\bm{x}_{0}}\\{\ell_{t}\\}$ $\displaystyle\simeq$
$\displaystyle
R\biggl{(}\ln\bigl{(}\sqrt{4Dt}/r_{0}\bigr{)}-\gamma/2\biggr{)},$ (43)
$\displaystyle\mathrm{var}_{\bm{x}_{0}}\\{\ell_{t}\\}$ $\displaystyle=$
$\displaystyle{\mathbb{E}}_{\bm{x}_{0}}\\{\ell_{t}^{2}\\}-\bigl{(}{\mathbb{E}}_{\bm{x}_{0}}\\{\ell_{t}\\}\bigr{)}^{2}$
(44) $\displaystyle\simeq$ $\displaystyle
R^{2}\left(\bigl{(}\ln\bigl{(}\sqrt{4Dt}/R\bigr{)}-\gamma/2\bigr{)}^{2}-\ln^{2}(r_{0}/R)-\frac{\pi^{2}}{12}\right).$

These expressions generalize the former results from [23] to an arbitrary
starting point.

##  Appendix B Numerical computation

Using Eqs. (33, 32), one can easily check that the integrals in Eqs. (3, 23)
converge logarithmically slowly near $z=0$. To improve the accuracy of
numerical computations, one can choose an appropriate $\varepsilon\ll 1$ and
split the integral into two parts, from $0$ to $\varepsilon$, and from
$\varepsilon$ to $\infty$. The second integral is then evaluated numerically,
whereas the contribution of the first integral, $I_{\varepsilon}$, can be
found approximately.

Let us first consider the computation of $Q(\ell,t|r_{0})$ in Eq. (23), for
which

$\displaystyle I_{\varepsilon}$ $\displaystyle=$
$\displaystyle\frac{4}{\pi}\int\limits_{0}^{\varepsilon}\frac{dz}{z}e^{-z^{2}Dt/R^{2}}{\rm
Re}\left(A(z,r)e^{-B(z)\ell/R}\right)$ $\displaystyle\approx$
$\displaystyle\frac{4}{\pi}\int\limits_{0}^{\varepsilon}\frac{dz}{z}{\rm
Re}\left(\frac{i}{2}\bigl{(}1-B(z)\ln r\bigr{)}e^{-B(z)\ell/R}\right)$
$\displaystyle\approx$
$\displaystyle-\frac{2}{\pi}\int\limits_{\ln(1/\varepsilon)}^{\infty}dy\,{\rm
Im}\left(\biggl{(}1-\frac{\ln
r}{y+b}\biggr{)}\exp\biggl{(}-\frac{c}{y+b}\biggr{)}\right),$

where $c=\ell/R$ and $b=\ln 2-\gamma-\frac{1}{2}\pi i$. Here we approximated
$e^{-z^{2}Dt/R^{2}}\approx 1$ for $0\leq z\leq\varepsilon$, and made the
change of variables: $y=\ln(1/z)$. Note that the above approximation naturally
imposes the constraint on $\varepsilon$: $\varepsilon^{2}\ll R^{2}/(Dt)$. In
particular, as $t$ increases, one has to choose smaller and smaller
$\varepsilon$.

The last integral can be evaluated as

$I_{\varepsilon}\approx 1+\frac{2}{\pi}\,{\rm
Im}\biggl{(}b_{\varepsilon}\,e^{-c/b_{\varepsilon}}-\bigl{(}c+\ln(r_{0}/R)\bigr{)}\,{\rm
Ei}_{1}\bigl{(}c/b_{\varepsilon}\bigr{)}\biggr{)},$ (45)

where ${\rm Ei}_{1}(z)=\int\nolimits_{1}^{\infty}dk\,e^{-zk}/k=\Gamma(0,z)$ is
the exponential integral, and $b_{\varepsilon}=b+\ln(1/\varepsilon)$. As a
consequence, we get

$Q(\ell,t|r_{0})=1-I_{\varepsilon}-\frac{4}{\pi}\int\limits_{\varepsilon}^{\infty}\frac{dz}{z}e^{-z^{2}Dt/R^{2}}\,{\rm
Re}\left(A(z,r_{0}/R)e^{-B(z)\ell/R}\right).$ (46)

A similar correction can be used for an accurate computation of
$\rho(\ell,t|r_{0})$:

$\displaystyle\rho(\ell,t|r_{0})$ $\displaystyle=$ $\displaystyle
S_{\infty}(t|r_{0})\delta(\ell)-I_{\varepsilon}^{(\rho)}$ $\displaystyle-$
$\displaystyle\frac{4}{\pi
R}\int\limits_{\varepsilon}^{\infty}\frac{dz}{z}e^{-z^{2}Dt/R^{2}}\,{\rm
Re}\left(A(z,r_{0}/R)B(z)e^{-B(z)\ell/R}\right),$

with

$\displaystyle I_{\varepsilon}^{(\rho)}$ $\displaystyle=$
$\displaystyle\frac{4}{\pi
R}\int\limits_{0}^{\varepsilon}\frac{dz}{z}e^{-z^{2}Dt/R^{2}}\,{\rm
Re}\left(A(z,r_{0}/R)B(z)e^{-B(z)\ell/R}\right)$ (48) $\displaystyle\approx$
$\displaystyle\frac{2}{\pi R}\,{\rm Im}\left({\rm
Ei}_{1}(c/b_{\varepsilon})-\frac{\ln(r_{0}/R)}{c}e^{-c/b_{\varepsilon}}\right).$

Note that there is no need for such a correction for computing
$U(\ell,t|r_{0})$ because $I_{\varepsilon}$ does not depend on time and thus
disappears after differentiating $Q(\ell,t|r_{0})$ with respect to $t$.

##  Appendix C Alternative representation

In [23, 18], alternative representations for the probability densities
$\rho(\ell,t|\bm{x}_{0})$ and $U(\ell,t|\bm{x}_{0})$ were developed in terms
of the eigenmodes of the Dirichlet-to-Neumann operator ${\mathcal{M}}_{p}$ for
an arbitrary Euclidean domain $\Omega$ with a smooth bounded boundary
$\partial\Omega$. For a given function $f$ on $\partial\Omega$, the operator
${\mathcal{M}}_{p}$ associates another function
$g=(\partial_{n}w)|_{\partial\Omega}$ on that boundary, where $w$ satisfies
the modified Helmholtz equation $(p-D\Delta)w=0$ in $\Omega$ with Dirichlet
boundary condition $w|_{\partial\Omega}=f$ [59, 60, 61, 62, 63]. For instance,
for diffusion outside the disk of radius $R$, the probability $Q(\ell,t|R)$
reads

$Q(\ell,t|R)={\mathcal{L}}_{p,t}^{-1}\biggl{\\{}\frac{1}{p}\exp\bigl{(}-\mu_{0}^{(p)}\ell\bigr{)}\biggr{\\}},$
(49)

where ${\mathcal{L}}_{p,t}^{-1}$ is the inverse Laplace transform with respect
to $p$,

$\mu_{0}^{(p)}=\sqrt{p/D}\frac{K_{1}(R\sqrt{p/D})}{K_{0}(R\sqrt{p/D})}$ (50)

is the smallest eigenvalue of the Dirichlet-to-Neumann operator
${\mathcal{M}}_{p}$ in the exterior of the disk, and $K_{\nu}(z)$ is the
modified Bessel function of the second kind [23]. The densities
$\rho(\ell,t|R)$ and $U(\ell,t|R)$ follow immediately by taking the
derivatives with respect to $\ell$ and $t$, respectively:

$\rho(\ell,t|R)={\mathcal{L}}_{p,t}^{-1}\biggl{\\{}\frac{\mu_{0}^{(p)}}{p}\exp(-\mu_{0}^{(p)}\ell)\biggr{\\}}.$
(51)

and

$U(\ell,t|R)={\mathcal{L}}_{p,t}^{-1}\biggl{\\{}\exp(-\mu_{0}^{(p)}\ell)\biggr{\\}}.$
(52)

In this Appendix, we briefly describe a direct derivation of Eq. (52) that
illustrates the meaning of the boundary local time as a proxy for the number
of encounters. Let us consider a regularized version of the problem, in which
the particle after hitting the circle at some point $(R,\varphi)$ (in polar
coordinates) is released from a bulk point $(R+a,\varphi)$, a distance $a$
above the circle (Fig. 1). The first-crossing time ${\mathcal{T}}_{\ell}$ can
then be represented as
${\mathcal{T}}_{\ell}=\tau_{1}+\tau_{2}+\ldots+\tau_{N}$, where $\tau_{i}$ is
the duration of bulk diffusion after $i$-th release (here we assumed that the
particle starts on the circle), and $N\approx\ell/a$ is the number of such
bulk explorations, see Eq. (2). The random variables $\tau_{i}$ are
independent and identically distributed, with a well-known moment-generating
function [26]:

${\mathbb{E}}\\{e^{-p\tau_{i}}\\}=\frac{K_{0}\bigl{(}(R+a)\sqrt{p/D}\bigr{)}}{K_{0}\bigl{(}R\sqrt{p/D}\bigr{)}}\,.$
(53)

As a consequence, one has

$\displaystyle{\mathbb{E}}\\{e^{-p{\mathcal{T}}_{\ell}}\\}$ $\displaystyle=$
$\displaystyle\biggl{(}{\mathbb{E}}\\{e^{-p\tau_{1}}\\}\biggr{)}^{N}$
$\displaystyle=$
$\displaystyle\exp\left(N\ln\left(\frac{K_{0}((R+a)\sqrt{p/D})}{K_{0}(R\sqrt{p/D})}\right)\right)$
$\displaystyle\simeq$
$\displaystyle\exp\left(\frac{\ell}{a}\ln\left(1+\frac{K^{\prime}_{0}(R\sqrt{p/D})}{K_{0}(R\sqrt{p/D})}a\sqrt{p/D}\right)\right)$
$\displaystyle\to$
$\displaystyle\exp\Biggl{(}-\ell\underbrace{\sqrt{p/D}\frac{K_{1}(R\sqrt{p/D})}{K_{0}(R\sqrt{p/D})}}_{=\mu_{0}^{(p)}}\Biggr{)}\qquad(a\to
0).$

As a consequence, the inverse Laplace transform of the right-hand side with
respect to $p$ yields $U(\ell,t|R)$ in Eq. (52).

##  Appendix D Concentric annulus

In order to quantify the relative contributions of far-reaching trajectories,
it is convenient to introduce an absorbing circle of radius $L$. In other
words, we are interested in diffusion inside a circular annulus,
$\Omega=\\{\bm{x}\in{\mathbb{R}}^{2}~{}:~{}R<|\bm{x}|<L\\}$. In this Appendix,
we adopt the spectral approach from [23, 18, 25, 24] to deduce the
distribution of the boundary local time on the inner circle with the
constraint of not hitting the outer circle. For brevity, we skip details and
only sketch the main formulas.

Due to rotational invariance, the eigenfunctions of the Dirichlet-to-Neumann
operator ${\mathcal{M}}_{p}$ associated to the inner circle are the Fourier
harmonics,

$v_{n}^{(p)}=\frac{e^{in\varphi}}{\sqrt{2\pi R}}\qquad(n\in{\mathbb{Z}}),$
(54)

whereas the associated eigenvalues are

$\mu_{n,L}^{(p)}=-\alpha\frac{K_{n}(\alpha L)I^{\prime}_{n}(\alpha
R)-I_{n}(\alpha L)K^{\prime}_{n}(\alpha R)}{K_{n}(\alpha L)I_{n}(\alpha
R)-I_{n}(\alpha L)K_{n}(\alpha R)}\,,$ (55)

with $\alpha=\sqrt{p/D}$ and the prime denoting the derivative with respect to
the argument (see similar computations in [23, 24]). The decomposition (26) of
the first-crossing time ${\mathcal{T}}_{\ell}$ implies that its probability
density is obtained by convolution of two densities, which in the Laplace
domain reads

$U_{L}(\ell,t|r_{0})={\mathcal{L}}_{p,t}^{-1}\biggl{\\{}\tilde{H}_{L}(p|r_{0})\,\exp(-\mu_{0,L}^{(p)}\ell)\biggr{\\}}.$
(56)

Here

$\tilde{H}_{L}(p|r_{0})=\frac{K_{0}(\alpha L)I_{0}(\alpha r_{0})-I_{0}(\alpha
L)K_{0}(\alpha r_{0})}{K_{0}(\alpha L)I_{0}(\alpha R)-I_{0}(\alpha
L)K_{0}(\alpha R)}$ (57)

is the Laplace transform of the probability density of the first-passage time
to the inner circle in the presence of an outer absorbing circle (expectedly,
this function vanishes on the outer circle, $r_{0}=L$, and is equal to $1$ on
the inner circle, $r_{0}=R$). In turn, $\exp(-\mu_{0,L}^{(p)}\ell)$ is the
Laplace transform of the probability density $U(\ell,t|R)$ (compare with Eq.
(52)). In the limit $L\to\infty$, the eigenvalue $\mu_{0,L}^{(p)}$ from Eq.
(55) tends to $\mu_{0}^{(p)}$ from Eq. (50), and one retrieves Eq. (52).
Similarly, one gets

$Q_{L}(\ell,t|r_{0})={\mathcal{L}}_{p,t}^{-1}\left\\{\tilde{H}_{L}(p|r_{0})\,\frac{\exp(-\mu_{0,L}^{(p)}\ell)}{p}\right\\}$
(58)

and

$\rho_{L}(\ell,t|r_{0})={\mathcal{L}}_{p,t}^{-1}\left\\{\tilde{H}_{L}(p|r_{0})\,\frac{\mu_{0,L}^{(p)}}{p}\exp(-\mu_{0,L}^{(p)}\ell)\right\\}.$
(59)

To get the long-time behavior, one can evaluate the above expressions in the
limit $p\to 0$; in particular,

$\lim\limits_{p\to
0}\tilde{H}_{L}(p|r_{0})=\frac{\ln(L/r_{0})}{\ln(L/R)}\,,\qquad\lim\limits_{p\to
0}\mu_{0}^{(p)}=\frac{1}{R\ln(L/R)}\,,$ (60)

from which

$Q_{L}(\ell,\infty|R)=\frac{\ln(L/r_{0})}{\ln(L/R)}\exp\biggl{(}-\frac{\ell/R}{\ln(L/R)}\biggr{)}.$
(61)

This is the probability that the threshold $\ell$ is crossed by the boundary
local time $\ell_{t}$ before the trajectory $\bm{X}_{t}$ hit the outer circle
of radius $L$. As $L\to\infty$, this probability tends to $1$, but this
approach is slow and controlled by $\ln(L/R)$.

## References

  * [1] Lévy P 1965 Processus Stochastiques et Mouvement Brownien (Paris: Gauthier-Villard).
  * [2] Ito K and McKean HP 1965 Diffusion Processes and Their Sample Paths (Berlin: Springer-Verlag).
  * [3] Freidlin M 1985 Functional Integration and Partial Differential Equations (Annals of Mathematics Studies, New Jersey: Princeton University Press, Princeton).
  * [4] Darling DA and Kac M 1957, Trans. Am. Math. Soc. 84 444-458.
  * [5] Ray D 1963 Illinois J. Math. 7 615-630.
  * [6] Knight FB 1963 Trans. Amer. Math. Soc. 109 56-86.
  * [7] Agmon N 1984 J. Chem. Phys. 81 3644-3647.
  * [8] Berezhkovskii AM, Zaloj V, and Agmon N 1998 Phys. Rev. E 57 3937-3947.
  * [9] Dhar A and Majumdar SN 1999 Phys. Rev. E 59 6413-6418.
  * [10] Yuste SB and Acedo L 2001 Phys. Rev. E 64 061107\.
  * [11] Godrèche C and Luck JM 2001, J. Stat. Phys. 104 489-524.
  * [12] Majumdar SN and Comtet A 2002 Phys. Rev. Lett. 89 060601\.
  * [13] Bénichou O, Coppey M, Klafter J, Moreau M, and Oshanin G 2003 J. Phys. A.: Math. Gen. 36 7225-7231.
  * [14] Condamin S, Bénichou O, and Moreau M 2005 Phys. Rev. E 72 016127\.
  * [15] Condamin S, Tejedor V, and Bénichou O 2007 Phys. Rev. E 76 050102R.
  * [16] Burov S and Barkai E 2007 Phys. Rev. Lett. 98 250601\.
  * [17] Burov S and Barkai E 2011 Phys. Rev. Lett. 107 170601\.
  * [18] Grebenkov DS 2020 Phys. Rev. Lett. 125 078102\.
  * [19] Grebenkov DS 2007 Phys. Rev. E 76 041139\.
  * [20] Saisho Y 1987 Probab. Theory Rel. Fields 74 455-477.
  * [21] Papanicolaou VG 1990 Probab. Theory Rel. Fields 87 27-77.
  * [22] Bass RF, Burdzy K, and Chen ZQ 2008 Proc. London Math. Soc. 96 273-311.
  * [23] Grebenkov DS 2019 Phys. Rev. E 100 062110\.
  * [24] Grebenkov DS 2020 Phys. Rev. E 102 032125\.
  * [25] Grebenkov DS 2020 J. Stat. Mech. 103205\.
  * [26] Redner S 2001 A Guide to First Passage Processes (Cambridge: Cambridge University press).
  * [27] Epstein CL and Schotland J 2008 SIAM Rev. 50 504-520.
  * [28] Carslaw HS and Jaeger JC 1959 Conduction of Heat in Solids, 2nd Ed. (Oxford: Clarendon).
  * [29] Grebenkov DS, Metzler R, and Oshanin G 2018 Phys. Chem. Chem. Phys. 20 16393-16401.
  * [30] Levitz PE, Zinsmeister M, Davidson P, Constantin D, and Poncelet O 2008, Phys. Rev. E 78 030102(R).
  * [31] Mörters P and Peres Y 2010 Brownian Motion (Cambridge Series in Statistical and Probabilistic Mathematics, Cambridge: Cambridge University Press).
  * [32] Smoluchowski M 1917 Z. Phys. Chem. 92U 129-168.
  * [33] Collins FC and Kimball GE 1949, J. Coll. Sci. 4 425-437.
  * [34] Sano H and Tachiya M 1979, J. Chem. Phys. 71 1276-1282.
  * [35] Shoup D and Szabo A 1982 Biophys. J. 40 33-39.
  * [36] Zwanzig R 1990 Proc. Natl. Acad. Sci. USA 87, 5856-5857.
  * [37] Sapoval B 1994 Phys. Rev. Lett. 73 3314-3317.
  * [38] Filoche M and Sapoval B 1999 Eur. Phys. J. B 9 755-763.
  * [39] Bénichou O, Moreau M, and Oshanin G 2000 Phys. Rev. E 61 3388-3406.
  * [40] Grebenkov DS, Filoche M, and Sapoval B 2003 Eur. Phys. J. B 36, 221-231.
  * [41] Berezhkovskii AM, Makhnovskii Y, Monine M, Zitserman V, and Shvartsman S 2004 J. Chem. Phys. 121 11390\.
  * [42] Grebenkov DS 2006 Partially Reflected Brownian Motion: A Stochastic Approach to Transport Phenomena, in “Focus on Probability Theory”, Ed. L. R. Velle, pp. 135-169 (New York: Nova Science Publishers).
  * [43] Grebenkov DS, Filoche M, and Sapoval B 2006, Phys. Rev. E 73 021103\.
  * [44] Bressloff PC, Earnshaw BA, and Ward MJ 2008, SIAM J. Appl. Math. 68 1223-1246.
  * [45] Reingruber J and Holcman D 2009 Phys. Rev. Lett. 103 148102\.
  * [46] Lawley SD and Keener JP 2015 SIAM J. Appl. Dyn. Sys. 14 1845-1867.
  * [47] Grebenkov DS and Oshanin G 2017 Phys. Chem. Chem. Phys. 19 2723-2739.
  * [48] Bernoff A, Lindsay A, and Schmidt D 2018 Multiscale Model. Simul. 16 1411-1447.
  * [49] Grebenkov DS 2019 J. Chem. Phys. 151 104108\.
  * [50] Tejedor V, Schad M, Bénichou O, Voituriez R, and Metzler R 2011 J. Phys. A.: Math. Theor. 44 395005\.
  * [51] Amitai A, Kupka I, and Holcman D 2012 Phys. Rev. Lett. 109 108302\.
  * [52] Bressloff PC and Newby J 2013 Rev. Mod. Phys. 85 135-196.
  * [53] Giuggioli L, Pérez-Becker S, and Sanders DP 2013, Phys. Rev. Lett. 110 058103\.
  * [54] Tzou JC, Xie S, and Kokolnikov T 2014, Phys. Rev. E 90 062138\.
  * [55] Agliari E, Blumen A, and Cassi D 2014 Phys. Rev. E 89 052147\.
  * [56] Agliari E, Cassi D, Cattivelli L, and Sartori F 2016, Phys. Rev. E 93 052111\.
  * [57] Peng J and Agliari E 2019 Phys. Rev. E 100 062310\.
  * [58] Le Vot F, Yuste SB, Abad E, and Grebenkov DS 2020 Phys. Rev. E 102 032118\.
  * [59] Arendt W, ter Elst AFM, Kennedy JB, and Sauter M 2014 J. Funct. Anal. 266 1757-1786.
  * [60] Daners D 2014 Positivity 18 235-256.
  * [61] Arendt W and ter Elst AFM 2015 Potential Anal. 43 313-340.
  * [62] Hassell A and Ivrii V 2017 J. Spectr. Theory 7 881-905.
  * [63] Girouard A and Polterovich I 2017 J. Spectr. Theory 7 321-359.


# Globally visible singularity in an astrophysical setup

Karim Mosani<EMAIL_ADDRESS>BITS Pilani K.K. Birla Goa Campus,
Sancoale, Goa-403726, India  Dipanjan Dey<EMAIL_ADDRESS>International Center for Cosmology, Charusat University, Anand 388421,
Gujarat, India  Pankaj S. Joshi<EMAIL_ADDRESS>International
Center for Cosmology, Charusat University, Anand 388421, Gujarat, India

###### Abstract

The global visibility of a singularity as an end state of the gravitational
collapse of a spherically symmetric pressureless cloud is investigated. We
show the existence of a non-zero measured set of parameters: the total mass
and the initial mean density of the collapsing cloud, giving rise to a
physically strong globally visible singularity as the end state for a fixed
velocity function. The existence of such a set indicates that such singularity
is stable under small perturbation in the initial data causing its existence.
This is true for marginally as well as non-marginally bound cases. The
possibility of the presence of such suitable parameters in the astrophysical
setup is then studied: $1)$ The singularities’ requirements at the center of
the M87 galaxy and at the center of our galaxy (SgrA*) to be globally visible
are discussed in terms of the initial size of the collapsing cloud forming
them, presuming that such singularities are formed due to gravitational
collapse. $2)$ The requirement for the primordial singularities formed due to
a collapsing configuration after getting detached from the background
universe, at the time of matter-dominated era just after the time of matter-
radiation equality, to be globally visible, is discussed. $3)$ The scenario of
the collapse of a neutron star after reaching a critical mass, which is
achieved by accreting the supernova ejecta expelled by its binary companion
core progenitor, is considered. The primary aim of this paper is to show that
globally visible singularities can form in astrophysical setups under
appropriate circumstances.

##  I Introduction

Observations about compact objects such as the shadow of the M87 galactic
center by the Event Horizon Telescope Akiyama , earlier observations such as
Sagittarius A* (Sgr A*) at our galactic center Schodel ; Ghez , and similar
compact objects at the center of other galaxies Kormendy  may hint that
singularities do exist in our universe. General relativity predicts that a
massive cloud under its gravitational influence collapses to form an
infinitely dense singularity, curving the spacetime infinitely Hawking .
Escape of timelike or null geodesics may or may not be possible from such
singularities. The weak cosmic censorship hypothesis, proposed by Penrose
Penrose , states that singularities cannot be visible to an asymptotic
observer. However, it can now be shown that a collapsing dust cloud can give a
globally visible singularity as its end state for certain mass profiles and
velocity profiles Deshingkar ; Mosani .

In addition to this, some theoretical predictions motivate us to consider the
existence of a naked singularity. One of them is that shadows are not unique
to black holes, but even naked singularities under certain circumstances can
cast a shadow Shaikh . Apart from this, the precession of the orbit of stars
around the compact object at the galactic center can also put light on the
causal structure of the central singularity. It has been found that while the
precession of particle orbit in the external Schwarzschild spacetime is in the
direction of the particle motion, the naked singularity spacetimes like the
Joshi-Malafarnia-Narayan (JMN) spacetime JMN  and the Janis-Newman-Winicour
(JNW) spacetime JNW  can cause the orbiting particle to precess in the
opposite direction of the particle’s motion under certain conditions
Bambhaniya ; AJoshi . The S2 star orbiting around the center of our galaxy may
also exhibit such behavior Dey . Some observational signatures to look into,
which help us in distinguishing from black holes and naked singularities, like
the accretion disk property, the Einstein ring formed due to the gravitational
lensing and the shadow, has been studied lately and can be found in Dey2 ;
Dey3 ; Shaikh2 ; Bambi ; Dey4 .

The existence of a singularity may represent the incomplete nature of the
classical theory of gravity Wheeler ; HawkingCUP ; Bergmann  and a quantum
version of gravity possibly resolves it. However, if a singularity exists, it
has its quantum aspects, which makes it more interesting to study. If such
singularity allows null geodesic to escape from it so that it reaches an
external observer without encountering trapped surfaces, then it could act as
an astrophysical lab to test the quantum theory of gravity, which is still not
understood JoshiCUP .

Mathematically, it is possible to achieve a globally naked singularity.
However, whether suitable configurations required for its existence are
possible astrophysically is our matter of concern, and we will discuss it in
this paper. The paper is arranged as follows: Section II consists of the
mathematical formalism of the gravitational collapse of a spherically
symmetric dust cloud and the collapse criteria to end in a globally visible
singularity. In section III, the astrophysically reasonable parameters, which
can give a globally visible singularity, are discussed. The collapse of the
cloud possibly forming the central singularity at the galactic center, the
collapse of the dark matter cloud in the primordial time, and the collapse of
an accreting neutron star are investigated. We end the paper with the
concluding remarks and open concerns in section IV.

##  II Gravitational collapse formalism

We consider the gravitational collapse of a spherically symmetric cloud of
dust. The spacetime of such dust cloud is described by Lemaitre-Tolman-Bondi
metric Lemaitre ; Tolman ; Bondi , which is given by

$ds^{2}=-c^{2}dt^{2}+\frac{R^{\prime}}{1+f}dr^{2}+R^{2}d\Omega^{2}.$ (1)

Here $R=R(t,r)$ is the physical radius of the collapsing cloud and $f(r)$ is
called the velocity function. The energy-momentum tensor in the comoving frame
is given by

$T^{\mu\nu}=\rho(t,r)U^{\mu}U^{\nu}.$ (2)

where $\rho(t,r)$ is the energy density of the cloud $U^{\mu}$ is the four
velocity. Using the Einstein’s field equations, we get the following relation:

$\frac{F^{\prime}}{R^{2}R^{\prime}}=8\pi\rho$ (3)

and

$\dot{F}=0,$ (4)

where

$F=\frac{c^{2}}{G}\left(\frac{R\dot{R}^{2}}{c^{2}}-fR\right).$ (5)

The superscripts dot and prime denotes the partial derivative with respect to
time and radial coordinate respectively. The term $F$ signifies the mass of
the cloud inside a shell of radial coordinate $r$ at time $t$. It is called
the Misner-Sharp mass function Misner . In case of a collapsing dust, the
Misner-Sharp mass function is a conserved quantity for a given radial
coordinate and hence is a function of only $r$ as seen in Eq.(4). Integrating
Eq.(5), gives

$t-t_{s}(r)=-\frac{R^{\frac{3}{2}}}{\sqrt{G}\sqrt{F}}\mathcal{G}\left(-\frac{c^{2}fR}{GF}\right),$
(6)

where $\mathcal{G}(y)$ is a positive, real and bound function given by

$\begin{split}&\mathcal{G}(y)=\left(\frac{\sin^{-1}\sqrt{y}}{y^{\frac{3}{2}}}-\frac{\sqrt{1-y}}{y}\right),\hskip
28.45274pt0<y\leq 1,\\\ &\mathcal{G}(y)=\frac{2}{3},\hskip
28.45274pty=0.\end{split}$ (7)

Now, we can rescale the physical radius using the coordinate freedom such that

$R(0,r)=r.$ (8)

This along with Eq.(6) gives

$t_{s}(r)=\frac{r^{\frac{3}{2}}}{\sqrt{G}\sqrt{F}}\mathcal{G}\left(-\frac{c^{2}fr}{GF}\right),$
(9)

where $t_{s}(r)$ is the time taken by a shell of radial coordinate $r$ to
collapse to a singularity. From Eq.(6), Eq.(9) and Eq.(7), one can write $R$
explicitly as a function of $r$ and $t$ as

$R(t,r)_{\textrm{NMB}}=\frac{5GF}{2c^{2}f}\bigg{(}1-\bigg{(}1-\frac{4fc^{2}}{5FG}\bigg{(}r^{\frac{3}{2}}-\frac{3c^{2}f}{10GF}r^{\frac{5}{2}}-\frac{3}{2}\sqrt{G}\sqrt{F}t\bigg{)}^{\frac{2}{3}}\bigg{)}^{\frac{1}{2}}\bigg{)},$
(10)

for $f\neq 0$. The subscript “NMB” stands for non-marginally bound collapse.
Here, we have used the Taylor expansion representation of $\mathcal{G}(y)$ up
to first order, which is given by

$\mathcal{G}(y)=\frac{2}{3}+\frac{y}{5}+o(y^{2}),\hskip 28.45274pt0<y<1.$ (11)

Eq.(10) is a good approximation only if the velocity function is such that
$|f|<<\frac{GF}{c^{2}r}$. Velocity functions with higher magnitude require
higher orders to be considered in the Taylor expansion of $\mathcal{G}(y)$.
For the case of marginally bound (MB) collapse, $R$ is expressed as

$R(t,r)_{\textrm{MB}}=\left(r^{\frac{3}{2}}-\frac{3}{2}\sqrt{G}\sqrt{F}t\right)^{\frac{2}{3}}.$
(12)

The condition for the apparent horizon curve can be obtained from the
following equality:

$g^{\mu\nu}R,_{\mu}R,_{\nu}=0.$ (13)

This gives the following relation needed to be satisfied for the apparent
horizon curve:

$\frac{\dot{R}^{2}}{c^{2}}=1+f.$ (14)

The above equation along with Eq.(5) gives

$F=\frac{c^{2}}{G}R.$ (15)

The event horizon evolves like a null geodesic wavefront and meets the
apparent horizon at the boundary of the collapsing cloud $r_{c}$. Hence, the
event horizon is described by the solution of the differential equation:

$\frac{dt}{dr}=\frac{1}{c}\frac{R^{\prime}}{\sqrt{1+f}},$ (16)

satisfying the condition:

$F(r_{c})=\frac{c^{2}}{G}R(t,r_{c}).$ (17)

We consider the Misner-Sharp mass function to be of the form

$F=F_{0}\left(\frac{r}{r_{c}}\right)^{3}+F_{3}\left(\frac{r}{r_{c}}\right)^{6}.$
(18)

According to Tipler Ellis ; Tipler , the volume element formed by three
independent Jacobi fields along timelike geodesics terminating at a physically
strong singularity should vanish. Using the sufficient criteria by Clarke and
Krolak Clarke  for a singularity to be Tipler strong, one can conclude that in
the case of a marginally bound collapse, the mass function as stated in
Eq.(18) maintains the strength of the singularity in the sense of Tipler Joshi
. Additionally, for $F_{3}<0$, such singularity is locally naked as well. In
the case of non-marginally bound collapse, the strength of the singularity and
its local causal property depends on the velocity function in addition to the
mass function Mosani  (Please refer to the Appendix for further discussion on
the local visibility of the singularity).

For the density to decrease as we move away from the center and smoothly
vanish at the boundary of the collapsing cloud ($\rho(0,r_{c})=0$), $F$ has to
be rewritten as

$F(r)=F_{0}\left(\frac{r}{r_{c}}\right)^{3}\left(1-\frac{1}{2}\left(\frac{r}{r_{c}}\right)^{3}\right).$
(19)

It should be noted that in the case of marginally bound pressureless
collapsing cloud with density vanishing smoothly at its boundary, the above
expression is the only possible mass function which corresponds to a Tipler
strong singularity as its end state (See the Appendix).

Substituting Eq.(10) and Eq.(19) in Eq.(17), we get

$t_{\textrm{EH}_{\textrm{NMB}}}(r_{c})=\frac{2\sqrt{2}}{3\sqrt{G}\sqrt{F_{0}}}r_{C}^{\frac{3}{2}}-\frac{2\sqrt{2}}{5}\frac{c^{2}f(r_{c})}{\left(GF_{0}\right)^{\frac{3}{2}}}-\frac{GF_{0}}{3c^{3}}\left(1-\frac{f(r_{c})}{5}\right)^{\frac{3}{2}}$
(20)

for $f\neq 0$. For $f=0$, substituting Eq.(12) and Eq.(19) in Eq.(17), we get

$t_{\textrm{EH}_{\textrm{MB}}}(r_{c})=\frac{2\sqrt{2}}{3\sqrt{F_{0}}\sqrt{G}}\left(r_{c}^{\frac{3}{2}}-\left(\frac{GF_{0}}{2c^{2}}\right)^{\frac{3}{2}}\right).$
(21)

Figure 1: Non-zero measured set of initial data giving rise to a globally
visible singularity (shaded region) and globally hidden singularity (unshaded
region) as an end state of a marginally bound collapse (a) and a non-
marginally bound collapse (b) of a pressureless inhomogeneous dust cloud
having the mass function as in Eq.(19). Here, $\bar{\rho_{0}}$ and
$\frac{F_{0}}{4}$ are the initial mean density and the total mass of the
collapsing cloud having units $\textrm{kg}/\textrm{m}^{3}$ and kg
respectively.

The event horizon curve $t_{EH}(r)$ is now the solution of the differential
Eq.(16) satisfying Eq.(20) in case of non-marginally bound collapse (Eq.(21)
in case of marginally bound case). At the time of formation of the central
singularity, the apparent horizon (the boundary of all trapped surfaces)
starts forming from the center. It evolves in the outward direction for an
inhomogeneous mass function as in Eq.(19) Mosani ; Joshi2 . Hence, any null
geodesic leaving the center after the formation of the singularity, due to the
collapse of the central shell, gets trapped by trapped surfaces. Therefore,
the event horizon forms either before or during the formation of the central
singularity. If it forms at the center before the formation of this central
singularity, any null geodesic escaping from the center after $t_{EH}(0)$,
having a positive tangent at $r=0$, will fall back to the singularity, thereby
making the singularity globally hidden. Hence, the necessary criteria for a
singularity to be globally visible is as follows:

$t_{s}(0)=t_{\textrm{EH}}(0).$ (22)

We note here that the singularity has to be a nodal point apart from
satisfying the above equality. This is to make sure that an entire family of
null geodesic escapes from the point $(t_{s}(0),0)$ in the $(t,R)$ plane so
that the central singularity remains visible to an asymptotic observer for an
infinite time Joshi2 . These escaping null geodesics are solutions of the
differential equation (16) starting from the point $(t_{s}(0),0)$ in the
$(t,R)$ plane, and reaching the boundary $r_{c}$ before the event horizon, for
which it has to satisfy the following inequality:

$F(r_{c})<\frac{c^{2}}{G}R(t,r_{c}).$ (23)

. We now use Eq.(22) to plot Fig.(1) on the basis of which we study the
astrophysical relevance of globally visible singularity in the following
section.

##  III Astrophysical relevance

It is evident that for a pressureless marginally bound collapsing cloud, with
smoothly vanishing density at its boundary, ending up in a Tipler strong
singularity, the global visibility or otherwise of the singularity depends on
two parameters: the total mass $\frac{F_{0}}{4}$ and the initial size $r_{c}$
of the collapsing cloud at the beginning of its collapse (or equivalently
$\frac{F_{0}}{4}$ and its initial mean density
$\bar{\rho_{0}}=\frac{3F_{0}}{8\pi r_{c}^{3}}$). This is because, in the case,
$f=0$, the equality or otherwise of the singularity curve Eq.(9) and the event
horizon curve (which is a solution of the differential Eq.(16) satisfying the
condition Eq.(17)), which decides the global visibility of the singularity,
depends on only $F_{0}$ and $r_{c}$. A similar situation may or may not occur
in the case of a non-marginally bound collapse. It is possible to show some
examples of the mass profile, values of $F_{0}$, and $\bar{\rho_{0}}$, which
gives rise to the globally visible singularity as an end state. However, do
such values arise in nature? Fig. (1) is a plot of $\textrm{log}_{10}F_{0}$
v/s $\textrm{log}_{10}\bar{\rho_{0}}$ (for marginally bound case as well as
non marginally bound case) in which the dotted line is the cut off for
$\bar{\rho_{0}}$ for a given $F_{0}$ having the property that any
$\bar{\rho_{0}}$ less than this cut off gives a globally visible singularity
for a given $F_{0}$. In the case of non-marginally bound collapse (Fig(1.b)),
the velocity function is taken such that mass and size remain the only
variables that control the global causal property of the singularity. Also, it
is small enough to avoid error due to the approximation of Taylor expansion of
Eq.(11) up to only the first order. Here, we study three types of
astrophysical scenarios giving rise to a singularity and investigate if they
can be observed by an asymptotic observer.

Figure 2: Causal structure of space-time singularity formed at the center of
the milky way galaxy (the compact object SgrA*). The initial unknown radius of
the collapsing core forming such singularity decides its end state. $r_{c}\leq
3.32\times 10^{11}\textrm{m}$ and $r_{c}\geq 3.33\times 10^{11}\textrm{m}$
ends up globally hidden and a globally visible singularity respectively.
Collapsing core is assumed to be spherical, pressureless and marginally bound.
Mass function is chosen such that the singularity is Tipler strong.

Figure 3: Causal structure of space-time singularity formed at the center of
the M87 galaxy. The initial unknown radius of the collapsing core forming such
singularity decides its end state. $r_{c}\leq 5.13\times 10^{14}\textrm{m}$
and $r_{c}\geq 5.14\times 10^{14}\textrm{m}$ ends up in a globally hidden and
a globally visible singularity respectively. Collapsing core is assumed to be
spherical, pressureless and marginally bound. Mass function is chosen such
that the singularity is Tipler strong.

###  III.1 The galactic center

We investigate the causal structure of the central singularity of SgrA*, which
is a compact object at the center of our galaxy. The mass of this object as
calculated by Gillessen et. al. Gillessen  is $4.31\times 10^{6}M_{\odot}$,
and its Schwarzschild radius is $\sim 1.26\times 10^{10}$ m. For the SgrA*
singularity to be globally visible, the collapsing spherical cloud forming
such singularity should have the initial radius greater than equal to $\sim
3.33\times 10^{11}$ m (in case of a marginally bound case), which translates
to roughly at least $26.42$ times the size of its Schwarzschild radius. A
lesser initial radius will make the singularity globally hidden. This
corresponds to having the initial mean density less than
$55.42\textrm{kg}/\textrm{m}^{3}$. Fig.(2) is the visual representation of
both the possible outcomes.

Similarly, the compact object at the center of the M87 galaxy has mass
$6.6\times 10^{9}M_{\odot}$ as calculated by the Event Horizon Telescope
collaboration Akiyama2 , and a Schwarzschild radius $\sim 1.9\times 10^{13}$
m. For the central singularity to be globally visible, the collapsing
spherical cloud forming such singularity should have the initial radius
greater than or equal to $\sim 5.14\times 10^{14}$ m (in case of marginally
bound collapse), which roughly translates to at least $27.05$ times its
Schwarzschild radius. A lesser initial radius will make the singularity
globally hidden. This corresponds to having the initial mean density less than
$2.31\times 10^{-5}\textrm{kg}/\textrm{m}^{3}$. Fig.(3) is the visual
representation of both the possible outcomes.

###  III.2 Singularities having primordial origin

Another astrophysically motivated entity that we consider here is the
primordial singularities. Primordial fluctuations in the very early universe
could give rise to such singularities Zeldovich ; Hawking2 . For a primordial
singularity to form due to gravitational collapse, we investigate the collapse
of matter in the era just after the matter-radiation equality, i.e., the start
of the matter-dominated era. The reason for considering this time as the
initiation of the collapse is that the temperature during this epoch drops
down to $9000$ K, and the dark matter configuration can now collapse. We know
that the epoch of matter-radiation equality occurred at redshift $z\sim 3000$
Frieman . For a $\Lambda\textrm{CDM}$ model, the deceleration to acceleration
transition of the universe is calculated to begin at $z\sim 0.72$ Farooq . The
relation between the density of the universe at these epochs depends on $z$ at
both these epochs as follows:

$\frac{\rho_{\textrm{MRE}}}{\rho_{\textrm{DA}}}=\left(\frac{1+z_{\textrm{MRE}}}{1+z_{\textrm{DA}}}\right)^{3}.$
(24)

This is because the era bounded below and above by the time of matter-
radiation equality and the time of deceleration to acceleration transition,
respectively, has matter as the dominant fluid, and the density of the
universe is dependent on the linear equation of state parameter $\omega$ of
the dominating fluid, as follows Liddle :

$\rho(\omega)\propto\left(1+z\right)^{3(1+\omega)}.$ (25)

Hence for $\omega=0$ (dust domination), Eq.(24) holds. If we consider the dark
energy candidate driving the present accelerated expansion as the cosmological
constant ($\omega=-1$), its density $\rho_{\Lambda}$ remains constant in time.
Therefore, we can say that

$\rho_{\textrm{DA}}=\rho_{\Lambda}.$ (26)

Using the above equality in Eq.(24), one can obtain the density of the
universe at the time of matter-radiation equality as

$\rho_{\textrm{MRE}}=3.092\times 10^{-17}\textrm{kg}/\textrm{m}^{3}.$ (27)

Here, we have used the present Hubble parameter, $H_{0}=67.4\textrm{km
s}^{-1}\textrm{Mpc}^{-1}$, and the density parameter of $\Lambda$,
$\Omega_{\Lambda}=0.686$, from the Planck 2018 results Planck  to get the
value of $\rho_{\Lambda}=5.82298\times 10^{-27}\textrm{kg}/\textrm{m}^{3}$
using the relation $\rho_{\Lambda}=\frac{3H_{0}^{2}\Omega_{\Lambda}}{8\pi G}$
Liddle . The overdensity configuration, formed due to primordial fluctuation
just after the epoch of matter-radiation-equality, detaches itself from the
background universe and starts contracting when the density contrast
$\delta\rho/\rho$ reaches order one Khlopov . Hence, at the start of the
collapse, one can assume the order of the initial mean density of the
collapsing configuration, which we denote by $\rho_{\textrm{CONFIG}}$, same as
the order of the density of the background universe, i.e.
$10^{-17}\,\textrm{kg}/\textrm{m}^{3}$. Therefore, we have

$\rho_{\textrm{CONFIG}}<\bar{\rho_{0}}\left(\frac{F_{0}}{4}\right)$ (28)

for all configuration having $\frac{F_{0}}{4}<10^{45}\,\textrm{kg}$ as is
evident from Fig.(1) in both marginally bound as well as non-marginally bound
case for very small velocity function. Hence, we can state that any
configuration with $\frac{\delta\rho}{\rho}$ having order one, and having mass
less than $10^{45}\,\textrm{kg}$ collapses to give a globally visible
singularity if the collapse initiates just after the time $t_{\textrm{MRE}}$.
It is worth mentioning that the overdense region collapses even if the density
contrast is of the order greater than one, in which case, the result may vary
accordingly.

###  III.3 Collapse of a matter accreting neutron star

We now consider a close binary system scenario, one of which is a star
undergoing supernova explosion and another is a neutron star, as investigated
by Rueda and Ruffini Reuda . We give a brief overview of this scenario.

The material expelled from the core progenitor is accreted by the companion
neutron star and can get captured if it falls in the region at a distance less
than

$R_{cap}=\frac{2GM_{NS}}{v_{rel}^{2}}$ (29)

from the center of the neutron star. Here, $M_{NS}$ is the mass of the neutron
star, and $v_{rel}$ is the velocity of the ejected particle relative to the
neutron star orbital motion, and is given by

$v_{rel}=\sqrt{v^{2}_{orb}+v^{2}_{ej}},$ (30)

where

$v_{orb}=\sqrt{\frac{G(M_{prog}+M_{NS})}{a}}.$ (31)

Here, $v_{ej}$ is the ejecta velocity, which at the start of the supernova
explosion has a value of the order of $10^{7}$ $ms^{-1}$, $v_{orb}$ is the
orbit velocity of the neutron star, $M_{prog}$ is the mass of the progenitor,
and $a$ is the binary separation. We neglect the neutron star magnetic field’s
effect by considering the magnetospheric radius $R_{m}$ comparatively smaller
than $R_{cap}$. This is possible for a suitably high rate of accretion since

$R_{m}=\frac{B^{2}R_{NS}^{6}}{(\dot{M}\sqrt{2GM_{NS}})^{\frac{2}{7}}},$ (32)

where $\dot{M}=\frac{dM}{dt}$ is the matter accretion rate of the supernova
ejecta by the neutron star, and $B$, $M_{NS}$ and $R_{NS}$ are the magnetic
field, mass and the radius of the neutron star respectively Toropina . The
mass accretion rate is given by Hoyle

$\dot{M}=\epsilon\pi\rho_{ej}R^{2}_{cap}v_{rel}=\epsilon\pi\rho_{ej}\frac{(2GM_{NS})^{2}}{v^{3}_{rel}}.$
(33)

Here, $0<\epsilon<1$, and depends on the medium in which the accretion process
takes place. $\rho_{ej}(t)=\frac{3M_{ej}(t)}{4\pi r_{ej}^{3}(t)}$ is the
density of the ejected material which decreases with time. The supernova
ejecta radius can be assumed to expand as $r_{ej}=at^{b}$, where $a$ and $b$
are constants Chevalier . Eq.(33) can now be integrated to get the magnitude
of the mass accreted by the neutron star. After some time, provided the
accretion rate is high,in order that the effect of the magnetic field is
negligible, the mass of the neutron star reaches the critical mass and further
collapses unhindered. Neutron star stable configuations have been studied
taking into account weak, strong, electromagnetic and gravitational
interaction in the framework of general relativity Belvedere . Tab.(1)
mentions the critical mass $M_{crit}$ and its corresponding radii for some
parameterizations of neutron star models studied in Belvedere . If we model
the unhindered collapse of the neutron star after achieving the critical
masses mentioned in Tab.1, by a spherical pressureless marginally bound cloud
collapse, having the mass function of the form of Eq.(19), then each model
(i.e. $NL3$ Lalazissis , $NL-SH$Sharma , $TM1$ Sugahara , and $TM2$ Hirata )
ends up in a singularity which is hidden globally. This is because the event
horizon forms before the singularity at $r=0$ in each case as seen using
equations (9), (16) and (17). However, the singularity is visible locally, as
explained in the appendix.

Model |  $M_{crit}$ ($\times M_{\odot}$) |  $R$ ($km$) |  $t_{s}(0)$ $(\times 10^{-5}s)$ |  $t_{EH}(0)$ $(\times 10^{-5}s)$
---|---|---|---|---
NL3 | $2.67$ | $12.33$ | $3.42937$ | $0.83442$
NL-SH | $2.68$ | $12.54$ | $3.51078$ | $0.91288$
TM1 | $2.58$ | $12.31$ | $3.48018$ | $0.98919$
TM2 | $2.82$ | $13.28$ | $3.72990$ | $1.00047$
Table 1: Neutron star critical masses and the corresponding radii are obtained
by fixing certain nuclear parameters, which includes the coupling constants
and the meson masses Belvedere . Each of these four models corresponds to one
such fixed-parameter set. For each model, in the case of marginally bound
unhindered gravitational collapse, the event horizon forms before the
formation of the singularity at $r=0$. Hence the singularity thus formed is
hidden globally. Figure 4: Neutron star critical masses and the corresponding
radii are obtained by fixing certain nuclear parameters which includes the
coupling constants and the meson masses Belvedere . Each of the four models:
$NL3$, $NL-SH$, $TM1$, and $TM2$ corresponds to one such fixed parameter set.
The Misner-Sharp mass function is chosen as
$F(r,R)=F_{0}\left(\frac{r}{r_{c}}\right)^{3}+F_{3}\left(\frac{r}{r_{c}}\right)^{6}+F_{R}\left(\frac{R}{r_{c}}\right)^{3}$.
We have considered a marginally bound collapse, i.e. $\mathcal{Y}=1$. The
total initial mass, i.e. the critical mass of the collapsing cloud is given by
$M_{crit}=(F_{0}+F_{3}+F_{R})/2$. The magnitudes of $F_{0}$ and $F_{3}$
depicted in the plot are considerably smaller than that of $F_{R}$. Hence,
$M_{crit}$ is approximately $F_{R}/2$. The region of parameter space
$(F_{0},F_{3})$ below (above) the curve gives a globally visible (hidden)
singularity at the end of the collapse of a neutron star after reaching the
critical mass by accreting the supernova ejecta from its binary companion core
progenitor.

One could argue that modeling the cloud formed due to matter accreting neutron
star by zero pressure is not justified since neutron star has non-zero
pressure, following a certain equation of state. To investigate the end state
of the collapsing cloud formed after the neutron star accretes the supernova
ejecta from its binary companion core progenitor, we rely on the work by
Giambo et al. Giambo ; Giambo2 . We give a brief overview of the formalism
presented in these articles.

The general spacetime metric of a spherical collapsing cloud in the comoving
coordinates is given by

$ds^{2}=-c^{2}\frac{\dot{R}^{2}}{\mathcal{H}}dt^{2}+\frac{R^{\prime
2}}{\mathcal{Y}}dr^{2}+R^{2}d\Omega^{2},$ (34)

where $R$, $\mathcal{Y}$ and $\mathcal{H}$ are functions of $t$ and $r$. The
above metric can be rewritten in the transformed comoving area-radial
coordinates $(r,R)$ as

$ds^{2}=-Adr^{2}-2BdRdr-\frac{1}{\mathcal{H}}dR^{2}+R^{2}d\Omega^{2},$ (35)

where $A$ and $B$ are functions of $r$ and $R$. $\mathcal{H}$ can also be
written as a function of $r$ and $R$ as

$\mathcal{H}(r,R)=\frac{F(r,R)}{R}+\mathcal{Y}(r,R)^{2}-1.$ (36)

The transformed metric is advantageous to obtain information about the global
causal structure of the singularity. A quantity $\Delta(r,R)$ is defined as

$\Delta=B^{2}-\frac{A}{\mathcal{H}}=\frac{R^{\prime
2}}{\mathcal{Y}\mathcal{H}}.$ (37)

Using the Einsteins’ field equation, one can get an integral expression of
$\sqrt{\Delta}$ as

$\sqrt{\Delta}=\int_{R}^{r}\frac{1}{\sqrt{\mathcal{Y}(r,R)}}\frac{\partial}{\partial
r}\left(\frac{1}{\mathcal{H}(r,R)}\right)dR+\frac{1}{\sqrt{\mathcal{Y}(r,r)\mathcal{H}(r,r)}}.$
(38)

The formalism is restricted to only those cases where the energy density
$\rho$ and matter density $\rho_{m}$ are related as

$\rho=\omega\rho_{m}c^{2},$ (39)

where

$\omega=\frac{E(r)}{2}\left(\frac{F_{r}(r,R)}{\sqrt{\mathcal{Y}(r,R)}}+\frac{F_{R}(r,R)R^{\prime}(r,R)}{\sqrt{\mathcal{Y}(r,R)}}\right).$
(40)

The energy density, the radial pressure, and the tangential pressure are
respectively given by

$\rho=\frac{c^{2}}{8\pi R^{2}}\left(\frac{F,_{r}}{R^{\prime}}+F,_{R}\right),$
(41) $p_{r}=-\frac{c^{2}F,_{R}}{8\pi R^{2}},$ (42)

and

$p_{t}=-\frac{c^{2}}{16\pi
RR^{\prime}}\left(F,_{rR}-\sqrt{\mathcal{Y}}F,_{r}\frac{\partial}{\partial
R}\left(\sqrt{\mathcal{Y}}\right)+F,_{RR}R^{\prime}\right).$ (43)

In a collapsing cloud, provided: $(1)$ the weak energy condition hold, $(2)$
the regularity conditions hold, $(3)$ the shell crossing singularity do not
occur, and $(4)$ the shell focusing singularity form in a finite comoving
time, the singularity formed as an end state of the gravitational collapse is
globally visible if $n=3$ and $\frac{\zeta}{\alpha}>\frac{26+15\sqrt{3}}{2}$,
where $\zeta$ and $n$ arise in the Taylor expansion of $\sqrt{\Delta}(r,0)$
around $r=0$ as

$\sqrt{\Delta}(r,0)=\zeta r^{n-1}+.....,$ (44)

and $2\alpha$ is the coefficient of $r^{3}$ in the expression of the Misner-
Sharp mass function $F$ in the $(r,R)$ coordinates. One can choose a suitable
$\mathcal{Y}$ and $F$, due to two degrees of freedom available to us, such
that the end state of a non-zero pressured collapsing cloud, satisfying the
weak energy condition and the regularity condition, is globally visible. Let
us consider the Misner-Sharp mass function of the form

$F(r,R)=F_{0}\left(\frac{r}{r_{c}}\right)^{3}+F_{3}\left(\frac{r}{r_{c}}\right)^{6}+F_{R}\left(\frac{R}{r_{c}}\right)^{3}.$
(45)

We set $\mathcal{Y}=1$, which corresponds to an acceleration-free marginally
bound case. The total mass of the collapsing cloud is $(F_{0}+F_{3}+F_{R})/2$.

It should be noted that for the collapsing fluid to satisfy the weak energy
condition, one of the necessary criteria is to have the mass function to be a
monotonically increasing function of the areal coordinate $R$ in the $(r,R)$
frame. Hence, one realizes that the formalism discussed in Giambo ; Giambo2
is restricted to studying the end state of only those collapsing fluid with
negative radial pressure, as seen in Eq.(42). It can be seen that for suitable
values taken by the components of the Misner-Sharp mass function, one gets a
globally visible end-state singularity.

##  IV Concluding Remarks

Some concluding points and concerns are discussed below:

  1. 1.

A non-zero measured set of initial parameters, which relate to the total mass
and the initial mean density of the collapsing cloud, is obtained, which leads
to a globally visible singularity as its end state. The Misner-Sharp mass
function is chosen such that the singularity formed is physically strong in
the sense of Tipler. Hence, it can be concluded that such globally visible
singularity is stable under small perturbations in the initial data
($\frac{F_{0}}{4}$ and $\bar{\rho_{0}}$) of the collapsing spherical cloud.

  2. 2.

Globally visible singularity is possible to achieve in principle; however,
here, we have discussed the possibility of its existence in an astrophysical
scenario by looking into the possibility of occurrence of suitable mass and
size required for its existence. The singularity at the center of the M87
galaxy may be globally visible if the initial radius of the collapsing cloud
forming it is more than $27.05$ times its Schwarzschild radius. Similarly, the
collapsing cloud’s initial radius forming the SgrA* singularity should be more
than $26.42$ times its Schwarzschild radius to end up being globally visible.
Here, we have modeled the collapsing cloud to be marginally bound, spherically
symmetric, and pressureless. The outcome of the causal structure of the
singularity may vary if these assumptions are dropped.

  3. 3.

If a density contrast $\frac{\delta\rho}{\rho}$ of order one is achieved in
the matter-dominated era just after the time of matter-radiation equality, the
configuration detaches from the background universe and starts contracting due
to gravity. A configuration having a total mass of the order less than
$10^{45}$ kg forms a physically strong (in the sense of Tipler) singularity,
which can be visible globally. It should be noted that we have not taken into
account the process of virialization, which is an essential phenomenon in
cosmology, and which could oppose the unhindered gravitational collapse of the
configuration during the primordial time, causing it to attain an equilibrium,
thereby avoiding the formation of the singularity at all.

  4. 4.

A neutron star, which is a part of the binary system, can reach a critical
mass by accreting the supernova ejecta of the companion exploding star, after
which it can collapse unhindered. We consider four such models having
different critical masses and their corresponding radii. Each of these models
is obtained by fixing six nuclear parameters, as discussed in Belvedere . The
final state of such accreting neutron star is found to be a singularity that
is globally hidden for all these models. One may claim that these act as
positive evidence for the hypothesis of weak cosmic censorship. However, a
loophole in this argument, which is of concern, is that the collapse formalism
used for our purpose assumes zero pressure; however, neutron stars may have
non-zero pressure satisfying the polytropic equation of state. To address this
problem, we have incorporated the formalism proposed by Giambo et. al. Giambo
; Giambo2 . However, the formalism works only for negative radial pressure if
one has to satisfy the weak energy condition. One could obtain a suitable mass
function $F(r,R)$ and the function $\mathcal{Y}(r,R)$ for the collapse of
supernova ejecta accreting neutron star after achieving the critical mass and
obtain a globally naked singularity as the end state of the collapse, as seen
in Fig.(4). This end state is stable under small perturbation in the initial
data (here $F_{0}$ and $F_{3}$). The formalism developed in Giambo ; Giambo2
works only for the subset of configurations that have negative radial
pressure. The formalism of finding out the global causal structure of the
singularity formed due to gravitational collapse is not developed enough to
incorporate arbitrary pressures and equation of state as demanded by the
neutron star’s complicated structure. A simplified model giving rise to both
the possibilities: hidden and visible, globally, indicates that even a more
realistic model may give rise to a globally visible singularity.

  5. 5.

An external observer encountering an escaping singular null geodesic may find
traces of quantum gravity encoded in it, making such singularities more
tempting to investigate. We have shown here that the formation of a globally
visible singularity may also arise in the astrophysical setup and is not
merely a mathematical artifact. One could therefore interpret that such
singularities, after all, may not be so elusive.

##  V Acknowledgement

KM would like to acknowledge the support of the Council of Scientific and
Industrial Research (CSIR, India, Ref: 09/919(0031)/2017-EMR-1) for funding
the work.

##  VI Appendix

Here, we discuss the local visibility of the singularity formed due to a
marginally as well as non- marginally bound collapse of the dust sphere, The
Taylor expansion of the time curve around $r=0$ is expressed as

$t(r,v)=t(0,v)+r\chi_{1}(v)+r^{2}\chi_{2}(v)+r^{3}\chi_{3}(v)+O(r^{4}),$ (46)

where

$\chi_{i}(v)=\frac{1}{i!}\frac{d^{i}t}{dr^{i}}\bigg{|}_{r=0}.$ (47)

Here $v$ is called the scaling function and is related to the physical radius
and the comoving radius as

$R(t,r)=rv(t,r).$ (48)

Coordinate freedom gives us the option of rescaling the physical radius as

$R(t_{i},r)=r,$ (49)

where $t_{i}$ is the initial time of the collapse. This is equivalent to
having $v(t_{i},r)=1$. First three $\chi_{i}$’s can be calculated by
integrating Eq.(5), and is expressed as Joshi3 ; Mosani2 :

$\chi_{1}(v)=-\frac{1}{2}\int^{1}_{v}\frac{\frac{G\mathcal{F}_{1}}{v}+f_{1}c^{2}}{\left(\frac{G\mathcal{F}_{0}}{v}+f_{0}c^{2}\right)^{\frac{3}{2}}}dv,$
(50)
$\chi_{2}(v)=\int^{1}_{v}\left[\frac{3}{8}\frac{\left(G\frac{\mathcal{F}_{1}}{v}+f_{1}c^{2}\right)^{2}}{\left(G\frac{\mathcal{F}_{0}}{v}+f_{0}c^{2}\right)^{\frac{5}{2}}}-\frac{1}{2}\frac{\frac{G\mathcal{F}_{2}}{v}+f_{2}c^{2}}{\left(\frac{G\mathcal{F}_{0}}{v}+f_{0}c^{2}\right)^{\frac{3}{2}}}\right]dv$
(51)

and

$\chi_{3}(v)=\int_{v}^{1}\frac{f_{1}}{\left(\frac{G\mathcal{F}_{0}}{v}+f_{0}\right)^{\frac{3}{2}}}\left(-\frac{5}{16}\left(\frac{f_{1}}{\frac{G\mathcal{F}_{0}}{v}+f_{0}c^{2}}\right)^{2}+\frac{3}{4}\left(\frac{\frac{G\mathcal{F}_{2}}{v}+f_{2}c^{2}}{\frac{G\mathcal{F}_{0}}{v}+f_{0}c^{2}}\right)\right)-\frac{1}{2}\frac{\left(\frac{G\mathcal{F}_{3}}{v}+f_{3}c^{2}\right)}{\left(G\frac{\mathcal{F}_{0}}{v}+f_{0}c^{2}\right)^{\frac{3}{2}}}dv.$
(52)

Here, $f_{i}$ and $\mathcal{F}_{i}$ are such that

$F=\sum_{i=0}^{\infty}\mathcal{F}_{i}r^{i+3}$ (53)

and

$f=\sum_{j=0}^{\infty}f_{j}r^{j+2}.$ (54)

near $r=0$. It was shown in Mosani  that for a strong singularity to be
locally naked, we should have $\chi_{1}=\chi_{2}=0$ and $\chi_{3}>0$ at $v=0$.
This condition is always satisfied for the case of marginally bound collapse
with the collapsing fluid having mass function as in Eq.(18) with $F_{3}<0$,
as is apparent from Eq.(50-52). The condition may or may not satisfy for a non
zero velocity function. However, for the velocity function
$f=10^{-5}\frac{GF}{c^{2}r}$, as chosen in Fig.(1b), the above mentioned
condition is satisfied. Hence, all the singularities arising in this article
are at least locally naked.

## References

  * (1) K. Akiyama et al. [Event Horizon Telescope Collabora- tion], Astrophys. J. 875, 1 (2019).
  * (2) R. Schodel et. al., Nature, 419, 694 (2002).
  * (3) A. M. Ghez, et. al., Astrophys. J., 689, 1044 (2008).
  * (4) J. Kormendy and L. C. Ho, ARAA 51, 511 (2013).
  * (5) S. W. Hawking and G. F. R. Ellis, The large scale structure of space-time, Cambridge University Press (1973).
  * (6) R. Penrose, Riv. Nuovo Cimento Soc. Ital. Fis. 1, 252 (1969). S.W. Hawking and G. F. R. Ellis, The Large Scale Structure of Spacetime (Cambridge University Press, Cambridge, United Kingdom, 1973).
  * (7) S. S. Deshingkar, S. Jhingan and P. S. Joshi, Gen. Relativ. Gravit. 30, 1477 (1998).
  * (8) K. Mosani, D. Dey and P. S. Joshi, arXiv:2003.07092v1
  * (9) R. Shaikh, P. Kocherlakota, R. Narayan, and P. S. Joshi, Mon. Not. Roy. Astron. Soc. 482, 52 (2019).
  * (10) P. S. Joshi, D. Malafarina and R. Narayan, Classical Quantum Gravity, 28, 23 (2011).
  * (11) A. I. Janis, E. T. Newman, and J. Winicour  Phys. Rev. Lett. 20, 878 (1968).
  * (12) P. Bambhaniya, A. B. Joshi, D. Dey, and P. S. Joshi, Phys. Rev. D 100, 124020 (2019).
  * (13) A. B. Joshi, P. Bambhaniya, D. Dey, and P. S. Joshi, arXiv: 1909.08873.
  * (14) D. Dey, P. S. Joshi, A. Joshi, and P. Bambhaniya, Int. J. Mod. Phys. D 28, no. 14, 1930024 (2019).
  * (15) D. Dey, K. Bhattacharya, and T. Sarkar, Phys. Rev. D 87, 103505 (2013).
  * (16) D. Dey, K. Bhattacharya, and T. Sarkar,  Phys. Rev. D 88, 083532 (2013).
  * (17) R. Shaikh and P. S. Joshi, Journal of Cosmology and Astroparticle Physics 10, 064 (2019).
  * (18) A. B. Abdikamalov, A. A. Abdujabbarov, D. Ayzenberg, D. Malafarina, C. Bambi, and B. Ahmedov,  Phys. Rev. D 100, 024014 (2019).
  * (19) K. Bhattacharya, D. Dey, A. Mazumdar, and T. Sarkar, Phys. Rev. D 101, 043005 (2020).
  * (20) J. A. Wheeler, in Relativity, groups and topology, Les Houches, eds. C. DeWitt, and B. DeWitt, Gordon and Breach, New York (1964).
  * (21) S. W. Hawking, In General Relativity: An Einstein Centenary Survey, eds S. W. Hawking and W. Israel. Cambridge: Cambridge University Press (1979).
  * (22) P. G. Bergmann, Some strangeness in the proportion, eds. H. Woolf, Addison-Wesley Reading, Massachusetts (1980).
  * (23) P. S. Joshi, Gravitational Collapse and Spacetime Singularities (Cambridge University Press, Cambridge, England, 2007).
  * (24) G Lemaitre, Ann. Soc. Sci. Bruxelles I A 53, 51 (1933).
  * (25) R. C. Tolman, Proc. Natl. Acad. Sci. USA 20, 410 (1934).
  * (26) H. Bondi, Mon. Not. Astron. Soc. 107, 343 (1947).
  * (27) C. W. Misner and D. H. Sharp, Phys. Rev. 136, B571 (1964).
  * (28) G.F.R. Ellis and B. G. Schmidt,  Gen. Relativ. Gravit.8, 915 (1977).
  * (29) F. J. Tipler, Phys. Lett. 64A, 8 (1977).
  * (30) C. J. S. Clarke and A. Krolak, J. Geo. Phys. 2, 127 (1986).
  * (31) P. S. Joshi, I. H. Dwivedi,Phys. Rev. D 47, 5357 (1993).
  * (32) P. S. Joshi, Global Aspects in Gravitation and Cosmology (Clendron Press, Oxford, 1993).
  * (33) P. S. Joshi, D. Malafarina, and R. V. Saraykar, Int. J. Mod. Phys. D. 21, 08, 1250066 (2012).
  * (34) K. Mosani, D. Dey, and P. S. Joshi, Phys. Rev. D 101, 044052 (2020).
  * (35) O. D. Toropina, M. M. Romanova, and R. V. E. Lovelace, Mon. Not. Roy. Astron. Soc., 420, 810 (2012).
  * (36) H. Bondi and F. Hoyle, Mon. Not. Roy. Astron. Soc., 104, 273 (1944).
  * (37) R. A. Chevalier, Astrophys. J 346, 847 (1989).
  * (38) R. Belvedere, D. Pugliese, J. A. Rueda, R. Ruffini, and S.-S. Xue, Nuclear Phys. A., 883, 1 (2012).
  * (39) S.Gillessen, F. Eisenhauer, S. Trippe, T. Alexander, R. Genzel, F. Martins, and T. Ott, The Astrophysical Journal, 692, 2 (2009).
  * (40) K. Akiyama et al. [Event Horizon Telescope Collabora- tion], The Astrophysical Journal Letters, 875, 1 (2019).
  * (41) Y. B. Zel’dovich and I. D. Novikov, Sov. Astron. 10, 602 (1967).
  * (42) S. Hawking, Mon. Not. R. Astron. Soc. 152, 75 (1971).
  * (43) J. A. Frieman, M. S. Turner, and D, Huterer, Annual Review of Astronomy and Astrophysics 46, 385 (2008).
  * (44) O. Farooq, F. Ranjeet Madiyar, S. Crandall, and B. Ratra, The Astrophysical Journal, 835, 1 (2016).
  * (45) A. Liddle, An introduction to modern cosmology, Wiley publication (2003).
  * (46) Planck Collaboration, arXiv:1807.06209.
  * (47) M. Y. Khlopov and A. G. Polnarev, Phys. Lett. B, 97,3-4 (1980).
  * (48) J. A. Rueda and R. Ruffini, Astrophys. J., 758, 1 (2012).
  * (49) G. A. Lalazissis, J. König, and P. Ring,  Phys. Rev. C 55, 540 (1997).
  * (50) M.M.Sharma, M.A.Nagarajan, and P.Ring, Phys. Lett. B, 312, 377 (1993).
  * (51) Y. Sugahara and H. Toki, Nuclear Phys. A 579, 557 (1994).
  * (52) D. Hirata, H. Toki, and I. Tanihata, Nuclear Phys. A, 589, 239 (1995).
  * (53) R. Giambò, F. Giannoni, G. Magli and P. Piccione, Communications in Mathematical Physics 235, 545–563 (2003).
  * (54) R. Giambò, Journal of Mathematical Physics 47, 022501 (2006).