††thanks: These authors contributed equally††thanks: These authors contributed
equally

# Optimally scrambling chiral spin-chain with effective black hole geometry

Aiden Daniel  School of Physics and Astronomy, University of Leeds, Leeds LS2
9JT, United Kingdom  Andrew Hallam  School of Physics and Astronomy,
University of Leeds, Leeds LS2 9JT, United Kingdom  Matthew D. Horner  School
of Physics and Astronomy, University of Leeds, Leeds LS2 9JT, United Kingdom
Aegiq Ltd., Cooper Buildings, Arundel Street, Sheffield, S1 2NS, United
Kingdom  Jiannis K. Pachos  School of Physics and Astronomy, University of
Leeds, Leeds LS2 9JT, United Kingdom

###### Abstract

There is currently significant interest in emulating the essential
characteristics of black holes, such as their Hawking radiation or their
optimal scrambling behavior, using condensed matter models. In this article,
we investigate a chiral spin-chain, whose mean field theory effectively
captures the behavior of Dirac fermions in the curved spacetime geometry of a
black hole. We find that within the region of the chain that describe the
interior of the black hole, strong correlations prevail giving rise to many-
body chaotic dynamics. Employing out-of-time-order correlations as a
diagnostic tool, we numerically compute the associated Lyapunov exponent.
Intriguingly, we observe a linear increase in the Lyapunov exponent with
temperature within the black hole’s interior at low temperatures, indicative
of optimal scrambling behavior. This contrasts with the quadratic temperature
dependence exhibited by the spin-chain on the region outside the black hole.
Our findings contribute to a deeper understanding of the interplay between
black hole geometry and quantum chaos, offering insights into fundamental
aspects of quantum gravity.

Introduction:– Black holes pose direct challenges to our understanding of
fundamental laws of nature. Central to these open questions is the black hole
information paradox, first articulated by Stephen Hawking in the 1970s [1, 2].
According to general relativity, the gravitational pull of black holes is so
intense that it creates a region known as the event horizon, beyond which
information appears to be irretrievably lost. However, the unitarity of
quantum mechanics suggests information cannot be destroyed, leading to the
question of what happens to the information of an object that falls into a
black hole. This apparent contradiction between quantum physics and general
relativity has given rise to intense theoretical investigations and remains
unresolved to this day.

Recent advancements in theoretical physics have provided new insights into
black hole dynamics, particularly through the investigation of quantum
information scrambling. Information scrambling refers to the rapid and
thorough mixing of information within a quantum system. It is thought that
black holes exhibit optimal scrambling behavior, leading to the rapid
thermalization of newly engulfed quantum information. Currently, toy models
that exhibit maximal scrambling, such as the $(0+1)$D Sachdev-Ye-Kitaev (SYK)
model [3, 4, 5, 6, 7, 8], are related to $(1+1)$D black holes only through the
AdS1+1/CFT0+1 duality [9], with direct black hole models that exhibit this
behavior still lacking.

In this study, we explore the quantum properties of $(1+1)$D black holes using
a recently introduced chiral spin-chain model [10, 11]. The mean field theory
limit of this model effectively describes Dirac fermions in a black hole
background geometry, which is similar to the semiclassical limit of quantum
gravity [12]. However, within the region of the chain representing the black
hole’s interior, the mean field theory description breaks down due to the
dominance of strong correlations. The natural question arises if this strongly
correlated region gives rise to optimal scrambling behaviour.

To probe the scrambling behaviour of the chiral spin-chain, we numerically
investigate its Out-of-Time-Order correlators (OTOCs). OTOCs are a special
class of quantum correlation functions that determine the Lyapunov exponent,
capable of diagnosing early-time chaotic behaviour [13, 14]. Our model
exhibits similar scrambling behavior as the SYK model, so we employ the
approach presented in [15] for the numerical analysis of the SYK model. By
focusing on the system’s Lyapunov exponent at low temperatures [9, 16, 17, 6]
we observe that it behaves linearly with temperature at the region of the
chain that describes the inside of the black hole, where the chiral
interactions are dominant. This signature of optimal scrambling contrasts with
the quadratic behaviour observed in the outside region of the black hole,
where the chiral interactions have a perturbative effect on top of an XY
coupling [18]. The functional dependence of the Lyapunov exponent on
temperature is analysed for various coupling regimes and system sizes showing
a robust behaviour and fast convergence to the expected thermodynamic values.
Therefore, our chiral model reveals an intricate interplay between black hole
geometry and quantum chaos behaviour as expected from a comprehensive quantum
gravity description.

Figure 1:  (a) The chiral spin-chain (1) with position dependent chiral
coupling $v$, while $u$ is constant. The chiral interaction, $\chi_{i}$, acts
on three successive spins and has alternating orientation. The fermionic sites
$a$ and $b$ represent the unit cell of the mean field theory (MFT), where the
spin of the Dirac field is encoded. The MFT gives the black hole (BH) geometry
with $v$ small on the left of the chain ($\frac{v}{2}<u$, outside BH) or large
on the right ($\frac{v}{2}>u$, inside BH), with the horizon positioned at
$\frac{v}{2}=u$. (b) The dispersion relation $E(p)$ obtained from the MFT
description of homogeneous periodic chains. In the low energy limit it
faithfully reproduces the behaviour of Dirac fermions in black hole geometry
given by (3). (c) The lightcones, reciprocal to the dispersion relation that
describe the effective metric (4) of the system for various values of $v$. (d)
The Lyapunov exponent of the fully interacting model exhibits quadratic
behaviour $\lambda\propto T^{2}$ (orange line) on the outside of the BH, where
chiral interactions act perturbatively. Maximal scrambling is observed on the
inside of the BH with $\lambda\propto T$. The orange line corresponds to the
predicted optimal scrambling behaviour with slope $2\pi v/2$. The constant
offsets tend to zero with system size. Plots are for $N=8$ and $u=1$.

The chiral spin-chain:– The model, first introduced in Refs. [10, 11], is
shown in Fig. 1(a). It describes an interacting chain of spin-$1/2$ particles
with Hamiltonian

$H=\frac{1}{2}\sum_{i=1}^{N}\left[-{u}\left(S^{x}_{i}S^{x}_{i+1}+S^{y}_{i}S^{y}_{i+1}\right)+\frac{v}{2}\boldsymbol{S}_{i}\cdot\boldsymbol{S}_{i+1}\times\boldsymbol{S}_{i+2}\right],$
(1)

where $u,v$ are taken here to be positive numbers and
$\boldsymbol{S}_{i}=({\sigma}^{x}_{i}/2,{\sigma}^{y}_{i}/2,{\sigma}^{z}_{i}/2)$
with $\sigma_{i}^{\alpha}$ ($\alpha=x,y,z$) is the $\alpha$-Pauli matrix of
the $i$th spin. Notably, this is the XY model with an additional three-spin
chirality term

$\chi_{i}=\boldsymbol{S}_{i}\cdot\boldsymbol{S}_{i+1}\times\boldsymbol{S}_{i+2},$
(2)

that introduces interactions [19, 20, 10, 11]. The enumeration of the sites in
the chiral interaction term causes the chirality to alternate along the chain,
as shown in Fig. 1(a). Unless otherwise stated we adopt open boundary
conditions.

As $v$ is increased the model undergoes a quantum phase transition from a
gapless XY-phase, to a gapless chirally-ordered phase, where the total
chirality $\langle\chi\rangle=\sum_{i}\langle\chi_{i}\rangle$ acts as an order
parameter. The critical point is located at $\frac{v}{2}\approx 1.12u$ [10,
11]. For $\frac{v}{2}\lessapprox 1.12u$, the ground state is in a free XY
phase with $\langle\chi\rangle=0$. Using standard bosonisation techniques [21,
22, 23] we found that the interactions were irrelevant and the low energy
physics is described by free fermions with renormalised Fermi velocities. On
the other hand, for $\frac{v}{2}\gtrapprox 1.12u$ the chiral interaction
dominates the XY term and the model transitions to a chiral phase with
$\langle\chi\rangle\neq 0$. In this phase bosonisation is more complicated due
to the system possessing two Fermi points. We find that in this phase, the
model does not remain at the free fermion point, revealing the dominance of
the interactions [11].

Black hole background geometry:– Intriguingly, the model has a geometric
interpretation in terms of black hole physics [10]. We first apply the Jordan-
Wigner transformation that maps the spins into fermions. The system then has a
unit cell of two fermion sites, $a$ and $b$, sitting at opposite ranks of the
triangular ladder shown in Fig. 1. By employing self-consistent mean field
theory (MFT), one can map the interacting spin model to a model of free
fermions on a lattice. We can investigate the behaviour of the model by taking
homogeneous coupling $v$ and adopting periodic boundary conditions to extract
the dispersion relation. This description faithfully captures the phase
diagram of the model, albeit with critical point at $\frac{v}{2}=u$. The
dispersion relation of the model, shown in Fig. 1(b), at low energies, i.e.,
in the continuum limit, can be faithfully reproduced by the Dirac action on a
fixed curved spacetime background

