††thanks: Current affiliation: University of Oxford, Parks Rd, Oxford OX1 3PJ,
United Kingdom††thanks: Current affiliation: PlanQC GmbH, Lichtenbergstr. 8,
85748 Garching, Germany

# Boltzmannian state counting for black hole entropy in Causal Set Theory

Vid Homšak<EMAIL_ADDRESS>Blackett Laboratory, Imperial College
London, Prince Consort Rd, London SW7 2AZ, United Kingdom  Stefano Veroni
<EMAIL_ADDRESS>Blackett Laboratory, Imperial College London, Prince
Consort Rd, London SW7 2AZ, United Kingdom

(August 31, 2024)

###### Abstract

This paper presents the first numerical study of black hole thermodynamics in
Causal Set Theory, focusing on the entropy of a Schwarzschild black hole as
embodied in the distribution of proposed horizon molecules. To simulate causal
sets we created a highly parallelized computational framework in C++ which
allowed for the generation of causal sets with over a million points, the
largest causal sets in a non-conformally flat spacetime to date. Our results
confirm that the horizon molecules model is consistent with the Bekenstein-
Hawking formula up to a dimensionless constant that can be interpreted as the
fundamental discreteness scale in the order of a Planck length. Furthermore,
the molecules are found to straddle the horizon of the black hole to within a
few Planck lengths, indicating that entropy lives on the surface of the black
hole. Finally, possible implications for the information paradox are drawn. In
particular, we show how the horizon molecules model could yield a finite black
hole temperature cut-off or even prevent full black hole evaporation.

##  I Introduction

Causal Set Theory (CST) is an approach to the problem of quantum gravity which
assumes that continuous spacetime emerges merely as a macroscopic
approximation to the fundamental discrete causal structure, a causal set. A
causal set can be imagined as a network, where the events embedded into
spacetime represent the nodes that connect if and only if they are causally
related.

If CST aims to be the approach to quantum gravity, it is imperative that it
gives predictions that agree with our current understanding of the universe
and works in the area of physics where quantum gravity would be required. A
perfect example of this is black hole thermodynamics, which arises from the
application of Quantum Field Theory on a curved spacetime background [1]. In
this paper, we numerically investigate the horizon molecules model for black
hole entropy. First developed in Ref. [2], this model assumes that the
microstructure of a black hole that gives rise to the Bekenstein-Hawking
entropy is formed by molecules of spacetime situated on its horizon.

An algorithm for the simulation of the Schwarzschild black hole in CST was
outlined in Ref. [3]. However, due to the computationally intensive nature of
the simulations required to study the behaviour of causal sets in curved
spacetime, no extensive numerical studies of CST in the area of black hole
thermodynamics exist until now. To our knowledge, this paper is the first
numerical application of causal sets to the area of black hole thermodynamics.
More generally, before this paper, Ref. [3] was the only numerical study of
CST in a non-conformally flat spacetime. This paper builds on that original
work, creating causal set black holes of cardinality many orders of magnitude
larger.

The manuscript is organised as follows. Section II provides an introduction to
the main disciplines involved in our study, namely black hole thermodynamics
and causal sets theory. In Sec. III, we explain how causal sets can be
simulated in Schwarzschild spacetime. Section IV introduces the general theory
of horizon molecules, explains their relation to black hole thermodynamics,
and outlines the specific characteristics of the ones used in our study.
Section V dives into the highly parallelized C++ simulation framework we
created to study Schwarzschild causal sets extensively for the first time ever
and describes how we tested it to make sure the simulations worked as
expected. Lastly, Sec. VI displays our numerical results, which are then
discussed in Sec. VII.

##  II Background

###  II.1 Black Hole Thermodynamics

A black hole is a region of spacetime where gravity is so strong that nothing
can escape it. The boundary of the black hole is the event horizon: any object
crossing the horizon will inevitably keep falling towards the centre, where a
singularity is present. This is a region of zero volume where curvature and
density diverge, and at which any object’s worldline - the object’s trajectory
in space and time - abruptly interrupts, as if ceasing to exist [4].

Black holes have been shown to be thermal objects endowed with a temperature
and entropy [5, 6]. However, as per the no-hair theorem, established via a
series of contributions [7, 8, 9, 10, 11, 12], black holes appear to be the
simplest objects in the universe, fully specified by only 3 observables: mass
$M$, charge $Q$ and angular momentum $J$. Black hole entropy, which is usually
associated with complexity, suggests the existence of hidden degrees of
freedom beyond those recognised by the no-hair theorem.

Entropy is a property of statistical mechanics and is usually related to the
distribution of microstates of molecules inside the object of interest. It is
an ongoing debate whether the entropy of a black hole counts microstates of
the interior of the black hole or whether it is associated to the horizon
[13]. In this paper we take the attitude that the entropy of a black hole is a
property of the horizon [14].

Heuristically, this is supported by the well-established Bekenstein-Hawking
formula [15, 6, 1]

$S_{BH}=\frac{A}{4\ell_{p}^{2}},$ (1)

which states that the entropy of a black hole scales with the area $A$ of the
event horizon, where $\ell_{p}=\sqrt{{G\hbar}/{c^{3}}}$ is the Planck length.

Since the horizon is just a boundary in empty spacetime, this suggests that
black hole entropy resides in the microstructure of spacetime itself. As
pioneered by Dou and Sorkin [2], in a discrete spacetime one can connect
“spacetime atoms” across the black hole horizon, under certain constraints, to
form horizon molecules, and then seek in these molecules the origin of the
entropy.

Assuming a large number $N_{mol}$ of independent molecules making up the
horizon, each occupying a microstate $n$ with probability $p_{n}$, one can
count them and look at their distribution to measure the Boltzmann entropy of
the event horizon

$S=-N_{mol}\sum_{n}p_{n}\ln{p_{n}}.$ (2)

Since the Bekenstein-Hawking formula states that entropy scales with the area,
we demand that the number of horizon molecules is proportional to the event
horizon area, i.e $N_{mol}\propto A$. This can be studied by simulating
causets embedded into Schwarzschild spacetime.

###  II.2 Causal Set Theory

This section gives an overview of CST, explains how causal sets can be
approximated by a spacetime continuum and provides the relevant nomenclature.

CST began in earnest by Bombelli et al. in 1987 [16]. A causal set (causet) is
a discrete partially ordered set of events $C$ with an order relation $\prec$.
If element $x$ precedes $y$ then $y$ lies in the causal future of $x$, i.e.
$x\prec y$. Here, $x$ is said to be an ancestor ($x$ lies in the past) of $y$
and $y$ is a descendant (lies in the future) of $x$. Formally, a set of events
$C$ with an order relation $\prec$ is defined to be a causal set if it obeys:

  1. 1.

Transitivity: $x\prec y\land y\prec z\Rightarrow x\prec z$.

  2. 2.

Acyclicity: $x\prec y_{1}\prec...\prec y_{N}\neq x$.

  3. 3.

Local Finiteness: $|\\{z|\ x\prec z\prec y\\}|$ is finite.

Transitivity and acyclicity together define the structure to be a partially
ordered set. Local finiteness asserts the discreteness of spacetime. One can
display a causal set using a Hasse diagram, where the causet takes the form of
a network-like structure in which the edges connect the pairs of events that
are directly causally related (see Fig. 1(a)). They represent links, the
irreducible relations between the parent $x$ and its child $y$.

((a))

((b))

\begin{overpic}[width=411.93767pt]{causet_and_BHlightcones.pdf}
\put(-2.0,75.0){(a)} \put(33.0,75.0){(b)} \end{overpic}

Figure 1:  (a) Causal sets are depicted as network of events. They obey
transitivity, acyclicity and local finiteness. (b) Simplistic depiction of a
Schwarzschild black hole in (2+1)D. Light cones bend as you get closer to the
black hole singularity. Once crossing the event horizon at Schwarzschild
radius $r_{S}$, nothing can escape from the black hole, not even light.

Belief in CST stems from the Hawking-King-McCarthy-Malament (HKMM) theorem
[17, 18], which states that continuous spacetime is fully described by its
causal structure up to a local conformal factor, i.e. the causal structure
specifies (9/10)th of the metric in 4-dimensional spacetime, only missing an
intrinsic measure of volume. The central conjecture of CST is that this also
holds for a discrete spacetime. Then, due to local finiteness, a continuum
spacetime volume would correspond to a finite number of causet elements which
can be counted to provide the missing measure of volume. Therefore, causal
sets are conjectured to be equivalent to spacetime.

Causet events need to be uniformly distributed across the spacetime volume
measure to satisfy the Number $\sim$ Volume correspondence. This is possible
with a Poissonian relation between the two [19, 20]. Then, the probability of
having $N$ elements in a region of volume $V$ at density $\rho$ is given by

$\text{Poiss}(N)=\frac{(\rho V)^{N}e^{-\rho V}}{N!}.$ (3)

This relation preserves Lorentz invariance, which is imperative for quantum
gravity, as no observations so far point towards Lorentz symmetry violation
[21, 22, 23].

Finally, for an extensive review of CST, we refer the reader to Ref. [24].

##  III Causal Sets in Schwarzschild Spacetime

This section describes the computational method used to simulate causal sets
in Schwarzschild spacetime. The process of causet embedding into a chosen
(3+1)D spacetime volume region is known as Poisson sprinkling, and is outlined
as follows:

  1. 1.

Pick the sprinkling density $\rho=\ell^{-4}$, where $\ell$ is the discreteness
scale, the average distance between causal set events in 4D spacetime.

  2. 2.

Distribute $N\sim\text{Poiss}(\rho V)$ points uniformly, such that the
probability of having $N$ points in a volume $V$ is described by Eq. 3.

  3. 3.

Connect the elements according to the causal structure of the spacetime
region.

###  III.1 Schwarzschild spacetime

In this paper, we study Schwarzschild black holes, the simplest type of black
holes. They are spherically symmetric, uncharged and have zero angular
momentum. They are only defined by their mass $M$. The curvature of the
spacetime increases as you move closer to the singularity, which causes the
light cones to bend towards the centre of the black hole (see Fig. 1(b)). At
distance $r_{S}=2M$ from the centre, the outer edge of the light cone becomes
vertical, therefore allowing information to flow only in the direction towards
the centre of the black hole. This implies that every event inside the event
horizon is causally unrelated to future events outside of it.

Schwarzschild spacetime in 3+1 dimensions is described by a line element
which, in Schwarzschild coordinates $(t_{S},r,\theta,\phi)$, takes the form of

$ds^{2}=-\left(1-\frac{2M}{r}\right)dt_{S}^{2}+\left(1-\frac{2M}{r}\right)^{-1}dr^{2}+r^{2}d\Omega^{2},$
(4)

