# OpenEMS: an open-source Package for Two-Stage Stochastic and Robust
Optimization for Ambulance Location and Routing with Applications to Austin-
Travis County EMS Data

Joshua Onga, David Kulpanowski b, Yangxinyu Xie c, Evdokia Nikolova a, Ngoc
Mai Tran d Correspondence: Joshua Ong Email<EMAIL_ADDRESS>aDepartment of
Electrical and Computer Engineering, University of Texas at Austin, TX 78712,
USA.; bDepartment of Emergency Medical Services, City of Austin; cDepartment
of Computer Science, University of Texas at Austin, TX 78712, USA.;
dDepartment of Mathematics, University of Texas at Austin, TX 78712, USA.

###### Abstract

Emergency Medical Systems (EMS) provide crucial pre-hospital care and
transportation. Faster EMS response time provides quicker pre-hospital care
and thus increases survival rate. We reduce response time by providing optimal
ambulance stationing and routing decisions by solving two stage stochastic and
robust linear programs. Although operational research on ambulance systems is
decades old, there is little open-source code and consistency in simulations.
We begin to bridge this gap by publishing OpenEMS, in collaboration with the
Austin-Travis County EMS (ATCEMS) in Texas, an end-to-end pipeline to optimize
ambulance strategic decisions. It includes data handling, optimization, and a
calibrated simulation. We hope this open source framework will foster future
research with and for EMS. Finally, we provide a detailed case study on the
city of Austin, Texas. We find that optimal stationing would increase response
time by 88.02 seconds. Further, we design optimal strategies in the case where
Austin EMS must permanently add or remove one ambulance from their fleet.

###### keywords:

OR in Health Systems; Emergency medical systems; Ambulance location; Ambulance
dispatching; Data-driven optimization

##  1 Introduction

The Emergency Medical Service (EMS) provides pre-hospital treatment and
transportation. Receiving fast pre-hospital treatment is vital in scenarios
such as cardiac arrest, stroke, and other severe trauma cases [29]. It is
generally agreed upon that these calls need pre-hospital treatment within 8-10
minutes [36, 12, 7, 25]. If a call is not answered with in this threshold, it
is considered a shortfall. To help EMS services achieve the best response time
the operational research community has widely studied how ambulance stationing
and routing can optimize response times. [7, 6, 27, 34, 11, 35, 39, 32, 22,
24].

Strategic decisions can also help EMS departments save money. White et al.[41]
found health care spending has dramatically increased in the US compared to
other countries. Further, Heightman [20] approximates that it costs half a
million dollars a year to staff one ambulance on a 24/7 basis. One way to
combat increased spending is to use current EMS resources more efficiently.
Optimization could help a city with a decreasing budget maintain performance
while removing an ambulance from its operation.

The operational research community commonly constructs the ambulance problem
into a two-stage optimization problem [7, 6, 27]. The first stage is the
ambulance location problem (ALP). This answers the question: where should we
station our ambulances around the city to achieve an optimal coverage?
Ambulances can be stationed around the city in locations such as hospitals,
fire stations, ect. In the ALP we want to maximize the coverage of the
ambulances to respond to calls, using historic call data as empirical
probability distributions.

###  1.1 Our Contributions

There have been many case studies on ambulance stationing and routing [34, 11,
35, 39, 32, 22, 24], but the methodologies regarding discrete-event
simulations vary widely. Some are theoretical and include no simulations, some
use private softwares, while only a few provide open-source tools [7]. In this
work, we develop OpenEMS, a modular software with open-source code consisting
of four well-separated parts: data handling, optimization, numerical
computations of discrete time events, and visualization. This allows future
research that innovates on one or several components to be directly compared
with previously published work. We hope that this increased transparency would
facilitate the translation from research results to field implementations by
EMS departments. Our main contributions are as follows:

  1. 1.

A software pipeline at https://github.com/joshua-ong/AmbulanceDeployment that
can be applied to any city. The software includes: discretizing a set area
into a grid, calibration, validation, optimization using state-of-the-art
(SOTA) methods from Bertsimas et al.[7], simulation, and visualizing
improvements compared to historical performance.

  2. 2.

A case study of Austin EMS, including 246,809 thousand calls within 1039
square miles, demonstrating the impact of optimal decision making.

The rest of the paper is organized as follows. Section 2 gives a brief
overview of ambulance dispatch research. Next, we include the necessary theory
of the optimization models in Section 3. We detail the simulation model in
Section 4. Next, we have a case study of the city of Austin in section 5.
Next, we include a socio-economic of Austin in section 6. Next, there is a
user code sample in section 7. We conclude in Section 8.

##  2 Past Work

There exist many case studies across various countries showing that EMS
systems’ heuristics are not optimal. Morohosi [34], Eaton et al.[14], and
Burkey et al.[11] use popular coverage linear programs to study Tokyo, Japan;
Austin, Texas; and 4 states in the southern USA respectively. Eaton et
al.[14], and Burkey et al.[11] also use coverage to examine social equality.
Nilsang et al.[35] included a case study in Bangkok, Thailand where they
utilize social media to design a covering model. They found they could reduce
the number of ambulance locations by 81.6 % without reducing coverage. Swalehe
et al.[39] included a case study in Eskisehir Province, Turkey using GIS
(Geogrpahic Information System) methods and found they could reduce the
response time from 10 minutes to 5 minutes for 77.6 % of ambulance demand
areas. McCormack et al.[32] included a case study in London, England using
genetic algorithms and testing the algorithm on a calibrated and validated
simulation model. Doener et al.[13] includes a broader case study of Austria
including coverage linear programs for ambulance systems but also other
related areas like inventory routing and disaster relief. Nonetheless, none of
these papers include open-source code.

These case studies can also answer important hypothetical questions. McCormack
et al.[32] used simulations to find how performance would decrease if there
were resource gaps that day. Ingolfsson et al.[22] tested the policy of a
’single start station’, where all ambulances start at the same station using
the ProModel software [19]. Although the policy of having all ambulances start
at the same station seems inefficient, it has the benefits of having pooled
resources in one location. Inakawa et al.[24] used simulation on the city of
Seto, Japan to find that adding a new location was twice as effective as
adding two ambulance stations. Aboueljinane et al.[3] used simulation in Val-
de-Marne department of France using the ARENA software [18]. In it they
recommended repositioning teams to increase coverage and found that it
improved response time by 3.5$\%$.

Implementing novel algorithms into an EMS department has proven to be
difficult. For example, in 1992, the London Ambulance Service (LAS) tried to
implement a computerized system including a Computerized Dispatch System (CAD)
(rather than human operators), and automatic vehicle stationing [15].
According to Adamu et al.[15] the software had a $\pounds$ 7.5 million budget,
but only lasted 9 days before it was closed due to failure and increased
mortality. The software was abandoned after 9 days as implementations did not
have time to learn and fail, because there are human lives at stake. It
further shows the significance of testing novel algorithms on confident
simulations.

A notable success began with Mason and Henderson’s [21] simulation model
BartSim, which later matured into a company [30] called Optima. Optima
currently has case studies in New Zealand, Australia, Denmark, United
Kingdoms, Canada, and the United States. However, Optima is a commercial
software designed for EMS operators. In contrast, OpenEMS is open-source and
designed to support researchers making the latest advances in the ambulance
allocation and routing problem with data validations while incurring minimal
technical setup and no charge.

Some of the earliest work began when ReVelle et al.[37] proposed the Maximum
Available Location Problem (MALP), inspired by the more general Maximum
Coverage Problem (MCP). In MALP they optimize the location of stations for
ambulances to maximize the coverage, where a caller is covered if there is an
ambulance that can reach the call within a fixed time. This is called the
ambulance location problem (ALP). Later Batta et al.[5] introduced Maximum
Expected Covering Location Problem (MEXCLP), which used an ambulance’s ”busy
fraction” as a Bernoulli random variable to approximate the expected coverage.
Other classical models include the Double Standard Model (DSM) by Gendreau et
al.[16]. These remain the important benchmarks in many case studies [7, 34,
14, 11].

Further research from the operational research community has since been
applied to solving the ambulance problems. Beraldi et al.[6] solved the ALP
with chance constraints to model routing ambulances to an unknown demand. For
more information on two stage stochastic formulation, we refer the reader to
[4] and [9]. A two-stage robust formulation was first implemented by Bertsimas
et al.[7] who applied a data-driven model to generate scenarios. A two-stage
distributional robust formulation was implemented by Liu et al.[27].

Other studies integrate a tool called the approximate dynamic programming
(ADP). Some of the first papers on applying ADP to the ambulance problem were
from Maxwell et al.[31] and Schmid et al.[38]. Later Jagtenberg et al.[23]
used ADP to show that the closest ambulance routing policy is not necessarily
optimal. For a more comprehensive overview, we refer readers to the review by
Aboueljinane et al.[2] and Tassone et al.[40].

##  3 Optimization

In this section, we lay out the theoretical formulation of two-stage
stochastic and robust formulations appeared in Bertsimas et al.[7] to solve
ALP and ARP. We adopt their notation in this paper. In a two-stage
optimization problem, there is the first stage here-and-now decision which
must be made under uncertainty and the wait-and-see variable that happens
afterwards.

In our context, at the first stage we must select how many ambulances are
stationed at each location. Given the set of ambulance stations $I$; let
$x_{i}$ be the number of ambulances at location $i\in I$. This is the first
stage decision, as we must station the ambulances before we know the
realization of the demand. In the second stage the calls are revealed and we
must determine the optimal routing of stationed ambulances to these calls. Let
$J$ be the set of regions, then each call can be assigned to a region $j\in
J$. Let $y_{ij}$ be the routing of ambulance stationed at $i$ to region $j$.
Let $\tau$ be the longest time it takes for an ambulance to serve a call and
return to their station. Then we can group calls within periods of $\tau$
where we assume an ambulance can serve one call in each period. A complete
table of notation is included in table 1.