$S_{\rm{MFT}}=\int\mathrm{d}^{1+1}x|e|\bar{\psi}(x)\left(ie_{a}^{\
\mu}\gamma^{a}\overset{\leftrightarrow}{\partial_{\mu}}\right)\psi(x),$ (3)

where $a=0,1$ are local inertial frame indices, $\mu=t,x$ are coordinate
indices; the spinor field is given in terms of the unit cell fermions as
$\psi(x)=(a(x),b(x))^{T}/\sqrt{|e|}$ as shown in Fig. 1(a);
$A\overset{\leftrightarrow}{\partial_{\mu}}B=\frac{1}{2}\left(A\partial_{\mu}B-(\partial_{\mu}A)B\right)$;
$\gamma^{a}=(\sigma^{z},-i\sigma^{x})$; and $|e|=\text{det}(e^{a}_{\ \mu})$.
The zweibein $e_{a}^{\ \mu}$ are related to the spacetime metric by
$g_{\mu\nu}=e^{a}_{\ \mu}e^{b}_{\ \nu}\eta_{ab}$, where
$\eta_{ab}=\text{diag}(1,-1)$ is the Minkowski metric and $g_{\mu\nu}$ is
given by [10]

$\mathrm{d}s^{2}=\left(1-\frac{v^{2}}{4u^{2}}\right)\mathrm{d}t^{2}-\frac{4v}{u^{2}}\mathrm{d}t\mathrm{d}x-\frac{16}{u^{2}}\mathrm{d}x^{2}.$
(4)

This is the Schwarzschild metric of a black hole expressed in Gullstrand-
Painlevé coordinates [24, 25, 26, 27] (see SM for more details), which has
also been observed in other synthetic black hole systems [28, 29, 30, 31].

We now take the coupling $v$ to be a function of position, $v(x)$, varying
monotonically from small to large values. If it is slowly-varying compared to
the lattice spacing, then the continuum description in terms of the Dirac
equation remains valid. In this case the event horizon is located at
$\frac{v}{2}=u$, where $\frac{v}{2}<u$ corresponds to the outside of the black
hole and $\frac{v}{2}>u$ corresponds to the inside, as shown in Fig. 1(c). In
Gullstrand-Painlevé coordinates the light cone tilts when approaching the
black hole, having eventually both light directions pointing towards its
centre inside the black hole, as shown in Fig. 1(c) (Right). We see that the
event horizon aligns well with the boundary between the two phases of the
spin-chain, where the chiral phase is inside the black hole and the XY phase
outside. Using the mean field description it was shown in Ref. [10] that a
free particle that passes through the phase boundary of our model emerges as a
thermal radiation with the Hawking temperature, similar to other models [32,
33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45].

Apart from reproducing the semiclassical behaviour of a black hole, our model
also exhibits a chaotic behaviour, as we shall see in the following. This can
be quantified by the Lyapunov exponent shown in Fig. 1(d). For small $v$ we
obtain a Lyapunov exponent $\lambda\propto T^{2}$, as expected from
perturbative interactions [46]. For large $v$, i.e. inside the black hole, the
numerically obtained Lyapunov exponent exhibits linear behaviour
$\lambda\propto T$ [47], with a slope which is in excellent agreement with the
predicted optimal scrambling behaviour. The constant offsets in Fig. 1(d) for
$v=1$ and $v=10$ go to zero with system size, as we shall see in the
following.

Quantum chaos inside black hole:– A natural question to ask is if this black
hole geometric analogy of Eq. (1) extends to the thermalising dynamics
expected at the interior of a black hole. In particular, we investigate
whether the spin model is chaotic for $\frac{v}{2}>u$ (with homogeneous $v$
along the chain) and, more importantly, whether it exhibits maximal
information scrambling as expected of a black hole.

One of the most effective methods for diagnosing the chaotic behaviour of a
many-body quantum system is to study its energy level statistics, provided all
relevant symmetries have been resolved. We consider the chain with periodic
boundary conditions which has translational symmetry, U$(1)$ symmetry and
global spin flip symmetry $X=\prod_{i}\sigma^{x}_{i}$. We restrict to the
symmetry sector with quantum numbers $k=0,z=0,x=+1$ of these symmetries,
respectively, and determine the eigenvalues, $\left\\{E_{n}\right\\}$ of (1).
We then take the set $\left\\{s_{n}\right\\}$, where $s_{n}=E_{n}-E_{n-1}$,
and evaluate [48]

$r_{n}=\frac{\min\\{s_{n},s_{n-1}\\}}{\max\\{s_{n},s_{n-1}\\}}.$ (5)

The average of this value $\langle r\rangle$ and the probability distribution
over all $r_{n}$, are shown in Fig. 2. For $v\neq 0$, we find Wigner-Dyson
statistics indicating that this model is chaotic, a characteristic that
becomes more prominent with system size. Notably, we find $\langle
r\rangle\approx 0.53$ [49], which corresponds to the GOE ensemble. This is
perhaps unexpected since $H$ possesses complex matrix elements. Despite this,
the model retains time-reversal symmetry due to satisfying the relation
$PH^{T}P=H$ where $P$ is the parity operator, reminiscent of Ref. [50]. We
also calculated the Spectral Form Factor which further probes the spectrum for
evidence of chaotic behaviour but did not observable any notable additional
behavior [51, 52, 53].

Figure 2:  (a) The average $r$-value (5) of the eigenspectrum of Hamiltonian
(1) for different values of $v$ and $N$ ($u=1$). (b) The probability
distribution of all r-values in the reduced symmetry sector at $N=20,v=4,u=1$.
In both cases we see the model matches the Gaussian orthogonal ensemble (GOE)
indicating a non-integrable model. Results were computed in the $k=0,z=0,x=+1$
symmetry sector.

Lyapunov exponent and optimal scrambling:– The average level spacing $\langle
r\rangle$ is a crude measure of the chaotic behaviour of a quantum system, in
the thermodynamic limit we expect Wigner-Dyson statistics for all models
except fine-tuned integrable systems. Therefore, we need a more precise
measure of the chaotic behaviour of the system inside the black hole to
determine whether the scrambling of quantum information is optimal, as it is
the case for the SYK model [15].

As a diagnostic tool we will employ the Lyapunov exponent, $\lambda$, that
quantifies the rate of thermalisation of a chaotic system [54]. In the quantum
mechanical framework, $\lambda$ is calculated using the decay in out-of-time-
order correlators (OTOCs) with respect to some local operator $O_{i}$ on site
$i$, as a function of time, $t$. We primarily focus on the regularised OTOC

$C(t)=\langle
O_{i}(t)\rho^{1/4}O_{j}(0)\rho^{1/4}O_{i}(t)\rho^{1/4}O_{j}(0)\rangle,$ (6)

where $\rho=\exp(-\beta H)/\mathcal{Z}$, with the partition function
$\mathcal{Z}=\textup{Tr}\exp(-\beta H)$ and $\beta=1/T$ is the inverse
temperature. We also scale such that $C(0)=1$.

The regularised version of OTOCs is suitable for investigating small
temperature behaviours and exhibits fast convergence even for small system
sizes [15]. This should be contrasted to the unregularized correlator which is
subject to stronger finite-size corrections at low temperatures [55]. Using
the regularised version, $C(t)$, we are restricted to exact diagonalisation
techniques. Further restrictions are placed on system size $N\leq 13$ due to
the need of time evolution. Fortunately, due to the fast convergence of our
model with system size, we find this to be sufficient for our study. Unless
otherwise stated, we take $O_{i}=S_{N/2}^{x}$ and $O_{j}=S_{N/2-2}^{x}$ for
$N$ even and $O_{i}=S_{(N+1)/2}^{x}$, $O_{j}=S_{(N-3)/2}^{x}$, for $N$ odd.
While the choice of $O=S^{x}$ restricts the use of system’s symmetries, it is
reminiscent of the Majorana fermionic operators in the SYK model and
successfully unearths the desired optimality behavior. We show below that the
regularised correlators we choose, allow us to faithfully extract the chaotic
behaviour of the model even with moderate system sizes.

With maximally quantum information scrambling models, one expects an
exponential decay in the OTOCs defined in Eq. (6) with an associated Lyapunov
exponent, $\lambda$, indicating the rate of decay. The numerical recipe we
employ here is identical to that presented in Ref. [15] in the context of SYK.
We find that the same method is effective at identifying the scrambling
behaviour of our model. Mirroring this study of the SYK model, we fit the
numerical data of Eq. (6) to the semiclassical function of the OTOC at low
temperatures

$C(t)=U\Big{(}\frac{1}{2},1,Ne^{-\lambda t}\Big{)}\sqrt{N}e^{-\lambda t/2},$
(7)

where $U$ is the Kummer’s confluent hypergeometric function and $\lambda$ is
the fitted Lyapunov exponent. In general, we expect $\lambda$ to depend on the
coupling $v$ and the temperature $T$ of the model, while we keep $u=1$.