where $d\Omega^{2}=d\theta^{2}+\sin^{2}{\theta}d\phi^{2}$. However,
Schwarzschild coordinates are problematic as they give a mathematical
singularity at $r=r_{S}$, hence we use Eddington-Finkelstein original 111The
term ”original” is used to distinguish them from the better known ingoing
(outgoing) coordinates with same name, often denoted as $u$ ($v$). It is
motivated by the fact that in the original papers, they use coordinates that
resemble those in our paper. (EFO) coordinates [25, 26]

$\left(t^{*}=t_{S}+2M\ln{\left|\frac{r}{2M}-1\right|},\,r,\,\theta,\,\phi\right),$
(5)

which are well-behaved everywhere but at the black hole singularity. This
allows us to rewrite the line element into

$\displaystyle ds^{2}=$
$\displaystyle-\left(1-\frac{2M}{r}\right)dt^{*2}+\frac{4M}{r}dt^{*}dr+\left(1+\frac{2M}{r}\right)dr^{2}$
$\displaystyle+r^{2}d\Omega^{2}.$ (6)

The volume element is $dV=r^{2}\sin{\theta}\>dt^{*}\>dr\>d\theta\>d\phi$ and
has the symmetry of a (3+1)D cylinder.

An important feature of the EFO coordinates is that, unlike other sets of
coordinates, the constant $t^{*}$ hypersurface is spacelike everywhere. This
allows us to naturally order events by their time $t^{*}$ coordinates without
further consideration. To prove this statement, let us define the hypersurface
$\Sigma$ as

$\Sigma(x^{\mu})=t^{*}=\mathrm{const.}$ (7)

and its normal vector $n^{\mu}$ as

$n^{\mu}=\partial^{\mu}\Sigma.$ (8)

Given the inner product

$n^{\mu}n_{\mu}=g^{\mu\nu}n_{\nu}n_{\mu},$ (9)

the hypersurface $\Sigma$ is null if Eq. (9) yields 0, is timelike if larger
than 0 and is spacelike if $n^{\mu}n_{\mu}<0$. Since the corresponding
covector is $n_{\mu}=(1,0,0,0)$ and the inverse metric in EFO coordinates
reads

$g^{\mu\nu}=\begin{pmatrix}-\left(1+\frac{2M}{r}\right)&\frac{2M}{r}&0&0\\\
\frac{2M}{r}&1-\frac{2M}{r}&0&0\\\ 0&0&\frac{1}{r^{2}}&0\\\
0&0&0&\frac{1}{r^{2}\sin^{2}{\theta}}\end{pmatrix},$ (10)

we thus obtain from Eq. (9) that $\Sigma$ is indeed spacelike as

$n^{\mu}n_{\mu}=-\left(\frac{2M}{r}+1\right)<0\;\forall r\in\mathbb{R^{+}}.$
(11)

###  III.2 Schwarzschild causal structure

Until now, all extensive CST numerical studies have been done in flat
Minkowski spacetime [24], where it is rather trivial to obtain the causal
relations amongst the events: to determine whether event
$E_{1}=(t_{1},\vec{x}_{1})$ is causally related to event
$E_{2}=(t_{2},\vec{x}_{2})$, one simply demands a null or positive spacetime
interval $\Delta s^{2}=(\Delta t)^{2}-|\Delta\vec{x}|^{2}$.

On the other hand, things become drastically more complicated in spacetimes
that are not conformally flat, such as the one of the Schwarzschild black
hole. In that case, integration of infinitesimal invariant distance $ds$ along
all possible paths from $E_{1}$ to $E_{2}$ is in principle required to
determine the causality, which is computationally very expensive.

\begin{overpic}[width=390.25534pt]{Lightcone_causality_better.pdf}
\put(-2.0,40.0){(a)} \put(54.0,40.0){(b)} \end{overpic} Figure 2:  Causal
relations in Schwarzschild spacetime in (1+1)D, outside (a) and inside (b) the
event horizon. Specifically, null geodesics departing from event $E_{1}$ are
shown and used to infer its causal relations with other events. We define
$t^{*}_{\gamma^{-}}$ as the time at which the earliest null geodesic departing
from $E_{1}$ reaches the spatial position of $E_{2}$. Similarly,
$t^{*}_{\gamma^{+}}$ is the latest. a) Outside the event horizon, only the
earliest $t^{*}_{\gamma}$ exists, and it is sufficient to check whether this
occurs before $t^{*}_{2}$, i.e if $t^{*}_{\gamma}\leq t^{*}_{2}$. b) Inside
the event horizon, one needs to check whether
$t^{*}_{2}\in[t^{*}_{\gamma^{-}},t^{*}_{\gamma^{+}}]$, correcting the method
given in Ref. [3].

However, a better algorithm was introduced by He and Rideout [3]. For a
general pair of events $E_{1}=(t_{1}^{*},r_{1},\theta_{1},\phi_{1})$ and
$E_{2}=(t_{2}^{*},r_{2},\theta_{2},\phi_{2})$, the method determines their
causal relation by checking whether $E_{2}$ lies inside the future light cone
of $E_{1}$. Since the boundaries of light cones are determined by null curves,
the method distils itself into finding the null geodesics emanating from
$E_{1}$, and comparing $t^{*}_{2}$ with the times $t^{*}_{\gamma}$ when null
geodesics reach the spatial position of $E_{2}$ (see Fig. 2).

However, Ref. [3] assumes that a line of constant spatial position is a
possible trajectory everywhere, which is wrong. Due to the bending of the
light cones, no object can remain at the same spatial position inside the
horizon. The corrected procedure is derived in the following subsections,
first for the simple (1+1) dimensional spacetime, then in dimensions up to
(3+1).

####  III.2.1 Causality in (1+1)D Schwarzschild spacetime

Consider two events $E_{1}=(t^{*}_{1},r_{1})$ and $E_{2}=(t^{*}_{2},r_{2})$.
The line element from Eq. (III.1) in 1+1 dimensions yields

$ds^{2}=-{dt^{*}}^{2}+dr^{2}+\frac{2M}{r}(dt^{*}+dr)^{2},$ (12)

which can be rewritten into

$ds^{2}=(dt^{*}+dr)\left[\left(\frac{2M}{r}-1\right)dt^{*}+\left(\frac{2M}{r}+1\right)dr\right].$
(13)

Setting $ds^{2}=0$ yields the two radial null geodesics (see Fig. 2). The
first one gives

$t^{*}_{2}-t^{*}_{1}=r_{1}-r_{2},$ (14)

describing an always ingoing null geodesic (left/lower straight line in Fig.
2) while the second one gives

$t^{*}_{2}-t^{*}_{1}=\left[r+4M\ln{(r-2M)}\right]_{r_{1}}^{r_{2}}.$ (15)

describing the other geodesic, which is ingoing inside the horizon and
outgoing outside of it (upper curved in Fig. 2).

Therefore, we are left with the following possible scenarios.

  1. 1.

If $E_{1}$ is outside of the event horizon, then the necessary and sufficient
condition is

$t^{*}_{2}-t^{*}_{1}\geq\begin{cases}r_{1}-r_{2},&\textrm{if}\;\;r_{2}\leq
r_{1}\\\
r_{2}-r_{1}+4M\ln{\left(\frac{r_{2}-2M}{r_{1}-2M}\right)}\;&\textrm{if}\;\;r_{2}>r_{1}\end{cases}.$
(16)

  2. 2.

If $E_{1}$ is inside the event horizon and $r_{2}>r_{1}$, they are trivially
causally unrelated.

  3. 3.

If both events are inside the event horizon and $r_{1}\geq r_{2}$, then the
two events are connected if and only if

$r_{1}-r_{2}\leq t_{2}^{*}-t^{*}_{1}\leq
r_{2}-r_{1}+4M\ln{\left(\frac{2M-r_{2}}{2M-r_{1}}\right)}.$ (17)

The third point corrects Ref. [3], where it is stated that “for both events
inside the horizon, the necessary and sufficient condition for their causal
relation is still $t^{*}_{2}\geq t^{*}_{1}+r_{1}-r_{2}$”, which is wrong as it
neglects the upper bound from Eq. (15).

####  III.2.2 Causality in (3+1)D Schwarzschild spacetime

Noting that the method always involves only two events
$E_{1}=(t_{1}^{*},r_{1},\theta_{1},\phi_{1})$ and
$E_{2}=(t_{2}^{*},r_{2},\theta_{2},\phi_{2})$ simultaneously, we can simplify
our calculations in (3+1)D EFO coordinates by rotating the spatial part of the
coordinate system to an equatorial plane containing both of the events, such
that $\theta^{\prime}_{1}=\theta^{\prime}_{2}=\pi/2$ and
$\phi^{\prime}_{1}=0$. As the equatorial plane is totally geodesic, we can
effectively neglect $\theta$ and rewrite their coordinates as

$\displaystyle E_{1}=\left(t_{1}^{*},r_{1},0\right),$ (18) $\displaystyle
E_{2}=\left(t_{2}^{*},r_{2},\arccos{\left(\sin{\theta_{1}}\sin{\theta_{2}}\cos{\Delta\phi}+\cos{\theta_{1}}\cos{\theta_{2}}\right)}\right),$
(19)

reducing the problem to 2+1 dimensions. Compared to the radial (1+1)D case, in
(3+1)D the procedure to determine the causal relations becomes more involved
since it is impossible to solve for null geodesics analytically.

Given an affine parameter $\lambda$, they satisfy

$\displaystyle\left(1-\frac{2M}{r}\right)\left(\frac{dr}{d\lambda}\right)^{2}-\left(1-\frac{2M}{r}\right)\left(\frac{dt^{*}}{d\lambda}\right)^{2}+r^{2}\left(\frac{d\phi}{d\lambda}\right)^{2}$
$\displaystyle=0,$ (20)

and since the line element in Sec. III.1 is independent of $t^{*}$ and $\phi$,
the following are conserved quantities

$\displaystyle E$
$\displaystyle=\left(1-\frac{2M}{r}\right)\frac{dt^{*}}{d\lambda},$ (21)
$\displaystyle L$ $\displaystyle=r^{2}\frac{d\phi}{d\lambda}.$ (22)

Moreover, we can rewrite the derivatives as

$\displaystyle\frac{dr}{d\lambda}$
$\displaystyle=\frac{dr}{d\phi}\frac{d\phi}{d\lambda}=\frac{L}{r^{2}}\frac{dr}{d\phi},$
(23) $\displaystyle\frac{dt^{*}}{d\lambda}$
$\displaystyle=\frac{dt^{*}}{d\phi}\frac{d\phi}{d\lambda}=\frac{L}{r^{2}}\frac{dt^{*}}{d\phi},$
(24)

to get rid of $\lambda$ and re-express $dr/d\phi$ in terms of $u=1/r$, which
is computationally advantageous as it leads to expressions with lower powers,
such that

$\frac{dr}{d\phi}=u^{-2}\frac{du}{d\phi}.$ (25)

Furthermore, we can define the ratio

$\eta=\frac{E}{L},$ (26)