Name | Description
---|---
$I$ | Set of ambulance locations
$J$ | Set of demand regions
$E$ | Set of feasible edges from $I$ to $J$
$M$ | Set of observed ambulance calls
$n$ | Total number of ambulances
$\tau$ | Total time (minutes) for an ambulance to serve a call and be available again
$y_{ij}$ | Number of ambulances routed from location $i$ to region $j$
$x_{i}$ | Number of ambulances at location $i$
$d_{j}$ | Number of ambulance calls from region $j$
$\bm{B}_{I}$ | Incidence matrix representing incoming edges
$\bm{B}_{J}$ | Incidence matrix representing outgoing edges
Table 1: Summary of notation for ambulance location and dispatch.

###  3.1 Stochastic Formulation

The stochastic formulation objective is to minimize the expectation of
shortfalls. To solve this problem we use Monte-Carlo sampling to construct a
discrete distribution $\hat{\mathbb{P}}(\hat{\bm{d}})$ where we assume each
sample of calls is equally likely. The stochastic formulation is written as:

$\text{minimize}\quad\frac{1}{M}\sum_{m=1}^{M}\bm{z}^{m}$ (1a) $\text{subject
to}\quad\displaystyle\bm{e}^{\top}\bm{x}\leq n$ (1b)
$\quad\bm{x}-\bm{B}_{I}\bm{y}^{m}\geq 0\>\>\forall m\in[M]$ (1c)
$\quad\bm{z}^{m}\geq\bm{d}^{m}-\bm{B}_{J}\bm{y}^{m}\>\>\forall m\in[M]$ (1d)
$\quad\bm{x}\in\mathbb{Z}_{+}^{I},\bm{y}^{m}\in\mathbb{Z}_{+}^{E},\bm{z}^{m}\in\mathbb{Z}_{+}^{J}\>\>\forall
m\in[M]$ (1e)

The objective 1a contains a slack variable that is bounded from below later.
Constraint 1b limits the number of ambulances stationed to be less than or
equal to the total number of ambulances. Constraint 1c iterates through each
Monte-Carlo sample $(m\in M)$, for each sample we limit a stationed ambulance
to only be sent once per time period, $\tau$. Constraint 1d makes the slack
variable greater than or equal to the shortfalls, which is equal to calls
minus calls answered. (Note this can be thought of as a maxflow problem of a
bipartite graph, where $x$ and $d$ are the vertices of the bipartite graph and
$y$ are the edge values.)

###  3.2 Robust Formulation

In the robust optimization formulation, we want to minimize the worst-case
scenarios. To do this we first model demand as a Poisson distribution,
$\gamma$. (Single is the demand of a single grid cell. Local is the demand of
adjacent grid cells. Regional is the demand within 10 minutes of a grid cell.
Global is the demand of all of Austin.) We then construct the uncertainty set
$\mathbb{D}(\alpha)$ as the $\alpha$ Value-at-Risk of these Poisson
distributions. Finally, we take the worst scenarios from and minimize over
them in the robust formulation like:

$\text{maximize}_{\bm{d}\in\mathbb{D}(\alpha)}\text{minimize}\quad\frac{1}{M}\sum_{m=1}^{M}\bm{d}^{m}-\bm{B}_{J}\bm{y}^{m}$
(2a) $\text{subject to}\quad\displaystyle\bm{e}^{\top}\bm{x}\leq n$ (2b)
$\quad\bm{x}-\bm{B}_{I}\bm{y}^{m}\geq 0\>\>\forall m\in[M]$ (2c)
$\quad\bm{z}^{m}\geq\bm{d}^{m}-\bm{B}_{J}\bm{y}^{m}\>\>\forall m\in[M]$ (2d)
$\quad\bm{x}\in\mathbb{Z}_{+}^{I},\bm{y}^{m}\in\mathbb{Z}_{+}^{E},\bm{z}^{m}\in\mathbb{Z}_{+}^{J}\>\>\forall
m\in[M]$ (2e)

where:

