# Spacetime continuity and quantum information loss

Michael R. R. Good1<EMAIL_ADDRESS>1 Nazarbayev University, Astana
010000, Kazakhstan

###### Abstract

Continuity across the shock wave of two regions in the metric during the
formation of a black hole can be relaxed in order to achieve information
preservation. A Planck scale sized spacetime discontinuity leads to unitarity
(a constant asymptotic entanglement entropy) by restricting the origin of
coordinates (moving mirror) to be timelike. Moreover, thermal equilibration
occurs and total evaporation energy emitted is finite.

black hole evaporation, information loss, remnants

##  I Introduction

In this note, the role of continuity and information loss (via the tortoise
coordinate, $r^{*}$), in understanding the crux of the phenomena of particle
creation from black holes Hawking:1974sw  is explored. In particular the
relationship to the simplified model of the moving mirror Davies:1976hi ;
Davies:1977yv ; Fulling:2018lez , is investigated because it has identical
Bogolubov coefficients Good:2016oey .

Uncompromising continuity is relaxed in favor of unitarity by a single
additional parameter generalization of the $r^{*}$ coordinate. An
understanding of the correspondence between the black hole and the moving
mirror in this new context is initiated Good:2016atu ; GTC ; Myrzakul:2018bhy
. We prioritize information preservation, and find finite evaporation energy,
thermal equilibrium, and analytical beta coefficients. Moreover, we determine
and assess a left-over remnant (e.g. Wilczek:1993jn ; remnants ).

Broadly motivating, delving into the ramifications of external effects on
quantum fields have led to interesting results on a wide variety of phenomena
all the way from, e.g. relativistic superfluidity super ; oxford  to the
creation of quantum vortexes in analogue spacetimes vortex ; Liberati:2018uev
. The fertile enterprise of QFT under external conditions (like the external
condition of a rotating black hole spacetime) has resulted in simple and
revealing relations, (e.g. springy  $2\pi T=g-k$, temperature of a Kerr black
hole).

Specifically motivating, the physical effect of the black hole origin is that
it amplifies quantum field fluctuations by reflecting virtual particles into
real ones. There is however, a problem: the black mirror Good:2016oey , which
is the aforementioned usual tortoise coordinate associated boundary condition,
extracts energy indefinitely and does not preserve unitarity.

We sketch a possible simple way out - ‘a timelike worldline for the origin’ -
which resolves these two problems (as probed in Good:2016atu ). It was
demonstrated recently GTC ; Myrzakul:2018bhy  that a timelike worldline
compels the introduction of a second parameter, $\xi$, the asymptotic sub-
light speed of the origin (in addition to mass, $M$), which generalizes the
tortoise coordinate and characterizes the dynamics, evaporative energy,
spectra, information, and continuity.

##  II Asymptotic Null Origin

It is first appropriate to survey the situation with information loss. This is
plotted as the usual Penrose diagram in Fig. (2). To that end, the matching
solution for the outside and inside of the black hole over the shock wave is
derived, following the notation of Wilczek Wilczek:1993jn  and in the textbook
of Fabbri-Navaro-Salas Fabbri:2005mw . The Regge-Wheeler coordinate, $r^{*}$,

$r^{*}\equiv r+2M\ln\left(\frac{r}{2M}-1\right),$ (1)

sets the dynamics. In the usual set up, the metric is matched on both sides of
the shock wave, $v_{0}$, inside and out, respectively:

$r(v_{0},U)=r(v_{0},u),$ (2)

where

$r(v_{0},U)=\frac{v_{0}-U}{2},\quad\quad\text{and}\quad\quad
r^{*}(v_{0},u)=\frac{v_{0}-u}{2}.$ (3)

Using Eq. (1), Eq. (2) and Eq. (3) to solve for the trajectory of the origin,
$u$, in null coordinates gives:

$u=U-4M\ln|\kappa U|.$ (4)

Here we define the surface gravity, $\kappa\equiv 1/4M$. As $u\to+\infty$, the
formation of a event horizon occurs at $v_{H}\equiv v_{0}-4M$, where for
simplicity and without loss of physical generality we set $v_{H}=0$ (shock
wave at $v_{0}=4M$). The result, Eq. (4), is the matching solution
Fabbri:2005mw ; Wilczek:1993jn ; GTC  required for the Schwarzschild geometry
(exterior) to the Minkowski geometry (interior) with a strict event horizon.
The regularity condition (the field must be zero at $r=0$) restricts the
behavior of the modes in the interior such that $U\leftrightarrow v$.
Converting to Cartesian coordinate $x$ by the definitions $u\equiv t(x)-x$ and
$v\equiv t(x)+x$ into Eq. (4), and solving for $t(x)$ gives the equation of
motion of the black mirror111The black mirror, which is not early-time thermal
(like the Carlitz-Willey trajectory Good:2012cp ; Carlitz:1986nh ), is often
called Omex after the transcendental Omega, $\Omega e^{\Omega}=1$, and
exponential Good:2016oey . It is also called the Black Hole Collapse (BHC)
trajectory Cong:2018vqx :

$t(x)=-x-4Me^{x/2M},$ (5)

which is studied222Likewise see MG14 proceedingsGood:2015jwa ;
Anderson:2015iga  and the 2nd LeCosPA Symposium Good:2016bsq . in Good:2016oey
. This is plotted as mirror in Fig. (3). The coordinates span $0<r<\infty$ and
$-\infty<x<\infty$. The transcendentally invertible trajectories of Eq. (4)
and Eq. (5) characterize a one-parameter system (the mass $M$) and determine
the late-time (thermo)dynamics which results in information loss. See Fig. (1)
for a visual of how modes traveling through the horizon are ‘lost’.

Figure 1: Left: The 1+1 dimensional Penrose diagram for the black hole
collapse mirror. The red line is the mirror, Eq. (5), computed with $\kappa=2$
for illustration. The dotted black line is the horizon, $v_{H}=0$. The blue
dashed line is the upper-half of the left future null-infinity Cauchy surface,
$\mathscr{I}_{L}^{+}$ (unlabeled), that acts as a singularity analog. Notice
the light green opaque region is where left-moving modes never reflect off the
mirror and get ‘trapped’, never to become right-movers. The red arrows are
Hawking radiation emitted by the mirror or they can be left-movers that
reflected off the mirror into right-movers reaching $\mathscr{I}^{+}$. Right:
The usual 3+1 dimensional Schwarzschild causal structure captured in a Penrose
diagram. The blue dashed line is the space-like singularity while the blue
solid line is $r=0$ where modes pass through but still hit the singularity.
The black dotted dashed line is the event horizon. The red arrows are Hawking
radiation or they can be left-movers that passed through $r=0$, ‘reflected’
and become right-movers. They will eventually escape to the observer at
$\mathscr{I}^{+}$. Notice the light green opaque region where left-movers are
lost into the singularity.

Figure 2: Left: The usual Penrose diagram for an evaporating Schwarzshild
black hole with a singularity (thick dashed horizontal line) and information
loss. The dotted dashed line is the event horizon. The outgoing blue arrow,
$\mathit{f}$, is the Hawking flux and the ingoing black arrow, $v_{0}$, is the
shock wave. Right: Penrose diagram for the 3+1 dimensional case where a black
hole remnant remains with infinite lifetime, for a similar diagram see
remnants . Notice the redshifted outgoing modes labeled as $\mathit{r}$ with a
red arrow, escape to null-infinity $\mathscr{I}^{+}$. Even though Hawking
evaporation stops, all field modes passing through the remnant undergo an
eternal late-time redshift after passing through and out of the remnant region
$\mathcal{R}$. The light green area labeled $\mathcal{P}$ is the Planckian
region where quantum geometric effects are important (discontinuity) and has a
non-zero, maximum, finite height centered at $r=0$. All field modes can
‘reflect’ off of $r=0$ and pass through the $\mathcal{R}$ region without loss
of information.

##  III Asymptotic Timelike Origin

