# Technologies and Computing Paradigms: Beyond Moore’s law?

Daniel Etiemble  Computer Science Laboratory (LISN)
University of Paris Saclay
Gif-sur-Yvette, France
<EMAIL_ADDRESS>
###### Abstract

As it is pretty sure that Moore’s law will end some day, questioning about the
post-Moore era is more than interesting. Similarly, looking for new computing
paradigms that could provide solutions is important. Revisiting the history of
digital electronics since the 60’s provide significant insights on the
conditions for the success of a new emerging technology to replace the
currently dominant one. Specifically, the past shows when constraints and «
walls » have contribute to evolution through improved technical techniques and
when they have provoked changes of technologies (evolution versus
breakthrough). The main criteria for a new technology or a new computing
paradigm is a significant performance improvement (at least one order of
magnitude). Cost, space requirement, power and scalability are the other
important parameters.

##  I Introduction

For many dozen of years, the evolution of digital electronics has been driven
by Moore’s law that states that the number of transistors in a circuit doubles
every N months, even if the process is slowing down, N increasing from 12 to
18 to 24. As a result of Moore’s law, new generations of CMOS/FinFET
technologies, called a technological node, have been regularly launched. The
different successive nodes are shown in Fig. 1.

However, everybody agrees on the fact that Moore’s law will end some time,
even if it is quite difficult to predict when. 1-nm node is often announced as
the last one. So, the big question is : What will be the future of computing
beyond Moore’s Law [1]. This question can be decomposed into several ones:

  * •

Will there be a new technology to continue the exponential performance
increase ?

  * •

Are there promising technologies that could replace CMOS/FinFET ?

  * •

Will there be new computing paradigms to overcome the lack of a replacement of
last node CMOS technologies ?

  * •

Does a Post-CMOS era mean a Post-binary era? [2]

Figure 1: CMOS technological nodes

Revisiting the history of digital electronics could provide significant
insights on the conditions for the success of a new emerging technology to
replace the currently dominant one. Specifically, the past shows when
constraints and “walls" have contribute to evolution through improved
technical techniques and when they have provoked changes of technologies.

##  II From transistors to Integrated Circuits with bipolar transistors.

While the concept of field-effect transistors was proposed in 1926, the first
working transistor was invented by Shockley in 1948, one year after Bardeen
and Brattain invented the point-contact transistor. Both were awarded the 1956
Nobel Prize in Physics for their discoveries. The bipolar transistor has been
used in different logic families, either as discrete transistors or integrated
circuit gates : RTL then DTL then TTL, ECL, I2L… ECL was the fastest family:
it was used in the famous CRAY-1 vector supercomputer.

###  II-A CDC6600 with discrete transistors

CDC 6600 was a mainframe used from December 1965 to May 1977 that is
considered as the first supercomputer. The clock frequency was 10-MHz. This
can be compared with the Intel 4004 (first microprocessor) 750 KHz clock
frequency launched in 1971. The logic circuitry of the CDC6600 is based on
resistor-transistor-logic (RTL) in which elementary transistors are used. It
should be mentioned that it included some architectural features that became
popular later: superscalar execution, multithreading between the central CPU
and the peripheral IO processors. Some details can be found in [3]

###  II-B Cray-1 with integrated circuits

The Cray-1 was launched in 1975 as the first implementation of the vector
processor design. The clock frequency was 80 MHz. The Cray-1 used only four
different IC types, an ECL dual 5-4 NOR gate (one 5-input, and one 4-input,
each with differential output) from Fairchild, another slower MECL 10K 5-4 NOR
gate used for address fanout, a 16×4-bit high speed (6 ns) static RAM (SRAM)
used for registers and a 1,024×1-bit 48 ns SRAM used for the main memory [4].

###  II-C End of bipolar technology as the dominant technology

The binary circuitries of I2L, TTL and ECL logic families are presented in
Fig. 2. While ECL uses an actual current generator, I2L and TTL use a pseudo
current generator. In both cases, the corresponding current flows towards
input or output for I2L and TTL, and towards left or right sides of the
differential pair for ECL. A significant current contributes to the static
power dissipation: it is not needed to consider dynamic power dissipation to
understand why bipolar digital circuits have disappeared as soon as MOS
technology has become mature and CMOS circuitry, that had no static power
dissipation at that period, was able to replace pMOS and nMOS circuitries.

Power dissipation was the key issue for switching from bipolar to MOS
technologies. Bipolar technologies were also far more complicated to integrate
than the MOS ones.

Figure 2: Bipolar circuitries

Bipolar technology coexisted for a while with MOS technology in a technology
called BiCMOS. It was supposed to combine the strengths of the two different
process technologies into a single chip: Bipolar transistors offer high speed
and gain, which are critical for high-frequency analog sections, whereas CMOS
technology excels for constructing simple, low-power logic gates. However,
combining two different processes was costly. As soon as MOS technologies were
able to provide buffers with sufficient fan-out, the BiCMOS technology was no
longer needed for digital circuits.

##  III MOS and CMOS technologies

###  III-A The first period of MOS technologies

The MOSFET transistor was invented by M. M. Atalla and D. Kahng at Bell Labs
in 1959, and first presented in 1960. By 1964, MOS chips had reached higher
transistor density and lower manufacturing costs than bipolar chips. MOS chips
further increased in complexity at a rate predicted by Moore’s law, leading to
large-scale integration (LSI) with hundreds of transistors on a single MOS
chip by the late 1960s.

###  III-B Moore’s law

MOS technologies obey an empirical law, which is a self-fulfilling prophecy,
stated in 1965 and known as Moore’s law [5]: the number of transistors
integrated on a chip doubles every N months. Fig. 3 presents the evolution for
DRAM memories, processors (MPU) and three types of read-only memories from
1959 to 2005. The growth rate decreases with years, from a doubling every 12
months to every 18 months (DRAM) and 24 months for microprocessors. The semi-
log scale clearly highlights the exponential growth of the number of
transistors per chip. Evolution after 2005 is presented in Fig. 4. Even if N
tends to increase, Moore’s law is still there in 2022.

Figure 3: Moore’s law presented in [7]

Figure 4: Moore’s law slowdown

###  III-C Scaling rules

Until mid-2000’, MOS technologies follow scaling rules known as Dennard
scaling [6]. With every technology generation:

  * •

Transistor dimensions could be scaled by –30% (0.7x).

  * •

Their area reduces by 50%, because area is length*width.

  * •

To keep the electric field constant, the voltage, V, is reduced by 30% (0.7x),
because voltage is field*length.

  * •

Circuit delays reduce by 30% (0.7x), because time is length*velocity.

  * •

The 30% reduction in delay allows an increase in operating frequency, f, by
about 40% (1.4x), because frequency varies as 1/delay.

  * •

The 30% reduction in all distances and related 50% drop in areas lead to a
decrease in capacitance, C, by 30% (0.7x), because capacitance varies as area
/distance.

  * •

Power consumption in turn decreases by 50%, because dynamical power is
proportional $V^{2}f$ (equation (1)

Therefore, with every new node, area halves and power consumption halves. As
the transistor density doubles, power consumption with twice the number of
transistors stays the same. The limits to this perfect scaling are discussed
in III-E2

###  III-D Two equations

We now present two equations that allow to understand the evolution of MOS
technologies and the performance of CPUs. First equation:

$P_{d}=P_{s}+P_{d}=V_{dd}*I_{leakage}+\alpha.\Sigma C_{i}*V_{dd}^{2}*F$ (1)

MOS power dissipation can be broken down into static and dynamic powers. For
dynamic power, $V_{dd}$ is the supply voltage, F is the clock frequency,
$\Sigma C_{i}$ is the combined gate and interconnection capacitances while
$\alpha$ is the average percentage of circuit capacitances that are switching
at each clock cycle: in other words, $\alpha$ is the activity factor of the
overall circuit.

Second equation:

$CPUtime=IC*CPI*Tc=\frac{IC}{IPC*F}$ (2)

  * •

IC is the instruction count of a program

  * •

CPI is the average Clock cycles Per Instruction and IPC = 1/CPI is the
Instruction count Per Clock cycle

  * •

Tc is the clock cycle time and F=1/Tc is the clock frequency

To reduce static power dissipation, a first technique consists in replacing
pMOS and nMOS circuits for which a resistor-like transistor is used as a load
by complementary n and p transistors (CMOS). The first MOS circuits used pMOS
transistors with 12V power supplies. pMOS circuits were replaced by nMOS
circuits that were widely used for computers in the 70’s. The first CMOS logic
families appeared in the late 60’s. At that time, CMOS circuits didn’t exhibit
any static power dissipation.

To reduce dynamic power dissipation, one technique consists in reducing
$V_{dd}$ as this power dissipation is proportional to $V_{dd}^{2}$. This
explains the evolution of power supplies for the successive nodes, as shown in
Fig.5.

Figure 5: Evolution of power supplies

###  III-E Free lunch

####  III-E1 Perfect scaling?

The scaling rules of CMOS technology allowed the successive technological
nodes (Fig. 1). These successive nodes leaded to a wonderful period sometimes
called “free lunch". It is illustrated by Fig. 6. From 1986 to 2003 (Intel 386
to Intel Pentium 4), there was a 200x increase of clock frequency, 200x
increase of transistor count, 11x reduction of transistor channel length and
1000x increase of performance according to Specint2000 benchmark.

Figure 6: Free lunch according to Intel

The performance increase mainly based on three factors:

  * •

Clock frequency (F): According to the equation (2), execution time is
inversely proportional to F. The scaling from one node to the next one reduces
the propagation delays and allow to increase clock frequency.

  * •

Instructions per Cycle (IPC). According to equation (2), execution time is
inversely proportional to IPC (instruction parallelism). This is achieved by
replacing scalar CPUs with a simple pipeline (IPC <1) by superscalar CPUs that
launch and execute several instructions per clock, with in-order or out-of
order execution. The basic scheme for out-of-order execution is presented in
Fig. 7

  * •

SIMD parallelism: SIMD execution of instructions is illustrated in Fig 8.
Arithmetic instructions compute simultaneously several sub-words, reducing the
instruction count for programs exhibiting data-parallelism.

Figure 7: Out-of-order execution in superscalar CPUs

Figure 8: SIMD execution of instructions