Figure 3: Out-of-time-ordered correlators, $C(t)$, and Lyapunov exponent,
$\lambda$, of the chiral spin-chain for various coupling regimes. (a) Coloured
dots show the numerically evaluated $C(t)$ given in Eq. (6) for various values
of $v$ (shown in the color bar), and $T=\infty$. Lines show the fit of Eq. (7)
with $\lambda=0.78v$, which improves for $v$ deep in the chiral phase. (b)
Using the same process and parameters as in (a), we compute the OTOCs, while
varying both $T$ and $v$ and extract $\lambda$ via fitting (7). Large values
of $\lambda$ are observed for large $v$ (chiral regime) and large temperatures
$T$ that probe the full spectrum of the Hamiltonian. Black dashed line
indicates the phase transition at ${v}/{2}\approx u$, where a clear change in
behaviour at low temperatures is witnessed. Both plots are computed with
$N=10$, $u=1$.

Fig. 3 (a) shows the behaviour of OTOCs for various values of coupling $v$
when $u=1$ and the system size is $N=10$. We verify that the OTOCs
exponentially decrease in time as seen by the fit (solid lines) in Fig. 3 (a)
with $\lambda=0.78v$. This exponential scrambling behaviour is present for
large $v$, that take the spin-chain in the chiral regime, i.e. inside the
black hole. These results remain largely unchanged as we vary the system size
demonstrating the fast convergence in the properties of the chiral spin-chain.

With the confirmed presence of exponential decrease in the OTOCs with time, we
have already substantiated the argument for the presence of exponential
scrambling characterized by the Lyapunov exponent. We now investigate the
change in the Lyapunov exponent with temperature where temperature defines the
average energy of the density matrix. Zero temperature corresponds to the
ground state while infinite temperature corresponds to a uniform superposition
of all eigenstates. Due to expected thermal behaviour of mid spectrum
eigenstates in chaotic models, one expects the Lyapunov exponent, $\lambda$,
to increase with $T$ to some maximal bound. The behaviour of $\lambda$ as a
function of the coupling $v$ and the temperature $T$ is given in Fig. 3 (b),
where the increase in $\lambda$ is observed with $v$ and $T$, as expected. We
also note a sudden change in the behaviour of $\lambda/v$ as $v$ crosses the
phase transition point $\frac{v}{2}\approx u$ at small temperatures, revealing
the corresponding dramatic change in the scrambling behaviour of the model.

Figure 4: The Lyapunov exponent of the chiral spin-chain simulating the black
hole. (a) Outside the black hole ($v=1$ and $u=1$) the Lyapunov exponent
exhibits at low temperatures a quadratic behavior. The observed even/odd
system size dependence is decreased when $N$ increases. (b) Inside the black
hole ($v=8$ and $u=1$) a linear behaviour in the Lyapunov exponent is
observed, indicating optimal scrambling. (c) The functional dependence of the
linear regime of $\lambda$ on the coupling $v$. When scaling both $\lambda$
and $T$ by $v$ we find that with increasing $v$, the scaled Lyapunov exponent
eventually flattens out ($N=8$). (Inset) Upon scaling both $\lambda$ and $T$
axis by $v$, the linear plots collapse on top of each other ($N=12$). (d) In
the linear regime we fit $\lambda=a(T-c)$ as shown in (b) and extract the
slope $a$ and the offset $c$. (d, Top) We plot $a$ scaled by $v/2$ and by
$2\pi$ against $1/N$ ($N$ from 6 to 13). We find the slope $\lambda/(v/2)$
tends towards $2\pi$ (dashed line) as $N$ increases. (d, Bottom) The constant
offset, $c$, is also extracted. The linear fit to the data (dashed line)
showing $c$ tending to 0 in the large $N$ limit.

It is predicted that the Lyapunov exponent, describing the scrambling
behaviour in quantum systems, satisfies the universal bound $\lambda/J\leq
2\pi T$, where $J$ is a characteristic coupling of the system [47]. The
quantum gravity description of black holes and their holographic duals, such
as the SYK model, are known to saturate this bound [56], when $\lambda$ is
normalised by appropriately chosen coupling $J$. This optimal scrambling
behaviour is analytically and numerically identified in the SYK model for low
temperatures [47]. We will now investigate whether the chiral spin-chain is
optimally scrambling or not. To achieve that we will quantitatively determine
the functional dependence of the Lyapunov exponent on temperature, $T$, for
various regimes of $v$ that correspond to the inside and outside of the black
hole. Again, we employ the same method seen in [15] when fitting with Eq. (7).

In the low temperature limit we expect that in the weakly-interacting regime,
described by $\frac{v}{2}<u$, the Lyapunov exponent will grow quadratically
with temperature. In Fig. 4(a) we observe that the Lyapunov exponent obtains
the quadratic behaviours as $T$ goes to zero, when $v$ is small, i.e. outside
the black hole. Note that odd system sizes have a zero offset, i.e.
$\lambda\to 0$ as $T\to 0$, while for even $N$ there is a non-zero offset.
This offset is a finite system size effect that tends to zero as $N$
increases.

If the strongly interacting regime, $\frac{v}{2}>u$, exhibits black hole
phenomena, then optimal scrambling is expected, witnessed by a linear growth
of the Lyapunov exponent, $\lambda/J=2\pi T$. In Fig. 4(b) we observe that,
similar to the SYK model [15], the numerically obtained Lyapunov exponent has
a liner dependence on temperature. Hence, we anticipate that the chiral spin-
chain exhibits optimal scrambling at the coupling regime that describes the
inside of a black hole. To quantify how accurately the quadratic and linear
behaviors are manifested in our system we module a fit in the data of the form
$\lambda=a(T^{b}-c)$ and study $b$ with system size. We find that $b$ takes
values 1 and 2 in the corresponding regimes to a very good accuracy (see SM
for more details). Furthermore, Fig. 4(c) shows that for large enough $v$, the
Lyapunov exponent remains more or less constant taking the value
$\lambda/v\approx 0.20$. This should be contrasted to other sub-optimal models
with a rate of chaos that is parametrically slower than the SYK model [57,
58].

We next investigate the slope of the linear behaviour exhibited by the
Lyapunov exponent when $v$ is large. In Fig. 4(b) we identify the linear
behavior $\lambda=a(T-c)$ for a range of temperatures between $T_{\rm{min}}$
and $T_{\rm{max}}$. Here, $T_{\rm{min}}$ depends on the discreteness of the
finite system and tends to zero as $N$ increases, while $T_{\rm{max}}$ depends
on the rest of the dynamics of the chiral model. The saturation of the
scrambling bound is achieved for a linear gradient $2\pi$ normalised by the
coupling of interactions, given in Hamiltonian (1) by $v/2$. In Fig. 4(d)
(Top) we see that the slope of $\lambda/(v/2)$ as a function of temperature,
$2a/v$, tends to $2\pi$ with increasing system size as expected from the
optimal scrambling behaviour. Moreover, we observe that the constant offset,
$c$, tends to zero with increasing system size, $N$, where the dashed line in
Fig. 4(d) (Bottom) is a linear fit to the data. Note that both $a$ and $c$
show a strong oscillatory behaviour as a function of $N$, indicating the
significance of the boundary effects for the system sizes we considered.
Hence, in the strong chiral regime of our simulator, i.e. inside the black
hole, we expect to have $\frac{\lambda}{v/2}=2\pi T$ in the thermodynamic
limit. Although, resolving the ambiguity associated with choosing the
appropriate energy scale of the model, $J=v/2$, needs a theoretical
investigation that goes beyond the scope of this article, we postulate that
this expression corresponds to optimal scrambling.

Our numerical investigation shows that as the coupling $v$ varies from small
to large values the spin-chain (1) undergoes a quantum phase transition. This
transition does not only changes its ground state properties from non-chiral
to chiral [10, 11], but its thermalisation properties change from weakly
scrambling to optimal scrambling at $\frac{v}{2}\approx u$, much in the same
way as in [18]. Notably, our system does not have random all-to-all
interactions as it is the case for the SYK model [15]. This locality and
uniformity facilitates the fast convergence of our numerical simulations with
system size to the expected behaviour. This is manifestly seen in Fig. 4,
where system size effects become less pronounced with increasing $N$.

Conclusions:– Our investigation establishes the connection between quantum
information scrambling and the behavior of chiral spin chains that encode
black hole spacetime geometry. Through numerical analysis of out-of-time-order
correlators, we have provided compelling evidence that at the coupling regime
representing the interior of the black hole information encoded in the chain
scrambles at an optimal rate.

Our findings open avenues for further exploration in several directions.
Firstly, a theoretical analysis to determine the Lyapunov exponent of our
chiral spin chain would complement our numerical findings, offering a deeper
understanding of its chaotic behavior. Additionally, investigating the quantum
phase transition at $\frac{v}{2}\approx u$, where both the ground state and
scrambling behavior undergo significant changes, presents an intriguing topic
for future research. Moreover, this initial investigation lends itself to
generalisations to $(2+1)$ or $(3+1)$ dimensional black holes following the
methodology presented in [59].