Prioritizing information preservation requires the existence of an asymptotic
comoving observer with the origin. This is a timelike restriction on the black
mirror: its maximum speed must always travel slower than light, even
asymptotically333Interestingly, the speed of light can be timelike approached
if the acceleration asymptotes to zero sufficiently fast, see
Good:2016reflectatlight .. The black hole origin, $r=0$, which corresponds to
the location of the moving mirror, will also be made to be timelike. Similar
to the mirror, the origin, expressed in null coordinates (outside retarded
time as a function of inside retarded time) will not asymptotically approach
null-future infinity, $\mathscr{I}^{+}$. This also means that the situation in
the black hole case resolves without $r=0$ attaching to a space-like
singularity (to do so would mean that $r=0$ is no longer timelike). The result
of a particular timelike trajectory for the origin is that the usual
evaporation process effectively ends with the formation of a remnant, depicted
in the Penrose diagram of Fig. (2). Notice that for the corresponding moving
mirror as pictured in the Penrose diagram, Fig. (1) or Fig. (3), the worldline
never finds itself located at null-future infinity. To make this quantitative,
and to raise the principle of unitarity to the most overarching priority of
the model, we modify Eq. (5) to a new trajectory Good:2016atu  for the origin
of coordinates (one that has an asymptotic drift):

$t_{\xi}(x)=-\frac{x}{\xi}-4Me^{x/2M\xi}.$ (6)

Here $\xi$ is the asymptotic coasting speed, $0<\xi<1$. This is plotted as
mirror in Fig. (3). The trajectory in null coordinates changes dramatically
from Eq. (4) to (using the regularity condition result $U\leftrightarrow v$):

$u_{\xi}=v-\frac{\xi}{\kappa}\ln\left[\frac{\epsilon}{2}\mathcal{W}\left(\frac{2}{\epsilon}e^{-\frac{2}{\epsilon}\kappa
v}\right)\right],$ (7)

where $\mathcal{W}$ is the Lambert function and we have introduced the
definition, $\xi\equiv 1-\epsilon$. When $\xi\to 1$, the trajectory Eq. (7) is
equivalent to Eq. (4), which demonstrates an operative formation of an
effective event horizon. Thus the timelike restriction (by the introduction of
a second parameter $\xi$ in addition to mass $M$) in Eq. (5), resulting in Eq.
(6), permits the genesis of an event horizon 444For all practical purposes. It
is an interesting question whether any type of horizon is formed during
gravitational collapse taking into account quantum effects see, e.g.
Mann:2018jcf . For horizonless models see Good:2017kjr ; Good:2017ddq ..

The origin, $r=0$, as far as the quantum field, $\psi$, is concerned, is a
perfectly reflecting boundary because $\psi(r\leq 0)=0$. Eq. (7) is the
timelike world line of the origin and re-tracing our steps from the matching
condition Eq. (4) to Eq. (1), we arrive at:

$r^{*}_{\xi}\equiv
r+2M\xi\ln\left[\frac{\epsilon}{2}\mathcal{W}\left(\frac{2}{\epsilon}e^{\frac{2}{\epsilon}[\frac{r}{2M}-1]}\right)\right].$
(8)

In the limit that $\xi\to 1$, Eq. (8) $\to$ Eq.(1), $r^{*}_{\xi}\to r^{*}$. No
singularity in Eq. (8) exists at $r=2M$. On the contrary, unlike Eq. (1),
using $\xi\neq 1$ and examining $r=2M$, the function takes on a finite value:
$r^{*}_{\xi}=2M\left[1-\xi\mathcal{W}(2/\epsilon)\right]$. The asymptotic
inertial drift, $\xi$, can be close to $\xi\approx 1$ for operative continuity
and still be strictly, $\xi<1$, for information preservation. At the cost of
classical continuity, we have purchased quantum purity.

Figure 3: Left: Penrose diagram for the case of the moving mirror (solid red
line) which has a one-to-one correspondence to the usual evaporating black
hole with information loss. The solid red line here is Eq. (5), computed with
$\kappa\equiv 1/4M=1$. Notice that the mirror goes light-like in the far
future. The dashed blue line is the singularity analog where left-movers
(opaque yellow region) get ‘trapped’ i.e. they fall into the black hole analog
never to return or be observed by $\mathscr{I}_{R}^{+}$. Right: Penrose
diagram for the case of the moving mirror (solid red line) which corresponds
to the remnant end-state with no information loss. Ultra-late time outgoing
reflected modes escape to null-infinity, $\mathscr{I}^{+}_{R}$, with a
redshift dependent on the final asymptotic speed of the mirror, $1>\xi>0$,
where in the diagram, $\xi=0.9$, for illustration. In this case, the solid red
line, Eq. (6), is calculated with $\kappa=1$ and is always time-like,
permitting all left-movers to reflect and become right-movers.

##  IV How Discontinuous?

The discontinuity is located in the Planckian region of Fig. (2), where one
expects quantum geometric effects to play an important role and spacetime
discontinuity to possibly occur. This discontinuity can be characterized at
late times by first considering the whole spacetime using the interior
coordinates Wilczek:1993jn ,

$\operatorname{d}\\!{s}^{2}=\begin{cases}-\operatorname{d}\\!{U}\operatorname{d}\\!{V},~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\text{for}~{}v\leq
v_{0},\\\
-f(u,v)f^{-1}(u,v_{0})\operatorname{d}\\!{U}\operatorname{d}\\!{V},~{}~{}\text{for}~{}v>v_{0}.\end{cases}$
(9)

The difference at late times, $u\to\infty$, of

$\delta s^{2}=1-f(u,v)f^{-1}(u,v_{0}),$ (10)

is a measure of the discontinuity. The difference, located at the shock wave,
$v=v_{0}$, utilizing Eq. (7) and its inverse, $v_{\xi}(u)$, is

$\delta
s^{2}|_{u\to\infty}=\lim_{u\to\infty}\left[1-\frac{dv_{\xi}}{du}\frac{du_{\xi}}{dv}\right]=2\frac{\epsilon-1}{\epsilon-2}\left(1-\frac{1}{W\left(\frac{2}{\epsilon}e^{-2/\epsilon}\right)+1}\right)$
(11)

This late-time measure helps clarify how tiny the discontinuity $\delta s^{2}$
can be when $\xi$ is ultra-relativistic,

$\delta s^{2}|_{u\to\infty}\approx 0,\quad\textrm{for}\;\xi\gtrapprox 0.9.$
(12)

as seen in Figure 4, which also includes a log plot of $\delta s^{2}$ as a
function of $\xi$.

Figure 4: Left: The discontinuity, Eq. (12), at late times over the shock wave
for ultra-relativistic speed $\xi\gtrapprox 0.9$ is approximately zero. Right:
A log plot of Eq. (12) to demonstrate the rapid fall off at high speeds. Here
$\epsilon\equiv 1-\xi$, and the final coast speed is greater than zero but
less than the speed of light, $0<\xi<1$.

##  V Restriction on Entropy

The limit of the von-Neumann entanglement entropy Chen:2017lum  from the
dynamics of Eq. (5) as a function of lab time $t$ is found from the rapidity
Good:2016atu , $\eta(t)=\tanh^{-1}\dot{z}(t)$, by using $6S(t)=\eta(t)$:

$\bar{S}\equiv\lim_{t\to\infty}S(t)=\lim_{t\to\infty}\frac{1}{6}\textrm{tanh}^{-1}\left[(\mathcal{W}[2e^{-2\kappa
t}]+1)^{-1}\right]=\infty.$ (13)

This is a consequence of the fact that not all field modes hit the origin and
information is lost (see e.g. the $\bar{S}=\infty$ in GoodS ). In contrast,
the asymptotically timelike trajectory of Eq. (6) has the limit,

$\bar{S}_{\xi}\equiv\lim_{t\to\infty}S_{\xi}(t)=\lim_{t\to\infty}\frac{1}{6}\textrm{tanh}^{-1}\left[\xi(\mathcal{W}[2e^{-2\kappa
t}]+1)^{-1}\right]=\frac{1}{6}\textrm{tanh}^{-1}(\xi)\neq\infty.$ (14)

The final asymptotic entropy, $\bar{S}_{\xi}$, is the final rapidity (an
additive quantity appropriate for an entropy). This is a result of the fact
that all the field modes intersect the timelike worldline origin. We will see
in the next sections that while $\bar{S}=\eta/6$ is finite, it will be
extremely large (e.g. relative to usual accelerator rapidities) in order reach
effective equilibrium ($\xi\approx 1$).