####  III-E2 Limits of perfect scaling

It can be noticed that limits of perfect scaling have appeared in this period.

  * •

Clock frequency: Increasing clock frequency is the simplest way to increase
performance. The successive nodes of CMOS technologies lead to a 1.4x decrease
in terms of gate delays. This is why there was a yearly 25% increase in clock
frequencies, from 740 kHz (Intel 4004) to 3 GHz (Intel Xeons with 45-nm
nodes). But power dissipation is proportional to clock frequency. And
increasing clock frequency contributes to the “Heat Wall" discussed in the
next subsection.

  * •

IPC: there are limits to reachable instruction parallelism, as illustrated in
Fig. 9. Intel (and AMD) superscalar CPUs break IA-32 and Intel64 CISC
instructions down into micro-operations (called µops), which are simpler
instructions similar to the RISC instructions of RISC instruction sets. For
Intel CPUs, from 1995 to 2013, the hardware resources (physical registers,
ROB, reservation stations) have grown significantly. However, this has only
resulted in the maximal µop count per cycle to increase from 3 to 4. The
additional resources are used to enlarge the number of µops considered for
simultaneous launching without changing the maximum number of µops per clock
(law of diminushing return).

Figure 9: Limits of instruction parallelism

Figure 10: Power density for successive nodes [8]

###  III-F The Heat Wall

####  III-F1 Clock Frequency

Without any change in micro-architectural features, moving from one node to
the next allowed a gain in frequency and a subsequent gain in performance even
without any CPI change. So, why clock frequencies are limited in the 5-6 GHz
range in 2022 while 45-nm nodes already operated at 3-GHz clock frequencies in
2008? With 2022 5-nm technologies, it is quite obvious than tens of GHz would
be possible. The answer is simple. Equation (1) tells us that the dynamic
power dissipation is proportional to the clock frequency and power density
exponentially increased from one generation to the next one. This is known as
the “heat wall” first coined in [8] and illustrated by Fig. 10. Consequently,
the clock frequencies must be limited just to keep the energy budget limits
compatible with a safe operation of the components and the typical IC packages
and cooling techniques. For higher frequencies, costlier cooling techniques
must be used, as shown by the IBM z14 CPU (5.2 GHz) that is water-cooled. The
“heat wall” is the main reason for the shift towards multicore architectures,
with more, simpler cores replacing single large cores.

####  III-F2 The shift towards multicores

Fig. 9 has outlined the limits of monoprocessors. Fig. 11 has been presented
by Intel. The figure compares a small core, a large core (4x larger than the
small one) and a multicore with four small cores. As a first approximation,
power consumption is proportional to the area, while performance is
proportional to $\sqrt{area}$. The 4-small core has 2x the performance of the
1-large core with the same power consumption. The 4-small core and the 1-small
core have the same performance/power ratio. Obviously, we don’t consider here
the issues of parallel programming.

According to equation (2), parallel processing is also a way to reduce IC. In
case of perfect speed-up, IC is divided by the number of cores!

Until mid-2000’s, the parallel machines have been confined to the mainframes
and the supercomputers. Then, after IBM, Intel and AMD launched the first bi-
cores, multicores have become mainstream and are now being used in the whole
spectrum of computing devices from the low-end to the high-end application
domains.

####  III-F3 Techniques to reduce power dissipation

Many different techniques have been used to reduce power dissipation, at
technological, circuit or architectural levels.

Old CMOS circuitry had no significant static power consumption. The $V_{dd}$
supply voltage has been reduced for a long time to the minimum value that
allows the transistors to operate with correct “on” and “off” states. Until
hitting the “heat wall”, CMOS power consumption was mainly due to the clock
frequency on one hand, and to $\alpha*\Sigma{C_{i}}$ that increases with the
number of transistors and interconnections that the successive technology
nodes allow. The situation has changed starting with 90-nm nodes due to
increased leakage currents. The current going through an “off” transistor is
several orders of magnitudes less than the current traversing an “on”
transistor. However, the “off” current is not zero. With more and more
transistors, the static power consumption becomes increasingly important and
can dominate the overall power consumption, as shown in Fig. 12. Intel tri-
gate (Fig. 13) shows how a technological improvement can reduce the leakage
current by one order of magnitude.

At circuit level, a chip can be decomposed into different domains. As power
dissipation is proportional to $V_{dd}$, different power domain can be defined
to adapt power supplies to the circuit needs. Power dissipation is also
proportional to clock frequency. All parts of the circuits do not have the
same speed requirements, leading to the decomposition into different clock
domains.

There are many other techniques that we do not present here.

Figure 11: Performance and Power of monocores and multicores

Figure 12: Static versus dynamic power dissipation

Figure 13: Intel trigate transistors

Figure 14: Domain decomposition

Figure 15: Evolution of cache hierarchies from 1986 to 2014

Figure 16: Computations versus Memory Accesses

###  III-G The memory wall

Processors and DRAM memories both use semiconductor technologies, but they are
different. Performance of DRAM access times and bandwidth depend on the DRAM
sizes. This situation has leaded to a situation coined as the “Memory Wall",
which is a growing gap between processor and memory performance.

####  III-G1 Caches

One popular technique consists in using different levels of cache, with a
growing complexity, as illustrated in Fig. 15 for microprocessors launched
between 1986 (Intel 386) and 2014 (IBM Power 8). The number of cache levels
increases from 1 to 4, only the last level being off-chip. For multicores with
shared memory organization (such as the commercial multicore multiprocessors),
cache coherency issues increase significantly. We do not consider here the
distributed memory architectures.

####  III-G2 Memory bandwiths