which is necessarily finite since $L=0$ describes the (1+1)D case with
$\Delta\phi=0$. Finally, inserting previous equations into Eq. (III.2.2), we
obtain the two differential equations

$\displaystyle\frac{d\phi}{du}$
$\displaystyle=\pm_{du}\left(\eta^{2}+2Mu^{3}-u^{2}\right)^{-1/2}$ (27)
$\displaystyle\frac{dt^{*}}{du}$
$\displaystyle=\frac{1}{u^{2}(2Mu-1)}\left[\frac{\pm_{du}\eta}{\sqrt{\eta^{2}+u^{2}(2Mu-1)}}+2Mu\right]$
(28)

where

$\pm_{du}=\mathrm{sign}(u_{2}-u_{1})$ (29)

depends on the order in which the two events were initially chosen.

Since the initial and final spatial boundary conditions $(u,\phi)$ are known,
we can fit Eq. (27) for $\eta^{2}$ to obtain $\pm\eta$ values and then use
them to integrate Eq. (28). This gives us two values $t^{*}_{+}$ and
$t^{*}_{-}$, corresponding to the times when each null geodesic reaches the
spatial coordinate of $E_{2}$.

In Appendix A, we prove explicitly the conditions for causality. We show that
in the case that $E_{1}$ is outside of the event horizon, only one of the
times $\\{t^{*}_{-},t^{*}_{+}\\}$ is physically meaningful, $t^{*}_{-}$. Then,
$t_{2}>t^{*}_{-}$ is the sufficient and necessary result to determine that
$E_{2}$ lies in the causal future of $E_{1}$, i.e that $E_{1}\prec E_{2}$. On
the other hand, if both points lie inside the horizon such that
$r_{S}>r_{1}>r_{2}$, $E_{1}\prec E_{2}$ if and only if
$t_{2}\in[t^{*}_{-},t^{*}_{+}]$. In fact, the bending of light cones (see Fig.
2), poses an upper time bound. Finally, we also show that the value of $\phi$
should be restricted to lie in $[0,2\pi]$.

####  III.2.3 Sufficient conditions for determining the causality

It is worth noting that a large number of event pairs in (3+1)D does not
require the general approach with integration and fitting for $\eta^{2}$ to
determine the causal structure of the Schwarzschild spacetime.

Instead, there exists a set of simple sufficient but not necessary conditions
that we check, for each pair of events, before applying the generic method
outlined above to determine whether the points are causally related or not.
This speeds up the procedure significantly since the checks are
computationally very inexpensive compared to the fitting and integration. To
get the full list and the derivations of the conditions, refer to Sec. 2.3 in
Ref. [3].

##  IV Horizon molecules

This section defines horizon molecules, determines the coefficient for the
scaling of number of molecules with the area of the black hole horizon and
calculates the curvature corrections in Schwarzschild spacetime.

In this paper, we investigated a new type of horizon molecule, the
$\Lambda_{n}$ molecule, which will be discussed in detail in the following
subsection.

Figure 3: Schematic representation of $\Lambda_{n}$ and Link horizon
molecules. They consist of connected points inside and outside the event
horizon $\mathcal{H}$, satisfying certain constraints. The spacelike
hypersurface $\Sigma:\>t^{*}=\mathrm{constant}$ defines the time of the
entropy measurement. $\Lambda_{n}$ molecules are defined by having an element
$k$ inside the horizon which is maximal with respect to $\Sigma$ and connects
to $n$ elements outside the horizon that only contain $k$ in their causal
future below $\Sigma$. A single $\Lambda_{n}$ molecule is equivalent to $n$
Links.

In general, a horizon molecule consists of one point inside the event horizon
$\mathcal{H}$ connected to one or more points outside it via links and lies
below some $t^{*}=constant$ hypersurface $\Sigma$. In EFO coordinates,
$\Sigma$ is spacelike everywhere as shown in Eq. (11), which allows us to
think of $t^{*}$ coordinate as time. Since entropy is only defined as a
function of time, $\Sigma$ signifies the moment when entropy is “measured”.
Therefore, one can count the number of molecules, study their distribution,
and obtain results for the proportionality factor for $N_{mol}\propto A$ in
Eq. (2).

On the other hand, scaling can be derived analytically for certain molecules.
Knowing that the probability that a region $\Delta R$ of spacetime with volume
$V$ contains $N$ events is given by the Poisson distribution as in Eq. (3), we
can calculate the expected number of horizon molecules and compare it with the
numerical result.

Analytic studies [27, 2] assume that at sufficiently high sprinkling density,
molecules are localised in space and time. This allows us to work in flat
approximation with Rindler coordinates, where the event horizon and the
constant time hypersurface are given by

$\mathcal{H}:t=x,\quad\Sigma:t=\mathrm{const.}$ (30)

This flat-space representation is used to display $\Lambda_{n}$ and Link
horizon molecules, described in the following subsection, in Fig. 3.

###  IV.1 $\Lambda_{n}$ and Link molecules

The $\Lambda_{n}$ molecule is given by a maximal element $k$ inside the event
horizon that lies in the causal future of $n$ mutually disconnected elements
$\\{p_{n}\\}$ outside of it (see Fig. 3), which can only contain $k$ in their
causal future below $\Sigma$. Unfortunately, we cannot derive the analytic
expression for the average expected number of $\Lambda_{n}$ molecules.

However, we can derive the analytical expression for Link molecules, which are
closely related to $\Lambda_{n}$. They are links crossing the horizon whose
past element is maximal but one with respect to $\Sigma$ 222Upper-case Link
denotes the horizon molecule, whereas lower-case link the nearest-neighbour
relation between causet events.. The analytic result for the average expected
number of Links $\langle N_{L}\rangle$ can then be compared with the numerical
result obtained by weighted counting of $\Lambda_{n}$ molecules, where each
$\Lambda_{n}$ contributes $n$ Links.

The expression for $\langle N_{L}\rangle$ has been computed by Barton et al.
[27]. Here we sketch the calculation. Given the discreteness density $\rho$,
the probability that a point $p$ outside the horizon in the region $\delta
R_{p}$ connects only to one maximal future point inside the horizon is given
by

$\displaystyle P_{L}(\delta R_{p})$ $\displaystyle=P(1\;\mathrm{in}\;\delta
R_{p})\cdot P(1\;\mathrm{in}\;I^{+}_{in}(p))$ (31)
$\displaystyle=\rho^{2}\delta V_{p}V^{+}_{in}(p)e^{-\rho V^{+}_{in}(p)}.$ (32)

$V^{+}(p)$ is the volume of the future light cone of $p$ capped by the
$\Sigma$ hypersurface (coloured orange in the $\Lambda_{1}$ schematic in Fig.
3) made of $V^{+}_{in}(p)$ and $V^{+}_{out}(p)$ partial volumes, inside and
outside the horizon respectively, and represents $I^{+}(p)$, the causal future
of $p$. We denote $\mathcal{J}$ as the intersection of the event horizon
$\mathcal{H}$ and hypersurface $\Sigma$, which represents the event horizon at
time $t_{\Sigma}$, and is just a point in 1+1 dimensions.

Integrating over all spacetime regions $\delta R_{p}$ in the past of the
horizon and hypersurface intersection $I^{-}(\mathcal{J})$, yields the
expected number of Links

$\langle
N_{L}\rangle=\rho^{2}\int_{I^{-}(\mathcal{J})}dV_{p}\;V^{+}_{in}(p)\;e^{-\rho
V^{+}_{in}(p)}.$ (33)

Noting that we are working with Rindler horizon in (3+1)D flat spacetime,
analytic expressions for the volumes can be found. The average expected number
of Links for the horizon area $A$ in the flat approximation is then given by

$\langle N_{L}^{flat}\rangle=\frac{\sqrt{3}}{10}\sqrt{\rho}A.$ (34)

In terms of the discreteness scale $\ell=\rho^{-1/4}$, this reads

$\langle
N_{L}^{flat}\rangle=a_{L}^{(0)}\frac{A}{\ell^{2}}=\frac{\sqrt{3}}{10}\frac{A}{\ell^{2}}.$
(35)

###  IV.2 Curvature corrections in Schwarzschild spacetime

Figure 4: The C++ simulation framework available on GitHub. It allows for an
efficient generation of causal sets in Minkowski and Schwarzschild spacetime.

In general, the average number of molecules, for any horizon $\mathcal{H}$, is
given by an infinite series

$\langle
N\rangle=a\frac{A}{\ell^{2}}=\left[a^{(0)}+a^{(1)}(\mathcal{H},\Sigma)\frac{\ell}{A}+\mathcal{O}\left(\frac{\ell^{2}}{A}\right)\right]\frac{A}{\ell^{2}}.$
(36)

To go beyond the flat spacetime approximation $a^{(0)}$ and estimate the 1st
order curvature correction $a^{(1)}$, which depends on the structure of
$\mathcal{H}$ and $\Sigma$, we follow the procedure described by Sec. 5 in
Barton et. al [27].

The general expression for the 1st order correction of Links in terms of the
trace $K$ of the extrinsic curvature of $\Sigma$, its component $K_{m}$
tangential to $\Sigma$ and orthogonal to $\mathcal{J}$, and the null expansion
$\vartheta$ of the horizon, all evaluated on $\mathcal{J}$ in 3+1 dimensions,
is given by

$\frac{a_{L}^{(1)}}{\sqrt{A}}\approx-(0.036K+0.088K_{m}+0.021\vartheta).$ (37)

In Appendix B, we compute $a_{L}^{(1)}$ for the specific case of the
Schwarzschild spacetime. This gives

$a_{L}^{(1)}\approx-0.0558.$ (38)

Therefore, the average expected number of Links in Schwarzschild spacetime is

$\langle
N_{L}^{BH}\rangle=a_{L}\frac{A}{\ell^{2}}=a_{L}^{(0)}\frac{A}{\ell^{2}}\left(1-0.322\frac{\ell}{\sqrt{A}}+\mathcal{O}(\ell^{2})\right).$
(39)

The correction becomes negligible when $A$ increases to the size of ordinary
black holes since $\ell\sim\ell_{p}$ is infinitesimal.

##  V Simulation Framework

This section describes the simulation framework we constructed to efficiently
simulate causal sets embedded into Schwarzschild spacetime and the tests we
implemented to ensure it works as intended.

((a))

((b))

\begin{overpic}[width=138.76157pt]{BHSprinkling_2D_projection_2.pdf}
\put(-206.0,200.0){(a)} \put(-206.0,97.0){(b)} \put(-102.0,97.0){(c)}
\put(-1.0,97.0){(d)} \end{overpic}

Figure 5: Main tests for causet generation in Schwarzschild spacetime. (a)
Test of the Poisson distribution of events by sampling the number of points in
regions of equal volume over 100 causet realizations with the cardinality of
$10\,000$ each. The first row refers to sprinkling into a filled cylinder,
whereas the second to a hollow cylinder. The columns refer to cutting regions
of spacetime limited in $t^{*}$, $r$, $\theta$ and $\phi$, respectively.
Points are indeed distributed following a Poissonian distribution. (b) 2D
projection of a very thin slice ($0.25^{\circ}$) of a 3D black hole causet.
(c) and (d) highlight some details of (b), showing that even when events are
extremely close to the null geodesics, they are only connected if and only if
inside the light cone, as it should be.