In the first case, Eq. (13), some left-mover field modes never hit the origin
and consequently never become right-movers. These modes are lost forever in
the black hole. However, in the second case, Eq. (14), an asymptotic approach
to time-like future infinity, $i^{+}$, means that all the left-movers hit the
origin and they all become right-movers, preserving information measured by an
observer at $\mathscr{I}^{+}$.

##  VI Restriction on Energy

Evaluation of the total energy emitted Good:2016atu ,

$E=\int_{-\infty}^{\infty}F(u)\;du,$ (15)

from the energy flux, $F(u)$, can be done with the Schwarzian derivative
Davies:1976hi ; Davies:1977yv ,

$-24\pi
F(u)=\\{p,u\\}\equiv\frac{p^{\prime\prime\prime}}{p^{\prime}}-\frac{3}{2}\left(\frac{p^{\prime\prime}}{p^{\prime}}\right)^{2},$
(16)

of the null-coordinate trajectory, $p(u)$ (the $v$ position of the origin as
function of $u$)- in our case, the function $p(u)=v_{\xi}(u)$, which is the
inverse of Eq. (7). The result is finite and analytic:

$E=\frac{\hbar\kappa}{24\pi c}\left(\gamma^{2}-\frac{\eta}{\xi}\right),$ (17)

where $\gamma\equiv 1/\sqrt{1-\xi^{2}}$ is the asymptotic coasting Lorentz
factor, $\eta\equiv\tanh^{-1}\xi$ is the asymptotic coasting rapidity, and
$\xi<1$ is the asymptotic coasting speed. [Notice the negative sign which is
correct, as opposed to the plus sign misprint in the published version of this
note]. When $\xi\approx 1$, near the speed of light, then
$\gamma^{2}\gg\eta/\xi$, and the second term can be neglected:

$E\sim\kappa\gamma^{2}.$ (18)

The energy diverges as the origin moves asymptotically null, (i.e. $\xi\to 1$
and $\gamma\to\infty$) behaving like the total divergent energy using the
inverse of Eq. (4), $v(u)$, there $\xi=1$. Consistent with conservation of
energy, Eq. (18) gives the total finite emission when the system has a chance
to equilibrate.

##  VII Mass of the Remnant

The energy of the spacetime is $M$ and the evaporation energy is $E$. The
difference between them is the mass, $m$, of the remnant:

$M-E=m,$ (19)

where for a steady-state system, $\xi\approx 1$, Eq. (18) holds,

$M-\frac{\gamma^{2}}{96\pi M}=m.$ (20)

The two parameters of the system, $\gamma$ and $M$, are simply related if the
mass of the remnant is so far less than the black hole, $m\ll M$, that it may
be neglected:

$\gamma^{2}=96\pi M^{2}.$ (21)

For perspective, a solar mass black hole, $M\sim 10^{38}$ Planck masses, via
$\cosh\eta=\gamma$, has final coasting origin rapidity $\eta\sim 91$. The
fastest particle ever measured, the OMG particle, had rapidity $\eta\sim 27$.

##  VIII Constant Energy Flux and Equilibrium Temperature

In the 1+1 dimensional black hole case, visually depicted in Fig. (5), a
temperature, $2\pi T=\kappa$, gives a constant energy flux,

$F=\frac{\kappa^{2}}{48\pi}=\frac{\pi}{12}T^{2},$ (22)

which emerges in the form of an extended plateau for assignment of high coast
speed, $\xi\approx 1$. A series expansion of the temperature,

$T_{\xi}=\sqrt{\frac{12}{\pi}F^{\textrm{max}}_{\xi}},$ (23)

as a function of maximum energy flux Good:2016atu ,

$F^{\textrm{max}}_{\xi}=\frac{\kappa^{2}}{48\pi}\left[1-3\sqrt[3]{6}\epsilon^{2/3}+O\left(\epsilon\right)\right]$
(24)

(where the radiation is closest to equilibrium) gives the temperature of the
black hole:

$T_{\xi}=\frac{\kappa}{2\pi}\left[1-3\left(\frac{3}{4}\right)^{1/3}\epsilon^{2/3}+O\left(\epsilon\right)\right],$
(25)

to lowest order in $\epsilon\equiv 1-\xi$. The higher order terms correspond
to the very small deviation due to sub-light speed drift, $\xi<1$, which are
safely ignored for small $\epsilon\approx 0$. The timelike modification of Eq.
(5) to Eq. (6) locks in a constant energy flux plateau and effective long-term
thermal equilibrium of Eq. (25) for $\xi\approx 1$. This corroborates the
stability of the steady-state (long-time Planckian distributed particle
radiation) due to the timelike worldline origin.

Figure 5: Left: The usual 1+1 Penrose diagram for a Schwarzschild black hole
with a singularity and information loss. The solid red line is the Cauchy
surface over which integration is done. $v_{H}$ is the event horizon, and
$v_{0}$ is the shock wave. The blue arrow, $f$, is Hawking flux. Right: The
1+1 Penrose diagram for the case of a remnant with infinite life time (solid
red lines), region $\mathcal{R}$. In this situation the remnant end-state
allows modes to pass through it (no information loss). Ultra-late time
outgoing modes escape to null-infinity, $\mathscr{I}^{+}$, with a redshift
dependent on the final asymptotic speed of the origin, $1>\xi>0$. The green
lines indicate the Planckian region $\mathcal{P}$ where non-trivial quantum
spacetime effects are important. The dotted dashed gray line is the location
of the original event horizon. The dashed black line, inside $\mathcal{P}$, is
the location of the original singularity, now hidden. The discontinuity occurs
in $\mathcal{P}$ over the shock wave $v_{0}$ like in the 3+1 case.

##  IX Particle Spectrum

The beta Bogolubov coefficient in the particle spectrum,

$\left\langle N_{\omega}\right\rangle\equiv\langle
0_{\textrm{in}}|N_{\omega}^{\textrm{out}}|0_{\textrm{in}}\rangle=\int_{0}^{\infty}|\beta_{\omega\omega^{\prime}}|^{2}d\omega^{\prime},$
(26)

can be computed Good:2016atu  via an integral with input of Eq. (6). The
result is analytically tractable:

$\beta_{\xi}=-\frac{\xi\sqrt{\omega\omega^{\prime}}}{2\pi\kappa\omega_{p}}\left(\frac{i\kappa}{\omega_{p}}\right)^{A}\Gamma(A),$
(27)

where $A\equiv\frac{i}{2\kappa}[(1+\xi)\omega+(1-\xi)\omega^{\prime}]$ and
$\omega_{p}\equiv\omega+\omega^{\prime}$. How is this Planckian? The integrand
of Eq. (26) is Planckian using Eq. (27) with $\xi\approx 1$,

$|\beta_{\omega\omega^{\prime}}(\xi\approx
1)|^{2}\approx\frac{\omega^{\prime}}{2\pi\kappa(\omega+\omega^{\prime})^{2}}\frac{1}{e^{2\pi\omega/\kappa}-1}\to\frac{1}{2\pi\kappa\omega^{\prime}}\frac{1}{e^{2\pi\omega/\kappa}-1}$
(28)

where in the second step, we took the high frequency approximation,
$\omega^{\prime}\gg\omega$. The pre-factor is an continuum artifact removed
with wave packets Hawking:1974sw . The result Eq. (28) is the usual Planck
distribution consistent with Eq. (25). Further consistency of the result can
be shown by using Eq. (27) to obtain, via a numerical computation, the total
energy,

$E\equiv\int_{0}^{\infty}\hbar\omega\;\left\langle
N_{\omega}\right\rangle\;d\omega\;$ (29)

and comparing this with the total energy, Eq. (17) (both the left and right
sides of the origin) from the stress tensor Good:2016atu ; Good:2013lca . The
results are in agreement.

##  X No Restriction on Total Particle Count: Soft Hair

It is good to mention a limitation and an extension of this approach. Assuming
unitarity from the start, this model regularizes the total emission energy,
the asymptotic entanglement entropy, and the asymptotic rapidity of the
origin. However, the total particle count is still infinite. That is, the
$\xi$-model is still impaired by an infrared divergence in the integrand,
$N_{\omega}$, when computing,