Importantly, unlike the SYK model, our model incorporates uniform local chiral
interactions, making it more experimentally accessible. Previous studies have
shown the feasibility of realizing chiral interactions in optical lattice
systems [60, 20]. This suggests the potential for experimental verification of
optimal scrambling behavior in laboratory settings. The plausible experimental
accessibility of our model opens new possibilities for studying quantum
gravity-inspired phenomena in controlled environments.

Acknowledgements:– We would like to thank Zlatko Papić and Cristian Voinea for
inspiring conversations. This work was supported by EPSRC with Grant No.
EP/R020612/1 and by the Leverhulme Trust Research Leadership Award
RL-2019-015.

## References

  * Hawking [1993] S. W. Hawking, Particle Creation by Black Holes, Advanced Series in Astrophysics and Cosmology 8, 85 (1993).
  * Wald [2001] R. M. Wald, The thermodynamics of black holes, Living Reviews in Relativity 4, 10.12942/lrr-2001-6 (2001).
  * Sachdev and Ye [1993] S. Sachdev and J. Ye, Gapless spin-fluid ground state in a random quantum Heisenberg magnet, Phys. Rev. Lett. 70, 3339 (1993).
  * Maldacena and Stanford [2016] J. Maldacena and D. Stanford, Remarks on the sachdev-ye-kitaev model, Phys. Rev. D 94, 106002 (2016).
  * Kitaev [2015] A. Kitaev, A simple model of quantum holography (2015), the Kavli Institute for Theoretical Physics.
  * Kitaev and Suh [2018] A. Kitaev and S. J. Suh, The soft mode in the Sachdev-Ye-Kitaev model and its gravity dual, Journal of High Energy Physics 2018, 10.1007/jhep05(2018)183 (2018).
  * Polchinski and Rosenhaus [2016] J. Polchinski and V. Rosenhaus, The spectrum in the Sachdev-Ye-Kitaev model, Journal of High Energy Physics 2016, 1–25 (2016).
  * Fu and Sachdev [2016] W. Fu and S. Sachdev, Numerical study of fermion and boson models with infinite-range random interactions, Phys. Rev. B 94, 035135 (2016).
  * Jensen [2016] K. Jensen, Chaos in ${\mathrm{AdS}}_{2}$ Holography, Phys. Rev. Lett. 117, 111601 (2016).
  * Horner  _et al._ [2023] M. D. Horner, A. Hallam, and J. K. Pachos, Chiral spin-chain interfaces exhibiting event-horizon physics, Phys. Rev. Lett. 130, 016701 (2023).
  * Forbes  _et al._ [2023] E. Forbes, M. D. Horner, A. Hallam, J. Barker, and J. K. Pachos, Exploring interacting chiral spin chains in terms of black hole physics, Phys. Rev. B 108, 245142 (2023).
  * Morsink and Mann [1991] S. M. Morsink and R. B. Mann, Black hole radiation of dirac particles in 1+1 dimensions, Classical and Quantum Gravity 8, 2257 (1991).
  * Shenker and Stanford [2014] S. H. Shenker and D. Stanford, Black holes and the butterfly effect, Journal of High Energy Physics 2014, 10.1007/jhep03(2014)067 (2014).
  * Hosur  _et al._ [2016] P. Hosur, X.-L. Qi, D. A. Roberts, and B. Yoshida, Chaos in quantum channels, Journal of High Energy Physics 2016, 10.1007/jhep02(2016)004 (2016).
  * Kobrin  _et al._ [2021] B. Kobrin, Z. Yang, G. D. Kahanamoku-Meyer, C. T. Olund, J. E. Moore, D. Stanford, and N. Y. Yao, Many-Body Chaos in the Sachdev-Ye-Kitaev Model, Physical Review Letters 126, 10.1103/physrevlett.126.030602 (2021).
  * Maldacena  _et al._ [2016a] J. Maldacena, D. Stanford, and Z. Yang, Conformal symmetry and its breaking in two-dimensional nearly anti-de Sitter space, Progress of Theoretical and Experimental Physics 2016, 12C104 (2016a), https://academic.oup.com/ptep/article-pdf/2016/12/12C104/9620935/ptw124.pdf .
  * Engelsöy  _et al._ [2016] J. Engelsöy, T. G. Mertens, and H. Verlinde, An investigation of AdS2 backreaction and holography, Journal of High Energy Physics 2016, 10.1007/jhep07(2016)139 (2016).
  * Banerjee and Altman [2017] S. Banerjee and E. Altman, Solvable model for a dynamical quantum phase transition from fast to slow scrambling, Physical Review B 95, 10.1103/physrevb.95.134302 (2017).
  * D’Cruz and Pachos [2005] C. D’Cruz and J. K. Pachos, Chiral phase from three-spin interactions in an optical lattice, Phys. Rev. A 72, 043608 (2005).
  * Tsomokos  _et al._ [2008] D. I. Tsomokos, J. J. García-Ripoll, N. R. Cooper, and J. K. Pachos, Chiral entanglement in triangular lattice models, Physical Review A 77, 10.1103/physreva.77.012106 (2008).
  * Giamarchi [2003] T. Giamarchi,  _Quantum Physics in One Dimension_ (Oxford University Press, 2003).
  * Miranda [2002] E. Miranda, Introduction to bosonization, Brazilian Journal of Physics 33, https://doi.org/10.1590/S0103-97332003000100002 (2002).
  * Aditya and Sen [2021] S. Aditya and D. Sen, Bosonization study of a generalized statistics model with four fermi points, Phys. Rev. B 103, 235162 (2021).
  * Volovik [2009] G. E. Volovik,  _The Universe in a Helium Droplet_ (Oxford University Press, 2009).
  * Volovik and Zhang [2017] G. Volovik and K. Zhang, Lifshitz Transitions, Type-II Dirac and Weyl Fermions, Event Horizon and All That, J. Low Temp. Phys. 189, 276–299 (2017).
  * Volovik and Huhtala [2002] G. Volovik and P. Huhtala, Fermionic microstates within the Painlevé-Gullstrand black hole, J. Exp. Theor. Phys. 94, 853–861 (2002).
  * Volovik [2016] G. Volovik, Black hole and Hawking radiation by type-II Weyl fermions, Jetp Lett. 104, 645–648 (2016).
  * Beule  _et al._ [2021] C. D. Beule, S. Groenendijk, T. Meng, and T. L. Schmidt, Artificial event horizons in Weyl semimetal heterostructures and their non-equilibrium signatures, SciPost Phys. 11, 095 (2021).
  * Morice  _et al._ [2021] C. Morice, A. G. Moghaddam, D. Chernyavsky, J. van Wezel, and J. van den Brink, Synthetic gravitational horizons in low-dimensional quantum matter, Phys. Rev. Res. 3, L022022 (2021).
  * Könye  _et al._ [2022] V. Könye, C. Morice, D. Chernyavsky, A. G. Moghaddam, J. van den Brink, and J. van Wezel, Horizon physics of quasi-one-dimensional tilted Weyl cones on a lattice, Phys. Rev. Res. 4, 033237 (2022).
  * Haller  _et al._ [2023] A. Haller, S. Hegde, C. Xu, C. D. Beule, T. L. Schmidt, and T. Meng, Black hole mirages: Electron lensing and Berry curvature effects in inhomogeneously tilted Weyl semimetals, SciPost Phys. 14, 119 (2023).
  * Parikh and Wilczek [2000] M. K. Parikh and F. Wilczek, Hawking radiation as tunneling, Phys. Rev. Lett. 85, 5042 (2000).
  * Yang  _et al._ [2020] R.-Q. Yang, H. Liu, S. Zhu, L. Luo, and R.-G. Cai, Simulating quantum field theory in curved spacetime with quantum many-body systems, Phys. Rev. Research 2, 023107 (2020).
  * Sabsovich  _et al._ [2022] D. Sabsovich, P. Wunderlich, V. Fleurov, D. I. Pikulin, R. Ilan, and T. Meng, Hawking fragmentation and Hawking attenuation in Weyl semimetals, Phys. Rev. Research 4, 013055 (2022).
  * Maertens  _et al._ [2024] D. Maertens, N. Bultinck, and K. Van Acoleyen, Hawking radiation on the lattice from floquet and local hamiltonian quench dynamics, Phys. Rev. B 109, 014309 (2024).
  * Mertens  _et al._ [2022] L. Mertens, A. G. Moghaddam, D. Chernyavsky, C. Morice, J. van den Brink, and J. van Wezel, Thermalization by a synthetic horizon, Phys. Rev. Res. 4, 043084 (2022).
  * Stone [2013] M. Stone, An analogue of Hawking radiation in the quantum Hall effect, Classical and Quantum Gravity 30, 085003 (2013).
  * Steinhauer [2016] J. Steinhauer, Observation of quantum Hawking radiation and its entanglement in an analogue black hole, Nature Physics 12, 959 (2016).
  * Kosior  _et al._ [2018] A. Kosior, M. Lewenstein, and A. Celi, Unruh effect for interacting particles with ultracold atoms, SciPost Phys. 5, 61 (2018).
  * Roldán-Molina  _et al._ [2017] A. Roldán-Molina, A. S. Nunez, and R. A. Duine, Magnonic black holes, Phys. Rev. Lett. 118, 061301 (2017).
  * Rodríguez-Laguna  _et al._ [2017] J. Rodríguez-Laguna, L. Tarruell, M. Lewenstein, and A. Celi, Synthetic Unruh effect in cold atoms, Phys. Rev. A 95, 013627 (2017).
  * Retzker  _et al._ [2008] A. Retzker, J. I. Cirac, M. B. Plenio, and B. Reznik, Methods for Detecting Acceleration Radiation in a Bose-Einstein Condensate, Phys. Rev. Lett. 101, 110402 (2008).
  * Guan  _et al._ [2017] S. Guan, Z.-M. Yu, Y. Liu, G.-B. Liu, L. Dong, Y. Lu, Y. Yao, and S. A. Yang, Artificial gravity field, astrophysical analogues, and topological phase transitions in strained topological semimetals, npj Quantum Materials 2, 23 (2017).
  * Liu  _et al._ [2020] H. Liu, J.-T. Sun, C. Song, H. Huang, F. Liu, and S. Meng, Fermionic analogue of high temperature hawking radiation in black phosphorus, Chinese Physics Letters 37, 067101 (2020).
  * Huang  _et al._ [2018] H. Huang, K.-H. Jin, and F. Liu, Black-hole horizon in the Dirac semimetal ${\mathrm{Zn}}_{2}{\mathrm{In}}_{2}{\mathrm{S}}_{5}$, Phys. Rev. B 98, 121110 (2018).
  * Kim  _et al._ [2020] J. Kim, X. Cao, and E. Altman, Scrambling versus relaxation in Fermi and non-Fermi liquids, Physical Review B 102, 10.1103/physrevb.102.085134 (2020).
  * Maldacena  _et al._ [2016b] J. Maldacena, S. H. Shenker, and D. Stanford, A bound on chaos, Journal of High Energy Physics 2016, 10.1007/jhep08(2016)106 (2016b).
  * Oganesyan and Huse [2007] V. Oganesyan and D. A. Huse, Localization of interacting fermions at high temperature, Phys. Rev. B 75, 155111 (2007).
  * Atas  _et al._ [2013] Y. Y. Atas, E. Bogomolny, O. Giraud, and G. Roux, Distribution of the ratio of consecutive level spacings in random matrix ensembles, Phys. Rev. Lett. 110, 084101 (2013).
  * Regnault and Nandkishore [2016] N. Regnault and R. Nandkishore, Floquet thermalization: Symmetries and random matrix ensembles, Phys. Rev. B 93, 104203 (2016).
  * Liu [2018] J. Liu, Spectral form factors and late time quantum chaos, Phys. Rev. D 98, 086026 (2018).
  * Gaikwad and Sinha [2019] A. Gaikwad and R. Sinha, Spectral form factor in non-gaussian random matrix theories, Phys. Rev. D 100, 026017 (2019).
  * Cotler  _et al._ [2017] J. S. Cotler, G. Gur-Ari, M. Hanada, J. Polchinski, P. Saad, S. H. Shenker, D. Stanford, A. Streicher, and M. Tezuka, Black holes and random matrices, Journal of High Energy Physics 2017, 10.1007/jhep05(2017)118 (2017).
  * Hallam  _et al._ [2019] A. Hallam, J. G. Morley, and A. G. Green, The Lyapunov spectra of quantum thermalisation, Nature Communications 10, 10.1038/s41467-019-10336-4 (2019).
  * Lantagne-Hurtubise  _et al._ [2020] E. Lantagne-Hurtubise, S. Plugge, O. Can, and M. Franz, Diagnosing quantum chaos in many-body systems using entanglement as a resource, Phys. Rev. Res. 2, 013254 (2020).
  * Shenker and Stanford [2015] S. H. Shenker and D. Stanford, Stringy effects in scrambling (2015), arXiv:1412.6087 [hep-th] .
  * Stanford [2016] D. Stanford, Many-body chaos at weak coupling, Journal of High Energy Physics 2016, 10.1007/jhep10(2016)009 (2016).
  * Chowdhury and Swingle [2017] D. Chowdhury and B. Swingle, Onset of many-body chaos in the $O(N)$ model, Phys. Rev. D 96, 065005 (2017).
  * Mohanta [2023] N. Mohanta, Chiral pair density wave as a precursor of the pseudogap in kagome superconductors, Phys. Rev. B 108, L220507 (2023).
  * D’Cruz and Pachos [2005] C. D’Cruz and J. K. Pachos, Chiral phase from three-spin interactions in an optical lattice, Physical Review A 72, 10.1103/physreva.72.043608 (2005).
  * Coleman [2015] P. Coleman,  _Introduction to Many-Body Physics_ (Cambridge University Press, 2015).
  * Nakahara [2003] M. Nakahara,  _Geometry, Topology and Physics_ , 2nd ed. (Taylor & Francis Group, Abingdon, England, 2003).
  * Bañados and Reyes [2016] M. Bañados and I. Reyes, A short review on noether’s theorems, gauge symmetries and boundary terms, International Journal of Modern Physics D 25, 1630021 (2016).
  * Dirac [2001] P. Dirac,  _Lectures on Quantum Mechanics_, Belfer Graduate School of Science, monograph series (Dover Publications, 2001).
  * Srednicki [2007] M. Srednicki,  _Quantum Field Theory_ (Cambridge University Press, 2007).
  * Barcelos-Neto  _et al._ [1987] J. Barcelos-Neto, A. Das, and W. Scherer, Canonical quantization of constrained systems, Acta Phys. Polon. B 18, 269 (1987).
  * Nielsen and Ninomiya [1981a] H. Nielsen and M. Ninomiya, A no-go theorem for regularizing chiral fermions, Physics Letters B 105, 219 (1981a).
  * Nielsen and Ninomiya [1981b] H. Nielsen and M. Ninomiya, Absence of neutrinos on a lattice: (II). Intuitive topological proof, Nuclear Physics B 193, 173 (1981b).
  * Fradkin [2013] E. Fradkin,  _Field theories of condensed matter physics_ (Cambridge University Press, 2013) pp. 106–108.
  * Golan and Stern [2018] O. Golan and A. Stern, Probing topological superconductors with emergent gravity, Physical Review B 98, 10.1103/physrevb.98.064503 (2018).