$\mathbb{D}(\alpha)=\\{\bm{d}\in\mathbb{Z}^{J}_{+}:d_{j}\leq
VaR(\gamma_{j}^{\text{single}}(\alpha))\>\>\forall j\in
J,VaR(\gamma_{j}^{\text{local}}(\alpha))\>\>\forall j\in J,\\\ d_{j}\leq
VaR(\gamma_{j}^{\text{regional}}(\alpha))\>\>\forall i\in I,d_{j}\leq
VaR(\gamma_{j}^{\text{global}}(\alpha)\\}$ (3)

The formulation is the same as the stochastic formulation except for the
objective. For more information we refer the reader to the original
formulation by Bertsimas et al. [7] on how to dualize 2e and make it tractable
with the column-constraint method. Zeng and Zhao [42] first designed the
column-constrained method for two-stage robust optimization. The column-
constraint method has the same idea as the Benders-style but with faster
heuristic implementation.

##  4 Simulation Model

In this section, we describe the different components of the simulation model
in the general setting. This includes the steps to replicate this on a
different case study. Notably, we measure the confidence of the simulation
response time so that an EMS department can understand the possible variance
in real-world applications. First, we explain the state space of a general
Emergency Medical System (EMS) and the correspondence of these events to the
simulation. Next,we explain how the simulation represents the real world by
discretizing it into a grid. We then verify the model by simulating historical
actions and comparing them with what happened to validate that the model works
in normal conditions. Finally, we explain the process of tuning the linear
program parameters for best results.

###  4.1 The state space for answering an emergency call

A general framework of the EMS procedure for answering an emergency call is
detailed in figure 1. The first state includes from the time the call was made
to when it was assigned to an ambulance, where an operator must collect enough
information like the seriousness of the injury and location. The time the
ambulance is en route to when it arrives at the call is labeled the response
time. This is arguably the most important metric for survival rate of the
patient. From the time the ambulance arrived to when it departs the scene
captures the time it took to administer pre-hospital treatment and load the
patient. The time departed from the scene to the hospital is drive time. The
ambulance is now available to serve another call, if there are no new calls,
it will return to its station.

Figure 1: A general state space of the EMS procedure for answering an
emergency call.

Our simulation model includes the same steps as the state space just
described. Each step is an event that is put into a priority queue with
respect to time. An overview of the code is included in pseudo-code 1. In the
beginning, all calls are input into the priority queue at their respective
times. The calls then propagate the other events following the state space. We
do not include the setup time, because it is only a fraction of the whole
ambulance travel time and is not impacted by our strategic decision. The
travel times are determined using Open Street Map, which was introduced by
Haklay and Weber [17]. The intervention time is modeled as a lognormal
distribution. They have been shown to effectively fit ambulance processes like
in Lam et al.[26] and McLay and Mayorga [33].

priorityQueue.enque(allEmergencyCalls)

while  _priorityQueue not empty_ do

event = priorityQueue.pop()

event.handle()

if  _event == newCall_ then

priorityQueue.push(callEnroute)

end if

if  _event == callEnroute_ then

priorityQueue.push(callArriveScene)

end if

if  _event == callArriveScene_ then

priorityQueue.push(callDepartScene)

end if

if  _event == callDepartScene_ then

priorityQueue.push(callArriveHospital)

end if

if  _event == callArriveHospital_ then

priorityQueue.push(AmbulanceAvailable)

end if

end while

Algorithm 1 An overview of the simulation model as pseudo code.

###  4.2 Grid construction and travel time

A grid construction is necessary for both solving the linear programs and
modeling the simulation. The main idea is to discretize space by representing
the continuous space with a grid, like in figure 11. If an event happens
somewhere within a grid cell, we assume it happens at the center of the grid
cell. Discretizing the ambulance space is common in both linear programming
theory [7] as well as other simulation models [32].

To construct the grid, first we inscribe the case study region into a
rectangle. Then we split the grid into uniformly sized rectangle members. Note
there is a tradeoff between how much resolution (number of grid cells) to
include and how much computation is required. If the grid is finer it has more
accurate travel times but longer computations and vice versa. Once the grid is
established, any call within a rectangle member will be considered to be at
the center. Then the travel time between any two points becomes the distance
of the members of the grid cells they belong to.

We compute travel times with Open Street Map (OSM). In general, we find
computed travel times are slower than the reported travel time. Brown et
al.[10] found that ambulances using lights and sirens traveled on average 105
seconds faster than ambulances without. This is to say that on average the
special privileges of an ambulance allow it to go faster by bypassing traffic.
Thus we recommend adjusting grid times to be more accurate. Regression will
depend on the data available and other apriori knowledge, but can be
relatively simple including temporal variables like day-of-week and hour or
spatial variables like the number of grid cells away the path is.

###  4.3 Verification

The verification step builds confidence that the model reflects the real
world. It evaluates how well the calibration step compares to real world
measurements. The verification model inputs historical actions and compares
the simulated times with the reported times. Because there is a lot of
randomness in the drive time to do with paths and traffic, we compare the
overall trends for $n$ = 1000 calls or about two weeks. We then look at the
mean response time error on average. To measure confidence, we then invoke
Multiple Replication Procedure (MRP) on multiple batches of $n$-calls and
measure the variance among batches.

###  4.4 Alpha Value

As suggested by Bertsimas et al.[7], the $\alpha$ hyperparameter is determined
experimentally for the robust linear program. This is done using cross-
validation. In step one we choose a list of alpha values, for example
$[0.1,.01,.001]$. In step two we split the dataset into a $k$-fold of training
and testing sets with an $80/20$ split. We optimize the robust optimization on
the training set and simulate it with our model on the testing set. We then
remove any alpha values that are saturated (every station includes exactly one
ambulance) and choose the alpha value which gives the best performance.

##  5 Case Study

###  5.1 Motivation on Austin EMS

In this section, we provide a case study of Austin, Texas, as an application
of our open-source OpenEMS. Austin-Travis County EMS (ATCEMS)
111https://www.austintexas.gov/content/ems-austin-travis-county has 37 full-
time ambulances that on average serve 338 calls per day. Data in this study
includes 246,809 calls from the 2019-2020 period. ATCEMS serves 2.2 million
residents within a 1,039 square mile space. Since Austin is one of the biggest
cities in the United States, it is an interesting and representative subject
to study.

###  5.2 Data and statistics

Here we give an outline of what publicly available data exists. The main data
entries included were the date of event, longitude, and latitude of call, and
reported ambulance response time. (There also exists other entries such as
reported ambulance assign time, reported ambulance on-site time, reported
ambulance hospital time, ect.) We then graphed all calls and constructed the
temporal distribution in figure 2. From this figure, we noticed that the most
consistently high demand calls are from 8:00-20:00 and from Monday-Friday. We
call these times peak hours because they have more consistently high call
demand than the offpeak hours. Further, we remove Saturday and Sunday because
the times of high demand are less consistent. For the following results, we
limit our experiments to these peak hours as they are the times which Austin
EMS system is most stressed. Offpeak hours tend to be slower and easier to
manage. We also plotted the data geospatially in 3. Like most cities, a
majority of calls happen in the downtown region. Austin also includes a
sprawling suburb area that has relatively fewer calls but takes a long time to
drive between. This makes Austin hard to solve using heuristic planning to
handle both the high-demand area of downtown while still being able to handle
the occasional suburban calls.

Figure 2: A heat map showing the frequency of calls with respect to day of
week and call hour. Figure 3: A heat map showing the frequency of calls with
respect to the spatial distribution.

###  5.3 Calibration via travel time adjustment

The calibration step is meant to calibrate the grid time to the reported
travel time. To do this we tried several regression models to fit our grid
time with the reported travel time. In general, we find that there is a low
correlation because while grid time is deterministic, the reported travel time
is dependent on many unknown factors such as traffic. We tried including
several other variables like time of day and the spatial distance from two
points but we found there were marginal benefits.

For this case study, we found that trimming the top and bottom 1 percent of
data was necessary to avoid outlier reported travel times like 0 minutes
(which is impossible) or calls that took more than an hour. Then we fitted the
log of grid time to the log of reported time. The $r^{2}$ value was $.1472$.
One possible reason is that the variance is mostly not dependent on the
average grid time, but on unobserved factors such as traffic. We computed the
adjusted grid times and compared them to the reported grid times in figure 4.
From this we observe that the adjusted times are squeezed with shorter
tailends than before. From the validation we argue this produces a fair mean
response time, but loses some resolution of extreme events as well as skewing
times around 10 minutes which makes it harder to calculate shortfalls.

Figure 4: An illustration of how the regression transforms the grid time.

###  5.4 Verification

The results of the verification are shown in figure 5. Each data point
represents the mean response time of 1000 calls or approximately 2 weeks’
worth of data. Before the grid time had an error of $-186.23\pm 70.13$
seconds, while with the adjusted grid time the error becomes $2.41\pm 35.31$.
From this, we see the calibration reduces the error. Thus we could expect the
simulation mean response time to be close to the actual events within a
quantified confidence.

Figure 5: This is a plot of 20 batches where each batch is one data point.

###  5.5 Alpha Value

Next, we determine the best $\alpha$ value experimentally with a 3-fold cross-
validation for the robust optimization. Results are in the appendix 5. The set
of optimal stationing was then tested on the same hold-out set. Interestingly,
a majority of the optimizations have the same optimal value. This is the case
except for when $\alpha$ becomes too big (.1) where it may not predict as many
scenarios and when $\alpha$ becomes too small (.0001) where it may predict for
too extreme scenarios. From this cross-validation we choose $\alpha=.01$ for
the following experiments.

###  5.6 Computational Setup

The code was primarily written in Julia [8] with auxiliary scripts in Python.
“Julia For Mathematical Programming” (JuMP) [28] was used to interface Gurobi
9.0 with Julia. The computations are done on a standard computer. The station
has 16 GB of RAM and a 2.5 GHz processor.

###  5.7 Experimental Setup

To compare with real data we use 40 ambulances, similar to Austin EMS during
the day. Simulations are run for 1000 calls. Then 12 batches of 1000 call
simulations are run to establish confidence using Multiple Replication
Procedure (MRP). (Note if a simulation is run on too few samples it is biased
to perform better because all of the ambulances are available in the
beginning. Then it will reach an equilibrium. If shortfalls become too
numerous in a certain region, they can plug the system and give a bias towards
worse performance.) The time it takes to serve is based on a lognormal
distribution with a mean 3.65 and variance 0.3. Calls are filtered to peak
hours to be more homogenous (see subsection 5.2). In this simulation, the
closest ambulance is always sent to the call. All parameters are summarized in
table 2.

Variable Name | Value
---|---
number of ambulances | 40
number of calls | 1000
number of batches | 12
serve time distribution | log-normal($\mu=$ 3.65, $\sigma=$ 0.3)
$\alpha$ | .01
Table 2: Summary of experimental parameters.

###  5.8 Numerical results

Results from the simulation can be summarized in table 3. We found the
stochastic optimization performed better than robust optimization. Although,
the robust optimization sometimes outperformed the stochastic optimization
when there were more ambulances. This may be because robust optimization is
more conservative, waiting for an anticipated scenario that does not happen in
these batches. Both results perform better than the reported call time though.
Namely, the stochastic optimization is expected to have a 88.02 second faster
mean response time than heuristic decision making.

The ambulance stationing was visualized in figure 6. The purple circles are
proportional to how many ambulances should be stationed at each location. We
note that the robust optimization will concentrate more ambulances in one
station. We believe this is to handle some anticipated extreme scenarios that
the robust optimization generated. Implementation of these policies is easy,
as they only need to change the ambulance stationings to the optimal ones we
calculated.

Batch Number | Reported MRT | Stochastic Mean Response Time | Robust MRT
---|---|---|---
Batch 1 | 8.021 | 6.145 | 6.706
Batch 2 | 8.418 | 6.22 | 6.747
Batch 3 | 7.306 | 6.224 | 6.76
Batch 4 | 8.539 | 6.239 | 6.81
Batch 5 | 7.385 | 6.185 | 6.769
Batch 6 | 7.799 | 6.2 | 6.831
Batch 7 | 7.512 | 6.24 | 6.787
Batch 8 | 7.276 | 6.134 | 6.71
Batch 9 | 7.891 | 6.264 | 6.856
Batch 10 | 7.064 | 6.231 | 6.893
Batch 11 | 7.423 | 6.229 | 6.704
Batch 12 | 7.495 | 6.21 | 6.828
Overall | 7.677 +/- 0.443 | 6.210 +/- 0.037 | 6.783 +/- 0.059
Table 3: The results of MRP batches on the simulation.

Figure 6: Ambulance location schemes for the Stochastic model (left) and
Robust ($\alpha$ = .01) mode (right) for $n=40$. The purple circle indicates
there is at least one ambulance located in the area. The radius of the circle
is proportional to the number of ambulances in the area.

###  5.9 Adding or removing ambulances.

Simulations can also help a city understand what to do when they have to
permanently add or remove an ambulance from their fleet. An ambulance may be
added if a city’s population is growing and the department has new funds to
employ a new ambulance engendering the question, where should they put it? The
opposite scenario is when a department budget is cut and they must remove an
ambulance, which one can they possibly remove to have the least impact? Thus
optimization can help a city with a small budget maintain performance while
removing an ambulance from operation.

We conduct the optimal stationing as the number of ambulances changes in
figure 7. As can be seen, both the adding and the removal of ambulances only
cause minor changes. This is good in the sense that if an ambulance must be
removed performance does not decrease that much. But disappointingly adding
ambulances also does not increase performance. These differences we found may
be magnified though, as the grid time adjustment we saw earlier tends to
”squeeze” the times together.

Figure 7: The effects of adding or removing ambulances with Austin EMS.

##  6 Understanding socio-economic factors

To test whether socio-economic factors play a role in influencing the EMS
response time across different neighborhoods across Austin, we applied the
following protocol:

Dataset For each census tract, we calculated the mean reported travel time of
all incidents located within the census tract as the dependent variable. As
for the independent variables, we included the minimum grid time and the
average grid time within each census tract. This was calculated in two steps.
First we calculated the minimum and average grid times for each grid cell.
Then we took the weighted average (based on the call frequency) of grid cells
within each census tract to represent them. Furthermore, we included 19 socio-
economic independent variables extracted from the 2018 CDC Social
Vulnerability Index dataset [1]. The variables included in this study are
listed in Table 6 in the appendix.

5-fold Cross Validation We applied $5$-fold cross-validation. For each fold,
we included $80\%$ as the training set and the remaining $20\%$ as the test
set. Before training predictive models, we normalized each variable of the
training set by subtracting the mean and then dividing by the standard
deviation. We then applied the same to the test set using the sample mean and
standard deviation from the training set.

Models We applied 5 models in total. The baseline model was simply the average
of all the census tracts’ mean reported travel time. Because we believe the
reported travel time is correlated with the grid time, we then fitted 3 linear
regressions models: the first model used minimum grid time as the sole
independent variable; the second model used average grid time as the sole
independent variable; the third model included minimum grid time and average
grid time as its independent variables. The last model is a LASSO (Least
Absolute Shrinkage and Selection Operator) regression model incorporating all
19 socio-economic factors besides minimum grid time and average grid time.

Average MSE on the test set To compare these 5 models, we calculated the mean
squared error (MSE) on the test set for each fold. Then we took the average
among these 5 folds. If the average MSE of the LASSO model incorporating
socio-economic factors is lower than all the other models, then we have
evidence to believe that socio-economic factors contribute to differences in
reported travel time among various census tracts. Otherwise, we believe that
the reported travel time solely depends on geo-location factors reflected by
the grid time and is not explicitly affected by socio-economic factors.

As reported in table 4, we found that the lowest average MSE was achieved by
the linear regression model fitted on the average grid time. Thus, we conclude
that socio-economic factors may not explicitly affect the travel time for
ambulances in Austin.

Model | Variable(s) | Average MSE
---|---|---
Mean in the train set | N/A | 9.7980
Linear Regression | min.station.time | 9.6429
Linear Regression | avg.station.time | 9.5944
Linear Regression | min.station.time + avg.station.time | 9.6126
Lasso | All 21 variables | 9.7782
Table 4: Comparing the 5 regression models.

##  7 User Guide

OpenEMS is open-source and published at https://github.com/joshua-
ong/AmbulanceDeployment. A user guide of the software package can be seen in
figure 8; an overview of the software variables are in table 7.

Figure 8: Main modules of the openEMS software package. Python modules are in
squares, Julia modules are in diamonds.

Variable |  Description
---|---
$|$ regions $|$ |  The number or regions is how many uniform discrete regions the continuous space is made into. This can be though of as the resolution of the grid.
grid |  An object containing the information of the discrete grid. This includes a time matrix between any two grid cells.
adjacent neighborhood |  Boolean matrix in $\\{0,1\\}^{|regions|\times|regions|}$. Where index $i,j$ is 1 if region $i$
is adjacent to region $j$.
coverage |  Boolean Matrix in $\\{0,1\\}^{|stations|\times|regions|}$ where index $i,j$ is 1 if station $i$ is within 10 minutes of region $j$.
LP train calls |  Matrix where each row contains how many calls happened each hour. So entry $i,j$ is the number of calls at hour $i$ and region $j$.
callibration train calls |  Vector of train calls including time and location.
test calls |  Vector of test calls including time and location.
callibration |  Object that takes in grid time and outputs the adjusted grid time.
lognormal$\\_$params |  This includes ($\mu$,$\sigma$) for a log normal distribution that are fit to the service time
$\textbf{x}^{*}$ |  Integer vector in $\mathbb{Z}^{|Stations|}$ where index $i$ is the optimal number of ambulances in station $i$.
ncalls |  number of calls to simulate
nbatches |  number of batches to simulate
List of variables for user guide. |

Open Street Map module creates the discrete grid described in section 4.2.
Inputs include the number of regions, as well as the minimum/maximum
longitudes/latitudes that bound the city of interest. The module returns the
grid object which can be saved as a .json file. This grid object includes both
the data about the grid cells as well as useful functions to handle that data.
For socioeconomic interest, it can also contain a mapping of census tracts to
grid points.

$\begin{array}[]{ll@{}ll}>\texttt{grid =
create$\\_$grid(num\\_regions,min$\\_$lattitude,max$\\_$lattitude}\\\
\texttt{,min$\\_$longitude,max$\\_$longitude) }\end{array}$

Data Pre-processing module collects the city EMS call data with headers for
date of call, latitude and longitude of the emergency call for training the
linear program and simulations. Further headers for reported travel time,
reported service time, and latitude and longitude of the ambulance at the time
of the call are needed for calibration and verification. The data is then
processed into an adjacent neighborhood matrix, coverage matrix, LP train call
matrix, calibration train call matrix, and test call matrix and saved. The
user may consider customizing training and testing sets. For example, if the
city’s data changed drastically because of the Covid-19 pandemic, they could
train before Covid-19 and test during Covid-19.

$\begin{array}[]{ll@{}ll}>\texttt{(adjacent$\\_$nbhd,coverage,LP$\\_$train$\\_$call,callibration$\\_$train\\_call,test\\_call)}\\\
\texttt{=data\\_pre$\\_$processing(city$\\_$data =
"path/to/city$\\_$data.csv",grid) }\end{array}$

Calibration module includes fitting a regression between grid times and
reported travel times to calibrate the model. There is an example in section
5.3. Datasets should be trained with the callibration_train_set. We included
linear and log-log linear regressions, but we suggest the reader look at other
regressions to find what works best for their dataset. The result is the
adjusted travel time.

Verification module measures the accuracy of the adjusted travel time. There
is an example in 5.4. The test calls vector includes the location and time of
the call, as well as the reported response time and longitude and latitude of
the ambulance at the time of the call. This is so we can verify the ambulance
reported travel path time with our grid. The calibration object takes grid
time and outputs adjusted grid time. The outputs are the differences in
adjusted grid time and reported time. This can be measured as a score or as a
plot.

$\begin{array}[]{ll@{}ll}>\texttt{verification(test\\_call,callibration =
log$\\_$callibration)}\end{array}$

Linear Program Solvers module solves the two stage stochastic and robust
programs detailed in section 3. Besides the pre-processed data, the other
input is the $\alpha$ value that determines how conservative the robust model
is. The function then returns the solved optimal stationing $x_{stoch^{*}}$
and $x_{robust^{*}}$.

$\begin{array}[]{ll@{}ll}>x_{stoch^{*}},x_{robust^{*}}=\texttt{generate$\\_$models($\alpha$
= .1,LPtrain$\\_$calls,adjacent$\\_$nbhd,}\\\
\texttt{coverage,callibration=log$\\_$callibration,grid)}\\\ \end{array}$

Simulation module performs numerical computations using the test calls,
optimal stationing, and adjusted grid time. There is an example in 5.8. The
inputs were all computed in previous sections. The output is the computed
response time. From this one can compute their metrics.

$\begin{array}[]{ll@{}ll}>\texttt{simulated\\_response\\_time =
simulation(x,test$\\_$calls,grid,n$\\_$calls}\\\ \texttt{=
1000,lognormal\\_params = ($\mu$,$\sigma$),callibration =
log$\\_$callibration)}\end{array}$

Results The results will require unique analysis for each case study. All of
the visualizations from this paper are included and can be reproduced.

##  8 Conclusion

This paper presents an open-source project OpenEMS, an end-to-end pipeline to
optimize ambulance strategic decisions. It includes two-stage stochastic and
robust programming to optimize ambulance stationing and routing. In this open-
source package, we include many of the details that prevent state of the art
research from being applied on a case study level. This includes constructing
a realistic grid using open street maps, calibrating these travel times to
account for ambulances’ special privileges, and verifying they work how we
think they should work. Further, we provide practical guidance to hypothetical
questions like the case where an EMS department has to add or remove one or
more ambulances from their fleet. This is a common question.

We then applied this software to a case study on the city of Austin, Texas. In
it, we found we could increase the mean response time by 88.02 seconds if the
stationing was changed. We also showed the stochastic optimization performed
better than the robust optimization. We also found the ideal locations to add
or remove an ambulance to increase or maintain performance respectively.

There is a lot of avenues for future work. Now that this code is open-source,
one could simply try some new mathematical program and simulate for new
results, for example, the distributionally robust optimization (DRO) that is
growing in popularity in the operation research community. Along the lines of
the simulation, we want to include more features such as sensitivity analysis
with the choice of the grid resolution. There are also more questions about
robust optimization; one main problem is the non-convex nature of the VaR
constraint in the formulation. This leads to some unexpected results. For
example, alpha = .1 and alpha = .01 could work fine, but this does not give
any guarantees on interpolations between those values.

## Funding details

This research is supported by the Good Systems Research Initiative, part of
University of Texas at Austin Bridging Barriers. Ngoc Mai Tran is supported by
NSF Grant DMS-2113468 and the NSF IFML 2019844 award to the University of
Texas at Austin. Evdokia Nikolova is supported by NSF CCF 1733832.

We also want to thank the Senior Project team from University of Texas at
Austin who helped code sections of this open-source project: Ethan Santoni-
Colvin, Michael Hilborn, Sudeep Narala, Xander Tedjo, Will Worthington

## Disclosure statement

The authors report there are no competing interests to declare.

##  9 References

## References

  * [1] Centers for disease control and prevention/ agency for toxic substances and disease registry/ geospatial research, analysis, and services program.  CDC/ATSDR Social Vulnerability Index Insert 2018 Database Texas, Accessed on 08/08/2021.
  * [2] L. Aboueljinane, E. Sahin, and Z. Jemai.  A review on simulation models applied to emergency medical service operations.  Computers and Industrial Engineering, 66(4):734–750, 2013.
  * [3] Lina Aboueljinane, Zied Jemaï, and E. Sahin.  Reducing ambulance response time using simulation: The case of val-de-marne department emergency medical service.  Proceedings Title: Proceedings of the 2012 Winter Simulation Conference (WSC), pages 1–12, 2012.
  * [4] Shabbir Ahmed.  Two-Stage Stochastic Integer Programming: A Brief Introduction.  American Cancer Society, 2011.
  * [5] Rajan Batta, June M. Dolan, and Nirup N. Krishnamurthy.  The maximal expected covering location problem: Revisited.  Transportation Science, 23(4):277–287, 1989.
  * [6] P. Beraldi, M.E. Bruni, and D. Conforti.  Designing robust emergency medical service via stochastic programming.  European Journal of Operational Research, 158(1):183–193, 2004\.
  * [7] Dimitris Bertsimas and Yeesian Ng.  Robust and stochastic formulations for ambulance deployment and dispatch.  European Journal of Operational Research, 279(2):557–571, 2019\.
  * [8] Jeff Bezanson, Alan Edelman, Stefan Karpinski, and Viral B. Shah.  Julia: A fresh approach to numerical computing, 2015.
  * [9] John R. Birge and François Louveaux.  Introduction to Stochastic Programming.  Springer New York, 2011.
  * [10] Lawrence H. Brown, Christa L. Whitney, Richard C. Hunt, Michael Addario, and Troy Hogue.  Do warning lights and sirens reduce ambulance response times?  Prehospital Emergency Care, 4(1):70–74, 2000.
  * [11] M.L. Burkey, J. Bhadury, and H.A. Eiselt.  A location-based comparison of health care services in four u.s. states with efficiency and equity.  Socio-Economic Planning Sciences, 46(2):157–163, 2012.  Modeling Public Sector Facility Location Problems.
  * [12] Charles Courtemanche, Andrew Friedson, Andrew P. Koller, and Daniel I. Rees.  The affordable care act and ambulance response times.  J. Health Econ., 67:102213, Sep 2019.
  * [13] Karl F. Doerner and Richard F. Hartl.  Health care logistics, emergency preparedness, and disaster relief: New challenges for routing problems with a focus on the austrian situation.  In Operations Research/Computer Science Interfaces, pages 527–550. Springer US, 2008.
  * [14] David J. Eaton, Mark S. Daskin, Dennis Simmons, Bill Bulloch, and Glen Jansma.  Determining emergency medical service vehicle deployment in austin, texas.  Interfaces, 15(1):96–108, feb 1985.
  * [15] Guy Fitzgerald and Nancy L. Russo.  The turnaround of the London Ambulance Service Computer-Aided Despatch system (LASCAD).  European Journal of Information Systems, 14(3):244–257, Sep 2005\.
  * [16] Michel Gendreau, Gilbert Laporte, and Frédéric Semet.  Solving an ambulance location model by tabu search.  Location Science, 5(2):75–88, 1997.
  * [17] M. Haklay and Patrick Weber.  Openstreetmap: User-generated street maps.  IEEE Pervasive Computing, 7:12–18, 2008.
  * [18] J.E. Hammann and N.A. Markovitch.  Introduction to arena [simulation software].  In Winter Simulation Conference Proceedings, 1995., pages 519–523, 1995.
  * [19] C.R. Harrell and R.N. Price.  Simulation modeling and optimization using promodel.  In 2000 Winter Simulation Conference Proceedings (Cat. No.00CH37165), volume 1, pages 197–202 vol.1, 2000.
  * [20] A. J. Heightman.  Resource overkill: we can do more for less.  JEMS, 34(3):16., Mar 2009.
  * [21] Shane G. Henderson and Andrew J. Mason.  Ambulance Service Planning: Simulation and Data Visualisation.  In Operations Research and Health Care, pages 77–102. Springer, Boston, MA, Boston, MA, USA, 2005.
  * [22] A Ingolfsson, E Erkut, and S Budge.  Simulation of single start station for edmonton ems.  Journal of the Operational Research Society, 54(7):736–746, 2003\.
  * [23] C. J. Jagtenberg, S. Bhulai, and R. D. van der Mei.  Dynamic ambulance dispatching: is the closest-idle policy always optimal?  Health Care Manag. Sci., 20(4):517–531, Dec 2017.
  * [24] Takehiro Furuta Atsuo Suzuki Keisuke Inakawa.  Effect of ambulance station locations and number of ambulances to the quality of the emergency service.  In Ninth International Symposium on Operations Research and Its Applications (ISORA’10) Chengdu-Jiuzhaigou, Jun 2010.  [Online; accessed 19. Jun. 2021].
  * [25] Sean Shao Wei Lam, Yee Sian Ng, Mohanavalli Rajagopal Lakshmanan, Yih Yng Ng, and Ong. Marcus Eng Hock.  Ambulance deployment under demand uncertainty.  Journal of Advanced Management Science, 2016.
  * [26] Sean Shao Wei Lam, Yee Sian Ng, Mohanavalli Rajagopal Lakshmanan, Yih Yng Ng, and Marcus Eng Hock Ong.  Ambulance deployment under demand uncertainty.  Journal of Advanced Management Science Vol, 4(3), 2016.
  * [27] Kanglin Liu, Qiaofeng Li, and Zhi-Hai Zhang.  Distributionally robust optimization of an emergency medical service station location and sizing problem with joint chance constraints.  Transportation Research Part B: Methodological, 119:79–101, 2019\.
  * [28] Miles Lubin and Iain Dunning.  Computing in Operations Research Using Julia.  INFORMS Journal on Computing, 27(2):238–248, May 2015.
  * [29] Mohammed-Najeeb Mahama, Ernest Kenu, Delia Akosua Bandoh, and Ahmed Nuhu Zakariah.  Emergency response time and pre-hospital trauma survival rate of the national ambulance service, Greater Accra (January – December 2014).  BMC Emerg. Med., 18(1):1–7, Dec 2018.
  * [30] Andrew James Mason.  Simulation and real-time optimised relocation for improving ambulance operations.  In Handbook of Healthcare Operations Management, pages 289–317. Springer New York, 2013.
  * [31] Matthew S. Maxwell, Mateo Restrepo, Shane G. Henderson, and Huseyin Topaloglu.  Approximate dynamic programming for ambulance redeployment.  INFORMS Journal on Computing, 22(2):266–281, may 2010.
  * [32] R.J. McCormack and G. Coates.  A simulation model to enable the optimization of ambulance fleet allocation and base station location for increased patient survival.  European journal of operational research., 247(1):294–309, May 2015\.
  * [33] Laura A. McLay and Maria E. Mayorga.  A model for optimally dispatching ambulances to emergency calls with classification errors in patient priorities.  IIE Transactions, 45(1):1–24, jan 2013.
  * [34] Hozumi Morohosi.  A case study of optimal ambulance location problems.  Operations Research and its Applications (The Sixth International Symposium, ISORA’08 Proceedings), 2008.
  * [35] Suriyaphong Nilsang, Chumpol Yuangyai, Chen-Yang Cheng, and Udom Janjarassuk.  Locating an ambulance base by using social media: a case study in bangkok.  Annals of Operations Research, 283(1-2):497–516, jun 2018.
  * [36] Peter T. Pons and Vincent J. Markovchick.  Eight minutes or less: does the ambulance response time guideline impact trauma patient outcome?1 1Selected Topics: Prehospital Care is coordinated by Peter T. Pons, MD, of Denver Health Medical Center, Denver, Colorado.  J. Emerg. Med., 23(1):43–48, Jul 2002.
  * [37] Charles ReVelle and Kathleen Hogan.  The maximum availability location problem.  Transportation Science, 23(3):192–200, aug 1989.
  * [38] Verena Schmid.  Solving the dynamic ambulance relocation and dispatching problem using approximate dynamic programming.  European Journal of Operational Research, 219(3):611–621, 2012\.  Feature Clusters.
  * [39] Masoud Swalehe and Semra Gunay Aktas.  Dynamic ambulance deployment to reduce ambulance response times using geographic information systems: A case study of odunpazari district of eskisehir province, turkey.  Procedia Environmental Sciences, 36:199–206, 2016.  International Conference on Geographies of Health and Living in Cities: Making Cities Healthy for All.
  * [40] Joseph Tassone and Salimur Choudhury.  A comprehensive survey on the ambulance routing and location problems, 2020.
  * [41] Chapin White.  Health care spending growth: How different is the united states from the rest of the OECD?  Health Affairs, 26(1):154–161, jan 2007.
  * [42] Bo Zeng and Long Zhao.  Solving two-stage robust optimization problems using a column-and-constraint generation method.  Operations Research Letters, 41(5):457–461, 2013.

##  Appendix A Additional figures for numerical results

This section presents additional tables and figures for the numerical
experiments.

Figure 9: Ambulance location schemes for the stochastic optimization (left column) and robust ($\alpha$ = .01) optimization (right column). Rows start with 30 ambulances, increase by 5 ambulances per row, and end with 50 ambulances. The purple circle indicates there is an ambulance station there. The radius of the circle is proportional to the number of ambulances in the area where the biggest circle is 5 and the smallest circle is 0. Figure 10: A scatter plot of all call data. Note they cannot represent density. Figure 11: A visualization of the discrete grid over Austin county. |  $\alpha$ = .1 |  $\alpha$= .05 |  $\alpha$ = .01 |  $\alpha$ = .001 |  $\alpha$ = .0001
---|---|---|---|---|---
cross validation 1 | 6.807 | 6.562 | 6.562 | 6.562 | 6.562
cross validation 2 | 6.806 | 6.562 | 6.562 | 6.562 | 7.190
cross validation 3 | 6.562 | 6.562 | 6.562 | 6.562 | 7.590
Table 5: Three-fold crossvalidation on 5 alpha values to find when $\alpha$
saturates.

##  Appendix B List of variables extracted from CDC SVI dataset

Variable | Description
---|---
E_TOTPOP | Population estimate
E_HU | Housing units estimate
E_HH | Households estimate
E_POV | Persons below poverty estimate
E_UNEMP | Civilian (age 16+) unemployed estimate
E_NOHSDP | Persons (age 25+) with no high school diploma estimate
E_AGE65 | Persons aged 65 and older estimate
E_AGE17 | Persons aged 17 and younger estimate
E_DISABL | Civilian noninstitutionalized population with a disability estimate
E_SNGPNT | Single parent household with children under 18 estimate
E_MINRTY | Minority (all persons except white, non-Hispanic) estimate
E_LIMENG | Persons (age 5+) who speak English ”less than well” estimate
E_MUNIT | Housing in structures with 10 or more units estimate
E_MOBILE | Mobile homes estimate
E_CROWD | At household level (occupied housing units), more people than rooms estimate
E_NOVEH | Households with no vehicle available estimate
E_GROUPQ | Persons in group quarters estimate
E_UNINSUR | Uninsured in the total civilian noninstitutionalized population estimate
E_DAYPOP | Estimated daytime population
Table 6: List of variables extracted from CDC SVI dataset


# Carbon pseudospheres and the BTZ black hole

A. Iorio<EMAIL_ADDRESS>Institute of Particle and Nuclear
Physics, Faculty of Mathematics and Physics, Charles University, V
Holešovičkách 2, 18000 Praha 8, Czech Republic.

###### Abstract

I first recall the uses of Dirac materials as table top realizations of high
energy physics scenarios. Then I point to a specific system that might
reproduce a massless BTZ black hole, where the key role is played by
hyperbolic carbon pseudospheres. Finally some considerations are offered on
the possibility to realize rotating black holes, along with some comments on
the future of the whole analog gravity enterprise.

##  I Introduction

The field of analogs has as noble father Richard Feynman who ignited the field
in a famous lecture titled “Electrostatic Analogs”, available in Feynman  (see
also twostories ). There he explains why analog systems do describe the same
physics, based on the thrilling hypothesis of more elementary constituents
than the ones we deem to be fundamental. Amazingly, when space itself is
included as an emergent phenomenon, these are also the conclusions of certain
completely independent arguments of contemporary quantum gravity bekenstein ;
scholtz ; carroll .

As for the field of gravity analogs, see Volovik:2003fe , the seminal paper is
that of Unruh of 1981, where he proposes to search for experimental signatures
of his and Hawking’s effects, in a fluid dynamical analog UnruhAnalog . Due to
our deeper understanding and experimental control of condensed matter systems,
it is now becoming increasingly popular to reproduce that and other aspects of
fundamental physics in analog systems. Examples include the Hawking phenomenon
in Bose–Einstein condensates Steinhauer:2015saa , the Weyl symmetry iorio  and
the related Hawking/Unruh phenomenon on graphene ioriolambiase1 ,
gravitational and axial anomalies in Weyl semimetals Gooth:2017mbd , and more
Ulf_LeonhardtPRL2019 . Actually, gravity analogs are not limited to condensed-
matter systems, as can be seen, e.g., by interpreting hadronization in heavy-
ion collisions as a consequence of the Unruh effect Castorina:2007eb ;
Castorina:2008gf .

Despite those impressive advances, there are still two milestones to reach.
One is to understand the epistemic role of analogs in fundamental hugh energy
physics, as not all theorists would agree that analogs are much more than mere
divertissements. In fact, experimental results obtained in analogs are not
used as feedbacks for the target theories they are analogs of (see, e.g.,
Dardashti2016 ; twostories ). Anothe milestone would be a reliable definition
of an analog BH entropy, or at least, of a QFT-like entanglement entropy that,
in the presence of horizons, might serve the scope of setting-up some form of
the second principle of BH thermodynamics.

Any progress in this direction would be truly important for the hep-th
research. Having some results there, we could eventually be able to address
the so-called information paradox, i.e., the apparent loss of information
during BH evaporation, a question that, most probably, cannot be entirely
solved via theoretical reasonings. See, e.g., Penrose1996 ; Hawking2004 ;
Almheiri2013 ; Hooft2016 ; Hooft2016_I ; Maldazena  for different points of
view.

On the analog side, theoretical work has shown over the years that black hole
physics can find an indirect realization in BEC systems tris , and thrilling
experimental evidences have confirmed this fact Steinhauer:2015saa . The
latter findings are often referred to as the first experimental examples of
the Hawking effect. Here we focus on the proposal of graphene as an analog of
high-energy fundamental physics111Inspired by those findings, fundamental
constituents of both matter and space have been proposed in scholtz ;
twostories ; smaldone . ioriolambiase1 ; iorio ; pabloStran ; ioriopaiswitten
; grapheneQFTreview ; reach the unreachable , based on the fact that its low-
energy excitations CastroNeto2009  are massless Dirac pseudo-relativistic
fermions (the matter fields $\psi$), propagating in a carbon two-dimensional
honeycomb lattice. The emergent (long-wave limit) description of the latter is
a surface (piece of spacetime described by the “emergent” metric
$g_{\mu\nu}$). Such a behavior is shared by a wide range of materials, ranging
from silicene and germanene through d-wave superconductors to topological
insulators wehling . Each of those materials has its own peculiarities, which
allow for further extensions of results obtained with graphene, and hence
permit to explore a wider range of the high-energy target systems. Let us now
give some details.

##  II Analog gravity on graphene

Graphene is a one-atom-thick allotrope of carbon, i.e. the closest in nature
to a 2-dimensional object. It was first theoretically speculated about wallace
, and, decades later, experimentally found geimnovoselovFIRST . Its honeycomb
lattice is made of two intertwined triangular sub-lattices. As is by now well
known, this structure is behind the description of its electronic properties
in terms of massless, (2+1)-dimensional, Dirac quasi-particles. If one
linearizes the tight-binding Hamiltonian around two Fermi points,
$\vec{k}^{D}_{\pm}=\left(\pm\frac{4\pi}{3\sqrt{3}\ell},0\right)$, then the
Hamiltonian becomes $H|_{\vec{k}_{\pm}}\simeq
v_{F}\sum_{\vec{p}}\left(\psi_{+}^{\dagger}\vec{\sigma}\cdot\vec{p}\;\psi_{+}-\psi_{-}^{\dagger}\vec{\sigma}^{*}\cdot\vec{p}\;\psi_{-}\right)$,
where $v_{F}=3\eta\ell/2\sim c/300$ is the Fermi velocity, $\psi_{\pm}$ are
two–component Dirac spinors, and $\vec{\sigma}\equiv(\sigma_{1},\sigma_{2})$,
$\vec{\sigma}^{*}\equiv(\sigma_{1},-\sigma_{2})$, with $\sigma_{i}$ the Pauli
matrices.

If one considers the linear regime only, the first scale is $E_{\ell}\sim
v_{F}/\ell\sim 4.2$eV. Notice that $E_{\ell}\sim 1.5\eta$, and that the
associated wavelength, $\lambda=2\pi/|\vec{p}|\simeq 2\pi v_{F}/E$, is
$2\pi\ell$. The electrons’ wavelength, at energies below $E_{\ell}$, is large
compared to the lattice length, $\lambda>2\pi\ell$. Those electrons see the
graphene sheet as a continuum. One Dirac point is enough, when only strain is
present (see, e.g., pabloStran ), and when certain approximations on the
curvature are valid ioriolambiase1 . The importance and relevance of the two
Dirac points for emergent hep-th descriptions has been discussed at length in
our work ioriopaiswitten , see also our recent tloop , where the focus though
is on torsion.

When only one Dirac point is necessary, the following Hamiltonian well
captures the physics of undeformed (planar and unstrained) graphene:
$H=-iv_{F}\int d^{2}x\;\psi^{\dagger}\vec{\sigma}\cdot\vec{\partial}\;\psi$,
where the two component spinor is, e.g., $\psi\equiv\psi_{+}$, we moved back
to configuration space, $\vec{p}\to-i\vec{\partial}$, and sums turned into
integrals because of the continuum limit. In various papers, we have exploited
this regime to a great extent, till the inclusion of curvature and torsion in
the geometric background. On the other hand, we also have investigated the
regimes beyond the linear one, where granular effects associated to the
lattice structure emerge, see GUP  and also the related GUPBTZ . When both
Dirac points are necessary, one needs to consider four component spinors
$\Psi\equiv\left(\begin{array}[]{c}\psi_{+}\\\ \psi_{-}\\\
\end{array}\right)$, and $4\times 4$ Dirac matrices
$\alpha^{i}=\left(\begin{array}[]{cc}\sigma^{i}&0\\\ 0&-{\sigma^{*}}^{i}\\\
\end{array}\right)$, $\beta=\left(\begin{array}[]{cc}\sigma^{3}&0\\\
0&\sigma^{3}\\\ \end{array}\right)$, $i=1,2$. These matrices satisfy all the
standard properties, see, e.g., grapheneQFTreview  and ioriopaiswitten . With
these, the Hamiltonian is $H=-iv_{F}\int
d^{2}x\left(\psi_{+}^{\dagger}\vec{\sigma}\cdot\vec{\partial}\;\psi_{+}-\psi_{-}^{\dagger}\vec{\sigma}^{*}\cdot\vec{\partial}\;\psi_{-}\right)=-iv_{F}\int
d^{2}x\;\bar{\Psi}\vec{\gamma}\cdot\vec{\partial}\;\Psi$.

In iorio  the goal was to identify the conditions for which graphene might
realize aspects of QFT in curved spacetime. Therefore, key issues had to be
faced, such as the proper inclusion of the time variable in a relativistic-
like description, and the role of the nontrivial vacua and their relation to
different quantization schemes for different observers. All of this finds its
synthesis in the Unruh or the Hawking effects ioriolambiase1 . Let us explain
here the main issues and the approximations made there.

Besides $E_{\ell}$, when we introduce curvature, we also have a second scale.
When this happens, $E_{\ell}$ is our “high energy regime”. This is so because
we ask the curvature to be small compared to a limiting maximal curvature,
$1/\ell^{2}$, otherwise: i) it would make no sense to consider a smooth
metric, and ii) $r<\ell$ (where $1/r^{2}$ measures the intrinsic curvature),
means that we should bend the very strong $\sigma$-bonds, an instance that
does not occur. Therefore, our second scale is $E_{r}\sim v_{F}/r$, with
$E_{r}=\ell/r\;E_{\ell}<E_{\ell}$. To have a quantitative handle on these
scales, let us take, e.g., $r\simeq 10\ell$ as a small radius of curvature
(high intrinsic curvature). To this corresponds an energy $E_{r}\sim 0.4$eV,
whereas, to $r\sim 1{\rm mm}\sim 10^{6}\ell$, corresponds $E_{r}\sim
0.6\mu$eV. The “high energy” to compare with is $E_{\ell}\sim 4$eV. When
energies are within $E_{r}$ (wavelengths comparable to $2\pi r$) the electrons
experience the global effects of curvature. That is to say that, at those
wavelengths, they can distinguish between a flat and a curved surface, and
between, e.g., a sphere and a pseudosphere. Therefore, whichever curvature
$r>\ell$ we consider, the effects of curvature are felt until the wavelength
becomes comparable to $2\pi\ell$. The formalism we have used, though, takes
into account all deformations of the geometric kind, with the exception of
torsion. Hence, this includes intrinsic curvature, and elastic strain of the
membrane (on the latter see pabloStran ), but our predicting power stops
before $E_{\ell}$, because there local effects (such as the actual structure
of the defects) play a role that must be taken into account into a QG type of
theory. On the latter the first steps were moved in GUP  (see also the related
GUPBTZ ).