The second technique consists in improving memory bandwidths. Fig. 16
presented by Texas Instruments shows the evolution of computing capability
(GFLOP/s)/memory bandwith (GWORD/s. More cores in multicores need more memory
bandwith, which increases less than the computing capabilities. Successive
generations of DDRAMs and HBMs try to improve bandwith. Improving memory
bandwith is a key issue. As an example, the successive generations of IBM
Power Systems shows a significant increase of maximal bandwith: 65 GB/s
(Power7 in 2010), 210 GB/s (Power8 in 2014), 650 GB/s (Power9 in 2020), 800
GB/s (Power10 in 2021).

####  III-G3 Processing in Memory

Computing with Memory suppresses the drawbacks of transfers between CPU and
Main Memory through a hierarchy of caches. This technique has been studied for
a while, and commercial devices are now proposed for applications such as
Inference for AI.

###  III-H Discussion on MOS technologies trends

The history of MOS technologies since the 60’s and the corresponding evolution
of electronics devices, including microprocessors and computer architecture
would need several books. In this paper, we have broadly sketched what we
believe to be the salient elements of this history.

One point is sure. The evolution we described is based on Moore’s law, i.e.,
an exponential function. As quoted by G. Moore in 2003, “No Exponential is
Forever: But “Forever” Can Be Delayed”. The same idea is developed by P.
Gardini, Intel Fellow and chairman of ITRS [7]

  * •

First lesson: “Predictors of Engineering Limits have Always been Proven Wrong
by The Right Improvements”.

  * •

Second lesson: “It Would be Wrong to Believe that the Right Fundamental Limits
Don’t Exist”.

There had been many examples on how the “first lesson" has been verified:

  * •

Reduction of $V_{dd}$

  * •

Trigate transistors to limit leakage currents

  * •

Switch from monoprocessors to multicores

  * •

Etc.

“Second lesson" outlines the limits of Moore’s law. However, it can be
observed than new technical advances are still announced:

  * •

In May 2021, IBM announced the first 2-nm chip technology [9] that allows to
integrate 50 billions transistors.

  * •

In December 2021, IBM and Samsung announced the VTFET that is a vertical FET
when the FinFET transistor is an horizontal one. Paper [10] presents some
features of this new transistor: “VTFET nanosheet and scaled finFET device
simulation results are compared at the same footprint and at an aggressive
sub-45nm gate pitch. VTFET nanosheets provides 2x performance of the scaled
finFET at equivalent power due to VTFET maintaining good electrostatics and
parasitics while finFET performance is impacted by severe scaling constraints.
Or VTFET could provide as much as 85% power reduction compared to the scaled
finFET architecture as compared at an equivalent frequency on the extrapolated
power-performance curves"

Whether these announces are evolution or breakthrough is a question. Will they
be an actual path to 1-nm node is an open question. However, this shows that
research of Right Improvements is continuing.

The exponential cost of IC manufacturing, including the costs of fabs, is
probably the feature that is the most illustrative of the approaching limits.
Better performance means new nodes. Fig. 17 shows the capability of the most
significant IC manufacturers to produce new nodes. By 2021, only TMS and
Samsung were able to provide the most performing nodes. The Intel troubles
with 10-nm nodes also outline the issue. There are technological issues, but
also financial issues. In 2021, TMSC dominated 55% of the market share of
semiconductor foundry companies. The second one (Samsung) had 17% market
share. The number of competitors is reducing.

Figure 17: Consequences of Fabs costs

##  IV Technological breakthrough: Quantum Computers

Quantum computers are an actual breakthrough compared to classical computers.
There are been a lot of announcements of significant performance gaps:

  * •

Google announced it has a quantum computer (Sycamore) that is 100 million
times faster than any classical computer in its lab. [13]

  * •

In November 2021, IBM announce a quantum computer with 127 qubits [14]

  * •

Quantum computers have also been developped in China [15]

Orders of magnitude performance improvement is the key point for the emergence
of quantum computers. At the same time, software limits are noticed [16].

###  IV-A Algorithms

Several quantum algorithms would solve specific problems, such as factoring
integers, exponentially faster than any known classical algorithm. However,
there are many others, such as playing chess, proving theorems, scheduling air
flights, for which quantum computers would suffer from the same algorithmic
limitations as classical computers. They would surpass conventional computers
only slightly [16]. This is why they would be used as coprocessors for the
exponentially faster applications.

###  IV-B Hardware

Quantum physics is totally different from classical physics, but operational
conditions for quantum devices are also completely different from the
conventional integrated circuits. We only list a first features of the
50-qubit IBM quantum computer:

  * •

The qubits are processed at 15 mK. Different stages operate at 4 K, 800 mK,
100 mK and 15 mK, very close to the absolute zero.

  * •

The quantum processor is located inside a shield to protect it from
electromagnetic radiation

  * •

The coaxial line between the first and second amplifying stages are made out
of superconductors.

  * •

Quantum amplifiers inside a magnetic shield capture and amplify processor
readout signal while minimizing noise.

Fig. 18 shows that this quantum computer has few similitudes with the
computing devices that we currently use. IBM and D-Wave quantum computers
operate at temperature close to 0o K. Another approach developed by Paskal
operates at room temperature by using neural atoms activated by lasers.
Quantum computing is a must and the evolution is quite fast. However, even
with more and more applications, it will remain a niche. Either with close to
0o K temperatures or with lasers, quantum computers or devices are truly
different of the whole world of personal and IoT devices.

Figure 18: The IBM 50-qubit computer

##  V Other technologies?

We now examine some technologies that are not actual potential candidates for
Post-Moore era. Some of them have features that make them useful among
currently used technologies. However, they have not the needed features in
terms of performance gap, integration density, scaling properties, costs, etc.
to be considered as candidate for Post Moore era.

###  V-A Magnetic technologies

Magnetic technology has been extensively used at the beginning of computing
history. From 1955 to 1975, magnetic-core memories were the main memory of
mainframes. These core memories were non-volatile. They were replaced by the
volatile semi-conductor DRAMs. Since the 90’s, different types of non volatile
magnetic RAMS have been developed and fabricated: Ferroelectric RAM (FRAM),
Magnetoresistive RAM (MRAM), Phase-Change memory (PCM)… Some techniques are
mature. A 1Gb MRAM has been launched by Everspin. Samsung announces MRAM-based
in-memory computing. These developments are very promising. Mature non-
volatile memories are perfect candidates to replace SRAM-based caches. In-
memory computing is still in the research phase. While these magnetic devices
find their place in the current electronics world, they do not seem to be a
global alternative to the main stream technologies.

###  V-B Optical computing

Optical computing has been studied for dozens of year. Paper [21] has got 217
paper citations. The conclusion of the abstract was: “Finally, the current
limitations and future needs of optical logic devices and digital optical
computing systems are outlined". Optical fibers and devices are widely used
for telecommunications. For example, they are the backbone of the Internet
infrastructure. But there has not been significant advances for computing.
Similar prospects and challenges of optical computing are described in the
blog [22]. It doesn’t mean that breakthroughs are impossible. In [23], a new
optical switch up to 1000x faster than transistors is presented. However, the
author Pavlos Lagoudakis, a physicist at the Skolkovo Institute of Science and
Technology in Moscow precise: “It took 40 years for the first electronic
transistor to enter a personal computer and the investment of many governments
and companies and thousands of researchers and engineers." He says: “It is
often misunderstood how long before a discovery in fundamental physics
research takes to enter the market." In other words, these advances are far
from being a solution for Moore end.

###  V-C CNTFET technologies

A lot of papers using 32-nm CNTFET technology have been published in the last
decade, especially for m-valued circuits. It turns out that this technology
lacks of integration density.

In 2012, IBM announced a breakthrough in nanotube computer chip fabrication
[17]. The circuit had 10,000 CNTFETs. The IBM researchers announced: “Carbon
nanotubes have the potential in the development of high-speed and power-
efficient logic applications. However, for such technologies to be viable, a
high density of semiconducting nanotubes must be placed at precise locations
on a substrate” and “This new placement technique is readily implemented,
involving common chemicals and processes, and provides a platform for future
CNTFET experimental studies". Other breakthroughs were announced in 2017 [18],
but it looks like IBM no longer focuses on CNTFET technology. In 2013, the
first carbon nanotube computer has been announced [19] by a Stanford group. It
was a significant advance for this technology. However, this 178 CNTFETs “one-
instruction-set computer" only operated at 1 KHz. The first commercial
microprocessor (Intel 4004) had 2300 transistors and ran at 780 KHz in 1971.
In 2019, a 16-bit RISC microprocessor has been built with 14,000 CNFET
transistors [20] by the same Stanford group. While this is a significant
advance for CNTFET technology, we may observe that the Intel 8086 CPU, which
was also a 16-bit microprocessor, has been launched in 1978 with 29,000
transistors, more than 40 years ago! In 2019, the largest transistor count in
a commercially available microprocessor was 39.54 billion MOSFETs, in AMD’s
Zen 2 fabricated using TSMC’s 7 nm FinFET semiconductor manufacturing process.
In 2020, the largest transistor count in a GPU (NVidia Ampere) was 54 billion
transistors with the same 7 nm process. The CNTFET 16-bit microprocessor
manufacturing process had 5 metal layers, while the number of metal layers in
nano-CMOS technologies ranges from 8 to 15, with a trade-off between
integration and cost. It does not seem that CNTFET technology can compete with
the most recent FinFET technologies.

##  VI New computing paradigms

Within the exponential increase of Integrated Circuits performance, new
computing paradigms have also appeared. The main feature is a significant
performance improvement.

###  VI-A Graphic Processing Units (GPUs)

Figure 19: GPU-CPU performance

GPUs have initially been designed for computer graphics. In the early 2000’s,
they started being used for scientific computation. The main reason is shown
in Fig. 19 presented by NVidia. An order of magnitude difference can be
observed both for Peak double precision FLOPS (floating point computation) and
Peak Memory Bandwidth. The cohabitation is not simple. CPUs and GPUs have two
different computing models. CPU uses large-grain parallelism while GPU uses
fine-grain parallelisme to exploit data parallelism. The programming models
are different and GPUs use specific application programming interfaces such as
CUDA (NVidia) or OpenCL. A GPU operates as a coprocessor. Tranfers between CPU
and GPU must be handled. While programming is more complex, the performance
advantage is such that 7 out of the top10 positions of the November 2021
TOP500 ranking were supercomputers using GPUs as accelerators.

###  VI-B FPGAs

Field Programmable Gate Arrays (FPGAs) appeared mid 80’s as a way to design
circuits by programming the interconnects of gates assembled in an array.
Thanks to Moore’s law, the transistor count of most efficient FPGAs is now
tens of billions. Some hardware features have been added to the array of logic
cells such as memory blocks, multipliers, hard microprocessors, I/O peripheral
circuits… and even AI tensor blocks. They constitute a cost effective approach
to get high performance devices without using the ASIC approach. They are
often part of Systems on Silicon (SoC). However, FPGAs belongs to the Moore
era: their future will depend on Moore’s law future.

###  VI-C Deep Neural Networks (DNNs)

Neural networks were first proposed in 1944 by W. McCullough and W. Pitts, two
University of Chicago researchers. Since that period, neural networks have
been continuously studied, but the breakthrough arrived in the early 2010’s
with Deep Neural Networks (DNNs). Fig. 20 shows a specific type of DNN called
CNN for image recognition. A first step consists in several levels of
operators (Convolution, ReLU, Pooling) to reduce the input data set. The
second step consists in several levels of fully connected layers to realize
the classification.

There are many types of applications using DNNs such as speech recognition,
image recognition, natural langage processing, etc. DNN is a very hot topics
in hardware design with:

  * •

New specific instructions in General Purpose Instruction Set such as Intel64

  * •

Specialized operators such as Tensor and Tensor cores in NVidia GPUs or Intel
FPGAs

  * •

Neuronal processors such as ARM Ethos, Google TPU, Intel Nirvana or Xilinx
Versal AI core

DNN networks are programmable (parameters of convolutions, pooling ; weights
of the neurons). However, this programming is dedicated to each application
and is different from programming general purpose microprocessors. DNN are
application-specific architectures.

While DNNs are a very successful approach for a lot of applications, their
success is based on Moore’s law, i.e. they profit from the huge number of
available transistors to improve performance.

Figure 20: A DNN for image recognition

###  VI-D Multivalued circuits

Multivalued circuits are sometimes presented as the solution for the post-CMOS
era [2]. Opposed to GPUs, FPGAs or DNNs, multivalued circuits are not a
solution. Assuming a perfect scaling between binary and multivalued circuits
performance, the gain would be $log(M)/log(2)$, i.e. x1.585 for ternary
circuits or x2 for binary circuits. However, comparing binary and multivalued
combinational circuits shows that the multivalued circuits are always more
complicated: more interconnects, more chip area, more power, etc [25]. There
is a fundamental reason: whatever is the electrical support of the different
values (voltage, current, charges) that is used, the different values are
totally ordered. At the mathematical level, it means that the Boolean algebra
(binary circuits) must be replaced by one variant of the Post algebras.

We illustrate the issue with a small example. Just consider a ternary system
in which 0<1<2\. All ternary designers use the binary functions NI and PI
presented in Table I. It means that the ternary variable is decoded into two
binary variables. After some binary computations, the binary result is encoded
as a ternary result. These ternary to binary decoding and binary to ternary
encoding means that the ternary combinational circuits are always more complex
than the binary ones that processes the same amount of information.

Any technology using a totally ordered set of values is faced to this issue.
There exists one successful multivalued technology: quantum devices and
computers because qubits have multiple unordered values.

Multivalued circuits are reduced to a very small niche [26]. In 2022, the only
significant use is for Flash Memories: as access time is not critical for
these devices used in USB memories, larger access times are compensated by
smaller chip areas.

TABLE I: NI and PI binary functions | NI | PI
---|---|---
0 | 2 | 2
1 | 0 | 2
2 | 0 | 0

##  VII Concluding remarks

While the end of Moore’s law is coming, we examined the lessons of the last 60
years to try to find a successor to the actual dominant semiconductor
technology. More precisely, we consider the required conditions for a
technology to replace the currently used one. We also consider the conditions
for a new computing paradigm to emerge, as a dominant one or as a complement
to the usual paradigms.

A first observation that a new technology should bring a significant
performance improvement, at least one or several orders of magnitude. Quantum
devices and computers are the perfect example. While they operate close to
0°K, while they look more like the ENIAC than to our desktops, laptops or
smartphones, they can solve specific problems can would need weeks or years of
computation with the best supercomputers. At the same time, they are efficient
for specific problems.

Without that huge performance gap, a new technology can be successful if it
has significant features such as cost, size, power and scaling properties. The
shift from bipolar technology to MOS technology corresponds to that situation.
The scaling rules that we described allows MOS technology, which was slower
than bipolar technology at the beginning to finally overcome bipolar
technology.

Within the Moore’s law era, the same remarks can be made for new computing
paradigms. GPUs performance advantage made them good candidates as
coprocessors of CPUs. They cannot replace CPUs. The emergence of Deep Neural
Networks is another good example. Based on a totally different computing
paradigm, they provide exceptional performance on a large range of specific
applications, without being a general purpose approach.

Performance, cost, size, power, scalability are the key features needed for a
Post Moore technology. Will such a technology be available at the end of the
Moore’s law ? It is difficult to predict the future! However, it seems even
more difficult now to find a good candidate.

## References

  * [1] J. Shalf, « The future of computing beyond Moore’s Law », Phil. Trans. R. Soc. A378:20190061. http://dx.doi.org/10.1098/rsta.2019.0061.
  * [2] Z. T. Sandhie, J. A. Patel, F. U. Ahmed, M. H. Chowdhury, “ Investigation of Multiple-valued Logic Technologies for Beyond-binary Era ", ACM Computing Surveys, January 2021 Article No.: 16 https://doi.org/10.1145/3431230
  * [3] D. Etiemble, “Coprocessors: failures and successes", https://arxiv.org/abs/1907.06948
  * [4] RM Russell, “The CRAY-1 Computer System," Comm. ACM, January 1978, pp. 63–72.
  * [5] G.E. Moore, “Cramming more components onto integrated circuits" Electronics Magazine. https://newsroom.intel.com/wp-content/uploads/sites/11/2018/05/moores-law-electronics.pdf
  * [6] R. H. Dennard, F. H. Gaensslen, H. -N. Yu, V. L. Rideout, E. Bassous and A. R. LeBlanc, “Design of ion-implanted MOSFET’s with very small physical dimensions," in IEEE Journal of Solid-State Circuits, vol. 9, no. 5, pp. 256-268, Oct. 1974, doi: 10.1109/JSSC.1974.1050511.
  * [7] P. Gargini, “The Roadmap to Success: 2013 ITRS Update", 2014 IEEE Seminars.
  * [8] F. Pollack, New Microarchitecture Challenges in the Coming Generation of CMOS Process Technologies, Micro32 conference key note, 1999
  * [9] “IBM Unveils World’s First 2 Nanometer Chip Technology, Opening a New Frontier for Semiconductors", https://newsroom.ibm.com/2021-05-06-IBM-Unveils-Worlds-First-2-Nanometer-Chip-Technology,-Opening-a-New-Frontier-for-Semiconductors
  * [10] https://research.ibm.com/blog/vtfet-semiconductor-architecture
  * [11] M. Tyson, “IBM and Samsung announce VTFET, “a breakthrough in semiconductor design”" Club 386, https://www.club386.com/ibm-and-samsung-announce-vtfet-a-breakthrough-in-semiconductor-design/
  * [12] Kerney, “The US semiconductor road to recovery: now or never", https://www.kearney.com/communications-media-technology/article/?/a/the-us-semiconductor-road-to-recovery-now-or-never
  * [13] E. Gibney, “Google published landmark Quantum supremacy claim", Nature, Vol 574 , 24 October 2019, https://media.nature.com/original/magazine-assets/d41586-019-03213-z/d41586-019-03213-z.pdf
  * [14] Eagle, “IBM Quantum breaks the 100-qubit processor barrier" https://research.ibm.com/blog/127-qubit-quantum-processor-eagle
  * [15] C.Q. Choi, “Two of World’s Biggest Quantum Computers Made in China: Quantum computers Zuchongzi and Jiuzhang 2.0 may both display quantum primacy over classical computers", IEEE Spectrum, 06 Nov 2021
  * [16] S. Aaronson, “The limits fo Quantum", Scientific American, pp. 62-69, March 2008.
  * [17] H. Park, A. Afzali, SJ. Han, et al. “High-density integration of carbon nanotubes via chemical self-assembly". Nature Nanotech 7, 787–791 (2012). https://doi.org/10.1038/nnano.2012.189
  * [18] https://www.ibm.com/blogs/research/2017/07/cnt-breakthroughs/
  * [19] M.M. Shulaker, G. Hills, N. Patil, H. Wei, H.Y Chen, H.S. Philip Wong and S. Mitra, “Carbon nanotube computer", Nature, Vol 501, 26 September 2013, doi: 10.1038/nature12502.
  * [20] G. Hills, C. Lau, A. Wright et al. “Modern microprocessor built from complementary carbon nanotube transistors", Nature 572, 595–602 (2019) doi:10.1038/s41586-019-1493-8
  * [21] A. A. Sawchuk and T. C. Strand, “Digital optical computing," in Proceedings of the IEEE, vol. 72, no. 7, pp. 758-779, July 1984, doi: 10.1109/PROC.1984.12937.
  * [22] “Optical Computing: Prospects and Challenges", https://www.findlight.net/blog/2019/02/01/optical-computing/
  * [23] C.Q. Choi, New Optical Switch up to 1000x Faster Than Transistors: “Optical accelerator” devices could one day soon turbocharge tailored applications", IEEE Spectrum, 15 Oct 2021
  * [24] SPRH, “Understanding Deep Learning: DNN, RNN, LSTM, CNN and R-CNN", https://medium.com/@sprhlabs/understanding-deep-learning-dnn-rnn-lstm-cnn-and-r-cnn-6602ed94dbff
  * [25] D. Etiemble, “Common Fallacies about Multivalued Circuits". Asian Journal of Research in Computer Science, 12(4), 67-83. https://doi.org/10.9734/ajrcos/2021/v12i430295
  * [26] D. Etiemble, “Why M-Valued Circuits are restricted to a Small Niche”, in Journal of Multiple Valued Logic and Soft Computing, Vol. 9, No1, 2003.


# General Relativistic magneto-hydrodynamical simulations of accretion flows
through traversable wormholes

Luciano Combi1,2, Huan Yang1,2,3, Eduardo Gutierrez4,5, Scott C. Noble6,
Gustavo E. Romero7,8, Manuela Campanelli9 1Perimeter Institute for Theoretical
Physics, Waterloo, Ontario, Canada, N2L 2Y5
2Department of Physics, University of Guelph, Guelph, Ontario, Canada, N1G 2W1
3 Department of Astronomy, Tsinghua University, Beijing 100084, China
4Institute for Gravitation and the Cosmos, The Pennsylvania State University,
University Park PA 16802, USA
5Department of Physics, The Pennsylvania State University, University Park PA
16802, USA
6Gravitational Astrophysics Lab, NASA Goddard Space Flight Center, Greenbelt,
MD 20771, USA
7Instituto Argentino de Radioastronomía (IAR, CCT La Plata, CONICET/CIC),
C.C.5, (1984) Villa Elisa, Buenos Aires, Argentina
8Facultad de Ciencias Astronómicas y Geofísicas, Universidad Nacional de La
Plata, Paseo del Bosque s/n, 1900 La Plata, Buenos Aires, Argentina
9Center for Computational Relativity and Gravitation, Rochester Institute of
Technology, Rochester, NY 14623, USA

###### Abstract

We present the first dynamical model of plasma accretion onto traversable
wormholes by performing General Relativistic magneto-hydrodynamical (GRMHD)
simulations of the flow on both sides of the wormhole. We evolve the ideal MHD
equations on a wormhole spacetime described by the spherically symmetric
Simpson–Visser metric. The disk is initialized on one side of the wormhole and
accretes onto the throat driven by the magneto-rotational instability (MRI).
We show that the inflowing plasma quickly settles in the throat and forms a
hot, rotating cloud. The wormhole cloud acts as an engine in which gas coming
from one side accumulates at the center, dissipates energy, and powers a
mildly relativistic thermal wind toward the other side. Our novel predictions
show that accreting wormholes behave very differently from black holes (BHs)
in astrophysical environments. In particular, one mouth presents outflows
without accretion signatures, contradicting the jet-disk symbiotic relation
that holds for black holes.

††preprint: APS/123-QED

##  I Introduction

A wormhole is a topological shortcut in spacetime connecting separate regions
of the Universe through a throat. These hypothetical, exotic objects provide a
rich conceptual framework to understand physics in curved spacetime and its
most extreme consequences [1]. A stable, traversable wormhole violates the
averaged null energy condition [2], which sets strong constraints on their
existence [3]; this usually requires non-standard physics like exotic matter
with extreme equations of state [4] or alternative theories of gravity[5, 6,
7, 8, 9]. Their radial and non-radial stability, on the other hand, depend on
specific wormhole configurations [10, 11, 12, 13, 14]

Despite the strong constraints on fundamental physics, the existence of
natural wormholes cannot be completely ruled out. Many properties of the
strong gravitational field near a BH horizon can be mimicked by the
horizonless spacetime of a wormhole outside the throat [15], including the
spectrum of thin disks surrounding one mouth [16, 17], and gravitational waves
[15, 18]; in other aspects, such as lensing of light rays, wormholes can show
distinctive features from normal BHs [19, 20]. All previous work on
astrophysical plasmas surrounding wormholes, to our knowledge, only considered
simple stationary models for matter, such as Novikov–Thorne disks surrounding
one mouth, and mostly focusing on one side of the throat.

In this work, we investigate the dynamics of magnetized plasma surrounding a
wormhole by performing GRMHD simulations of an accretion disk in traversable
wormhole spacetimes. To this purpose, we use the GRMHD code Harm3d which is
capable of evolving the ideal MHD equations on an arbitrary metric with
arbitrary coordinates. We use the Simpson–Visser metric as a simple model of a
two-way traversable wormhole with no rotation, assumed to be sustained by
exotic matter, and similar to a Schwarzschild BH from the outside. We analyze
the plasma dynamics on both sides and across the throat for different throat
sizes. We use geometrical units ($c=G=1$) and a $+2$ signature.

##  II Wormhole metric

We are interested in analyzing accretion onto a simple traversable wormhole
model. We focus on spherically-symmetric wormholes, which are completely
characterized by a shape function (which defines the geometry of the throat)
and a redshift factor (which defines the acceleration of stationary frames)
[1]. For this purpose, we consider the Simpson–Visser metric [21], a simple
extension of the Schwarzschild metric:

$ds^{2}=-f(r)dt^{2}+f(r)^{-1}dr^{2}+(r^{2}+\ell^{2})d\Omega^{2},$ (1)

with $f(r):=\Big{(}1-{2M}/{\sqrt{r^{2}+\ell^{2}}}\Big{)}$, and
$d\Omega^{2}=d\theta^{2}+\sin(\theta)^{2}d\phi^{2}$. The solution
parameterized with $\ell$ represents a two-way wormhole for $\ell>2M$, a
black-bounce solution for $\ell<2M$, and a Scwharzschild BH for $\ell=0$. For
$\ell>0$, the solution is regular at $r=0$, and for $\ell>2M$ it can also be
extended to negative values of $r$. Although the coordinate
$r\in[-\infty,+\infty]$ does not represent the proper size, it maps out the
entire spacetime. The wormhole solution has a traversable throat at $r=0$,
which is a time-like spherical surface of proper radius $\ell$. We define
mouth A as the positive entrance to the throat, $r>0$, and mouth B as the
negative side with $r<0$. The spacetime on both mouths is symmetric and
attractive on each side. The innermost stable orbit (ISCO) of the spacetime is
located at $r_{\rm ISCO}=6M\sqrt{1-(\ell/6M)^{2}}$ and the photon ring is at
$r_{\rm ph}=3M\sqrt{1-(\ell/3M)^{2}}$ [22]. In General Relativity, this
solution requires exotic matter near the throat; we assume that the field
supporting the wormhole is not dynamically affected by the accreting plasma,
nor vice-versa.

##  III Simulation set-up

We evolve the GRMHD equations using the flux-conservative code Harm3d [23],
which is a 3D MPI-parallelized, fully covariant, modernized version of the
public code HARM [24, 25]. The code allows us to simulate an accretion disk on
an arbitrary spacetime with arbitrary coordinates [26]. For this application,
we build a spherical grid symmetric on each side of the throat with more
resolution close to the throat. The uniform-spaced numerical coordinates
$(x_{1},x_{2},x_{3})$ are mapped to physical coordinates as
$r(x_{1})=c_{r}\sinh(x/c_{r})$, where we fix $c_{r}=1-\log{2}/\log{R_{\rm
out}}$, with $R_{\rm out}=600$, and
$\\{\theta(x_{2}),\phi(x_{3})\\}=\\{x_{2},x_{3}\\}$. We use a resolution of
$n_{r}\times n_{\theta}\times n_{\phi}=460\times 128\times 128$, with $240$
radial cells on each side of the wormhole. The radial extension of the grid is
$r\in[-600M,600M]$ and we impose outflow conditions on the outer boundaries.
We cover the full $(\theta,\phi)$ sphere with our grid and we impose
transmissive boundary conditions on the angular boundaries [27].

We assume the wormhole is accreting gas through mouth A from a hot thick
torus, as found in advection-dominated accretion flows (ADAF), either at low
accretion rates, such as M87, or at high accretion rates, such as some X-ray
binaries and tidal disruption events. Mouth B of the wormhole is initially set
to be gas-free. The equation of state of the plasma is $u=p/(\Gamma-1)$, where
$u$ is the internal energy and $p$ is the pressure. The adiabatic index is
chosen as $\Gamma=4/3$, consistent with a radiation-dominated flow. We do not
include any radiative cooling; all heat generated in dynamic processes, e.g.
turbulence, shocks, and magnetic heating, increases the gas entropy.

The initial data consists of a disk in hydrostatic equilibrium surrounding
mouth A with a power-law angular momentum profile following the procedure in
Refs. [28, 29]. We set the inner radius of the torus at $r=25M$ and the
maximum pressure of the disk at $r=35M$. Notice that at $r\approx 10M$, the
difference between the wormhole metric (Eq. 1) and the Schwarzschild metric is
less than $0.1\%$ for $\ell=2.1M$. In addition, we initialize a weak poloidal
magnetic field with a strength set by ensuring the ratio of the volume-
averaged gas and magnetic pressures is $100$, and small perturbations to the
internal energy to accelerate the magneto-rotational (MRI) instability and the
plasma’s development into non-axisymmetric turbulence. In mouth B, we set all
MHD variables to atmosphere values. We simulate two spacetime configurations
by setting the smaller throat size $\ell=2.1M$ and the larger throat size
$\ell=3M$. For $\ell=2.1M$, the ISCO is at $5.2M$ and there is a photon ring
at $\approx 2.1M$; for $\ell=3M$, the ISCO is at $\sim 5.6M$ and the spacetime
does not have a photon ring.

We focus our discussions in the next sections on the model with the smaller
throat $\ell=2.1M$. We will compare with the bigger throat in the last
section.

Figure 1: Mass contained in the throat ($|r|<5M$), mass flux and magnetic flux
measured at $r=\pm 5M$ as a function of time for our fiducial case $\ell=2.1$.
Figure 2: Surface-average properties of the accretion throat as a function of
radius, averaged in time between $t\in[10,20]\times 10^{4}M$ for the case
$\ell=2.1M$. We show the rest-mass density of the whole flow,
$\langle\rho\rangle$, rest-mass density of the bound fluid, $\langle\rho_{\rm
bound}\rangle$, where integrate only fluid elements with $-hu_{t}<1$), and
density-weighted temperature, $\langle T\rangle$, magnetization,
$\langle\sigma\rangle$ and entropy, $\langle s\rangle$. We plot the density of
a spherically symmetric stationary wind, $\rho_{\rm WS}\propto r^{-2}$.