##  Appendix A Lattice representation of Dirac field in black hole background

###  A.1 Mean field approximation

The system we investigate in the main text is the one-dimensional
spin-$\frac{1}{2}$ chain with the Hamiltonian

$H=\frac{1}{2}\sum_{n}\left[-{u}\left(S^{x}_{n}S^{x}_{n+1}+S^{y}_{n}S^{y}_{n+1}\right)+\frac{v}{2}\chi_{n}\right],$
(8)

where the spin chirality operator is given by

$\chi_{n}=\boldsymbol{S}_{n}\cdot\left(\boldsymbol{S}_{n+1}\times\boldsymbol{S}_{n+2}\right),$
(9)

where
$\boldsymbol{S}_{n}=\frac{1}{2}(\sigma^{x}_{n},\sigma^{y}_{n},\sigma^{z}_{n})$
is the spin vector, where $\sigma_{n}^{a}$ is the $a$-Pauli matrix acting on
the $n$th lattice site for $a\in\\{x,y,z\\}$, and the $u,v\in\mathbb{R}$ are
couplings with dimensions of energy. In this supplemental material, we assume
we are working in the thermodynamic limit.

First we transform from spin operators to Pauli operators. In terms of the
Pauli operators, the Hamiltonian is given by

$\displaystyle H=\sum_{n}\Big{[}$
$\displaystyle-\frac{u}{8}\left(\sigma^{x}_{n}\sigma^{x}_{n+1}+\sigma^{y}_{n}\sigma^{y}_{n+1}\right)$
(10)
$\displaystyle+\frac{v}{32}\epsilon_{abc}\sigma^{a}_{n}\sigma^{b}_{n+1}\sigma^{c}_{n+2}\Big{]}$

where the repeated Latin indices $a,b,c$ are summed over in the chirality
term. We now introduce the ladder operators
$\sigma_{n}^{\pm}=(\sigma^{x}_{n}\pm i\sigma^{y}_{n})/2$ and the Jordan-Wigner
transformation defined by [61]

$\displaystyle\sigma^{+}_{n}$
$\displaystyle=\exp\left(-i\pi\sum_{m<n}c_{m}^{\dagger}c_{m}\right)c^{\dagger}_{n}$
(11) $\displaystyle\sigma^{-}_{n}$
$\displaystyle=\exp\left(i\pi\sum_{m<n}c_{m}^{\dagger}c_{m}\right)c_{n}$ (12)
$\displaystyle\sigma^{z}_{n}$ $\displaystyle=2c^{\dagger}_{n}c_{n}-1$ (13)