The intrinsic curvature is taken here as produced by disclination defects,
that are customarily described in elasticity theory (see, e.g., Kleinert ), by
the (smooth) derivative of the (non-continuous) SO(2)-valued rotational angle
$\partial_{i}{\omega}\equiv{\omega_{i}}$, where $i=1,2$ is a curved spatial
index. The corresponding (spatial) Riemann curvature tensor is easily
obtained,
${R^{ij}}_{kl}=\epsilon^{ij}\epsilon_{kl}\epsilon^{mn}\partial_{m}\omega_{n}=\epsilon^{ij}\epsilon_{lk}2{\cal
K}$, where $\cal K$ is the Gaussian (intrinsic) curvature of the surface. In
our approach we have included time, although the metric we adopted is

$g^{\rm graphene}_{\mu\nu}=\left(\begin{array}[]{cc}1&0\quad 0\\\
\begin{array}[]{c}0\\\ 0\end{array}&g_{ij}\\\ \end{array}\right)\;,$ (1)

i.e., the curvature is all in the spatial part, and $\partial_{t}g_{ij}=0$.
Since the time dimension is included, the SO(2)-valued (abelian) disclination
field has to be lifted-up to a SO(1,2)-valued (non-abelian) disclination
field, ${\omega_{\mu}}^{a}$, $a=0,1,2$, with
$\omega_{\mu}^{\;a}=e^{b}_{\mu}\omega_{b}^{\;a}$ and the expression
$\omega_{a}^{\;d}=\frac{1}{2}\epsilon^{bcd}\left(e_{\mu
a}\partial_{b}E_{c}^{\mu}+e_{\mu b}\partial_{a}E_{c}^{\mu}+e_{\mu
c}\partial_{b}E_{a}^{\mu}\right)$, gives the relation between the disclination
field and the metric (dreibein). All the information about intrinsic curvature
does not change. For instance, the Riemann curvature tensor,
${R^{\lambda}}_{\mu\nu\rho}$, has only one independent component, proportional
to $\cal K$ (see iorio ). When only curvature is important, the long
wavelength/small energy electronic properties of graphene, are well described
by the action ${\cal A}=iv_{F}\int
d^{3}x\sqrt{g}\;\bar{\Psi}\gamma^{\mu}(\partial_{\mu}+\Omega_{\mu})\Psi$, with
$\Omega_{\mu}\equiv{\omega_{\mu}}^{a}J_{a}$, and $J_{a}$ are the generators of
SO(1,2), the local Lorentz transformations in this lower-dimensional setting.
In ioriopaiswitten  we have discussed at length this action within the Witten
approach witten3dgravity  to Poincaré ($ISO(2,1)$) or (A)dS gravity as gauge
theory, and within the USUSY approach, see also susyZanelli1 , and especially
the recent u-susy-graphene .