Figure 3: Top: $\theta$-averaged quantities at $t=20\times 10^{4}M$ projected
onto an embedding diagram of the wormhole for $\ell=2.1M$. The center of the
throat exhibits a hot overdensity— a wormhole cloud. Bottom: Meridional
($x=0$) rest-mass density plot at mouth A (left) and mouth B (right) with
magnetic field lines.

##  IV Results

###  IV.1 Accretion from disk to throat

The onset of turbulence by the MRI triggers accretion of an initially
equilibrated torus onto mouth A (Fig. 1). The gas initially flows freely
through the throat ($r=0$) and emerges from mouth B. Because the infalling gas
is bound on average ($-hu_{t}<1$, where $h$ is the specific enthalpy), most of
it precipitates back into the throat. This fall-back material then shocks
against the inflowing gas, drastically increasing the temperature and density
in the throat by several orders of magnitude with respect to the initial
atmosphere values, as shown in the mass evolution integrated over $|r|<5M$
around the throat (Fig. 1). Since $r=0$ is a saddle point of the effective
gravitational potential, the gas sloshing in the throat dissipates its kinetic
energy and accumulates there. The surface-averaged rest-mass density across
the throat $\langle\rho\rangle:={\int\rho d\mathcal{A}}/{\int d\mathcal{A}}$,
with $d\mathcal{A}:=\sqrt{-g}d\theta d\phi$, grows $\sim 4$ orders of
magnitude above that of the disk (Fig. 2). Similarly, the temperature in the
throat increases by $\sim 3$ orders of magnitude with respect to the disk. In
the top panel of Fig. 3, we show an embedding diagram of the wormhole for
different $\theta-$averaged quantities illustrating how matter accumulates in
the throat.