It is important to note that other than the work done in Ref. [3], which laid
out the foundations, there exist no extensive studies of Schwarzschild causal
sets. Particularly this is because the procedure used to determine the causal
structure of Schwarzschild spacetime, as described in Sec. III, is extremely
computationally expensive since for a generic pair of events minimization and
integration need to be performed. Based on our experience, generating
Schwarzschild causets is approximately 200 times slower compared to flat
spacetime, which explains why other causet numerical studies did not attempt
it.

While we initially used the causal set package written by Minz [28], we
realised that Python is too slow and decided to create our simulation
framework in C++, due to its speed. The package is publicly available on
GitHub and allows the user to simulate causal sets in regions of different
shapes in both Minkowski and Schwarzschild spacetime for the dimensions
$D=2,3,4$. It furthermore contains functions to analyse horizon molecules,
estimate dimensions of causal sets, plot them and more. It is the most
efficient causet package publicly available to our knowledge.

In general, during the causet generation, our program only saves events’
coordinates and the causal matrix $C$, where $C_{ij}=1$ if $i\prec j$ and
$C_{ij}=0$. The framework also allows for storing the sets of points in the
past and/or future sets and corresponding sets of past/future links for each
element as they could be preferable for certain applications. This drastically
speeds up the process while simultaneously reducing memory usage. The latter
in the end turned out to be an important bottleneck of the simulation, since
the realization of a causet with $N=10^{6}$ required up to 4 TB of RAM.

Our code uses CPU parallelization as much as possible to speed up the
simulations. Importantly, we used dynamical (unordered) parallelization by
OpenMP, which allows for threads that have already completed their tasks to
take on new ones, not requiring to wait for the slowest thread to finish. In
our simulations, we achieved maximum causet cardinality $N=1,011,423$, by far
the largest recorded causal set in curved space, which allowed us to obtain
sufficient statistics. Overall, the estimated speed-up of the causal set
generation in the parallelized C++ framework on 256 cores using Imperial
College’s HPC cluster compared to the non-parallelized Python setup on a
personal laptop, using the same generation algorithms, comes out to be
approximately a factor of 15000.

###  V.1 Testing the simulation framework

It is crucial to ensure that the simulation framework works as intended.
Therefore, we implemented various checks to test the validity of our
simulations. This section describes the checks relevant to the purposes of
this article, whereas other testing of the framework is described in Appendix
D.

As per Sec. III, the embedding of a causet in spacetime relies on two
ingredients: a Poissonian correspondence between the number of events and
spacetime volume, and a correct derivation of causal relations. To assess
whether events were correctly distributed, we divided the simulated embedding
into regions of equal spacetime volumes and counted the number of events in
each region. The distribution of regions with a specific number of events
indeed reproduced a Poisson distribution, both for a full and hollow cylinder
of spacetime (see Fig. 5(a)).

Due to the novelty of our simulations, we lacked numerical results to
benchmark our code with. In particular, to test the correctness of the causal
relation. However, we could visually assess the correctness of (1+1)D
simulations by checking that causal relation only occur within light cones,
for which analytic solutions are given by Eq. (16) and Eq. (17). In the same
way, we could test (2+1)D simulations in the 2-dimensional limit. This was
done by sprinkling causets in a very thin slice of a cylinder and projecting
the results on a plane (see Fig. 5(b)). Moreover, the number density scales
with $r$, as expected for (2+1)D spacetime. Note that this also assesses the
correctness of the 4D simulations, as the methods used to simulate black hole
causets in 4 and 3 dimensions are exactly the same, due to the symmetry about
$\theta=\pi/2$.

##  VI Results

This section describes the results we obtained on the entropy and the
properties of the horizon molecules.

We simulated causets in Schwarzschild spacetime with density
$\rho=\ell^{-4}=5000$ and cardinality up to $\sim 10^{6}$ events to count the
number of $\Lambda_{n}$ molecules and study their distribution. Each data
point in the following plots is averaged over 200 realizations.

The $\Sigma$ hypersurface was set at $t^{*}=0$. We sprinkled into a (3+1)D
hollow cylinder with points lying in $t^{*}\in[-4,0]$ and
$r\in[r_{S}-3,\,r_{S}+3]$ in terms of discreteness units $\ell$. These bounds
were chosen close enough to the $\Sigma$ and $\mathcal{H}$ hypersurfaces to
allow for an efficient simulation of large black holes, but far enough to
avoid affecting the behaviour of the molecules. All points that are part of
$\Lambda_{n}$ molecules consistently lie far from the chosen boundaries (see
Fig. 6(a)).

The number of $\Lambda_{n}$ shows clear proportionality with the horizon area
(see Fig. 6(b)), giving the Pearson correlation coefficients $r\approx 1$ for
$n<8$, at which point we lack enough statistics. The number of Link molecules,
which is a weighted sum of $\Lambda_{n}$, also scales linearly with $A$ (see
Fig. 6(c)). Therefore, the proportionality constants we obtain are:

  * •

Neglecting curvature corrections, the proportionality coefficient is estimated
as

$a^{(0)}_{L[\text{1}]}=0.1757\pm 0.0004.$

  * •

Setting the curvature correction $a_{L}^{(1)}$ to the theoretical value
$a^{(1)}_{L[\text{th}]}=-0.0558$, $a^{(0)}$ is estimated as

$a^{(0)}_{L[\text{2}]}=0.1761\pm 0.0004$

.

  * •

Fitting both $a^{(0)}$ and $a^{(1)}$ gives

$a^{(0)}_{L[\text{3}]}=0.173\pm 0.001,\qquad a^{(1)}_{L[\text{3}]}=0.4\pm
0.2.$ (40)

Recall that the analytic result from Barton et al. [27] is
$a^{(0)}_{L[\text{th}]}=\sqrt{3}/10\approx 0.1732$. This perfectly agrees with
the numerical fit when accounting for the 1st order curvature correction.
Oppositely, we note our best numerical estimate for $a^{(1)}_{L}$ differs from
the analytic estimate both in order of magnitude and sign. This will be
discussed later. Finally, non trivially, Fig. 6(c) also shows that the
relative uncertainty in the number of Links $\sigma_{L}/\langle N_{L}\rangle$
is small and further decreases with the size of the black hole.

\begin{overpic}[width=212.47617pt]{AAA_Rho5000_lambdas.pdf}
\put(-100.0,73.0){(a)} \put(0.0,73.0){(b)} \end{overpic}

\begin{overpic}[width=182.1196pt]{AAA_n_lambda_prob_distribution_expx_logy_withcloud.pdf}
\put(-136.0,88.0){(c)} \put(-2.0,88.0){(d)} \end{overpic}

((a))

((b))

((c))

((d))

Figure 6:  Numerical results of our simulations. a) Averages of the maximal
separation of elements forming horizon molecules from the relevant
hypersurfaces, with $1\sigma$ and $3\sigma$ deviations highlighted. Above, the
maximum temporal separation $\Delta t_{\mathrm{max}}$ from $\Sigma:t^{*}=0$.
Below, the maximal spatial separation $\Delta r_{\mathrm{max}}$ between the
horizon and the furthest elements in $\Lambda_{n}$ molecules, both from the
inside and outside. b) Number of $\Lambda_{n}$ molecules $\langle
N_{\Lambda_{n}}\rangle$ scaling linearly with the horizon area $A$. The
largest molecule’s size we found was 11. c) Average umber of Links $\langle
N_{L}\rangle$ (black) as a function of the horizon area $A$, and relative
standard deviation (blue crosses). The azure and green dashed lines represent
two analytical models: the analytical flat spacetime approximation from Eq.
(35), and the scaling with curvature correction from Eq. (39), respectively,
both with $a^{0}_{L}=\sqrt{3}/10$. The red line is our best fit of the
curvature-corrected expansion, specified in Eq. 40. Note,
$A_{\ell}=A/\ell^{2}$. The first two models lie within the uncertainties but
are consistently lower, as evident in the inset. The red line is a perfect
fit. d) The distribution of $\Lambda_{n}$ is well described by a falling
exponential function $p_{n}\propto e^{-\chi n}$, where $\chi=1.53\pm 0.03$.
The yellow cloud around the fit shows the $1\sigma$ deviation. Likely due to
the finiteness of our simulations, the number of molecules of size larger than
7, for which poor statistics was found, is slightly lower than the fit.

We also studied the distribution of $\Lambda_{n}$ molecules (see Fig. 6(d)).
Disregarding data for $n\geq 8$ which has insufficient statistics due to the
finiteness of our simulations, we discover that we can model the distribution
very well by a falling exponential

$p_{n}=\left(e^{\chi}-1\right)e^{-n\chi},\quad\chi=1.53\pm 0.03.$ (41)

Using the relation between the number of Links and $\Lambda_{n}$ molecules and
the aforementioned probability distribution, it is easy to infer that

$\langle N_{\Lambda}\rangle=\left(1-e^{-\chi}\right)\langle N_{L}\rangle.$
(42)

Inserting the numerical scaling factor $a^{(0)}_{L[\text{3}]}=0.173\pm 0.001$
into Eq. (39) and using the probability distribution of $\Lambda_{n}$
molecules, the Boltzmann entropy given by Eq. (2) yields

$S_{\Lambda}=(0.091\pm 0.002)\>\frac{A}{\ell^{2}}.$ (43)

Equating Eq. (43) with the black hole entropy $S_{BH}$ from Eq. (1) we can
give a physical meaning to the discreteness scale $\ell$. In terms of Plank
length $\ell_{p}$, we obtain

$\ell=(0.603\pm 0.003)\>\ell_{p}.$ (44)

Inverting, $G\hbar=(2.75\pm 0.03)c^{3}\ell^{2}$.

##  VII Discussion

This section discusses the obtained results and situates the paper into the
landscape of black hole thermodynamics in Quantum Gravity.

###  VII.1 Numerical results

We confirmed that the entropy of the Schwarzschild black hole modelled by the
horizon molecules unequivocally scales linearly with the horizon area for
$A\gg\ell^{2}$, in agreement with the well-established Bekenstein-Hawking
formula. This is a very important result for two reasons. It numerically
demonstrates for the first time that CST can yield correct predictions in
curved (non-conformally flat) spacetime and that it can be used to study black
hole thermodynamics.

The analytic estimate by Barton et al. [27] for $a_{L}^{0}$ perfectly matches
our numerical fit. However, the analytic and numerical estimates for the 1st
order curvature correction disagree. The analytic model lies within our
uncertainties but is systematically lower than our averages. We note that a
discrepancy was also found by Ref. [29] with the 1st order correction due to
horizon expansion in their study of a dynamical black hole in flat spacetime.