$N\equiv\int_{0}^{\infty}\left\langle
N_{\omega}\right\rangle\;d\omega\;=\int_{0}^{\infty}\int_{0}^{\infty}|\beta_{\omega\omega^{\prime}}|^{2}d\omega^{\prime}\;d\omega,$
(30)

which is the total particle count. Even after the energy has effectively
stopped being radiated, an extreme red-shift occurs. This end-of-time Doppler
shift suffered by the modes which continue to pass through the origin marks
the existence of the remnant, which is characterized by this divergence,

$\textrm{asymptotic drift}:\quad N\to\infty.$ (31)

The modes reflect off a constant, $\xi>0$, boundary condition (albeit one with
asymptotically zero acceleration). The massless scalar particles carry zero-
energy, e.g. H . Generally, this divergence is remedied by abandoning
asymptotically coasting trajectories and replacing them with asymptotically
static trajectories (the final asymptotic velocity would be zero, $\xi=0$),
which leave no evidence of any remnant (zero Doppler shift),

$\textrm{asymptotic static}:\quad N\to\textrm{finite}.$ (32)

An interesting extension of this model would be a construction of such an
asymptotically static dynamic which also has thermal emission and corresponds
to Eq. (5) or Eq. (6). This line of work is in progress.

##  XI Conclusion

Quantity | Continuity $\xi=1$ | Unitarity $\xi\neq 1$
---|---|---
$r^{*}$ | $r^{*}\equiv r+2M\ln\left(\frac{r}{2M}-1\right)$ | $r^{*}_{\xi}\equiv r+2M\xi\ln\left[\frac{\epsilon}{2}\mathcal{W}\left(\frac{2}{\epsilon}e^{\frac{2}{\epsilon}[\frac{r}{2M}-1]}\right)\right]$
$t(x)$ | $t=-x-4Me^{x/2M}$ | $t_{\xi}=-x/\xi-4Me^{x/2M\xi}$
$u$ | $u=v-4M\ln|\kappa v|$ | $u_{\xi}=v-4M\xi\ln\left[\frac{\epsilon}{2}\mathcal{W}\left(\frac{2}{\epsilon}e^{-\frac{2}{\epsilon}\kappa v}\right)\right]$
$\bar{S}$ | $\infty$ | $\bar{S}_{\xi}=\eta/6$
$E$ | $\infty$ | $E_{\xi}=\frac{\kappa}{24\pi}\left(\gamma^{2}-\frac{\eta}{\xi}\right)$
$T$ | $T=\kappa/2\pi$ | $T_{\xi}=\kappa/2\pi+\mathcal{O}(\epsilon),~{}~{}\textrm{for}~{}~{}\epsilon\approx 0$
$\beta$ | $\beta=-\frac{\sqrt{\omega\omega^{\prime}}}{2\pi\kappa\omega_{p}}\left(\frac{i\kappa}{\omega_{p}}\right)^{\frac{i\omega}{\kappa}}\Gamma(\frac{i\omega}{\kappa})$ | $\beta_{\xi}=-\frac{\xi\sqrt{\omega\omega^{\prime}}}{2\pi\kappa\omega_{p}}\left(\frac{i\kappa}{\omega_{p}}\right)^{A}\Gamma(A)$
$N$ | $\infty$ | $\infty$

This note postcedes the work of Good:2016atu ; GTC ; Myrzakul:2018bhy . The
evaporation model prioritizes unitarity over continuity, and one consequence
is that a finite total energy is emitted. Here two analytic answers are found:
the beta Bogolubov coefficient and the finite total energy. The dynamics are
described by the asymptotically timelike worldine, Eq. (7), and the two-
parameter Regge-Wheeler coordinate, Eq. (8).

Uncompromising continuity across the shock wave in the metric results in
information loss. At the Planck scale, at least, continuity is less certain.
High precision in effective continuity is permitted with very fast sub-light
drifting speeds. This relaxation results in preserved information, finite
energy and radiated Planckian distributed particles with constant energy flux
at equilibrium temperature.

## Acknowledgments

M.R.R.G. thanks Yen Chin Ong, Aizhan Myrzakul, and Khalykbek Yelshibekov.
M.R.R.G. also thanks Daniele Malafarina for clarifying points about the
Penrose diagrams. The Ministry of Education and Science of the Republic of
Kazakhstan are also acknowledged.

## References

  * (1) S. W. Hawking, Commun. Math. Phys. 43 199 (1975).
  * (2) S. A. Fulling, P. C. W. Davies, Proc. Roy. Soc. Lond. A 348 393 (1976).
  * (3) P. C. W. Davies, S. A. Fulling, Proc. Roy. Soc. Lond. A 356 237 (1977).
  * (4) S. A. Fulling and J. H. Wilson, Festschrift for Wolfgang Schleich (2018), [1805.01013].
  * (5) M. R. R. Good, P. R. Anderson, C. R. Evans, Phys. Rev. D 94 065010 (2016), [1605.06635].
  * (6) M. R. R. Good, K. Yelshibekov, Y. C. Ong, JHEP 1703 013 (2017), [1611.00809].
  * (7) M. R. R. Good, Y. C. Ong, A. Myrzakul, K. Yelshibekov, Gen Relativ Gravit (2019) 51: 92, [1801.08020].
  * (8) A. Myrzakul, M. R. R. Good, MG15 (2018), [1807.10627].
  * (9) F. Wilczek, IASSNS-HEP-93-012 (1993), [9302096].
  * (10) P. Chen, Y. C. Ong and D. h. Yeom, Phys. Rept. 603, 1 (2015), [1412.8366].
  * (11) C. Xiong, M. R. R. Good, Y. Guo, X. Liu, K. Huang, Phys. Rev. D 90, 125019 (2014), [1408.0779].
  * (12) O. Waldron and R. A. Van Gorder, Phys. Scripta 92, no. 10, 105001 (2017).
  * (13) M. R. R. Good, C. Xiong, A. J. K. Chua, K. Huang,  New J. Phys. 18, no. 11, 113018 (2016), [1407.5760].
  * (14) S. Liberati, S. Schuster, G. Tricella and M. Visser, [1802.04785].
  * (15) M. R. R. Good, Y. C. Ong,  Phys. Rev. D 91, 044031 (2015), [1412.5432].
  * (16) A. Fabbri, J. Navarro-Salas, London, UK: Imp. Coll. Pr., (2005).
  * (17) M. R. R. Good,  Int. J. Mod. Phys. A 28 1350008 (2013), [1205.0881].
  * (18) R. D. Carlitz, R. S. Willey,  Phys. Rev. D 36, 2327 (1987).
  * (19) W. Cong, E. Tjoa, and R. B. Mann, JHEP 1906 021 (2019), [1810.07359 ].
  * (20) P. R. Anderson, M. R. R. Good, C. R. Evans, MG14 (2015), [1507.03489].
  * (21) M. R. R. Good, P. R. Anderson, C. R. Evans, MG14 (2015), [1507.05048].
  * (22) M. R. R. Good, 2nd LeCoSPA Meeting (2017), [1602.00683].
  * (23) M. R. R. Good, World Scientific, Singapore (2018), [1612.02459].
  * (24) R. B. Mann, I. Nagle, D. R. Terno, Nucl. Phys. B 936, 19 (2018), [1801.01981].
  * (25) M. R. R. Good, E. V. Linder,  Phys. Rev. D 96, 125010 (2017), [1707.03670].
  * (26) M. R. R. Good, E. V. Linder,  Phys. Rev. D 97, 065006 (2018), [1711.09922].
  * (27) M. R. R. Good, E. V. Linder,  Phys. Rev. D 99, 025009 (2019), [1807.08632].
  * (28) P. Chen and D. H. Yeom,  Phys. Rev. D 96, 025016 (2017), [1704.08613].
  * (29) M. R. R. Good, P. R. Anderson, C. R. Evans, Phys. Rev. D 88 025023 (2013), [1303.6756].
  * (30) S. W. Hawking, M. J. Perry and A. Strominger,  Phys. Rev. Lett. 116, no. 23, 231301 (2016), [1601.00921].


# Universal Hyperuniform Organization of Cellular Structures in Leaf Vein
Networks