where $c_{n}$ are a set of fermionic modes obeying the anti-commutation
relations $\\{c_{n},c_{m}\\}=\\{c_{m}^{\dagger},c_{n}^{\dagger}\\}=0$ and
$\\{c_{n},c_{m}^{\dagger}\\}=\delta_{mn}$. After expressing the Hamiltonian in
terms of $\sigma_{n}^{\pm}$ and $\sigma^{z}_{n}$ and then applying the Jordan-
Wigner transformation, we arrive at

$\begin{split}H=\frac{1}{4}\sum_{n}\bigg{[}&-uc_{n}^{\dagger}c_{n+1}-\frac{iv}{4}c^{\dagger}_{n}c_{n+2}\\\
&+\frac{iv}{4}\left(c_{n}^{\dagger}c_{n+1}\sigma^{z}_{n+2}+c^{\dagger}_{n+1}c_{n+2}\sigma^{z}_{n}\right)\bigg{]}+\text{H.c.},\end{split}$
(14)

where for convenience we have left $\sigma^{z}_{n}$ alone under the assumption
that it represents the Jordan-Wigner transformation of Eq. (13). We see that
the model is intrinsically interacting as the fermionic Hamiltonian contains
quartic terms which arise from terms like
$c_{n}^{\dagger}c_{n+1}\sigma^{z}_{n+2}$ after explicitly substituting in Eq.
(13).

To analyse the behaviour of the interacting model, we apply mean field theory
(MFT) to transform the Hamiltonian into an effective quadratic Hamiltonian
which can be analytically diagonalised. MFT defines the fluctuation of an
operator $A$ as $\delta A=A-\langle A\rangle$, where $\langle A\rangle$ is the
expectation value of the operator $A$ with respect to the mean field ground
state $|\Omega\rangle$. For a product of two operators we have

$AB=\langle A\rangle B+A\langle B\rangle-\langle A\rangle\langle
B\rangle+\delta A\delta B,$ (15)

where the second order in fluctuations can be ignored. Applying this to the
interacting terms of Eq. (14) where we always consider $\sigma^{z}_{n}$ as one
of the operators in the product of Eq. (14), so replace
$\sigma^{z}_{n}\rightarrow\langle\sigma^{z}_{n}\rangle\equiv Z$ and
$c_{n}^{\dagger}c_{n+1}\rightarrow\langle
c_{n}^{\dagger}c_{n+1}\rangle\equiv\alpha$, where we have assumed
translational invariance. These expectation values are done with respect to
the ground state of the mean field Hamiltonian. The Hamiltonian becomes

$\begin{split}H_{\mathrm{MF}}(\alpha,Z)&=\frac{1}{4}\sum_{n}\bigg{[}-\left(u-\frac{iv}{2}Z\right)c^{\dagger}_{n}c_{n+1}\\\
&-\frac{iv}{4}c^{\dagger}_{n}c_{n+2}+4\mu
c^{\dagger}_{n}c_{n}\bigg{]}+E_{0}+\mathrm{H.c.},\end{split}$ (16)

where $\mu=v\mathrm{Im}(\alpha)/4$ is an effective chemical potential
controlling the number of particles in the ground state,
$E_{0}=v(Z-1)\mathrm{Im}(\alpha)/8$ is a constant energy shift.

Let $|\Omega(\alpha,Z)\rangle$ be the ground state of the Hamiltonian of Eq.
(16). Self consistency requires

$\displaystyle\langle\Omega(\alpha,Z)|\sigma^{z}_{n}|\Omega(\alpha,Z)\rangle$
$\displaystyle=Z$ (17)
$\displaystyle\langle\Omega(\alpha,Z)|c^{\dagger}_{n}c_{n+1}|\Omega(\alpha,Z)\rangle$
$\displaystyle=\alpha$ (18)

for all $n$. While these two equations have many solutions, we can single one
out on physical grounds: the fully interacting Hamiltonian of Eq. (14) has
particle-hole symmetry, $[H,U]=0$, where $U$ is the particle-hole
transformation with $Uc_{n}U^{\dagger}=(-1)^{n}c_{n}^{\dagger}$ and
$Uc_{n}^{\dagger}U^{\dagger}=(-1)^{n}c_{n}$. This symmetry implies that
$\langle c_{n}^{\dagger}c_{n}\rangle=1/2$ and $\langle
c_{n}^{\dagger}c_{n+1}\rangle\in\mathbb{R}$ in the ground state. If we require
the MFT to retain the particle-hole symmetry, then these conditions imply that
$Z=\mathrm{Im}(\alpha)=0$, and the MFT Hamiltonian becomes

$H_{\text{MF}}=\frac{1}{4}\sum_{n}\left(-uc^{\dagger}_{n}c_{n+1}-\frac{iv}{4}c^{\dagger}_{n}c_{n+2}\right)+\text{H.c.}.$
(19)

It was shown in Ref. [10] that this mean field limit faithfully describes the
second order phase transition exhibited by the full spin model.

###  A.2 The Dirac equation on curved spacetime

We now briefly introduce the Dirac field on a curved spacetime which we shall
use in the next section of the supplemental material. Suppose we have an
$N+1$-dimensional spacetime with metric $g_{\mu\nu}$ and a set of veilbein
$\\{e_{a}^{\ \mu}\\}$ and their inverses $\\{e_{a}^{\ \mu}\\}$ which are
related to the metric via

$g_{\mu\nu}=e^{a}_{\ \mu}e^{b}_{\ \nu}\eta_{ab},\quad g^{\mu\nu}=e_{a}^{\
\mu}e_{b}^{\ \nu}\eta^{ab},$ (20)

where $\eta_{ab}=\mathrm{diag}(1,-1,-1,\ldots,-1)$ is the Minkowski metric.
The veilbein and their inverses also obey

$e_{a}^{\ \mu}e^{b}_{\ \mu}=\delta^{a}_{b},\quad e_{a}^{\ \mu}e^{a}_{\
\nu}=\delta^{\mu}_{\nu}.$ (21)

The veilbein $\\{e_{a}^{\ \mu}\\}$ are a set of vector fields that form an
orthonormal basis at every point in some patch of $M$. The Latin indices
$a,b=0,1,\ldots$ refer to the orthonormal frame indices, whilst the Greek
indices $\mu,\nu=t,x,\ldots$ refer to the coordinate indices.

Using this, we can introduce spinor fields on a curved spacetime as a field on
$M$ which transforms as a spinor under local Lorentz transformations
(transformations that act on Latin indices) and as a scalar under coordinate
transformations (transformations that act on the Greek indices). The locally
Lorentz invariant and coordinate invariant action for spinor field $\psi$ of
mass $m$ on an $N+1$ dimensional spacetime $M$ with metric $g_{\mu\nu}$ is
given by [62]

$\displaystyle S$
$\displaystyle=\frac{i}{2}\int_{M}\mathrm{d}^{N+1}x|e|\left(\bar{\psi}\gamma^{\mu}D_{\mu}\psi-\overline{D_{\mu}\psi}\gamma^{\mu}\psi+2im\bar{\psi}\psi\right)$
(22) $\displaystyle\equiv\int_{M}\mathrm{d}^{N+1}x\mathcal{L},$

where the gamma matrices $\\{\gamma^{\mu}\equiv e_{a}^{\mu}\gamma^{a}\\}$ are
the curved space gamma matrices which obey the Clifford algebra
$\\{\gamma^{\mu},\gamma^{\nu}\\}=2g^{\mu\nu}$ and are related to the local
flat Minkowski space gamma matrices $\\{\gamma^{a}\\}$ which obey the flat
space Clifford algebra $\\{\gamma^{a},\gamma^{b}\\}=2\eta^{ab}$. The Dirac
adjoint is defined as $\bar{\psi}=\gamma^{\dagger}\gamma^{0}$ where
$\gamma^{0}$ is the flat space gamma matrix. We also have $|e|=\det e^{a}_{\
\mu}=\sqrt{-g}$. The covariant derivative of spinors $D_{\mu}$ is defined via

$D_{\mu}\psi=\partial_{\mu}\psi+\Omega_{\mu}\psi$ (23)

where $\Omega_{\mu}$ is the spin connection related to the connection of $M$
via

$\Omega_{\mu}=\frac{i}{2}\Omega_{ab\mu}\Sigma^{ab},\quad\Sigma^{ab}=\frac{i}{4}[\gamma^{a},\gamma^{b}]$
(24)

and $\Omega_{ab\mu}$ are the components of the connection. For more details,
see Ref. [62].

In this study we are interested in the $(1+1)$D spacetimes, in which case the
spin connection vanishes from the symmetrised action. To see this, we can
substitute in the covariant derivative explicitly