On mouth A, the over-density forming on the throat acts as a hard-surface
inner boundary condition to the accretion disk, restructuring the flow in the
first $t\approx 10^{4}M$; after this period, the mass flux, $\dot{M}:=-\int
u^{r}\rho d\mathcal{A}$, on mouth A reaches a new inflow equilibrium (at
$r<15M$) with an average value of $\langle\dot{M}\rangle_{A}\approx 2\times
10^{-4}$ (Fig. 4I) The mass flux on mouth B settles into an outflow steady-
state with an average value of $\langle\dot{M}\rangle_{B}\approx 3\times
10^{-5}$. The imbalance in mass flux across mouths A and B indicates that part
of the matter is deposited in the throat, consistent with the linear growth of
mass (see Fig. 2).

Figure 4: We plot different quantities for $\ell=2.1M$ in thick lines and
$\ell=3M$ in dashed lines. (I): Time-averaged mass flux (accretion on $r>0$
and outflow on $r<0$) as a function of radius for simulations with $\ell=2.1M$
($\dot{M}_{0}=2\times 10^{-4}$) and $\ell=3M$ ($\dot{M}_{0}=0.75\times
10^{-4}$). (II): Specific internal energy time-averaged. (III): Specific
angular momentum (blue lines) and angular velocity (purple lines) as a
function of radius, averaged on the sphere and in time between
$t\in[10,20]\times 10^{4}M$, weighted with rest-mass density

###  IV.2 An engine for thermal winds

The over-density at the throat forms a hot, rotating bound structure extending
to both sides with an average radius of $|r_{\rm cloud}|\approx 5M$ as shown
in Figs. 3 (bottom panels) and 2. This “wormhole cloud” is differentially
rotating with an average angular velocity,
$\langle\Omega\rangle_{\rho}:=\langle u^{\phi}/u^{t}\rangle_{\rho}$ peaking
around $r\approx 5M$, while decreasing towards the center (Fig. 4III). Outside
the cloud, $\langle\Omega\rangle_{\rho}$ follows a Keplerian profile $\propto
r^{-3/2}$, whereas the interior obeys $d\langle\Omega\rangle_{\rho}/dr>0$ (a
$d\langle\Omega\rangle_{\rho}/dr<0$ profile would decelerate the core and not
be sustainable), similar to neutron stars remnants in binary neutron star
mergers [30]. At the throat, $\Omega$ flattens and transitions to a sub-
Keplerian profile on mouth B. The average specific angular momentum $\langle
l\rangle_{\rho}:=\langle-u_{\phi}/u_{t}\rangle_{\rho}$ in mouth A follows a
Keplerian profile in the bulk of the disk (Fig. 4III), starts decreasing
inside the ISCO ($r=5.2M$), and peaks again at the center of the throat. This
can be understood as the effect of a net positive torque exerted by the outer
envelope of the cloud toward the center as consistent with
$d\langle\Omega\rangle_{\rho}/dr<0$ in the interior.

Figure 5: Equatorial ($z=0$) rest-mass density plot at mouth B for
$\ell=2.1M$. The wind has $m=1$ density modes triggered by oscillations in the
wormhole cloud.