Yuan Liu  Department of Physics, Qufu Normal University, Qufu, Shandong
Province, China 273165  Duyu Chen  Materials Research Laboratory, University
of California, Santa Barbara, California 93106, United States  Yang Jiao
<EMAIL_ADDRESS>Materials Science and Engineering, Arizona State
University, Tempe, AZ 85287  Department of Physics, Arizona State University,
Tempe, AZ 85287  Jianxiang Tian<EMAIL_ADDRESS>Department of Physics,
Qufu Normal University, Qufu, Shandong Province, China 273165

###### Abstract

Leaf vein network is a hierarchical vascular system that transports water and
nutrients to the leaf cells. The thick primary veins form a branched network,
while the secondary veins develop closed circuits forming a well-defined
cellular structure. Through extensive analysis of a variety of distinct leaf
species, we discover that the apparently disordered cellular structures of the
secondary vein networks exhibit a universal hyperuniform organization and
possess a hidden order on large scales. Disorder hyperuniform (DHU) systems
lack conventional long-range order, yet they completely suppress normalized
large-scale density fluctuations like crystals. Specifically, we find that the
distributions of the geometric centers associated with the vein network loops
possess a vanishing static structure factor in the zero-wavenumber limit,
i.e., $S(k)\sim k^{\alpha}$, where $\alpha\approx 0.64$, providing an example
of class III hyperuniformity in biology. This hyperuniform organization leads
to superior efficiency of diffusive transport, as evidenced by the much faster
convergence of the time-dependent spreadability $\mathcal{S}(t)$ to its long-
time asymptotic limit, compared to that of other uncorrelated or correlated
disordered but non-hyperuniform organizations. Our results also have
implications for the discovery and design of novel disordered network
materials with optimal transport properties.

Recently, it was discovered that apparently disordered organizations in
biological systems across scales, such distributions of photoreceptor cones in
avian retina Jiao et al. (2014), self-organization of immune cells Mayer et
al. (2015), and large-scale patterns of vegetation Ge (2023), can possess an
exotic hidden order called “hyperuniformity”. Disorder hyperuniform (DHU)
systems lack conventional long-range order as in an amorphous material, yet
they completely suppress normalized infinite-wavelength density fluctuations
like crystals Torquato and Stillinger (2003); Torquato (2018). In this sense,
DHU systems can be considered to possess a “hidden order” in between that of a
perfect crystal and a completely random system (e.g., an ideal gas). DHU is
characterized by a local number variance $\sigma_{N}^{2}(R)$ associated with a
spherical window of radius $R$ that grows more slowly than the window volume
(e.g., with scaling $R^{d}$ in $d$-dimensional Euclidean space) in the
large-$R$ limit Torquato and Stillinger (2003); Torquato (2018), i.e.,

$\lim_{r\rightarrow\infty}\sigma_{N}^{2}(R)/R^{d}=0.$ (1)

This is equivalently manifested as the vanishing static structure factor in
the infinite-wavelength (or zero-wavenumber) limit, i.e.,

$\lim_{|{\bf k}|\rightarrow 0}S({\bf k})=0,$ (2)

where ${\bf k}$ is the wavenumber and $S({\bf k})$ is related to the pair-
correlation function $g_{2}({\bf r})$ via $S({\bf k})=1+\rho\int e^{-i{\bf
k}\cdot{\bf r}}[g_{2}({\bf r})-1]d{\bf r}$ and $\rho=N/V$ is the number
density of the system. For statistically isotropic systems, the structure
factor only depends on the wavenumber $k=|{\bf k}|$. The small-$k$ scaling
behavior of $S(k)$, i.e.,

$S({k})\sim k^{\alpha}$ (3)

determines the large-$R$ asymptotic behavior of $\sigma_{N}^{2}(R)$, based on
which all DHU systems can be categorized into three classes:
$\sigma_{N}^{2}(R)\sim R^{d-1}$ for $\alpha>1$ (class I);
$\sigma_{N}^{2}(R)\sim R^{d-1}\ln(R)$ for $\alpha=1$ (class II); and
$\sigma_{N}^{2}(R)\sim R^{d-\alpha}$ for $0<\alpha<1$ (class III), where
$\alpha$ is the hyperuniformity exponent Torquato (2018).

Besides the aforementioned DHU examples in biological systems, a wide spectrum
of equilibrium Torquato et al. (2015); Batten et al. (2009); Jiao (2022) and
non-equilibrium Gabrielli et al. (2002); Donev et al. (2005); Zachary et al.
(2011a, b); Yuan et al. (2021) many-body and material systems, in both
classical Kurita and Weeks (2011); Hunter and Weeks (2012); Dreyfus et al.
(2015); Chen et al. (2023a); Zhang et al. (2023); Chremos and Douglas (2018)
and quantum mechanical Torquato et al. (2008); Feynman and Cohen (1956);
Gerasimenko et al. (2019); Sakai et al. (2022); Rumi et al. (2019); Sánchez et
al. (2019); Zheng et al. (2020); Chen et al. (2021, 2022, 2023b) varieties,
have been identified to possess the property of hyperuniformity. Notably,
certain driven non-equilibrium systems Hexner and Levine (2015); Jack et al.
(2015); Weijs et al. (2015); Salvalaglio et al. (2020); Nizam et al. (2021);
Zheng et al. (2023), active-particle fluids Lei et al. (2019); Lei and Ni
(2019); Huang et al. (2021); Zhang and Snezhko (2022); Oppenheimer et al.
(2022), and dynamic random organizing systems Hexner and Levine (2017); Hexner
et al. (2017); Weijs and Bartolo (2017); Wilken et al. (2022), which have been
used as models for, e.g., self-organization of biological cells, school of
fish or flock of birds, are also found to exhibit emergent hyperuniform
behaviors. Novel DHU materials have been engineered that can possess superior
properties compared to their crystalline counterpart, such as high-degree of
isotropy and robustness against defects Florescu et al. (2009); Klatt et al.
(2022); Zhang et al. (2016); Chen and Torquato (2018); Torquato (2021); Xu et
al. (2017); Leseur et al. (2016); Yu (2023); Shi et al. (2023); Torquato and
Chen (2018).

Although it is not completely understood how DHU emerges in many of the
aforementioned systems, the attempt for the system to achieve certain
functional optimality under constraints seems to be a crucial underlying
mechanism Ge (2023). For example, in the case of avian photoreceptors, the
cones of the same type favor maximal mutual separation and formation of
ordered triangular-lattice arrangement that maximizes light sampling
efficiency. However, this ideal arrangement is frustrated by the size
polydispersity, shape irregularity and incompatible number ratios of five cone
types. Therefore, the resulting DHU arrangement is one that optimizes light
sampling under these constraints Jiao et al. (2014). Similarly, it was found
the DHU pattern of arid and semi-arid vegetation ecosystem is likely resulted
from optimal utilization of water resources under constraints Ge (2023). These
exciting discoveries naturally lead to the question: Are there any other
biological systems that have developed hyperuniform organizations for optimal
functionality under constraints?

$\begin{array}[]{c}\\\
\includegraphics[width=210.30595pt]{fig1.png}\end{array}$

Figure 1: Cellular structures in leaf vein network. (a) Optical image of air-
dried ficus religiosa leaf containing a hierarchical vein network. (b) A
portion of looped network formed by the secondary veins in between two primary
vein branches. The scale bar is 0.5 cm. (c) Cellular structures extracted from
the looped network in (b). (d) Point configuration associated with the
geometrical centers of the cells in (c).

In this letter, we report the discovery of universal hyperuniform organization
of the cellular structures of vein networks. Leaf veins are present throughout
the leaf lamina, providing mechanical support of the leaf structure Scarpella
et al. (2010); Biedroń and Banasiak (2018). Another key functionality of leaf
veins is to transport water, minerals and nutrients to the leaf cells, which
are crucial for their metabolism. As shown in Fig. 1, these veins form a
hierarchical network: On the macroscopic scale, the thick primary veins
organize into a branched network, which is very efficient to mechanically
support the leaf structure, see Fig. 1(a). On the other hand, in between the
main branches, the secondary thinner veins form looped networks Katifori and
Magnasco (2012), which contain many cellular structures, see Fig. 1(b). Our
detailed analysis below indicates the distribution of the geometrical centers
associated with these cellular structures, for all of leaf species analyzed
here, possesses a high degree of hyperuniformity, which is characterized by a
vanishing structure factor in the zero-wavenumber limit, i.e., $S(k)\sim
k^{\alpha}$ with $\alpha\approx 0.64$, providing an example of class III
hyperuniformity in biology. Importantly, we show the DHU organization of these
cellular structures of vein networks leads to superior transport behaviors,
compared to other uncorrelated and correlated disordered organization models,
which is crucial to the functionality of the veins.