$\displaystyle\mathcal{L}$
$\displaystyle=\frac{i}{2}|e|\left(\bar{\psi}\gamma^{\mu}D_{\mu}\psi-\overline{D_{\mu}\psi}\gamma^{\mu}\psi\right)$
(25)
$\displaystyle=\frac{i}{2}|e|\left(\bar{\psi}\gamma^{\mu}\partial_{\mu}\psi-\partial_{\mu}\bar{\psi}\gamma^{\mu}\psi+\bar{\psi}\\{\gamma^{\mu},\Omega_{\mu}\\}\psi\right)$

However, in $(1+1)$D, we have
$\Omega_{\mu}\propto[\gamma^{0},\gamma^{1}]\propto\gamma^{3}$, where
$\gamma^{3}$ is the $(1+1)$D analogue of $\gamma^{5}$. As
$\\{\gamma^{\mu},\gamma^{3}\\}=0$ for all $\mu$, then the spin connection
vanishes and we arrive at the Lagrangian [62]

$\mathcal{L}=|e|\bar{\psi}\gamma^{\mu}\overset{\leftrightarrow}{\partial_{\mu}}\psi,$
(26)

where
$A\overset{\leftrightarrow}{\partial_{\mu}}B=\frac{1}{2}\left(A\partial_{\mu}B-(\partial_{\mu}A)B\right)$.

Throughout this we assume that we are working with time-independent metrics,
so that the vector $\xi=\partial_{t}$ is a time-like Killing vector which
obeys $\mathcal{L}_{\xi}g=0$. In order to perform canonical quantisation of
this theory, we introduce the canonical momentum $\pi_{a}(x)$ of the field and
apply the canonical commutation relations. We have

$\pi=\frac{\partial\mathcal{L}}{\partial\dot{\psi}}=\frac{i}{2}|e|\bar{\psi}\gamma^{t}$
(27)

The equal time canonical Poisson bracket reads

$\\{\psi_{\alpha}(t,x),\pi_{\beta}(t,y)\\}=i\delta(x-y)\delta_{\alpha\beta}$
(28)

where the indices $\alpha,\beta$ here refer to the spinor indices. This
implies that the spinor field obeys the commutation relation

$\\{\psi_{\alpha}(t,x),\psi_{\beta}^{\dagger}(t,y)\\}=\frac{(\gamma^{0}\gamma^{t})_{\alpha\beta}^{-1}\delta(x-y)}{|e|}.$
(29)

Note that the factor of $1/2$ is missing despite it being present in the
canonical momentum $\pi$. This is because the canonical momentum defines a
constraint on phase space, as $\pi$ is linearly related to $\psi^{\dagger}$,
which means we must employ the machinery of Dirac brackets instead of Poisson
brackets to quantise this theory. It is the Dirac bracket that we upgrade to
an anti-commutator using canonical quantisation. More information can be found
in Ref. [63, 64, 65, 66]. When taking this into account, the factor of $1/2$
vanishes.

###  A.3 Relativistic limit

Figure 5: To reveal the relativistic behaviour, we introduce a two-site unit
cell by bicolouring the lattice.

To make the link with relativity, the lattice sites are now labelled as
alternating between sub-lattices $A$ and $B$ by introducing a two-site unit
cell, as shown in Fig. 5. The mean field Hamiltonian of Eq. (19) can be re-
parameterised as

$\displaystyle H_{\mathrm{MF}}=\frac{1}{4}\sum_{n}\Big{[}$ $\displaystyle-
ua^{\dagger}_{n}(b_{n}+b_{n-1})$ (30)
$\displaystyle-\frac{iv}{4}(a_{n}^{\dagger}a_{n+1}+b^{\dagger}_{n}b_{n+1})\Big{]}+\text{H.c.},$

where the fermionic modes $a_{n}$ and $b_{n}$ belong to sublattice $A$ and
$B$, respectively, where now $n$ labels the unit cells. These modes obey the
commutation relations
$\\{a_{n},a^{\dagger}_{m}\\}=\\{b_{n},b_{m}^{\dagger}\\}=\delta_{nm}$, while
all mixed anti-commutators vanish. The index $n$ now labels the unit cells. We
Fourier transform the fermionic modes as

$a_{n}=\frac{1}{\sqrt{N}}\sum_{p\in\mathrm{B.Z.}}e^{ipan}a_{p},$ (31)

and similarly for $b_{n}$, where $N$ is the number of unit cells in the
system, $a$ is the unit cell spacing, and $\mathrm{B.Z.}=[0,2\pi/a)$ is the
Brillouin zone and $p=2m\pi/L$ for $L=Na$ and integer $m$. The Fourier
transformed Hamiltonian becomes

$H_{\mathrm{MF}}=\sum_{p\in\mathrm{B.Z.}}\chi^{\dagger}_{p}h(p)\chi_{p},\quad
h(p)=\begin{pmatrix}g(p)&f(p)\\\ f^{*}(p)&g(p)\end{pmatrix},$ (32)

where the two-component spinor is defined as
$\chi_{p}=(a_{p},b_{p})^{\mathrm{T}}$ and the functions are given by

$f(p)=-\frac{u}{4}(1+e^{-iap}),\quad g(p)=\frac{v}{8}\sin(ap).$ (33)

The dispersion relation is given by the eigenvalues of the single-particle
Hamiltonian $h(p)$ which yields

$\displaystyle E(p)$ $\displaystyle=g(p)\pm|f(p)|$ (34)
$\displaystyle=\frac{v}{8}\sin(ap)\pm\frac{u}{4}\sqrt{2+2\cos(ap)}.$

In Fig. 1 of the main text, it is seen that the parameter $v$ has the effect
of tilting the cones as it increases and one band becomes flat at the Fermi
point when $|v|/2=|u|$.

The Fermi points $\\{p_{i}\\}$, defined as the points for which $E(p_{i})=0$,
are found at

$p_{0}=\frac{\pi}{a},\quad
p_{\pm}=\pm\frac{1}{a}\arccos\left(1-\frac{8u^{2}}{v^{2}}\right).$ (35)

The roots $p_{\pm}$ only exist if the argument of $\arccos$ is in the range
$[-1,1]$ which implies $|v|/2\geq|u|$ for these to appear in the dispersion.
Therefore, if $|v|/2\leq|u|$, the only Fermi point is located at $p_{0}=\pi/a$
which is where the Dirac cone is located, as shown in Fig. 1 of the main text.
When the cone over-tilts, so when $|v|/2\geq|u|$, then the additional zero-
energy crossings at $p_{\pm}$ appear which is due to the Nielsen-Ninomiya
theorem which states that the number of left- and right-movers must be equal
[67, 68].

The continuum limit of a lattice model is an effective theory obtained by
letting the lattice spacing $a\rightarrow 0$ in such a way that the Fermi
velocity remains fixed. In this process, only the linear portion of the
dispersion relation, near the Fermi points, is relevant as the non-linear
portion of the dispersion goes off to infinite momentum. Therefore, the
continuum limit is equivalent to restricting ourselves to a small
neighbourhood of the Fermi points in momentum space. We outline this below.
See Ref. [69] for more detail.

First we Taylor expand the single-particle Hamiltonian $h(p)$ about the Fermi
point $p_{0}$ to first order in momentum which yields

$h(p_{0}+p)=\frac{1}{4}\left(u\sigma^{y}p-\frac{v}{2}\mathbb{I}p\right)\equiv
e_{a}^{\ i}\alpha^{a}p_{i},$ (36)

The coefficients in the second equality are defined as $e^{\ x}_{0}=-v/8,e^{\
x}_{1}=u/4$, where we have absorbed a factor of $a$ into the couplings as
$au\rightarrow u$ and $av\rightarrow v$, and the Dirac matrices are
$\alpha^{0}=\mathbb{I},\alpha^{1}=\sigma^{y}$. We then project the Hamiltonian
of Eq. (30) into a small region of momentum space centred on $p_{0}$ by
truncating the summation with a cutoff $\Lambda=O(1/a)$ as

$\displaystyle H$
$\displaystyle\approx\sum_{|q|<\Lambda}\chi^{\dagger}_{p_{0}+q}h(p_{0}+q)\chi_{p_{0}+q}$
(37) $\displaystyle=\sum_{|q|<\Lambda}\chi^{\dagger}(q)e_{a}^{\
i}\alpha^{a}p_{i}\chi(q)$

where we have defined the new momentum space fields
$\chi(q)\equiv\chi_{p_{0}+q}$, where $q$ measures the distance from the Fermi
point $p_{0}$.

We also truncate the discrete Fourier transform for the lattice fermions from
Eq. (31) as

$\displaystyle a_{n}$
$\displaystyle\approx\frac{1}{\sqrt{N}}\sum_{|q|<\Lambda}e^{i(p_{0}+q)an}a_{p_{0}+q}$
(38)
$\displaystyle=e^{ip_{0}an}\frac{1}{\sqrt{N}}\sum_{|q|<\Lambda}e^{iqan}a(q)$
$\displaystyle\equiv e^{ip_{0}an}a(n),$