The cloud receives matter from the accretion disk in mouth A and launches a
net outflow from mouth B. Because the flow coming from mouth A is originally
bound on average and the spacetime is symmetric on both sides of the throat, a
mechanism to release energy must be in place to power the outflow. Indeed,
this extra energy comes from the rest-mass density deposited in the throat,
shown in the imbalance of mass flux across the throat (Fig. 4I). The wormhole
cloud functions as an engine: matter deposited in the throat on mouth A is
compressed, loses kinetic energy, and heats up the cloud. The rest-mass
density in the throat increases, and a portion of the fluid in mouth B gains
thermal energy in the process, producing a wind (Fig. 3, bottom right).

###  IV.3 Properties of the wormhole wind

The outflow on mouth B has an asymptotic velocity of
$v_{\infty}=\sqrt{1-1/(-hu_{t})^{2}}\approx 0.2c$ becoming unbound at
$r\approx-8\>M$ as shown in Fig. 4II. Accreted material from mouth A passing
through the throat increases its temperature and entropy due to shocks (see
Fig. 3, top panel, and Fig.2); the entropic outflow on mouth B then cools down
as it expands (Fig. 2). Non-axisymmetric density modes in the wormhole cloud
endow this wind with an $m=1$ spiral pattern (Fig. 5). The wind on mouth B has
a sub-Keplerian-specific angular momentum profile, which remains almost
constant as the wind expands. The surface-averaged properties of the wormhole
wind are described fairly well by a spherically-symmetric stationary wind
(dotted-dashed lines in Fig. 2) that decays as $\rho_{\rm WS}\propto r^{-2}$
for large distances as the wind reaches a constant terminal velocity (see
Appendix). The efficiency of the wormhole wind can be defined from the jump of
mass flux (Fig. 4I), as
$\eta:=\langle\dot{M}\rangle_{B}/\langle\dot{M}\rangle_{A}\times 100\sim
20\%$. .

###  IV.4 Magnetic field in the throat

Magnetic fields are necessary to produce turbulent MRI-driven stresses in the
disk, triggering accretion from mouth A and advecting magnetic flux into the
throat. The magnetic field strength increases in the center of the throat as
it is compressed by gravity, see upper left diagram in Fig. 3. The poloidal
magnetic flux (using Gaussian units), $\Phi:=(1/2)\int|B^{r}|d\Omega$, reaches
a steady-state (Fig. 1) with a mass-flux normalized average value of
$\sqrt{4\pi}\Phi/\sqrt{\dot{M}}\approx 15-20$ on mouth A. Part of this
magnetic flux is trapped in the throat and part of it emerges on the other
side; indeed, on mouth B we observe a magnetic flux that is on average a
factor of $4-5$ less than the flux on mouth A. The spacetime is non-rotating
and thus the system does not produce strong Poynting luminosities by a
Blandford–Znajeck type process, which depends on $\Phi$. If the system has an
ergosphere, e.g. as in the Teo wormhole metric [31], these results show that
one side (mouth A) would likely produce a Poynting outflow — due to magnetic
reconnection in the current sheets [32]— with a larger luminosity than the
other side. The global topology of the field is conserved across $r=0$ but
with an inverse polarity (Fig. 3, bottom right). Winds launched from mouth B
are weakly magnetized, with average values of $\sigma=(b^{2}/2)/\rho\approx
10^{-3}$, see Fig. 2.

###  IV.5 Comparison between different throat sizes

In the simulation with throat’s size $\ell=3M$, a similar wind-emitting cloud
is formed, albeit with some noticeable differences. The cloud is bigger for
$\ell=3M$ and the wind becomes marginally unbound further away, at a typical
radius of $r\sim 14M$ (see Fig. 4II). The conversion of gravitational energy
to kinetic/thermal energy is expected to be less efficient in a wormhole with
a bigger throat because of the shallower effective potential. However, it is
easier to unbind a fluid element that is trapped in the potential if the
throat is bigger. We find that the efficiency for emitting winds in the
simulation with $\ell=3M$ is actually smaller, $\eta\sim 5\%$ (see Fig. 4I).

In $\ell=3M$, the flow does not form a strong shock on the throat and does not
heat up as much, transitioning from mouth A to B smoothly, see Fig. 6. In
$\ell=2.1M$, because of the strong heating across the throat, the wind is
launched quasi-isotropically on mouth B (Fig. 3, bottom right), while in
$\ell=3\>M$ the wind is colder and propagates closer to the equatorial plane,
i.e. with a smaller scale height, see Fig.7.

Figure 6: Surface-avareged rest-mass density, entropy, and temperature for
$\ell=2.1$ (thick lines) and $\ell=3M$ (dashed lines).

Figure 7: Meridional ($x=0$) rest-mass density plot at mouth A (top) and mouth
B (bottom) for $\ell=3M$.

Even though the outflow is less massive than in the $\ell=2.1M$ case, we find
that asymptotic outflow velocities are larger for $\ell=3M$ with
$v_{\infty}\sim 0.3c$, as shown in the specific internal energy in Fig. 4.
This difference is explained by the larger rotational energy of the outflow in
$\ell=3M$. Indeed, in the $\ell=2.1M$ case, matter with high specific angular
momentum is trapped in the throat, while in the $\ell=3M$ case, most of the
angular momentum is carried by the wind, as shown in Fig. 4III, dashed lines.
Notice that the rotational profile for $\ell=3M$ is symmetric with respect to
the throat, contrary to $\ell=2.1M$.

##  V Astrophysical consequences of accreting wormholes

In our simulations, we assumed that the total mass of the accreting gas is
negligible compared to the mass of the wormhole. The predicted linear growth
of the density in the throat, however, implies that the gas can eventually
become dynamically important, e.g. triggering the collapse of the wormhole
into a BH. Converting code units to physical units (see Appendix) we find that
$\rho_{\rm th}\sim 1$ gr cm-3 $(\dot{M}/\dot{M}_{\rm Edd})(M_{\odot}/M)$, and
temperatures of $T_{\rm th}\sim 10^{11}K$ (which might be sufficient to ignite
hydrogen fusion for stellar-mass systems). The mass of the gas in the throat
then would evolve as

$M_{\rm th}\sim\frac{4\pi}{3}r_{\rm th}^{3}\rho_{\rm th}\sim
10^{-6}M_{\odot}\frac{\dot{M}}{\dot{M}_{\rm
Edd}}\frac{M_{\odot}}{M}\frac{t}{\rm years},$ (2)

which can reach $M_{\rm th}\sim M$ in less than a Hubble time for stellar-mass
systems. This implies that astrophysical wormholes might have a finite life
span.

Different from an accreting BH, a wormhole does not have an absorbing
boundary. Stationary accretion models on wormholes have focused on the
properties of the disk emission near the ISCO [16, 17, 33] and ray-tracing
assuming the light can pass through the throat [22]. They ignore, however, the
self-interaction of the flow in the throat, which changes the global
properties of the system and its EM emission. On mouth B, the initially
optically thick outflow expands and eventually becomes optically thin,
producing thermal emission.

In our accretion model, all heat is retained in the flow, and for high-
accretion rates, the gas would be optically thick to radiation. As the
wormhole wind on mouth B expands and adiabatically cools, the plasma will
eventually become optically thin. Assuming a stationary mildly-relativistic
wind, where $u^{r}=\Gamma v^{r}\approx 0.2$ as extracted from the simulation,
the photosphere radius, $r_{\rm ph}$, is placed where the optical depth
becomes one

$\tau=\int^{\infty}_{r_{\rm ph}}\kappa_{\rm e}\rho_{\rm B}(r)dr=1,$ (3)

where $\rho_{\rm B}$ is the density at mouth B ($r<0$). We then find that

$r_{\rm ph}\sim\frac{\dot{M}_{\rm wind}\kappa_{\rm e}}{4\pi v^{r}}\approx
15\>\frac{GM}{c^{2}}\times\Big{(}\frac{\dot{M}}{\dot{M}_{\rm edd}}\Big{)},$
(4)

where $\dot{M}_{\rm wind}=0.25\dot{M}$ is the mass-flux of the wind, $\dot{M}$
is the accretion rate onto mouth A, and $\kappa_{\rm e}$ is the electron
scattering opacity. At the photosphere, we assume that the effective
temperature of thermalized photons is set by the total thermal energy:

$u=\rho\frac{k_{\rm B}}{m_{\rm p}(\Gamma-1)}T_{\rm eff}+aT_{\rm eff}^{4},$ (5)

where $a$ is the radiation constant. We find that, at the photosphere, $T_{\rm
eff}\sim 10^{7}K$, which peaks on the far UV-X-rays for stellar-mass systems;
here we plugged the numerical value of $u$ and $\rho$ from our simulation,
$M=M_{\odot}$, and $\dot{M}=\dot{M}_{\rm edd}$.

The formation of a hot overdensity pouring on both sides of the throat will
drastically change the electromagnetic radiation compared to an accreting BH.
In particular, on mouth A —the accreting side— there will be no shadow because
the cloud size is larger than the light ring, contrary to previous
expectations [34].

##  VI Conclusions

We have performed the first 3D GRMHD simulations of accreting traversable
wormholes. We found that the magnetized plasma accumulates in the throat and
forms a rotating hot cloud which is $\sim 10^{4}$ denser than the disk. The
wormhole cloud acts as an engine, taking gas from the accretion disk on one
side, depositing part of it in the throat, and launching a mildly relativistic
thermal wind on the other side with an efficiency of $\sim 20\%$ for a throat
of size $\ell=2.1M$. These results entail novel predictions for its EM
radiation: the hot rotating cloud in the throat will shine beyond the light
ring of the spacetime and thus there will be no shadow, even if the outer
metric is similar to a BH. We estimate that the cloud wind will produce
thermal emission peaking in the UV/X-rays for stellar mass systems. Moreover,
we show that the hot cloud grows linearly in time and could eventually trigger
the collapse of the wormhole into a BH.

Accumulation of gas in the throat occurs because the throat is a minimum
stable point in the gravitational potential, a feature in most traversable
wormholes; we then expect our main predictions to qualitatively hold in other
realizations of spherically symmetric wormholes. We notice, however, that the
wind efficiency depends on the generated heat in the throat; for accretion
flows that are able to cool efficiently, matter will clump in the throat and
the wind could be less powerful. Different initial conditions and the
inclusion of spin might alter some aspects of the system and can be studied
with the techniques introduced here (c.f. Refs.[35, 36] for GRMHDs onto other
exotic objects). Our work shows that simulations are crucial to producing
realistic models of accreting wormholes