Figure 1:  The hyperbolic pseudosphere for $a=1$, $C=1$. Here $\rho_{\rm
min}=1$ and $\rho_{\rm max}\simeq 1.4142$.

Within this approach, a nontrivial $g_{tt}$ in (1), hence a clean nontrivial
general relativistic effect (recall that $g_{tt}\sim U_{grav}$) can only
happen if specific symmetries and set-ups map the lab system into the wanted
one, one being the local Weyl symmetry. This produced a measurable predictions
of a Hawking/Unruh effect, for certain specific shapes. It was ofund that
surfaces of constant Gaussian curvature, and among them, those of negative
curvature (that necessarily have singular boundaries, see grapheneQFTreview
and icrystals ) are key. The above lead to the proposal of a variety of set-
ups, especially three key spacetimes with horizon: the Rindler, the de Sitter
and the BTZ black hole BTZ1992 .

##  III Carbon pseudospheres and the BTZ black hole

Let us write the BTZ black hole metric as in GUPBTZ

$ds_{BTZ}^{2}=f(r)^{2}c^{2}dt^{2}-f(r)^{-2}dr^{2}-r^{2}(d\phi+N^{\phi}cdt)^{2}$
(2)

where

$f^{2}(r)=-\frac{8GM}{c^{2}}-\Lambda
r^{2}+\frac{16G^{2}J^{2}}{c^{4}\,r^{2}}\,,\quad\quad\quad
N^{\phi}=-\frac{4GJ}{c^{2}\,r^{2}}\,,$ (3)