Sample preparation and image processing. To begin with, we collected and
analyzed a large number of air-dried leaf samples of 6 distinct species
including ficus religiosa (FR), ficus caulocarpa (FC), ficus microcarpa (FM),
smilax indica (SI), populus rotundifolia (PR), and yulania denudate (YD),
commonly found in middle and northern part of Asia, Europe and North America.
We took high-resolution optical images of the samples (see Fig. 1(a) for an
example of ficus religiosa leaf), and use Matlab for image correction and
enhancement to highlight the vein network (Fig. 1(b)). We employed the OpenCV
package in Python to identify the “cells” formed by the veins based on two
criteria: (i) a cellular region is enclosed by a closed loop of veins; and
(ii) a cellular region contains an open-end branched vein in the center. In
particular, such open-end branched structures play a similar role of
capillaries, releasing the minerals and nutrients to individual leaf cells
Scarpella et al. (2010). An example of the resulting cellular structure is
shown in Fig. 1(c), based which the geometrical center of each cell is
extracted to obtain a point configuration (with $N$ varying from $\sim 100$ to
$\sim 400$) for subsequent analysis, see Fig. 1(d). The average number density
of derived point configurations is $\rho=125$ cm-2.

$\begin{array}[]{c}\\\
\includegraphics[width=214.63977pt]{fig2.png}\end{array}$

Figure 2: Statistics of the point configurations with $\rho=125$ cm-2 derived
from the cellular structures of the secondary vein networks of ficus religiosa
leaf, including the static structure factor $S(k)$ (a) and (b), the number
variance $\sigma_{N}^{2}(R)$ (c), as well as the pair-correlation function
$g_{2}(r)$ (d), from 20 independent realizations with $N\approx 200$ randomly
selected from the leaf samples. The quantity $a$ is a characteristic length
scale, e.g., the average nearest neighbor distance in the system. The
corresponding statistics of point configurations for a Poisson process and
random sequential addition (RSA) process with the same number density are also
shown for purpose of comparison.

Hyperuniform organization of cellular structures of vein networks. We obtain
statistics of the point configurations derived from the cellular structures of
the secondary vein networks, including the static structure factor $S(k)$, the
number variance $\sigma_{N}^{2}(R)$, as well as the pair-correlation function
$g_{2}(r)$, by averaging over 20 independent realizations randomly selected
from the leaf samples. To better understand the unique hyperuniform
organizations of these cellular structures, we also consider two benchmark
systems at the same number density: (i) the Poisson distribution of points and
(ii) the distribution of disk centers generated via a random sequential
addition (RSA) process. The former represents a totally uncorrelated
disordered system, while the latter represents a correlated system through the
mutual exclusion effects of the disks.

Figure 2(a) and (b) respectively shows $S(k)$ in the linear and log scale for
the three systems. As expected, the uncorrelated Poisson system possesses a
constant $S(k)=1$, while the RSA system shows a strong suppression of
scattering for small $k$, but is non-hyperuniform. On the other hand, the
point configurations derived from the cellular structures of the leaf vein
network exhibit a high degree of hyperuniformity, characterized by a vanishing
$S(k)$ in the zero-$k$ limit with a hyperuniformity index Torquato (2018)
$H=S(0)/S(k_{p})\approx 10^{-3}$ (where $k_{p}$ is the wavevector associated
the first and highest peak of the structure factor) and a hyperuniformity
exponent $\alpha\approx 0.64$ (illustrated as the dashed line), suggesting
these cellular structures are class III hyperuniform. To understand the
boundary effects, we also carry out a system-size study and confirm that the
same small-$k$ behavior and hyperuniformity exponent can be robustly obtained
in systems with distinct $N$ from $N\sim 150$ to $N\sim 400$ (see Fig. S2 in
Supporting Information), with the latter corresponds to the largest patches of
secondary veins found in between the two main branches.

Figure 2(c) shows the number variance $\sigma_{N}^{2}(R)$ of the three systems
in the log scale. It can be clearly seen that the Poisson and RSA systems
follow the scaling $\sigma_{N}^{2}(R)\sim R^{2}$, while the leaf vein system
possesses the scaling $\sigma_{N}^{2}(R)\sim R^{\beta}$, where $\beta\approx
1.35$, which is consistent with the value estimated based on the
hyperuniformity exponent from the $S(k)$ analysis, i.e.,
$\beta=2-\alpha\approx 1.36$. Figure 2(d) shows $g_{2}(r)$, characterizing the
short-range correlations in the three different systems. Again, the Poisson
system possesses a $g_{2}(r)=1$, indicating its totally random nature.
$g_{2}(r)$ of the RSA system possesses a clear exclusion region, corresponding
to the diameter of the congruent circular disks, which is immediately followed
by a sharp peak, resulted from the nearly contacting particle configurations
in the RSA packing. On the hand other, $g_{2}(r)$ of the vein networks is
reminiscent of that of packing of hard disks with a size dispersity, which is
characterized by a well-defined exclusion region, followed by a smooth
increase and a broad peak. Similar $g_{2}(r)$’s were also observed for the
distribution of the photoreceptors in avian retina Jiao et al. (2014).

$\begin{array}[]{c}\\\
\includegraphics[width=214.63977pt]{fig3.png}\end{array}$

Figure 3: Excess spreadability $\mathcal{S}(\infty)-\mathcal{S}(t)$ of the DHU
leaf vein networks and different structural models. (a) Comparison of
$\mathcal{S}(\infty)-\mathcal{S}(t)$ directly computed from vein network
imaging data and a model structure consisting of the Voronoi network
associated with point configuration derived from the cellular vein networks,
as shown in the inset. (b) Comparison of $\mathcal{S}(\infty)-\mathcal{S}(t)$
of three model structures derived from the Poisson (blue), RSA (red) and vein-
network (green) point configurations by placing congruent circular disks at
the points, as shown in (c). Inset of (b) illustrates the diffusion of solutes
initially confined in certain sub-regions of the system.

Superior transport properties induced by hyperuniformity. As noted above, one
main function of the secondary vein network is to transport water and minerals
to the leaf cells. As can be seen in Fig. 1b, there is an open-end branched
structure in the middle of each “cellular region” of the secondary vein
network, from which the minerals are released and diffused to the individual
leaf cells in this region. Here we employ the time-dependent diffusion
spreadability $\mathcal{S}(t)$ (or simply “spreadability” for short) Torquato
(2021) to quantify the efficiency of diffusive transport supported by the DHU
cellular structures. Specifically, one assumes a fixed amount of solutes are
initially confined within a subset of regions $\mathcal{V}_{2}$ in the system
(e.g., within the veins), and then allows the solutes to start diffusing into
the remaining regions $\mathcal{V}_{1}$ at $t=0$, assuming the diffusivity $D$
is identical for both phases (see inset of Fig. 3b for illustration).
$\mathcal{S}(t)$ is defined as the fraction of total solutes present in
$\mathcal{V}_{1}$ as a function of time $t$, which is a measure of the
spreadability of diffusion information as a function of time. A related
quantity is the excess spreadability $\mathcal{S}(\infty)-\mathcal{S}(t)$,
i.e.,

$\mathcal{S}(\infty)-\mathcal{S}(t)=\frac{1}{(2\pi)^{2}\phi_{2}}\int_{\mathbb{R}^{2}}\tilde{\chi}_{V}({\bf
k})\exp[-k^{2}Dt]d{\bf k},$ (4)

where $\phi_{2}$ is the volume fraction of $\mathcal{V}_{2}$, and
$\tilde{\chi}_{V}({\bf k})$ is the spectral density function associated with
$\mathcal{V}_{2}$, which is the Fourier transform of the corresponding
autocovariance function $\chi_{V}({\bf r})=S^{(2)}_{2}({\bf r})-\phi_{2}^{2}$
and $S^{(2)}_{2}({\bf r})$ is the two-point correlation function of
$\mathcal{V}_{2}$ Shi et al. (2023). In particular, the large-$t$ behavior of
$\mathcal{S}(\infty)-\mathcal{S}(t)$ quantifies how fast the solutes can
achieves the asymptotic uniform distribution in the system, reflecting the
efficiency of diffusive transport in the system.