which defines a slowly-varying field $a(n)$, and similarly for $b_{n}$ which
is related to $b(n)$ analogously. We see that on the subspace near the ground
state, the fermionic operators $a_{n}$ and $b_{n}$ consist of a slowly-varying
field $a(n)$ and $b(n)$ respectively, with a high-frequency oscillation
$e^{ip_{0}an}$ on top [69].

We then take the limit that $a\rightarrow 0$ in such a way that the rescaled
couplings $u$ and $v$ remain finite (equivalently the Fermi velocity remains
fixed) and $Na=L$ remains constant. We also define $na\rightarrow x$ which we
must remember is the unit cell coordinate. The cutoff
$\Lambda\rightarrow\infty$ additionally, so the summation is from $\pm\infty$.
Performing this limit, real space becomes a continuum and the envelope
functions $a(n)$ become

$a(x)=\lim_{a\rightarrow
0}\frac{a(n)}{\sqrt{a}}=\frac{1}{\sqrt{L}}\sum_{q\in\mathrm{B.Z.}}e^{iqx}a(q)$
(39)

and similarly for $b(x)$, where now the Brillouin zone has extended to
infinity as $\mathrm{B.Z.}=[-\infty,\infty]$ with $p=2m\pi/L$ for
$m\in\mathbb{Z}$, where the re-scaling by $1/\sqrt{a}$ ensures that the limits
exist and the commutation relations become continuum commutation relations. If
we define the two-component spinor field $\chi(x)=(a(x),b(x))^{T}$, we see
that this is related to the momentum space fields derived in Eq. (36) by a
Fourier transform as

$\chi(q)=\frac{1}{\sqrt{L}}\int_{0}^{L}\mathrm{d}xe^{-iqx}\chi(x).$ (40)

With this result in hand, we are now able to transform the truncated
Hamiltonian of Eq. (37) back into real space, arriving at the Hamiltonian

$\displaystyle H$
$\displaystyle=\int_{\mathbb{R}}\mathrm{d}x\chi^{\dagger}(x)\left(-ie_{a}^{\
i}\alpha^{a}\overset{\leftrightarrow}{\partial_{i}}\right)\chi(x)$ (41)
$\displaystyle\equiv\int_{\mathbb{R}}\mathrm{d}x\mathcal{H},$

with
$A\overset{\leftrightarrow}{\partial_{\mu}}B=\frac{1}{2}\left(A\partial_{\mu}B-(\partial_{\mu}A)B\right)$
and the Dirac $\alpha^{a}=(\mathbb{I},\sigma^{y})$ and $\beta=\sigma^{z}$. We
have ignored the overall factor of $1/4$ here.

This is a Hamiltonian for the slowly-varying envelope function $\chi(x)$. The
corresponding action of this theory is given by

$\displaystyle S$ $\displaystyle=\int
d^{1+1}x\left(i\chi^{\dagger}\overset{\leftrightarrow}{\partial_{t}}\chi-\mathcal{H}\right)$
(42) $\displaystyle=\int d^{1+1}xi\bar{\chi}e_{a}^{\ \mu}\gamma^{a}\chi$

where we have defined $\bar{\chi}=\chi^{\dagger}\gamma^{0}$ and the gammas are
related to the alpha and beta matrices via $\gamma^{0}=\beta$ and
$\gamma^{i}=\beta\alpha^{i}$, where $\gamma^{0}=\sigma^{z}$ and
$\gamma^{1}=-i\sigma^{x}$ which obey the anti-commutation relations
$\\{\gamma^{a},\gamma^{b}\\}=2\eta^{ab}$, with
$\eta^{ab}=\mathrm{diag}(1,-1)$. The coefficients are given by

$e_{a}^{\ \mu}=\begin{pmatrix}1&-v/8\\\ 0&u/4\end{pmatrix},\quad e^{a}_{\
\mu}=\begin{pmatrix}1&v/(2u)\\\ 0&4/u\end{pmatrix}$ (43)

If we assume that the couplings are upgraded to slowly-varying functions as
$u\rightarrow u(x)$ and $v\rightarrow v(x)$, then the continuum limit is still
a good approximation.

This action looks very similar to the action of a Dirac field on a
$(1+1)$-dimensional curved spacetime, except two subtle differences. The first
is that the integration measure is missing the factor of $|e|$, so it is the
flat space volume element. The second is that the fields obey flat space
commutation relations

$\\{\chi_{\alpha}(x),\chi_{\beta}(y)\\}=\delta_{\alpha\beta}\delta(x-y)$ (44)

which can be obtained using Eqs. (39) and (40) and the fact that the component
momentum space modes obey $\\{a(p),a^{\dagger}(q)\\}=\delta_{pq}$, and
similarly for $b(p)$. Therefore, the theory in its current form describes a
generalised Dirac action on a flat $(1+1)$-dimensional space with space-
dependent coefficients.

In order to interpret this theory as a curved space theory, we introduce a new
field

$\psi=\frac{\chi}{\sqrt{|e|}}$ (45)

then the fields obey the curved space commutation relations

$\\{\psi_{\alpha}(x),\psi_{\beta}(y)\\}=\frac{\delta_{\alpha\beta}\delta(x-y)}{|e|}$
(46)

agreeing precisely with the general commutation relations of Eq. (29) using
the veilbein of Eq. (43). The action of Eq. (22) re-expressed in terms of the
field $\psi$ is precisely the Dirac action for a spinor on a spacetime with
veilbein given in Eq. (43). This veilbein corresponds to the metric

$\mathrm{d}s^{2}=\left(1-\frac{v^{2}}{4u^{2}}\right)\mathrm{d}t^{2}-\frac{4v}{u^{2}}\mathrm{d}t\mathrm{d}x-\frac{16}{u^{2}}\mathrm{d}x^{2}.$
(47)

If the variables $u$ and $v$ are upgraded to slowly-varying functions of
space, then the preceding calculation is still valid and the event horizon is
located at the point $x_{h}$, where $|v(x_{h})|/2=|u(x_{h})|$. In the small
region in which $v$ is a slowly-varying functions of $x$, the coupling of
different momentum modes will be small and can be ignored to a good
approximation, leaving the diagonal terms $a^{\dagger}_{p}a_{p}$ only. This is
quite standard to do in lattice model where the continuum is described by a
Dirac equation [34, 69, 70]. For coordinate dependent coefficients, this is
the Gullstrand-Painlevé metric [24] also know as the acoustic metric which is
the Schwarzschild metric of a $(1+1)$D black hole expressed in Gullstrand-
Painlevé coordinates. This metric is referred to here as an internal metric of
the model as it depends upon the internal couplings of the Hamiltonian and not
the physical geometry of the lattice. In addition, this is a fixed classical
background metric and the quantum fields have no back-reaction on the metric.
Quite remarkably, the phase boundary between the regions for $|v|/2<u$ and
$|v|/2>u$ can be interpreted as the inner and outer regions of the black hole,
where the phase boundary $|u|=|v|/2$ aligns with the event horizon.

In order to transform the metric of Eq. (47) into the Schwarzschild metric, a
coordinate transformation defined as $(t,x)\mapsto(\tau,x)$ is used, where

$\tau(t,x)=t-\int_{x_{0}}^{x}\mathrm{d}z\frac{V(z)}{U^{2}-V^{2}(z)},$ (48)

where we have absorbed some factors into the coupling as $u/4\rightarrow U$
and $v/8\rightarrow V$. This maps the metric to

$\mathrm{d}s^{2}=\left(1-\frac{V^{2}(x)}{U^{2}(x)}\right)\mathrm{d}\tau^{2}-\frac{1}{U^{2}(x)\left(1-\frac{V^{2}(x)}{U^{2}(x)}\right)}\mathrm{d}x^{2},$
(49)

which takes the general form of the Schwarzschild metric, where the horizon is
at the location where the metric becomes singular at $U=V$ which is equivalent
to $u=v/2$.

##  Appendix B Fitting the exponent of the Laypunov vs temperature

Figure 6: For different values of $v$, we fit $\lambda=(a(T^{b}-c)$ and
extract the exponent $b$, plotting this against $1/N$. We see in XY phase
($v=1$), this value tends towards 2 (red dashed line) while in the chiral
phase, $v>7.5$, this value tends towards 1 (black dashed line).

In Fig. 4, we present the Lyaponov exponent vs temperature in two regimes: in
the XY phase ($v=1$) where the growth is quadratic, and deep in the chiral
phase ($v=8$) where the growth is shown to be linear. In this section, we
explicity show that these fitting forms are suitable in the two limits by
additionally parameterising the fit of $\lambda$ vs $T$ to $a(T^{b}-c)$ and
plotting $b$ against $N$. The results are shown in Fig. 6 where we see in XY
phase, $b\approx 2$, while in the chiral phase $b\approx 1$, as one would
hope. It is important to note that the results of the fit fluctuate dependent
on the choice of $T_{\rm{min}}$ and $T_{\rm{max}}$ in the fitting window. To
achieve this, different values have been taken for different $N$ and $v$, so
one must be cautious during the fitting procedure. These results, nonetheless,
demonstrate that with a suitable choice of fitting window, a quadratic and
linear fit are reasonable to make.