with $M$ the mass, $\Lambda\equiv-1/\ell^{2}<0$ the negative cosmological
constant and $J$ the angular momentum. Horizons are located at the positive
zeros of $f(r)$

$r^{2}_{\pm}=\frac{4GM\ell^{2}}{c^{2}}\left[1\pm\left(1-\frac{J^{2}}{\ell^{2}M^{2}}\right)^{1/2}\right]\,.$
(4)

When $M>0$ and $|J|\leq M\ell$, we have a black hole, with $r_{+}$ a genuine
event horizon, and $r_{-}$ a Cauchy horizon.

###  III.1 The massless black hole

In GUPBTZ , see also ioriolambiase1 , the extremal case of $M\to 0$

$ds_{0}^{2}=(r/\ell)^{2}c^{2}dt^{2}-(r/\ell)^{-2}dr^{2}-r^{2}d\phi^{2}\;,$ (5)

was put into (conformal) correspondence to a graphene analog realization of
this scenario, when shaped in a very specific manner. There it is shown that
the most natural choice is to identify $\ell$ with $\ell_{L}$, the lattice
spacing. The shape that is necessary to realize is that of the hyperbolic
pseudosphere $\Sigma_{\rm HYP}$, see the figures, whose line element is

$dl_{\rm HYP}^{2}=du^{2}+C^{2}\cosh^{2}(u/a)d\phi^{2}\,,$ (6)