###### Acknowledgements.

We thank Daniel Siegel for reading the draft and providing insightful
comments. This research was enabled in part by support provided by SciNet
(www.scinethpc.ca) and Compute Canada (www.computecanada.ca). LC is a CITA
National fellow and acknowledges the support of the Natural Sciences and
Engineering Research Council of Canada (NSERC), funding reference
DIS-2022-568580. H. Y. is supported by the Natural Sciences and Engineering
Research Council of Canada and in part by Perimeter Institute for Theoretical
Physics. Research at Perimeter Institute is supported in part by the
Government of Canada through the Department of Innovation, Science and
Economic Development Canada and by the Province of Ontario through the
Ministry of Colleges and Universities. G.E.R. acknowledges financial support
from the State Agency for Research of the Spanish Ministry of Science and
Innovation under grants PID2019-105510GB-C31AEI/ 10.13039/501100011033/ and
PID2022-136828NB-C41/AEI10.13039/501100011033, and by “ERDF: A way of making
Europe”, by the ”European Union”, and through the “Unit of Excellence María de
Maeztu 2020-2023” award to the Institute of Cosmos Sciences
(CEX2019-000918-M). Additional support came from PIP 0554 (CONICET).

##  VII Appendix

###  VII.1 Properties of the Simpson–Visser wormhole

The spacetime of a spherically symmetric wormhole can be characterized by its
redshift function, $\varphi(r)$, and shape form $b(r)$ [1]. The redshift
factor determines the acceleration of the stationary observer at a fixed $r$
as well as the time of flight from nearby regions as seen by an asymptotic
observer; the shape factor $b(r)$, determines the embedding diagram of the
wormhole. The two functions set the tidal deformations in the throat which
will affect the fluid.

We can analyze time-like geodesics in this wormhole spacetime taking into
account the motion equation for the radial velocity:

$E=\dot{r}^{2}+V_{\rm eff}(r)$ (6)

where $E$ is the conserved energy and the effective potential is:

$V_{\rm
eff}(r):=\Big{(}1-\frac{2\>M}{\sqrt{r^{2}+\ell^{2}}}\Big{)}\Big{(}\frac{L^{2}}{(r^{2}+\ell^{2})}+1\Big{)}$
(7) $\equiv V^{\rm Schw}_{\rm
eff}(r)+\frac{(M/\ell)}{(r/\ell)^{3}}+\mathcal{O}(r^{-4}).$ (8)

Figure 8: Top: effective potential for fixed specific angular momentum $l=3M$
and different throat sizes. Bottom: effective potential for fixed throat size
$\ell=2.1M$ and different values of specific angular momentum.

From this effective potential, we can immediately see that the particle motion
is very similar to a Schwarzschild BH outside the ISCO. A remarkable feature
of the wormhole spacetime is that there is a valley in the potential at $r=0$.
Let us consider the zero specific angular momentum case. In a black hole, a
particle released at rest from a finite distance falls directly into the
horizon. In the case of the wormhole, a particle that is released at a finite
radius $r_{0}$ will fall into the throat, emerge from the other side and go
back, oscillating between $\pm r_{0}$. If a particle is released from
infinity, it will reach infinity on the other side if there is no dissipation
of energy; if the particle loses energy in the process, then it will get
easily trapped in the throat oscillating back and forth until it settles in
$r=0$ where it can rest.

To eject matter to the other side, the fluid must acquire some additional
kinetic energy from MHD processes or thermal energy. In our simulations, the
incoming fluid dissipates kinetic into thermal energy, which is transferred to
an outer portion of the fluid that will escape. If the throat is smaller
(smaller $\ell$), the effective gravitational potential is deeper and the gas
will be able to gain and dissipate more kinetic energy from the gravitational
field, launching a more efficient wind. On the other hand, a smaller throat
also means that the gas is more trapped in the gravitational field, acting
against the wind launching. However, because gas thermal pressure increases
the size of the cloud, the outer parts of the bound structure will gain
sufficient energy to escape as they receive thermal energy from the interior.

###  VII.2 Embedded diagram

We construct the embedding diagram following Refs. [1]. We embed a 2D spatial
sheet of the wormhole metric fixing $\theta=\pi/2$ and $t=0$, with line
element $ds^{2}=f(t)dr^{2}+(r^{2}+\ell^{2})d\theta$, into an Euclidean
3-dimensional spacetime with metric
$ds_{E^{3}}^{2}=dz^{2}+dR^{2}+R^{2}d\varphi$ in cylindrical coordinates.
Assuming a surface $z=z(r)$ and $R=R(r)$, we have:

$ds_{E^{3}}^{2}=\Big{(}\Big{[}\frac{dR}{dr}\Big{]}^{2}+\Big{[}\frac{dz}{dr}\Big{]}^{2}\Big{)}dr^{2}+R(r)^{2}d\varphi.$
(9)

Imposing $ds^{2}=ds_{E^{3}}^{2}$, we have $R(r)^{2}=(r^{2}+\ell^{2})$ and the
following equation for $z$:

$\frac{dz}{dr}=\sqrt{\frac{1}{f(r)}-\frac{r^{2}}{r^{2}+\ell^{2}}}$ (10)

which can be solved numerically to give a surface $z(r)$.

###  VII.3 Spherically-symmetric wind solution for wormholes

Let us consider a radial and stationary ideal fluid with four-velocity
$u^{a}=u^{t}(r)\partial_{t}+u^{r}(r)\partial_{r}$ in our wormhole spacetime
(Eq. 1). Considering the constraint $u^{a}u_{b}=-1$, together with the ideal
gas equation of state $p=(\Gamma-1)u$, our unknown fluid variables are given
by $u^{r}(r)$, the enthalpy, $h(r)=1+(u+p)/\rho\equiv 1+\Gamma u/\rho$, and
the rest-mass density $\rho(r)$. We thus need three equations to solve the
system [37].

We start using the mass continuity equation, $\nabla_{a}(\rho u^{a})=0$,
obtaining

$4\pi\rho u^{r}(r^{2}+\ell^{2})=\dot{M}.$ (11)

A basic feature of a wind solution that has reached a terminal constant
velocity $u^{r}=u_{\infty}$ is that the density decays as $\rho\propto r^{-2}$
at large separations. In a BH spacetime, the density necessarily diverges as
it approaches the horizon; in the case of the wormhole, the density and the
velocity can remain finite in the throat.

We now consider the conservation equation $\nabla_{a}T^{ab}=0$, which can be
separated into a momentum equation (projected to $h_{ab}$) and an energy
equation (projected to $u^{a}$). Since the fluid and spacetime are spherically
symmetric and stationary, the momentum equation is equivalent to the Bernoulli
equation:

$-hu_{t}=-h\sqrt{(u^{r})^{2}+f(r)}=E.$ (12)

The remaining equation is the energy equation $u^{a}\nabla_{a}T^{ab}=0$, which
can be rewritten as the first law of thermodynamics for an adiabatic fluid:

$\frac{d(\rho+u)}{d\rho}=h$ (13)

which gives $u=K\rho^{\Gamma}$. With these equations, we find a solution in
terms of $\dot{M}$, $E$, and K. For a wind solution, we set $\rho_{\infty}=0$,
and thus $h_{\infty}=1$ and $E=\sqrt{(u^{r}_{\infty})^{2}+1}$. On the other
hand $u^{r}(r)=\sqrt{E^{2}/h(r)^{2}-f(r)}$, and then:

$4\pi\rho(r)(r^{2}+\ell^{2})\sqrt{\frac{E^{2}}{1+\Gamma
K\rho(r)^{\Gamma-1}}-f(r)}=\dot{M}.$ (14)

We can analyze the solution in terms of its value at the critical point when
the velocity is equal to the sound speed; taking the differentials from the
Bernoulli and mass conservation equations, we have:

$\frac{du^{r}}{u^{r}}\Big{(}c_{\rm
s}^{2}-\frac{(u^{r})^{2}}{(u^{r})^{2}+f(r)}\Big{)}+$
$\frac{dr}{r}\Big{(}2r^{2}\frac{c_{\rm
s}^{2}}{r^{2}+\ell^{2}}-\frac{rf^{\prime}(r)}{2((u^{r})^{2}+f(r))}\Big{)}=0,$

where we define the relativistic sound speed $c_{\rm s}^{2}=(\rho/h)\partial
h/\partial\rho$. If the solution reaches the critical point at some radius
$r_{\rm c}$ where $v^{r}_{\rm c}=u^{r}_{\rm c}/((u^{r}_{\rm c})^{2}+f(r_{\rm
c}))=c_{\rm s,c}$, then to maintain regularity, the other parenthesis in the
previous equation should be zero, i.e.:

$2\frac{r_{c}^{2}}{r_{c}^{2}+\ell^{2}}c_{\rm s,c}=r_{\rm
c}\frac{f^{\prime}(r_{\rm c})}{u_{t}^{2}},$ (15)

which means

$u^{r}_{\rm c}=\frac{f^{\prime}(r_{\rm c})(r_{\rm c}^{2}+\ell^{2})}{4r_{\rm
c}}=\frac{M}{2\sqrt{r_{\rm c}^{2}+\ell^{2}}},$ (16)

and

$c_{\rm s,c}=\frac{u^{r}_{\rm c}}{1-3u^{r}_{\rm c}}.$ (17)

With these equations, we can check if the solution reaches or not a critical
point. It is interesting to check what happens if the critical point is on the
throat ($r=0$). In that case, $c_{\rm s,c}=M/(2\ell-3M)<1$, while $u^{r}_{\rm
c}=M/(2\ell)<1/4$. In turn, this implies that the density in the throat must
be $\rho(0)=(\dot{M}/M)/(2\pi\ell)$. We leave a complete analysis of the wind
properties of a wormhole for future work.

###  VII.4 Dependence on resolution

In this section, we show that the bulk properties of gas in the simulation are
robust under different resolutions and predict similar behavior for the
inflow-outflow balance across the throat. We investigate three different
resolutions: $180\times 96\times 96$ (LR), $224\times 112\times 112$ (MR), and
$480\times 128\times 128$ (HR) for radial, poloidal, and azimuthal cells. We
evolve the same wormhole accretion disk presented in the main text with
$\ell=3M$ for each resolution for $t\sim 10^{4}M$.