Figure 3a shows the excess spreadability $\mathcal{S}(\infty)-\mathcal{S}(t)$
of the DHU vein networks directly computed by obtaining the spectral density
function $\tilde{\chi}_{V}({\bf k})$ of the “vein phase” (i.e., the black
regions) from the imaging data. This assumes that the solutes are carried and
released from all veins. Also shown is $\mathcal{S}(\infty)-\mathcal{S}(t)$ of
a model system consisting of the Voronoi network associated with point
configuration derived from the cellular vein networks, including congruent
disks centered at those points. The size of the disks is chosen such that
their area is equivalent to the average area of the open-end branched
structures in the middle of the cellular region in the network. It can be seen
this simple model captures very well both short-time and long-time behavior of
$\mathcal{S}(\infty)-\mathcal{S}(t)$ for the actual vein networks. Notably,
the DHU vein network possesses a long-$t$ asymptotic scaling
$\mathcal{S}(\infty)-\mathcal{S}(t)\sim 1/t^{\gamma}$, where $\gamma\approx
1.38$. This is to contrast the scaling $\mathcal{S}(\infty)-\mathcal{S}(t)\sim
1/t$ for a general nonhyperuniform disordered system Torquato (2021).

To further illustrate the efficiency of diffusive transport in the DHU
cellular structures, we compare in Fig. 3(b)
$\mathcal{S}(\infty)-\mathcal{S}(t)$ of three model structures derived from
the Poisson, RSA and vein-network point configurations by placing congruent
circular disks at the points, see Fig. 3(c). Again, the size of the disks is
chosen such that their area is equivalent to the average area of the open-end
branched structures, and the models assume the solutes are only carried and
released by these branched structures. For the DHU system, we have
$\tilde{\chi}_{V}({k})=\rho\tilde{m}^{2}(k)S(k)$ for the disk regions, where
$\tilde{m}(k)$ is the Fourier transform of the shape function of the disks,
and $\tilde{m}(0)$ is a constant. Thus, the small-$k$ scaling
$\tilde{\chi}_{V}({k})\sim S(k)\sim k^{\alpha}$ leads to
$\mathcal{S}(\infty)-\mathcal{S}(t)\sim 1/t^{(1+\alpha/2)}$. The scaling
analysis from the numerical data yields $\sim 1/t^{1.33}$, which is consistent
with the calculation of the exponent from $\alpha\approx 0.64$. This indicates
the uniform solute distribution can be achieved much faster in the DHU system
than in other two disordered system, with the latter two characterized by a
scaling $\mathcal{S}(\infty)-\mathcal{S}(t)\sim 1/t$ (illustrated as dashed
lines in Fig. 3(b)). These analyses clearly indicate the superior efficiency
of diffusive transport behaviors associated with the DHU vein networks,
compared to other disordered organizations. In this sense, the DHU pattern can
be considered as a type of optimal organization of secondary vein networks for
transport properties under constraints induced by the primary vein branches.

$\begin{array}[]{c}\\\
\includegraphics[width=205.97214pt]{fig4.png}\end{array}$

Figure 4: Statistics of the point configurations with a rescaled density
$\rho=125$ cm-2 associated with distinct leaf species including ficus
religiosa (FR), ficus caulocarpa (FC), ficus microcarpa (FM), smilax indica
(SI), populus rotundifolia (PR), and yulania denudate (YD). The images of leaf
samples except for FR are shown in (a). The scale bars for the samples are:
7.8 cm for FC; 2.6 cm for FM; 1.6 cm for SI; 3.9 cm for YD; and 2.1 com for
PR. The statistics including the static structure factor $S(k)$ (a) and (b),
the number variance $\sigma_{N}^{2}(R)$ (c), as well as the pair-correlation
function $g_{2}(r)$ (d) indicate universal hyperuniform organizations in the
secondary vein networks of these different leaf species.

Universal hyperuniformity in vein networks across different leaf species. Last
but not least, we present the analyses on the leaf samples of other species
including ficus caulocarpa (FC), ficus microcarpa (FM), smilax indica (SI),
populus rotundifolia (PR), and yulania denudate (YD), as shown in Fig. 4(a).
The statistics including $S(k)$, $\sigma_{N}^{2}(R)$ and $g_{2}(r)$ for all
six leaf species studied here, when rescaled to the same number density, all
virtually collapse on to the corresponding universal curves, as shown in Fig.
4(b)-(e). Although the small-$k$ behaviors of $S(k)$ for the six systems
exhibit some small discrepancy, possibly due to the finite system size
effects, the estimated hyperuniformity exponent $\alpha$ for all systems agree
reasonably well with one another, with the numerical values falling into the
narrow range between [0.626, 0.648]. These results indicate the secondary vein
networks for all six leaf species studied here exhibit a universal
hyperuniform organization on large scales, characterized by $S(k)\sim
k^{\alpha}$ with $\alpha\approx 0.64$ which belongs to Class III
hyperuniformity. In addition, the different leaf species also possess very
similar local structures when properly rescaled, manifested as the collapsed
pair-correlation functions.

How such universal DHU organizations emerge during the development and
regulation of the vein network? It is well established that the leaf vascular
system develops in response to the flow and localization of the growth
regulator auxin as continuous strands of conducting tissues arranged in
repetitive spatial patterns, which eventually forms the hierarchical vein
networks Scarpella et al. (2010); Biedroń and Banasiak (2018). Recently,
successful substrate-depletion type reaction-diffusion models have been
developed, which have qualitatively reproduced the growth and branching
patterns of the leaf vascular system, including the closed circuits of the
veins Fujita and Mochizuki (2006a, b). Interestingly, it has been shown that a
subset of Turing patterns resulted from such reaction-diffusion equations can
evolve toward the unusual long-range property of hyperuniformity Ge (2023),
shedding lights on the emergence of hyperuniformity in the development of
these vein networks, which we will explore in future work.

In summary, we have discovered a remarkable universal hyperuniform
organization of the cellular structures of secondary leaf vein networks across
a wide spectrum of tree species. We have shown that such DHU organization
leads to superior efficiency of diffusive transport in the system, which is
associated with a crucial function of the vein networks for carrying and
transporting water and minerals to the leaf cells. In this sense, the DHU
pattern represents an optimal organization of secondary vein networks for
transport under constraints induced by the primary vein branches. These
results could inspire the design of disordered network materials, as such the
Voronoi networks derived from the DHU patterns investigated here with
desirable transport properties.