with $C=\ell_{L}$, $u$ the longitudinal coordinate, $\phi\in[0,2\pi]$ and the
constant negative Gaussian curvature given by $K=-1/a^{2}<0$. In terms of the
radial coordinate

$\rho(u)=C\cosh(u/a)\,,$ (7)

there is a singular boundary at the largest circle of radius $\rho_{\rm
max}=\sqrt{a^{2}+C^{2}}\equiv\rho_{Hh}$, where $C$ is the radius of the
smallest throat, $\rho_{\rm min}=C$, cf. Figs.1, 2 and 3.

Figure 2:  The hyperbolic pseudosphere for $a=1$, $C=1/10$. Here $\rho_{\rm
min}=C=1/10$ and $\rho_{\rm max}\simeq 1.005$.

Indeed, setting $J=0$, and easing a little the notation, $8G/c^{2}=1$

$ds^{2}_{BTZ}=\left(r^{2}/C^{2}-M\right)ds^{2}_{\rm HYP}\;,$ (8)

where $ds^{2}_{\rm HYP}$ is the line element of $\Sigma_{\rm
HYP}\times\mathbf{R}$. Therefore, all the relevant quantities of the BTZ black
hole are given in terms of measurable quantities

$\Lambda\equiv-1/\ell_{L}^{2}\quad,\quad M\equiv\ell_{L}^{2}/a^{2}\quad,\quad
r_{+}\equiv\ell_{L}^{2}/a\;.$ (9)

Before moving to the discussion of the $J\neq 0$ case, let us recall how the
event horizon, $r_{+}$ relates to the singular boundary of the hyperbolic
pseudosphere ioriolambiase1 ; GUPBTZ

$r_{Hh}\equiv r(u_{Hh})=r_{+}\coth\left({\rm
arccosh}\left(\sqrt{1+a^{2}/\ell_{L}^{2}}\right)\right)\;.$ (10)

In the limit of small $\ell_{L}/a$, these two horizons coincide. That is also
the limit where $M\to 0$, and, accordingly $r_{+}\to 0$, i.e. the announced
zero mass black hole. In the figures we show three cases of such pseudosphere,
for increasing closeness of $r_{Hh}$ to $r_{+}$.

There are three reasons to still call this a black hole.

First of all, it is continuously connected to the spectrum of BTZ black holes,
$M\geq 0$. In fact, unlike the $M<0$ cases, which are states under the black
hole threshold, and represent point-particles in $AdS_{3}$, the $M=0$ case is
not a point-particle. One way of looking at it is that in the flat limit
$\Lambda\to 0$ (i.e. $\ell\to\infty$), where point particles survive, the case
$M=0$ does not exist as a solution. Thus $M=0$ belongs more to the black holes
segment/set ($M>0$) than to the point-particles ($M<0$). On the point
particles see DJtH .

Second, the metric (5) has a horizon, even though it is not an event horizon.
Indeed, the massless BTZ is locally equivalent to $AdS_{3}$ in Poincaré
coordinates. If one decompactifies the angular coordinate, one finds that the
BTZ metric with $M=0$ is the $AdS_{3}$ Poincaré patch, since the “Poincaré
horizon” coincides with $r=0$. The Poincaré horizon is the $r=0$ of the $M=0$
BTZ in BTZ coordinates. It is written in coordinates $z=1/r$ sometimes. In
that case, obviously, the Poincaré horizon is $z=\infty$ while the boundary of
the space is $z=0$. It is well known, so it is not a result per se, but have a
look, e.g., at braga .

A third argument (related to the first) in favor of calling it a black hole is
that, from the point of view of holography, this solution ($M=0$) corresponds
to a well-defined state of the dual theory (although it is a state of $T=0$).
About the states of the dual theory, one should consider the link between
three-dimensional gravity and two-dimensional Liouville field theory, see,
e.g. the nice review modave , and also wittenM=0 , where the “black hole
threshold” is precisely the gap between $AdS_{3}$ and the $M=0$ BTZ.

###  III.2 The rotating black hole?

It is an intriguing and challenging question whether we could reproduce with
graphene a rotating black hole, $J\neq 0$, that is whether we could find a
specific configuration that can be associated to the spacetime described by
(2). One crucial observation is that such configuration needs to include a
nontrivial behavior along the time direction, on top of a nontrivial behavior
along the space direction. Of course, we are referring here to the important
role played in that metric by

$g_{t\phi}=2r^{2}N^{\phi}=-J\,.$ (11)

We have, in fact, at least two choices. One choice would be to insist with the
strategy above illustrated, and look for a configuration that is conformal to
the full BTZ metric, rather than conformal only to the $J=0$ case, see (8). In
this case, one should probably give up the possibility of having both, the
inner and outer horizons, $r_{\pm}$, in favor of the most important inner
horizon only. Given that now