Moreover, we were able to “count” the Planck length. If causal sets were the
underlying structure from which the spacetime continuum emerges, and if
correspondingly CST was the fundamental theory, then it would have no free
parameters. This would imply that Planck length $\ell_{p}$ is some multiple of
the discreteness length $\ell$, which is the fundamental length of size $1$ in
the fundamental units. As generally anticipated, we found that
$\ell\sim\ell_{p}$, which is where a lot of general physics breaks down.

It is interesting to note the sizes of the black holes we simulated. The
largest ones had an area of $A=42000\>\ell^{2}$ and since $\ell\sim\ell_{p}$
this corresponds to mass $M=\sqrt{A/16\pi}\sim 30$ in Planck units, i.e. $\sim
10^{-7}$ kg. For comparison, the masses of the observed black holes usually
range from $10^{30}$ kg to $10^{38}$ kg, yet, our tiny black holes yielded
such good results.

Furthermore, we were intrigued by the fact that the falling exponential
distribution for $\Lambda_{n}$ molecules had an exponent $\chi=1.53\pm 0.03$,
very close to $3/2$ (see Fig. 6(d)). If we treat $\Lambda_{n}$ molecules as
thermal molecules that exist in $n$ different microstates, we can write their
thermal energy distribution as $p_{n}\propto e^{-E_{n}/k_{B}T}$. Comparing
this with Eq. (41), we find that

$E_{n}\approx\frac{3n}{2}k_{B}T,$ (45)

which is the energy of $n$-atomic molecules in 3D space. Since $\Lambda_{n}$
molecules exist in 3D space at a given moment of time, we hoped to be able to
generalise this expression, such that Eq. (45) would have been just a specific
case for $f=3$ degrees of freedom. However, simulations of causets in (2+1)D
spacetime showed that $\chi\sim 1.5$ in (2+1)D as well (see App. E).
Therefore, it is false that $\chi\sim f/2$ in all dimensions, and
$\Lambda_{n}$ molecules cannot be interpreted as thermal molecules.

###  VII.2 Holographic principle and information paradox

Hawking argued that black hole radiation only depends on its external degrees
of freedom [17]: mass, charge and angular momentum. By emitting radiation,
black holes can evaporate away. In this case, all that is left, the radiation,
thus retains only information about those three features. All other
information about previously absorbed objects is lost. This breaks the
fundamental law of conservation of quantum information [17, 30], and
constitutes the black hole information paradox.

There are two possible solutions: either black hole evaporation stops at the
Planck scale, at which Hawking’s derivation breaks down, or radiation contains
more information than assumed. Recently, different suggestions [31, 32, 33,
34] have been made in support of the latter. These follow, albeit in different
ways, from the holographic principle, which states that a complete description
of a volume of space is given by the degrees of freedom of its surface [35,
36]. Horizon molecules are distributed over the horizon’s surface, within an
infinitesimal shell of thickness $\sim 4\ell$. They form the horizon. As they
encode the entropy of the black hole, they might provide a relation to the
holographic principle and the information paradox. This consideration requires
further studies.

###  VII.3 Planckian black holes

Black holes form from the gravitational collapse of massive stars. They form
with a mass of at least the order of the mass of the Sun [37], hence at most a
Hawking temperature of order $\sim 10^{-8}$K. This is 100 million times colder
than the Cosmic Microwave Background (CMB) [38]. Therefore, the energy black
holes absorb from CMB surpasses the losses from their own radiation, making it
impossible for current black holes to evaporate into smaller masses.

However, the CMB will drop. Therefore, black holes will eventually be allowed
to evaporate. Furthermore, it was suggested that primordial black holes of
Planckian size could have formed in the initial stages of the universe via
different processes [39]. Hawking’s work breaks down at Planckian scales,
begging the question of how a black hole of Planckian size behaves. In this
section, we provide an answer within the theory of horizon molecules. Given
that the molecular model might break down at such small scales, especially
given the current lack of quantum effects in the picture, the results of this
subsection need to be taken cum grano salis. Nevertheless, we believe they
might offer interesting insights for further research.

To answer in the context of our molecular approach to black hole entropy,
consider that, in general, the average number of horizon molecules is a sum of
the form

$\langle
N\rangle=\sum_{i=0}^{\infty}a^{(i)}\left(\frac{\sqrt{A}}{\ell}\right)^{2-i},$
(46)

where $(i)$ labels the order of the curvature correction - with $i=0$ being
the flat spacetime approximation and $a^{(i)}$ a real constant for all $i$.
This implies an entropy

$S=\sum_{i=0}^{\infty}\sigma^{(i)}\left(\frac{\sqrt{A}}{\ell}\right)^{2-i},$
(47)

where, if the molecules have distribution $\\{p_{n}\\}$, it yields that
$\sigma^{(i)}=-a^{(i)}\sum_{n}p_{n}\ln{p_{n}}$. The leading term gives the
standard black hole entropy, hence $\ell^{2}=4\sigma^{(0)}$ and

$S=\frac{A}{4}+\sum_{i=1}^{\infty}\sigma^{(i)}\left(\frac{\sqrt{A}}{\ell}\right)^{2-i}.$
(48)

Thus, for a Schwarzschild black hole $A=16\pi M^{2}$ we can write

$S=4\pi
M^{2}+\sum_{i=1}^{\infty}\varsigma^{(i)}\left(\frac{M}{\ell}\right)^{2-i},$
(49)

where $\varsigma^{(i)}=(4\sqrt{\pi})^{2-i}\sigma^{(i)}$. Applying the mass-
energy equivalence $U=M$, the 1st law of thermodynamics gives

$\frac{1}{T}=\frac{\partial S}{\partial M}.$ (50)

Hence, by rearranging we obtain

$\frac{1}{T}=8\pi
M+\frac{\varsigma^{(1)}}{\ell}-\frac{1}{\ell}\sum_{i=3}^{\infty}(i-2)\varsigma^{(i)}\left(\frac{\ell}{M}\right)^{i-1}.$
(51)

The first term dominates for large $M\gg\ell$, in which case this results in
the Hawking temperature. The second represents a correction of order unity. In
the limit $M\rightarrow 0$, the second term alone could give a cutoff at
approximately the Planck temperature $T_{p}\sim 10^{32}K$. A similar cutoff is
predicted by loop quantum gravity too, albeit due to different mechanisms
[40].

Before investigating the effects of the third contribution, let’s outline the
constraints on the real constants $\varsigma^{(i)}$. Define
${\mu}=\ell/M\in(0,\mu_{d}]$, where $\mu_{d}$ could be infinite as $M$
vanishes. We require $\forall\mu\in(0,\mu_{d}]$:

  1. (i)

Entropy Positiveness. From Eq. (49), relabelling the summation index with
$j=i-2$ and taking out the first two elements,

$\quad
S=\frac{4\pi\ell^{2}}{\mu^{2}}+\frac{\varsigma^{(1)}}{\mu}+\varsigma^{(2)}+\sum_{j=1}^{\infty}\varsigma^{(j+2)}\mu^{j}\;\geq
0.$ (52)

  2. (ii)

Entropy Finiteness. From Eq. (52),

$\sum_{j=1}^{\infty}\varsigma^{(j+2)}\mu^{j}\;\;\mathrm{is\;\;finite\;}.$ (53)

  3. (iii)

Temperature Positiveness. From Eq. (51),

$\frac{1}{T}=\frac{8\pi\ell}{\mu}+\frac{\varsigma^{(1)}}{\ell}-\frac{\mu}{\ell}\;\sum_{j=1}^{\infty}j\;\varsigma^{(j+2)}\mu^{j}\;\geq
0.$ (54)

A type of molecule could perhaps be envisaged such that black holes are not
allowed to fully evaporate. This would be possible if the following condition
is also satisfied:

  1. (iv)

Temperature Loss. From Eq. (54), $T\rightarrow 0^{+}$ for a non-zero mass if
there exists a finite $\mu_{d}$ such that

$\sum_{j=1}^{\infty}j\;\varsigma^{(j+2)}\mu_{d}^{j}\;\rightarrow-\infty\;$
(55)

Horizon molecules with real constants $\varsigma^{(i)}$ satisfying conditions
(i-iv) would give an always positive finite entropy and a temperature which is
always positive and finite but for $\mu=\mu_{d}$, at which it vanishes.
Physically, the black hole would radiate less and less, asymptotically
approaching a minimum mass $M_{\mathrm{min}}=\ell/\mu_{d}$, thus preventing
full evaporation. The black hole would reach an equilibrium, at a finite
temperature, when the emission matches the ingoing radiation from the CMB: the
black hole would be a tiny dark fragment wandering in the universe.

This could resolve the information paradox. The information related to the
absorbed material could be stored inside the horizon, forever. Furthermore, as
black holes can form where none was previously present, a mechanism preventing
their complete evaporation could break charge-parity-time-reversal (CPT)
symmetry [17].

Can such a set of real constants $\varsigma^{(i)}$ exist? It can. As an
example, consider the Riemann zeta function

$\zeta(x)=\sum_{i=1}^{\infty}i^{-x}.$ (56)

Assume that $\varsigma^{(i+2)}=ai^{-x}\;\;\forall i\in\mathbf{N}$, where $a$
and $x$ are real constants. This leaves $\varsigma^{(1)}$ and
$\varsigma^{(2)}$ momentarily undefined.

Then, the sums in (ii) and (iv) are $a\mathrm{Li}_{x}(\mu)$ and
$a\mathrm{Li}_{x-1}(\mu_{d})$ respectively, where

$\mathrm{Li}_{x}(\mu)=\sum_{j=1}^{\infty}j^{-x}\mu^{j}$ (57)

is the polylogarithm [41]. $\mathrm{Li}_{x}(\mu)$, for real $x$ and $\mu$,
converges for $\mu<1$, or for $\mu=1\land x>1$. Then (ii) converges for
$\mu_{d}<1\;\forall x$ or for $\mu_{d}=1\land x>1$. (iv) diverges for
$\mu_{d}>1\;\forall x$ or for $\mu_{d}=1\land x\leq 2$. For it to be negative,
$a<0$ is needed. Therefore, (ii) and (iv) are both satisfied for
$\mu_{d}=1,\;x\in(1,2]$ and $a<0$.

These satisfy (i) and (iii) as well. Taking $\varsigma^{(1)}$ and
$\varsigma^{(2)}$ to be positive, the only negative contribution in (i) comes
from the series, which converges since (ii) is satisfied. With $|a|$ small
enough, this negative term is smaller than others, hence (i) is satisfied.
Finally, as $a<0$ and $\mu>0$, the series in (iii) is necessarily positive,
hence (iii) is satisfied.