## References

  * Jiao et al. (2014) Y. Jiao, T. Lau, H. Hatzikirou, M. Meyer-Hermann, J. C. Corbo, and S. Torquato, Phys. Rev. E 89, 022721 (2014).
  * Mayer et al. (2015) A. Mayer, V. Balasubramanian, T. Mora, and A. M. Walczak, Proc. Natl. Acad. Sci. USA 112, 5950 (2015).
  * Ge (2023) Z. Ge, Proc. Natl. Acad. Sci. USA 120, e2306514120 (2023).
  * Torquato and Stillinger (2003) S. Torquato and F. H. Stillinger, Phys. Rev. E 68, 041113 (2003).
  * Torquato (2018) S. Torquato, Phys. Rep. 745, 1 (2018).
  * Torquato et al. (2015) S. Torquato, G. Zhang, and F. H. Stillinger, Phys. Rev. X 5, 021020 (2015).
  * Batten et al. (2009) R. D. Batten, F. H. Stillinger, and S. Torquato, Phys. Rev. Lett. 103, 050602 (2009).
  * Jiao (2022) Y. Jiao, Physica A: Statistical Mechanics and its Applications 585, 126435 (2022).
  * Gabrielli et al. (2002) A. Gabrielli, M. Joyce, and F. S. Labini, Phys. Rev. D 65, 083523 (2002).
  * Donev et al. (2005) A. Donev, F. H. Stillinger, and S. Torquato, Phys. Rev. Lett. 95, 090604 (2005).
  * Zachary et al. (2011a) C. E. Zachary, Y. Jiao, and S. Torquato, Phys. Rev. Lett. 106, 178001 (2011a).
  * Zachary et al. (2011b) C. E. Zachary, Y. Jiao, and S. Torquato, Physical Review E 83, 051309 (2011b).
  * Yuan et al. (2021) Y. Yuan, Y. Jiao, Y. Wang, and S. Li, Physical Review Research 3, 033084 (2021).
  * Kurita and Weeks (2011) R. Kurita and E. R. Weeks, Phys. Rev. E 84, 030401 (2011).
  * Hunter and Weeks (2012) G. L. Hunter and E. R. Weeks, Rep. Prog. Phys. 75, 066501 (2012).
  * Dreyfus et al. (2015) R. Dreyfus, Y. Xu, T. Still, L. A. Hough, A. G. Yodh, and S. Torquato, Phys. Rev. E 91, 012302 (2015).
  * Chen et al. (2023a) D. Chen, X. Jiang, D. Wang, H. Zhuang, and Y. Jiao, Acta Materialia 246, 118678 (2023a).
  * Zhang et al. (2023) H. Zhang, X. Wang, J. Zhang, H.-B. Yu, and J. F. Douglas, arXiv preprint arXiv:2302.01429 (2023).
  * Chremos and Douglas (2018) A. Chremos and J. F. Douglas, Phys. Rev. Lett. 121, 258002 (2018).
  * Torquato et al. (2008) S. Torquato, A. Scardicchio, and C. E. Zachary, J. Stat. Mech.: Theory Exp. p. P11019 (2008).
  * Feynman and Cohen (1956) R. P. Feynman and M. Cohen, Phys. Rev. 102, 1189 (1956).
  * Gerasimenko et al. (2019) Y. A. Gerasimenko, I. Vaskivskyi, M. Litskevich, J. Ravnik, J. Vodeb, M. Diego, V. Kabanov, and D. Mihailovic, Nat. Mater. 18, 1078 (2019).
  * Sakai et al. (2022) S. Sakai, R. Arita, and T. Ohtsuki, arXiv preprint arXiv:2207.09698 (2022).
  * Rumi et al. (2019) G. Rumi, J. A. Sánchez, F. Elías, R. C. Maldonado, J. Puig, N. R. C. Bolecek, G. Nieva, M. Konczykowski, Y. Fasano, and A. B. Kolton, Phys. Rev. Res. 1, 033057 (2019).
  * Sánchez et al. (2019) J. A. Sánchez, R. C. Maldonado, N. R. C. Bolecek, G. Rumi, P. Pedrazzini, M. I. Dolz, G. Nieva, C. J. van der Beek, M. Konczykowski, C. D. Dewhurst, et al., Commun. Phys. 2, 1 (2019).
  * Zheng et al. (2020) Y. Zheng, L. Liu, H. Nan, Z.-X. Shen, G. Zhang, D. Chen, L. He, W. Xu, M. Chen, Y. Jiao, et al., Sci. Adv. 6, eaba0826 (2020).
  * Chen et al. (2021) D. Chen, Y. Zheng, L. Liu, G. Zhang, M. Chen, Y. Jiao, and H. Zhuang, Proc. Natl. Acad. Sci. U.S.A. 118, e2016862118 (2021).
  * Chen et al. (2022) D. Chen, Y. Liu, H. Zhuang, M. Chen, and Y. Jiao, Physical Review B 106, 235427 (2022).
  * Chen et al. (2023b) D. Chen, H. Zhuang, M. Chen, P. Y. Huang, V. Vlcek, and Y. Jiao, Applied Physics Reviews 10 (2023b).
  * Hexner and Levine (2015) D. Hexner and D. Levine, Phys. Rev. Lett. 114, 110602 (2015).
  * Jack et al. (2015) R. L. Jack, I. R. Thompson, and P. Sollich, Phys. Rev. Lett. 114, 060601 (2015).
  * Weijs et al. (2015) J. H. Weijs, R. Jeanneret, R. Dreyfus, and D. Bartolo, Phys. Rev. Lett. 115, 108301 (2015).
  * Salvalaglio et al. (2020) M. Salvalaglio, M. Bouabdellaoui, M. Bollani, A. Benali, L. Favre, J.-B. Claude, J. Wenger, P. de Anna, F. Intonti, A. Voigt, et al., Physical Review Letters 125, 126101 (2020).
  * Nizam et al. (2021) Ü. S. Nizam, G. Makey, M. Barbier, S. S. Kahraman, E. Demir, E. E. Shafigh, S. Galioglu, D. Vahabli, S. Hüsnügil, M. H. Güneş, et al., Journal of Physics: Condensed Matter 33, 304002 (2021).
  * Zheng et al. (2023) Y. Zheng, M. A. Klatt, and H. Löwen, arXiv preprint arXiv:2310.03107 (2023).
  * Lei et al. (2019) Q.-L. Lei, M. P. Ciamarra, and R. Ni, Sci. Adv. 5, eaau7423 (2019).
  * Lei and Ni (2019) Q.-L. Lei and R. Ni, Proceedings of the National Academy of Sciences 116, 22983 (2019).
  * Huang et al. (2021) M. Huang, W. Hu, S. Yang, Q.-X. Liu, and H. Zhang, Proceedings of the National Academy of Sciences 118, e2100493118 (2021).
  * Zhang and Snezhko (2022) B. Zhang and A. Snezhko, Physical Review Letters 128, 218002 (2022).
  * Oppenheimer et al. (2022) N. Oppenheimer, D. B. Stein, M. Y. B. Zion, and M. J. Shelley, Nature communications 13, 804 (2022).
  * Hexner and Levine (2017) D. Hexner and D. Levine, Physical review letters 118, 020601 (2017).
  * Hexner et al. (2017) D. Hexner, P. M. Chaikin, and D. Levine, Proceedings of the National Academy of Sciences 114, 4294 (2017).
  * Weijs and Bartolo (2017) J. H. Weijs and D. Bartolo, Physical review letters 119, 048002 (2017).
  * Wilken et al. (2022) S. Wilken, A. Z. Guo, D. Levine, and P. M. Chaikin, arXiv preprint arXiv:2212.09913 (2022).
  * Florescu et al. (2009) M. Florescu, S. Torquato, and P. J. Steinhardt, Proc. Natl. Acad. Sci. U.S.A. 106, 20658 (2009).
  * Klatt et al. (2022) M. A. Klatt, P. J. Steinhardt, and S. Torquato, Proceedings of the National Academy of Sciences 119, e2213633119 (2022).
  * Zhang et al. (2016) G. Zhang, F. H. Stillinger, and S. Torquato, J. Chem. Phys. 145, 244109 (2016).
  * Chen and Torquato (2018) D. Chen and S. Torquato, Acta Mater. 142, 152 (2018).
  * Torquato (2021) S. Torquato, Physical Review E 104, 054102 (2021).
  * Xu et al. (2017) Y. Xu, S. Chen, P. Chen, W. Xu, and Y. Jiao, Phys. Rev. E 96, 043301 (2017).
  * Leseur et al. (2016) O. Leseur, R. Pierrat, and R. Carminati, Optica 3, 763 (2016).
  * Yu (2023) S. Yu, Nature Computational Science 3, 128 (2023).
  * Shi et al. (2023) W. Shi, D. Keeney, D. Chen, Y. Jiao, and S. Torquato, Physical Review E 108, 045306 (2023).
  * Torquato and Chen (2018) S. Torquato and D. Chen, Multifunct. Mater. 1, 015001 (2018).
  * Scarpella et al. (2010) E. Scarpella, M. Barkoulas, and M. Tsiantis, Cold Spring Harbor perspectives in biology 2, a001511 (2010).
  * Biedroń and Banasiak (2018) M. Biedroń and A. Banasiak, Plant Cell Reports 37, 1215 (2018).
  * Katifori and Magnasco (2012) E. Katifori and M. O. Magnasco, PloS one 7, e37994 (2012).
  * Fujita and Mochizuki (2006a) H. Fujita and A. Mochizuki, Developmental Dynamics 235, 2710 (2006a).
  * Fujita and Mochizuki (2006b) H. Fujita and A. Mochizuki, J Theoret Biol 241, 541 (2006b).