$r_{\pm}=\frac{1}{2}\frac{\ell_{L}^{4}}{a^{2}}\left(1\pm(1-a^{2}/\ell_{L}^{4}J^{2})^{1/2}\right)\,,$
(12)

we loose spherical symmetry, hence we should expect a deformed hyperbolic
pseudosphere. In that case, the parameter $J$ will have to be related to such
a geometric deformation, that might as well be a static one. Therefore, one
should not expect, in principle, $J$ to be a true angular momentum, like the
one stemming from a spinning pseudosphere, but $J$ can rather be just a
geometric measure of the deviation of actual boundary from the singular circle
of the previous discussion, of ray $r_{Hh}$. At this point, though the road
becomes quite steep. In fact, such a deformed surface might serve well the
scope of visually reproducing the deformed horizons, at least the inner one,
but then would probably fail to be a surface of constant Gaussian curvature.
In that case, the strategy adopted in the previous discussion cannot be
applied, because we would not be guaranteed to be in a conformally flat space
time iorio .

Therefore, why not trying a completely different road, that is to actually act
on the time components of (1), hence to construct the wanted metric directly?
For this to work, we could shape the membrane along the radial direction in
such a way to reproduce $1/f^{2}(r)$. This alone would clearly produce a flat
metric, as the curvature would only be in one direction. At this point it is
crucial to remember that the metric we are discussing, is that experienced by
the conductivity electron of the material, henceforth, at least in principle,
this could be achieved by letting such electrons interact with a suitably
fine-tuned external electromagnetic field PabloLaser , in such a way that
$U_{em}$ mimics the $U_{grav}$ that is in $g_{tt}$. Therefore, we need to act
with a space-dependent external electromagnetic field in such a way that the
$g_{tt}$ component is equal to $f^{2}(r)$, and we are nearly done.

Figure 3:  The hyperbolic pseudosphere for $a=1$, $C=1/100$. Here $\rho_{\rm
min}=C=1/100$ and $\rho_{\rm max}\simeq 1.00005$.

In fact, the question left to answer is: what about the reproduction of the
crucial $g_{t\phi}$? This would probably be the most delicate step, because we
need here a feedback between spatial action and temporal action. One way to go
could be to use the electromagnetic field to induce elastic strain of the
membrane. If that happens, such a strain could produce a gauge field
$A_{\mu}^{strain}$, as customary, see, e.g., pabloStran , that can be summed
up to $A_{\mu}^{em}$, to obtain an overall term that mixes space and time
component, that is the spatial feedback on the metric obtained acting on the
time components, and viceversa. Notice that, while the interaction withe the
external field need be very fine tuned, in order to reproduce exactly the
wanted $f^{2}(r)$, the feedback on the strain can be quite generic as that
component is actually defining the $J$, as seen in (11), hence can be
arbitrary.

All the above is very fascinating, but need a thoroughly investigation of all
details, theoretical and experimental, that are beyond the scope of this
essay, and the focus on ongoing research PabloLaser .

##  IV Conclusions

The exciting and rapidly evolving field of analog gravity is facing a new era.
The interest is shifting from the reproduction of the kinematical aspects of
the Hawking/Unruh phenomenon, that has reached a climax of precision and
accuracy, to the realization of some form of dynamics, a very challenging
problem. The most important phenomenon to study is black-hole evaporation, the
most prominent dynamical phenomenon of quantum gravity, with its plethora of
open fundamental issues, such as the possibility that information is not
preserved in this process, etc. Having recalled here how Dirac materials lend
themselves to realize crucial aspects of black hole physics, we believe that
the search for the realizations of such dynamical will benefit from Dirac
materials, such as graphene, that showed already to be a powerful and
versatile analog of both quantum fields on curved spaces and quantum gravity.

## Acknowledgements

The author is indebted to Gaston Giribet and Jorge Zanelli for discussions on
the special status of the massless BTZ black hole. He gladly acknowledges
support from Charles University Research Center (UNCE/SCI/013) and from the
grant SVV No. 260576.

## References

  * (1) R. Feynman, et al., 2006 The Feynman Lectures on Physics (Pearson/Addison-Wesley).
  * (2) A. Iorio, J. Phys.: Conf. Series 1275 (2019) 012013 [arXiv:1902.07096].
  * (3) J. D. Bekenstein, Phys. Rev. D 23 (1981) 287; Phys. Rev. E 89 (2014) 1; Sci. Am. 289 (2003) 58.
  * (4) G. Acquaviva, A. Iorio, M. Scholtz, Ann. Phys. 387 (2017) 317.
  * (5) N . Bao, S. M. Carroll, A. Singh, Internat. J. Modern Phys. D 26 (12) (2017) 1743013.
  * (6) G. E. Volovik, The Universe in a helium droplet, Int. Ser. Monogr. Phys. 117 (2006) 1; C. Barceló, S. Liberati, M. Visser, Liv. Rev. Rel. 14 (2011) 3.
  * (7) W.G. Unruh, Phys. Rev. Lett. 46 (1981) 1351.
  * (8) J. Steinhauer et al., Nature 569 (2019) 688; J. Steinhauer, Nature Phys. 12 (2016) 959.
  * (9) A. Iorio, Ann. Phys. 326 (2011) 1334.
  * (10) A. Iorio, G. Lambiase, Phys. Let. B 716 (2012) 334; Phys. Rev. D 90 (2014) 025006.
  * (11) J. Gooth, et al, Nature 547 (2017) 324.
  * (12) U. Leonhardt et al, Phys. Rev. Lett. 122 (2019) 010404.
  * (13) P. Castorina, D. Kharzeev, H. Satz, Eur. Phys. J. C 52 (2007) 187.
  * (14) P. Castorina, A. Iorio, H. Satz, Int. J. Mod. Phys. E 24 (2015) 1550056; P. Castorina, D. Grumiller, A. Iorio, Phys. Rev. D 77 (2008) 124034.
  * (15) R. Dardashti, K. P. Thebault, E. Winsberg, Brit. J. Phil. Science 68 (2015) 55.
  * (16) R. Penrose, Gen. Rel. Grav. 28 (1996) 581.
  * (17) S. Hawking, in General relativity and gravitation, Proc. GR17, Dublin (2004) pp. 56.
  * (18) J. Polchinski, et al., J. High Energy Phys. 2013 (2013) 62.
  * (19) R. B. Mann, Black Holes: Thermodynamics, Information, and Firewalls (Springer, Berlin, 2015).
  * (20) G. ’t Hooft, (2016), arXiv:1612.08640v1.
  * (21) J. Maldacena, L. Susskind, Fortsch. Phys. 61 (2013) 781.
  * (22) R. Balbinot, A. Fabbri, S. Fagnocchi, A. Recati, I. Carusotto, Phys. Rev. A 78 (2008) 021603(R).
  * (23) G. Acquaviva, A. Iorio, L. Smaldone, Phys. Rev. D 102 (2020) 106002.
  * (24) A. Iorio, P. Pais, Phys. Rev. D 92 (2015) 125005.
  * (25) A. Iorio, P. Pais, Ann. Phys. 398 (2018) 265.
  * (26) A. Iorio, Int. J. Mod. Phys. D 24 5 (2015) 1530013.
  * (27) A. Iorio, Frontiers in Materials 1 (2015) 36.
  * (28) A. H. Castro Neto, F. Guinea, N. M. R. Peres, K. S. Novoselov, A. K. Geim, Rev. Mod. Phys. 81 (2009) 109.
  * (29) T. O Wehling, A. M. Black-Schaffer, A. V. Balatsky, Adv. Phys. 76 (2014) 1.
  * (30) P. R. Wallace, Phys. Rev. 71 (1947) 622; G. W. Semenoff, Phys. Rev. Lett. 53 (1984) 2449.
  * (31) K. S. Novoselov, A. K. Geim, et al, Science 306 (2004) 666.
  * (32) M. Ciappina, A. Iorio, P. Pais, A. Zampeli, Phys. Rev. D 101 (2020) 036021.
  * (33) A. Iorio, P. Pais, I. A. Elmashad, A. F. Ali, Mir Faizal, L. I. Abou-Salem, Int. J. Mod. Phys. D 27 (2018) 1850080.
  * (34) A. Iorio, G. Lambiase, P. Pais, F. Scardigli, Phys. Rev. D 101 (2020) 105002.
  * (35) H. Kleinert, Gauge fields in condensed matter, Vol II, World Scientific (Singapore) 1989; M. O. Katanaev, I. V. Volovich, Ann. Phys. 216 (1992) 1.
  * (36) E. Witten, Nucl. Phys. B 311 (1988) 46.
  * (37) J. Zanelli, et al., J. High Energy Phys. 1204 (2012) 058; J. High Energy Phys. 85 (2016) 201.
  * (38) L. Andrianopoli, et al., J. High Energy Phys. 2001 (2020) 084.
  * (39) A. Iorio, J. Phys.: Conf. Ser. 442 (2013) 012056; et al., J. Phys.: Cond. Matt. 28 (2016) 13LT01.
  * (40) M. Bañados, C. Teitelboim, J. Zanelli, Phys. Rev. Lett. 69 (1992) 1849.
  * (41) A. Iorio, G. Lambiase, G. Vitiello, Ann. Phys. 309 (2004) 151.
  * (42) M. Cvetic and G. Gibbons, Ann. Phys. 327 (2012) 2617.
  * (43) S. Deser, R. Jackiw, G. ’t Hooft, Three-dimensional Einstein gravity: Dynamics of flat space, Annals of Physics, 152 (1984) 220.
  * (44) C. A. Bayona and N. R. F. Braga, Gen. Rel. Grav. 39 (2007) 1367.
  * (45) L. Donnay, PoS Modave2015 (2016) 001 and references therein.
  * (46) J.-M. Schlenker, E. Witten, No Ensemble Averaging Below the Black Hole Threshold, arXiv:2202.01372 [hep-th].
  * (47) A. Iorio, P. Pais, Laser-graphene interaction and the improvement of gravity analogs, in preparation.