Therefore, it is possible to find real constants $\varsigma^{(i)}$ satisfying
conditions (i-iv). Whether a molecule with such constants exists remains an
open question. It must also be noted that, at Planckian scales, the average
number of molecules is small, thus, the validity of Eq. (47) in this limit is
questionable in the first place. Moreover, as noted at the beginning of the
subsection, the molecular model itself is vague at such small scales.
Nevertheless, this result hints at the potential of CST to express new physics
and provides insights for further research. As computing curvature corrections
of arbitrary order seems an overwhelmingly hard task, numerical studies on the
behaviour of $\partial S/\partial M$ in the small $M$ limit are suggested.

##  VIII Conclusion

This paper shows for the first time a numerical study of the entropy of a
Schwarzschild black hole in the framework of Causal Set Theory, an approach to
quantum gravity which only adopts the assumptions of

  1. (I)

spacetime discreteness,

  2. (II)

transitive and acyclic causal relations,

  3. (III)

Poissonian correspondence between discrete and continuum.

We simulated causal sets in Schwarzschild spacetime following the procedure
outlined in section III.2.2, which corrects the algorithms provided by the
original paper [3]. In order to efficiently simulate causal sets we created a
highly parallelized computational framework in C++, described in Sec. V, which
allowed us to generate causets in non-conformally flat Schwarzschild spacetime
orders of magnitude larger than any previous simulation.

We showed for the first time that the minimal assumptions of Causal Set Theory
are enough to provide a molecular model of black hole entropy where the
horizon molecules

  1. (i)

are localised close to the horizon within $2\ell$ distance,

  2. (ii)

give a black hole entropy proportional to the horizon area in the
$A\gg\ell^{2}$ limit which agrees with the Bekenstein-Hawking formula,

  3. (iii)

imply a discreteness length $\ell$ of Planckian order $l_{p}\sim 10^{-35}$m.

Furthermore, we showed that the numerical scaling of Link molecules agrees
perfectly with the analytical result obtained by Barton et al. [27] in the
flat spacetime approximation. Oppositely, the 1st order curvature corrections,
which we analytically computed in section IV, were found to disagree with our
numerical estimate. Specifically, the analytical result lies within the
uncertainties of data points, but it is systematically lower. However, we note
this discrepancy does not impact the aforementioned conclusions of our study,
especially for black holes of physically relevant sizes, where the curvature
correction is negligible.

Lastly, we argued that the horizon molecule approach may yield a finite
temperature cut-off, and even prevent full black hole evaporation for black
holes of Planckian size. Whether a horizon molecule can be envisaged with
these properties remains an open challenge. Moreover, we warn that these
results need to be taken with caution, as the validity of the considered
molecular model is likely to break down at sizes close to the Planck scale,
especially since quantum effects are yet to be introduced. Nevertheless, these
conclusions are intriguing and may offer valuable insights.

There still lies a myriad of exciting unanswered questions about black hole
thermodynamics and causal sets in non-conformally flat spacetimes. Given the
simulation framework now publicly available on GitHub, we hope to have opened
the door for future causal set studies in curved spacetime, possibly
discretising more complex black holes, such as Riessner-Nordstrom and Kerr,
investigating other entropy models [42, 43], and exploring other curvature
features [44, 45].

###### Acknowledgements.

We would like to thank Prof Fay Dowker, who developed the idea of the
Boltzmannian state counting in causal sets theory, for her guidance, the fun
discussions, and the extensive feedback on the manuscript. We would also like
to thank Hinchi Lee, Narin Tananitaporn, Rohith Mukkapati and Olivia Ho for
motivating us and whose ideas contributed to the success of our project.

##  Appendix A Considerations on causality in (3+1)D

In this appendix, we give mathematical proof of the final statements of
Section III.2.2. There we concluded that we can fit

$\frac{d\phi}{du}=\pm_{du}\left(\eta^{2}+2Mu^{3}-u^{2}\right)^{-1/2}$ (58)

against the spatial coordinates of the events to obtain $\pm\eta$, and then
use

$\frac{dt^{*}}{du}=\frac{1}{u^{2}(2Mu-1)}\left[\frac{\pm_{du}\eta}{\sqrt{\eta^{2}+u^{2}(2Mu-1)}}+2Mu\right]$
(59)

to find time bounds on causality, where

$\pm_{du}=\mathrm{sign}(u_{2}-u_{1}).$ (60)

Making the signs of $\eta$ and $dr$ explicit, and re-converting the above
equations in terms of $r=1/u$, we have

$d\phi=\frac{\pm_{dr}}{r^{2}}\left[\eta^{2}-\left(1-\frac{2M}{r}\right)r^{-2}\right]^{-1/2}dr,$
(61)

and

$dt^{*}=\frac{|dr|}{\frac{2M}{r}-1}\cdot\begin{cases}\frac{|\eta|}{\sqrt{\eta^{2}+(\frac{2M}{r}-1)r^{-2}}}-(\pm_{dr}\frac{2M}{r}),&\text{$\eta>0$}\\\
\\\
\frac{-|\eta|}{\sqrt{\eta^{2}+(\frac{2M}{r}-1)r^{-2}}}-(\pm_{dr}\frac{2M}{r}),&\text{$\eta<0$}\end{cases}$
(62)

Next, we show which bounds are to be considered based on where the events lie
with respect to the horizon, and show that when fitting Eq. (58) we require
$\phi\in[0,2\pi)$.

###  A.1 Both points inside the horizon

When both points are inside the horizon, $2M>r$. As expected, our equations
force any object inside the horizon to always keep falling towards the
singularity. In fact, consider a movement opposite to the fall, with
$r_{2}>r_{1}$, hence $dr>0$. We impose $dt^{*}>0$, as $t_{2}^{*}>t_{1}^{*}$.
From Eq. (62), as $2M>r$, this requires

$\frac{\eta}{\sqrt{\eta^{2}+(2M/r-1)r^{-2}}}>\frac{2M}{r}.$ (63)

Trivially, this has no solution for $\eta<0$. For positive $\eta=|\eta|$,
substituting $2M/r=1+\delta$, with $\delta>0$, and rearranging, this yields

$|\eta|>(1+\delta)|\eta|\sqrt{1+\frac{\delta}{\eta^{2}r^{2}}},$ (64)

which also has no solution: objects inside the horizon keep falling in.

Therefore, we must have $dr<0$, hence $\mathrm{sign}(r_{2}-r_{1})<0$. Then,
Eq. (62) shows that the negative and positive values of $\eta$ give,
respectively, the smallest and largest times $t_{-}^{*}$ and $t_{+}^{*}$ of
intersection with the spatial coordinates of $E_{2}$. Therefore, the events
are related if $t_{-}^{*}\leq t_{2}^{*}\leq t_{+}^{*}$, similarly to the 2D
case. He and Rideout [3] do not mention that both signs have to be considered,
again neglecting the upper bound.

###  A.2 One or both points outside the horizon

If both points are outside the horizon, substituting $2M/r=\epsilon<1$ in Eq.
(62)

$dt^{*}=\frac{|dr|}{1-\epsilon}\cdot\begin{cases}\frac{-|\eta|}{\sqrt{\eta^{2}-(1-\epsilon)r^{-2}}}+\mathrm{sign}(r_{2}-r_{1})\epsilon&\text{if
$\eta>0$}\\\ \\\
\frac{+|\eta|}{\sqrt{\eta^{2}-(1-\epsilon)r^{-2}}}+\mathrm{sign}(r_{2}-r_{1})\epsilon&\text{if
$\eta<0$}\end{cases}.$ (65)

As $\epsilon<1$, we see that

$\frac{|\eta|}{\sqrt{\eta^{2}-(1-\epsilon)r^{-2}}}>1>\epsilon.$ (66)

Then, $\eta>0$ gives a negative $dt^{*}$ for both values of
$\mathrm{sign}(r_{2}-r_{1})$. This is useless for us, as we impose
$t_{2}^{*}>t_{1}^{*}$ i.e. $dt^{*}>0$. Therefore, we use the negative $\eta$.
This yields the earliest geodesic reaching $E_{2}$’s spatial coordinates.
Thus, $E_{1}$ and $E_{2}$ are connected if $t_{2}^{*}>t^{*}$, where $t^{*}$ is
the solution to the differential equation.

As in 2D, this result also applies to the case when only $E_{1}$ is outside
the horizon, because outside the horizon an object can remain at rest for an
arbitrarily long period of time, hence no upper bound applies.

###  A.3 The Optimal $\phi$

Null geodesics might travel an angular distance $\phi_{2}+2k\pi$, matching the
spatial coordinates of the events by wrapping around the singularity $k$
times. Fitting Eq. (61), a set $\\{\eta^{2}_{k}\\}$ of suitable parameters
exists, one for each number $k$ of wrappings. We show that the best $\eta^{2}$
corresponds to the 0 wrappings case, i.e. fitting Eq. (61) to $\phi_{2}$.

First note that with increasing $k$, as $|\Delta r|$ stays the same whereas
$|\Delta\phi|$ increases, a larger $|d\phi/dr|$ is required. This implies,
from Eq. (61), $|\eta_{k}|$ decreases with $k$. Then, let’s write Eq. (62) as

$dt^{*}=\left[\;a^{-1}(r)\mathrm{sign}(\eta)f\left(|\eta|;r\right)+b(r)\;\right]|dr|,$
(67)

where $a(r)=(2M/r-1)$ and

$f\left(|\eta|;r\right)=\frac{|\eta|}{\sqrt{\eta^{2}+r^{-2}a(r)}}=\frac{1}{\sqrt{1+r^{-2}a(r)/\eta^{2}}}.$
(68)

For both events inside the horizon, $a(r)>0$, hence $f$ monotonically
increases with $|\eta|$, and so does $a(r)f$. Then the largest $|\eta|$ yields
the earliest ($\eta<0$) and latest ($\eta>0$) null geodesics, thus setting the
boundaries of the light cone. This is given by $k=0$ wrappings.

For both events outside the horizon, $a(r)<0$, hence $f$ monotonically
decreases with $|\eta|$, hence $a(r)f$ monotonically increases. Thus, the
largest $|\eta|$ possible -corresponding to 0 wrappings, again sets the
earliest geodesic. Recall, outside the horizon $\eta<0$ is the only relevant
solution.

##  Appendix B Curvature corrections in Schwarzschild spacetime

In Sec. IV.1, we derived the analytic scaling for the number of Link molecules
with horizon area in the flat spacetime approximation. In this appendix, we
derive the 1st order curvature correction. Therefore, Eq. (35) is the
0${}^{\text{th}}$ order approximation to the full expression in Schwarzschild
spacetime given by

$\langle
N_{L}^{BH}\rangle=a_{L}\frac{A}{\ell^{2}}=\left(a_{L}^{(0)}+a_{L}^{(1)}\frac{\ell}{\sqrt{A}}+\mathcal{O}\left(\frac{\ell^{2}}{A}\right)\right)\frac{A}{\ell^{2}}.$
(69)

Barton et. al [27] outlined a general procedure to derive the
1${}^{\text{st}}$ order curvature correction term $a_{L}^{(1)}$ for arbitrary
horizon. Here we derive it for the Schwarzschild spacetime. The final result
of this derivation was shown in Eq. 39.