The quality of the resolution to produce and maintain a turbulent state can be
determined by the MRI quality factor, $Q$. Given the fastest growing
wavelength of the linear MRI defined as

$\lambda^{r,\theta}=\frac{2\pi}{\Omega\sqrt{\rho h+b^{2}}}b^{r,\theta}.$ (18)

the quality factor is defined by the number of cells resolving this scale,
$Q^{r,\theta}=\lambda^{r,\theta}/dx^{r,\theta}$. In our HR runs, the quality
factor in different directions is an average $Q^{r,\theta}\sim 10-15$, which
is a good number to maintain a turbulent state, as discussed in Ref. [38]. In
Figure 9 we plot a snapshot of $Q^{\theta}$ for our fiducial run for
$\ell=2.1$. Our LR run resolves only marginally the MRI $Q\sim 5$ but it is
enough to sustain accretion through the duration of the simulation.

Figure 9: Meridional cut of the MRI quality factor for the $\theta$ direction
for the $\ell\sim 2.1$ fiducial run. Figure 10: Surface-avareged rest-mass
density, entropy, and temperature for $\ell=3M$ as a function of radius for
different resolutions. Figure 11: Time evolution of the mass contained in the
throat for $\ell=3M$ at different resolutions.

The time-averaged properties of the flow across the throat change only
slightly with resolution, as shown in Figure 10. In Figure. 11, we plot the
mass on the throat as a function of time for each resolution. The mass flux is
a bit higher for higher resolution, likely due to a better-resolved MRI in the
accretion flow.

###  VII.5 Transformation from code units to physical units

Because the simulations are scale-free, the physical density is determined by
choosing the mass of the spacetime, $M$, and the physical accretion rate,
$\dot{M}$. These can be related to the quantities in code units using the
conservation equation. Although we can choose any value for these two
parameters, our thick-disk simulations are consistent with high or low
accretion rates (with respect to the Eddington rate), which result in
advection-dominated flows. Given the density in code units, $\rho_{\rm CU}$,
and the accretion rate in code units, $\dot{M}_{\rm CU}$, taken from the
simulation as a time average and radial average in $r\in[0,20M]$), the
physical density is

$\rho=\frac{\rho_{\rm CU}}{\dot{M}_{\rm
CU}}\frac{\dot{M}}{c}\Big{(}\frac{GM}{c^{2}}\Big{)}^{-2}$ $=2\times
10^{-3}\textrm{gr cm}^{-3}\frac{\rho_{\rm CU}}{\dot{M}_{\rm
CU}}\Big{(}\frac{\dot{M}}{\dot{M}_{\rm
edd}}\Big{)}\Big{(}\frac{M_{\odot}}{M}\Big{)},$ (19)

and the internal energy is given by:

$u=\frac{\rho_{\rm CU}T_{\rm CU}}{(\Gamma-1)\dot{M}_{\rm
CU}}c^{2}\frac{\dot{M}}{c}\Big{(}\frac{GM}{c^{2}}\Big{)}^{-2},$ (20) $=2\times
10^{18}\textrm{erg cm}^{-3}\frac{\rho_{\rm CU}T_{\rm
CU}}{(\Gamma-1)\dot{M}_{\rm CU}}\Big{(}\frac{\dot{M}}{\dot{M}_{\rm
edd}}\Big{)}\Big{(}\frac{M_{\odot}}{M}\Big{)}$ (21)

Because we have not included any type of cooling, the scale of the gas
temperature is determined by the composition of the material, which we assume
to be electrons. The physical temperature is then calculated from $\rho_{\rm
CU}$ and $u_{\rm CU}$ through the equation of state as:

$T=c^{2}m_{\rm p}k_{\rm B}^{-1}T_{\rm CU}=(\Gamma-1)\frac{u_{\rm
CU}}{\rho_{\rm CU}}\times 10^{12}K,$ (22)

where $c$ is the speed of light, $m_{\rm p}$is the mass of the proton, and
$k_{\rm B}$ is the Boltzmann constant.

## References

  * Morris and Thorne [1988] M. S. Morris and K. S. Thorne, American Journal of Physics 56, 395 (1988).
  * Visser [1995] M. Visser, Woodbury (1995).
  * Lobo [2017] F. S. Lobo,  _Wormholes, warp drives and energy conditions_ , Vol. 189 (Springer, 2017).
  * Visser [1990] M. Visser, Physical Review D 41, 1116 (1990).
  * Visser [1989] M. Visser, Physical Review D 39, 3182 (1989).
  * Maldacena and Milekhin [2021] J. Maldacena and A. Milekhin, Physical Review D 103, 066007 (2021).
  * Blázquez-Salcedo  _et al._ [2021] J. L. Blázquez-Salcedo, C. Knoll, and E. Radu, Physical Review Letters 126, 101102 (2021).
  * Konoplya and Zhidenko [2022] R. Konoplya and A. Zhidenko, Physical Review Letters 128, 091104 (2022).
  * Barcelo and Visser [2000] C. Barcelo and M. Visser, Classical and Quantum Gravity 17, 3843 (2000).
  * Poisson and Visser [1995] E. Poisson and M. Visser, Phys. Rev. D 52, 7318 (1995), arXiv:gr-qc/9506083 .
  * Cardoso and Pani [2019] V. Cardoso and P. Pani, Living Rev. Rel. 22, 4 (2019), arXiv:1904.05363 [gr-qc] .
  * Dias and Lemos [2010] G. A. S. Dias and J. P. S. Lemos, Phys. Rev. D 82, 084023 (2010), arXiv:1008.3376 [gr-qc] .
  * Armendariz-Picon [2002] C. Armendariz-Picon, Phys. Rev. D 65, 104010 (2002), arXiv:gr-qc/0201027 .
  * Yang  _et al._ [2023] H. Yang, B. Bonga, and Z. Pan, Phys. Rev. Lett. 130, 011402 (2023), arXiv:2207.13754 [gr-qc] .
  * Damour and Solodukhin [2007] T. Damour and S. N. Solodukhin, Physical Review D 76, 024016 (2007).
  * Bambi and Stojkovic [2021] C. Bambi and D. Stojkovic, Universe 7, 136 (2021).
  * Harko  _et al._ [2008] T. Harko, Z. Kovács, and F. S. Lobo, Physical Review D 78, 084005 (2008).
  * Cardoso  _et al._ [2016] V. Cardoso, E. Franzin, and P. Pani, Physical review letters 116, 171101 (2016).
  * Safonova  _et al._ [2001] M. Safonova, D. F. Torres, and G. E. Romero, Physical Review D 65, 023001 (2001).
  * Torres  _et al._ [1998] D. F. Torres, G. E. Romero, and L. A. Anchordoqui, Physical Review D 58, 123001 (1998).
  * Simpson and Visser [2019] A. Simpson and M. Visser, Journal of Cosmology and Astroparticle Physics 2019 (02), 042.
  * Bambhaniya  _et al._ [2022] P. Bambhaniya, K. Jusufi, P. S. Joshi, _et al._ , Physical Review D 105, 023021 (2022).
  * Noble  _et al._ [2009] S. C. Noble, J. H. Krolik, and J. F. Hawley, The Astrophysical Journal 692, 411 (2009).
  * Gammie  _et al._ [2003] C. F. Gammie, J. C. McKinney, and G. Tóth, ApJ 589, 444 (2003), arXiv:astro-ph/0301509 [astro-ph] .
  * Noble  _et al._ [2006] S. C. Noble, C. F. Gammie, J. C. McKinney, and L. Del Zanna, The Astrophysical Journal 641, 626 (2006).
  * Zilhao and Noble [2014] M. Zilhao and S. C. Noble, Classical and Quantum Gravity 31, 065013 (2014).
  * Armengol  _et al._ [2022] F. G. L. Armengol, Z. B. Etienne, S. C. Noble, B. J. Kelly, L. R. Werneck, B. Drachler, M. Campanelli, F. Cipolletta, Y. Zlochower, A. Murguia-Berthier, _et al._ , Physical Review D 106, 083015 (2022).
  * Chakrabarti [1985] S. K. Chakrabarti, The Astrophysical Journal 288, 1 (1985).
  * De Villiers  _et al._ [2003] J.-P. De Villiers, J. F. Hawley, and J. H. Krolik, The Astrophysical Journal 599, 1238 (2003).
  * Kastaun  _et al._ [2013] W. Kastaun, F. Galeazzi, D. Alic, L. Rezzolla, and J. A. Font, Phys. Rev. D 88, 021501 (2013), arXiv:1301.7348 [gr-qc] .
  * Teo [1998] E. Teo, Physical Review D 58, 024014 (1998).
  * East and Yang [2018] W. E. East and H. Yang, Phys. Rev. D 98, 023008 (2018), arXiv:1805.05952 [astro-ph.HE] .
  * Harko  _et al._ [2009] T. Harko, Z. Kovacs, and F. S. Lobo, Physical Review D 79, 064001 (2009).
  * Akiyama  _et al._ [2022] K. Akiyama, A. Alberdi, W. Alef, J. C. Algaba, R. Anantua, K. Asada, R. Azulay, U. Bach, A.-K. Baczko, D. Ball, _et al._ , The Astrophysical Journal Letters 930, L17 (2022).
  * Kocherlakota  _et al._ [2023] P. Kocherlakota, R. Narayan, K. Chatterjee, A. Cruz-Osorio, and Y. Mizuno, The Astrophysical Journal Letters 956, L11 (2023).
  * Chatterjee  _et al._ [2023] K. Chatterjee, P. Kocherlakota, Z. Younsi, and R. Narayan, arXiv preprint arXiv:2310.20040 (2023).
  * Bahamonde and Jamil [2015] S. Bahamonde and M. Jamil, The European Physical Journal C 75, 1 (2015).
  * Porth  _et al._ [2019] O. Porth, K. Chatterjee, R. Narayan, C. F. Gammie, Y. Mizuno, P. Anninos, J. G. Baker, M. Bugli, C.-k. Chan, J. Davelaar, _et al._ , The Astrophysical Journal Supplement Series 243, 26 (2019).