###  B.1 Preliminaries

The horizon $\mathcal{H}$ of the Schwarzschild black hole and the spacelike
hypersurface $\Sigma$ are in Eddington-Finkelstein-original coordinates
defined as

$\mathcal{H}(x^{\mu}):r=r_{S}=2M,\quad\quad\Sigma(x^{\mu}):t^{*}=0.$ (70)

We denote the intersection of the two hypersurfaces as $\mathcal{J}$ and it
represents the event horizon at time $t^{*}=0$. Hence, evaluating tensors at
$\mathcal{J}$ implies setting $r=r_{S}=2M$ and $t^{*}=0$. Moreover, for
completeness the metric $g_{\mu\nu}$ and its inverse $g^{\mu\nu}$ in EFO
coordinates are

$g_{\mu\nu}=\begin{pmatrix}-\left(1-\frac{2M}{r}\right)&\frac{2M}{r}&0&0\\\
\frac{2M}{r}&\left(1+\frac{2M}{r}\right)&0&0\\\ 0&0&{r^{2}}&0\\\
0&0&0&{r^{2}\sin^{2}{\theta}}\end{pmatrix},$ (71)

and

$g^{\mu\nu}=\begin{pmatrix}-\left(1+\frac{2M}{r}\right)&\frac{2M}{r}&0&0\\\
\frac{2M}{r}&\left(1-\frac{2M}{r}\right)&0&0\\\ 0&0&\frac{1}{r^{2}}&0\\\
0&0&0&\frac{1}{r^{2}\sin^{2}{\theta}}\end{pmatrix}.$ (72)

The 1${}^{\text{st}}$ order corrections must depend on mutually independent
invariant scalars of dimension $\ell^{-1}$, evaluated at $\mathcal{J}$ and
depend on the geometry of the problem at hand. The only such quantities are
$K$, the trace of the extrinsic curvature of $\Sigma$, its component $K_{m}$
tangential to $\Sigma$ an orthogonal to $\mathcal{J}$, and the null expansion
$\vartheta$ of the horizon. It can be shown that the 1${}^{\text{st}}$ order
correction for links in 3+1 dimension is [27]

$\frac{a_{L}^{(1)}}{\sqrt{A}}\approx-(0.036K+0.088K_{m}+0.021\vartheta).$ (73)

###  B.2 Determining the required geometrical objects

To determine the invariant scalars we first require finding the basic
geometrical objects from which they are constructed. The first one is the
metric, given by Eq. (71).

Then we also need $n^{\mu}$ and $k^{\mu}$, normal vectors to the hypersurfaces
$\Sigma$ and $\mathcal{H}$, respectively. Vector $n^{\mu}$ is future-pointing
($n^{0}>0$), normalised as $n^{\mu}n_{\mu}=-1$ and is given by

$n_{\mu}=a\>\partial_{\mu}\Sigma=a(1,0,0,0),$ (74)

Correspondingly, we have

$n^{\mu}=g^{\mu\nu}n_{\nu}=a\left(-1-\frac{2M}{r},\frac{2M}{r},0,0\right).$
(75)

The previous conditions yield

$a=\frac{-1}{\sqrt{\frac{2M}{r}+1}}.$ (76)

Similarly, $k^{\mu}$ is also future-pointing and is given by

$k_{\mu}=b\>\partial_{\mu}\mathcal{H}=b(0,1,0,0).$ (77)

Since it is normalised such that $[n^{\mu}k_{\mu}]_{\mathcal{J}}=-1/\sqrt{2}$,
where $[\,\cdot\,]_{\mathcal{J}}$ denotes that the quantity is evaluated at
$\mathcal{J}$, this immediately sets $b=1$.

Lastly, we require the vector $m^{\mu}=\sqrt{2}k^{\mu}-n^{\mu}$, which is
tangent to $\Sigma$ and orthogonal to $\mathcal{J}$.

###  B.3 Local geometric invariants

Next, we introduce the “projector” tensors

$h^{\alpha}_{\;\;\beta}=\delta^{\alpha}_{\beta}+n^{\alpha}n_{\beta}$ (78)

and

$\sigma^{\alpha}_{\;\;\beta}=h^{\alpha}_{\;\;\beta}-m^{\alpha}m_{\beta}.$ (79)

Noting that extrinsic curvature is defined as

$K_{\mu\nu}=(\nabla_{\beta}n_{\alpha})h^{\alpha}_{\;\;\mu}h^{\beta}_{\;\;\nu},$
(80)

we can therefore determine the mutually independent invariant scalars
evaluated on $\mathcal{J}$ that are required to obtain the curvature
correction. The relevant Christoffel symbols needed for the computation are
given in App. C. The trace $K$ of the extrinsic curvature is given by

$K(\mathcal{J})=[g^{\mu\nu}K_{\mu\nu}]_{\mathcal{J}},$ (81)

its component along the $m$ direction by

$K_{m}(\mathcal{J})=[K_{\mu\nu}m^{\mu}m^{\nu}]_{\mathcal{J}},$ (82)

and the null expansion of the horizon as

$\vartheta=[(\nabla_{\beta}k_{\alpha})\sigma^{\alpha\beta}]_{\mathcal{J}}.$
(83)

After determining the required Christoffel symbols in EFO coordinates and
other required calculations, we obtain that $\vartheta=0$, as expected since
the Schwarzschild black hole is static. The extrinsic curvature is found to be

$[K_{\mu\nu}]_{\mathcal{J}}=\frac{1}{32}\begin{pmatrix}\frac{3\sqrt{2}}{M}&\frac{6\sqrt{2}}{M}&0&0\\\
\frac{6\sqrt{2}}{M}&\frac{12\sqrt{2}}{M}&0&0\\\ 0&0&-32\sqrt{2}M&0\\\
0&0&0&-32\sqrt{2}M\sin^{2}{\theta}\end{pmatrix},$ (84)

its trace is given by

$K=-\frac{5\sqrt{2}}{16}\frac{1}{M}\approx-\frac{0.4419}{M},$ (85)

and the component along $m^{\mu}$ is

$K_{m}=\frac{3\sqrt{2}}{16}\frac{1}{M}\approx\frac{0.2651}{M}.$ (86)

Inserting the obtained $K$, $K_{m}$ and $\vartheta$ into Eq. (37) and noting
that $\sqrt{A}=4\sqrt{\pi}M$ yields the 1${}^{\text{st}}$ order curvature
correction term

$a_{L}^{(1)}\approx-0.0558.$ (87)

Then, the average expected number of Links in Schwarzschild spacetime, as a
function of the horizon area $A$ in terms of the discreteness scale $\ell$, is
given by

$\langle
N_{L}^{BH}\rangle=a_{L}\frac{A}{\ell^{2}}=a_{L}^{(0)}\frac{A}{\ell^{2}}\left(1-0.322\frac{\ell}{\sqrt{A}}+\mathcal{O}(\ell^{2})\right),$
(88)

where $a_{L}^{(0)}=\sqrt{3}/10$ is the $0^{\text{th}}$ order flat spacetime
term. The correction term becomes negligible when $A$ increases to the size of
ordinary black holes since $\ell\sim\ell_{p}$ is infinitesimal.

((a))

((b))

((c))

Figure 7: Sprinkled causet elements in (2+1)D Minkowski spacetime lie within
the shaded spacetime regions. The simulation framework allows for sprinkling
into ball, bicone and cylinder shapes (left to right) in $D\in\\{2,3,4\\}$
spacetime.

##  Appendix C Christoffel symbols in the Eddington-Finkelstein Original
coordinates

We computed the Christoffel symbols of the EFO coordinates introduced in
Section III.1 using SymPy and EinsteinPy libraries [46, 47]. The non-zero ones
are

$\displaystyle\Gamma^{t^{*}}_{{t^{*}}{t^{*}}}$
$\displaystyle=\frac{2M^{2}}{r^{3}},$
$\displaystyle\Gamma^{t^{*}}_{{t^{*}}r}=\Gamma^{t^{*}}_{r{t^{*}}}$
$\displaystyle=\frac{M(2M+r)}{r^{3}},$ $\displaystyle\Gamma^{t^{*}}_{rr}$
$\displaystyle=\frac{2M(M+r)}{r^{3}},$
$\displaystyle\Gamma^{t^{*}}_{\theta\theta}$ $\displaystyle=-2M,$
$\displaystyle\Gamma^{t^{*}}_{\phi\phi}$
$\displaystyle=-2M\;\sin^{2}{\theta},$
$\displaystyle\Gamma^{r}_{{t^{*}}{t^{*}}}$
$\displaystyle=\frac{-M(2M-r)}{r^{3}},$
$\displaystyle\Gamma^{r}_{{t^{*}}r}=\Gamma^{r}_{r{t^{*}}}$
$\displaystyle=-\frac{2M^{2}}{r^{3}},$ $\displaystyle\Gamma^{r}_{rr}$
$\displaystyle=\frac{-M(2M+r)}{r^{3}},$
$\displaystyle\Gamma^{r}_{\theta\theta}$ $\displaystyle=2M-r,$
$\displaystyle\Gamma^{r}_{\phi\phi}$ $\displaystyle=(2M-r)\;\sin^{2}{\theta},$
$\displaystyle\Gamma^{\theta}_{r\theta}=\Gamma^{\theta}_{\theta r}$
$\displaystyle=\frac{1}{r},$ $\displaystyle\Gamma^{\theta}_{\phi\phi}$
$\displaystyle=-\sin{\theta}\cos{\theta},$
$\displaystyle\Gamma^{\phi}_{r\phi}=\Gamma^{\phi}_{\phi r}$
$\displaystyle=\frac{1}{r},$
$\displaystyle\Gamma^{\phi}_{\theta\phi}=\Gamma^{\phi}_{\phi\theta}$
$\displaystyle=\frac{1}{\tan{\theta}}.$

##  Appendix D Other tests of the causal set simulation framework

This appendix describes other more basic and less relevant tests we performed,
besides the ones from Sec. V, to ensure our simulation framework works as
intended for studying Schwarzschild causets.

Plotting the sprinkling region as a shaded shape, we confirmed that sprinkling
into different shapes works as intended and that causets only lie within the
required boundaries (see Fig. 7).

To further check that the causal structure in Minkowski spacetime works
correctly we implemented a Myrheim-Mayer dimension estimator [48, 49], which
converges to the dimension of the spacetime that the causet is embedded in, as
the size $N$ increases. Simulating causal sets up to $N=32768$, sprinkling
into bicones in the dimensions $D=2,3,4$, we see that the dimension estimate
agrees with the true value within the uncertainty, giving substantial evidence
that our simulation framework works as intended (see Fig. 8).

Figure 8: Myrheim-Mayer dimension estimates converge to the dimensions of the
Minkowski spacetime for $D=2,3,4$. Each data point is averaged over 200
causets.

((a))

((b))

\begin{overpic}[width=177.78578pt]{AAA_n_lambda3D_prob_distribution_expx_logy_withcloud.pdf}
\put(-136.0,88.0){(a)} \put(-2.0,88.0){(b)} \end{overpic}

Figure 9: (a) Counting of Link molecules from the simulations of causets in
(3+1)D Schwarzschild spacetime, with simulations boundaries at distance
$6\ell$ from relevant hypersurfaces $\Sigma$ and $\mathcal{H}$. Each datapoint
corresponds to 100 simulations. The red line, with values from Eq. 89, is our
best fit. The azure and green lines are the analytical models: respectively,
the flat-spacetime approximation and the 1st-order curvature expansion. (b)
$\Lambda_{n}$ molecules distribution in (2+1)D spacetime, from the simulation
of 50 causets for each of 9 different area values in the range
$A\in[500,4500]\ell^{2}$.

##  Appendix E Other numerical results

In this appendix, we provide additional numerical results to the ones in Sec.
VI.

In the first place, we simulated causets in (3+1)D Schwarzschild spacetime in
a region with bounds $t^{*}\in[-6,0],\,r\in[r_{S}-6,\,r_{S}+6]$ in terms of
discreteness units $\ell$, looser than the ones used for Sec. VI.

The obtained results (see Fig. 9(a)) are equivalent to those in Sec. VI.
Firstly, the analytical model lies within uncertainty, but is consistently
lower than datapoints. Secondly, the best fit is

$a^{(0)}_{L[\text{3}]}=0.173\pm 0.005,\qquad a^{(1)}_{L[\text{3}]}=0.4\pm
0.4,$ (89)

i.e. the same as in Eq. 40, but with larger uncertainties. This provides an
additional confirmation of the validity of the simulations in the main text.

Furthermore, we performed simulations in (2+1) dimensions, to investigate how
the probability distribution $p_{n}$ of the horizon molecules varies in other
dimensions. The sprinkled region was bounded as the (3+1)D one in the main
text: $t^{*}\in[-4,0],\,r\in[r_{S}-3,\,r_{S}+3]$ in terms of discreteness
units $\ell$. We simulated 50 causets for each of 9 different equally-spaced
horizon areas in the range $A\in[500,4500]\ell^{2}$, with average
cardinalities ranging roughly from $50,000$ to $450,000$. We found that the
molecules obey the same distribution as in (3+1)D,
$p_{n}=\left(e^{\chi}-1\right)e^{-n\chi}$, with approximately the same
coefficient $\chi=1.50\pm 0.03$ (see Fig. 9(b)). Therefore, the hypothesis
outlined in Sec. VII.1 that horizon molecules could satisfy some thermodynamic
distribution is to be discarded.

## References

  * Hawking [1975] S. W. Hawking, Particle creation by black holes, Commun. Math. Phys. 43, 199–220 (1975).
  * Dou and Sorkin [2003] D. Dou and R. D. Sorkin, Black-hole entropy as causal links, Foundations of Physics 33, 279 (2003).
  * He and Rideout [2009] S. He and D. Rideout, A causal set black hole, Class. Quantum Grav. 26, 125015 (2009).
  * Hawking and Penrose [1970] S. W. Hawking and R. Penrose, The singularities of gravitational collapse and cosmology, Proc. R. Soc. Lond. A 314, 529 (1970).
  * Bardeen  _et al._ [1973] J. M. Bardeen, B. Carter, and S. W. Hawking, The four laws of black hole mechanics, Commun. Math. Phys. 31, 161 (1973).
  * Bekenstein [1973a] J. D. Bekenstein, Black holes and entropy, Phys. Rev. D 7, 2333 (1973a).
  * Israel [1967] W. Israel, Event horizons in static vacuum space-times, Phys. Rev. 164, 1776 (1967).
  * Israel [1968] W. Israel, Event horizons in static electrovac space-times, Commun. Math. Phys. 8, 245 (1968).
  * Carter [1971] B. Carter, Axisymmetric black hole has only two degrees of freedom, Phys. Rev. Lett. 26, 331 (1971).
  * Wald [1971] R. M. Wald, Final states of gravitational collapse, Phys. Rev. Lett. 26, 1653 (1971).
  * Robinson [1975] D. C. Robinson, Uniqueness of the kerr black hole, Phys. Rev. Lett. 34, 905 (1975).
  * Mazur [1982] P. O. Mazur, Proof of uniqueness of the kerr-newman black hole solution, J. Phys. A: Math. Gen. 15, 3173 (1982).
  * Jacobson  _et al._ [2005] T. Jacobson, D. Marolf, and C. Rovelli, Black hole entropy: Inside or out?, International Journal of Theoretical Physics 44, 1807–1837 (2005).
  * Bombelli  _et al._ [1986] L. Bombelli, R. K. Koul, J. Lee, and R. D. Sorkin, Quantum source of entropy for black holes, Phys. Rev. D 34, 373 (1986).
  * Bekenstein [1973b] J. D. Bekenstein, Black holes and the second law, Lettere al Nuovo Cimento 4, 737 (1973b).
  * Bombelli  _et al._ [1987] L. Bombelli, J. Lee, D. Meyer, and R. D. Sorkin, Space-time as a causal set, Phys. Rev. Lett. 59, 521 (1987).
  * Hawking [1976] S. W. Hawking, Breakdown of predictability in gravitational collapse, Phys. Rev. D 14, 2460 (1976).
  * Malament [1977] D. Malament, The class of continuous timelike curves determines the topology of spacetime, J. Math. Phys. 18, 1399– (1977).
  * Bombelli [1987] L. Bombelli,  _Space-time as a causal set_ , Ph.D. thesis, Syracuse University (1987).
  * Saravani and Aslanbeigi [2014] M. Saravani and S. Aslanbeigi, On the causal set–continuum correspondence, Classical and Quantum Gravity 31, 205013 (2014).
  * Dowker [2005] F. Dowker, Causal sets and the deep structure of spacetime, in  _100 Years of Relativity - Space-time Structure: Einstein and Beyond_, edited by A. Ashtekar (World Scientific, [Online], 2005) pp. 445–464.
  * Bombelli  _et al._ [2009] L. Bombelli, J. Henson, and R. D. Sorkin, Discreteness without symmetry breaking: a theorem, Mod. Phys. Lett. A 24, 2579 (2009).
  * Stoyan  _et al._ [ap 3] D. Stoyan, S. Chiu, W. Kendall, and J. Mecke,  _Stochastic geometry and its applications._, 3rd ed. (Wiley, Chichester, 2013. Chap. 3.) Chap. 3, pp. 35–63.
  * Surya [2019] S. Surya, The causal set approach to quantum gravity, Living Reviews in Relativity 22, 10.1007/s41114-019-0023-1 (2019).
  * Eddington [1924] A. Eddington, A comparison of whitehead’s and einstein’s formulae, Nature 113, 192 (1924).
  * Finkelstein [1958] D. Finkelstein, Past-future asymmetry of the gravitational field of a point particle, Phys. Rev. 110, 965 (1958).
  * Barton  _et al._ [2019] C. Barton, A. Counsell, F. Dowker, D. S. W. Gould, G. Taylor, and I. Jubb, Horizon molecules in causal set theory, Phys. Rev. D 100, 126008 (2019).
  * Minz [2020] C. Minz, Python-causets (2020).
  * Tananitaporn [2023] N. Tananitaporn, Black hole entropy in causal sets theory using horizon molecules (2023), MSci project report. Imperial College London.
  * Banks  _et al._ [1984] T. Banks, L. Susskind, and M. E. Peskin, Difficulties for the evolution of pure states into mixed states, Nucl. Phys. B 244, 125 (1984).
  * Lowe and Thorlacius [1999] D. A. Lowe and L. Thorlacius, Ads-cft and the information paradox, Phys. Rev. D 60, 104012 (1999).
  * Hawking [2005] S. W. Hawking, Information loss in black holes, Phys. Rev. D 72, 084013 (2005).
  * Hawking [2015] S. W. Hawking, The information paradox for black holes, arXiv 10.48550/arXiv.1509.01147 (2015).
  * Hawking  _et al._ [2016] S. W. Hawking, M. J. Perry, and A. Strominger, Soft hair on black holes, Phys. Rev. Lett. 116, 231301 (2016).
  * ’t Hooft [1993] G. ’t Hooft, Dimensional reduction in quantum gravity, arXiv 10.48550/arXiv.gr-qc/9310026 (1993).
  * Susskind [1995] L. Susskind, The world as a hologram, J. Math. Phys. 36, 6377 (1995).
  * Thompson  _et al._ [2019] T. A. Thompson, C. S. Kochanek, K. Z. Stanek, C. Badenes, R. S. Post, T. Jayasinghe, D. W. Latham, A. Bieryla, G. A. Esquerdo, P. Berlin, _et al._ , A noninteracting low-mass black hole-giant star binary system, Science 366, 637 (2019).
  * Fixsen [2009] D. J. Fixsen, The temperature of the cosmic microwave background, The Astrophysical Journal 707, 916 (2009).
  * Hawking [1971] S. Hawking, Gravitationally collapsed objects of very low mass, Mon. Not. R. Astron. Soc. 152, 75 (1971).
  * Röken [2013] C. Röken, On the nature of black holes in loop quantum gravity, Class. Quantum Grav. 30, 015005 (2013).
  * Lewin [ap 7] L. Lewin,  _Polylogarithms and associated functions_ , 1st ed. (North-Holland Publishing Co., New York, 1981. Chap. 7.).
  * Benincasa [2013] D. M. T. Benincasa,  _The action of a causal set_ , Ph.D. thesis, Imperial College London (2013), chap. 6.
  * Sorkin and Yazdi [2018] R. D. Sorkin and Y. K. Yazdi, Entanglement entropy in causal set theory, Class. Quantum Grav. 35, 074004 (2018).
  * Benincasa and Dowker [2010] D. M. T. Benincasa and F. Dowker, Scalar curvature of a causal set, Phys. Rev. Lett. 104, 181301 (2010).
  * Roy  _et al._ [2013] M. Roy, D. Sinha, and S. Surya, Discrete geometry of a small causal diamond, Phys. Rev. D 87, 044046 (2013).
  * Meurer  _et al._ [2017] A. Meurer, C. P. Smith, M. Paprocki, O. Certik, S. B. Kirpichev, M. Rocklin, A. Kumar, S. Ivanov, J. K. Moore, S. Singh, _et al._ , Sympy: symbolic computing in python, PeerJ Computer Science 3, 10.7717/peerj-cs.103 (2017).
  * EinsteinPy Development Team [2021] EinsteinPy Development Team,  _EinsteinPy: Python library for General Relativity_ (2021).
  * Myrheim [1978] J. Myrheim, Statistical geometry (1978).
  * Meyer [1988] D. A. Meyer,  _The dimension of causal sets_ , Ph.D. thesis, MIT (1988).