\section{Introduction}

Interfaces of quantum systems offer a fertile environment for rich and exotic physics to emerge that often cannot be met without the support of bulk systems. For example, domain walls between fractional quantum Hall states can give rise to parafermions, anyons with non-Abelian statistics \cite{Clarke}, while domain walls between 2D Heisenberg models can give rise to deconfined fractional excitations \cite{Batista}. Moreover, higher order topological phases support gapless edge states at boundary defects \cite{Benalcazar,Langbehn,Ezawa} and non-unitary conformal field theories can emerge at the boundary of interacting field theories \cite{Soderberg}. The complexity of interfaces, especially for interacting systems, is so high that simple and effective modelling is instrumental to obtaining a qualitative and quantitative understanding of their behaviour.

Here we investigate the interface between two spin-$1/2$ chains that are in different chiral phases. We consider the XY model Hamiltonian supplemented by a three-spin chiral operator, i.e., our system is intrinsically interacting. Such chiral systems are of much interest as they exhibit a rich spectrum of quantum correlations \cite{Pachos1} and can give rise to skyrmionic configurations \cite{Tikhonov}. We find that such a chiral interface can be effectively modelled by a black hole event horizon. Specifically, we demonstrate that the mean field theory of the system can be modelled as a Dirac fermion in the curved background geometry of a black hole, where the event horizon is positioned at the interface between the two phases. The chiral phase is identified with the inside of the black hole, whilst the non-chiral phase is identified with the asymptotically flat spacetime exterior to the black hole. 

To establish the validity of the mean field theory that gives rise to the back hole description and further understand the behaviour of the interacting chiral system, we model it both numerically with Matrix Product State (MPS) methods and analytically by bosonising the system and using the theory of Luttinger liquids. In this way, we show that the horizon description is quantitatively and quantitatively faithful to the one provided by mean field theory. We further demonstrate that the interface is between two conformal field theories with central charge $c=2$ inside the black hole and $c=1$ outside, thus identifying the change in the fermionic degrees of freedom across the interface. 

To test the faithfulness of the description of our system with black hole physics, we use the mean field description to investigate the time evolution following a quench that propagates through the horizon. We take an interface between two opposite chiralities that models a black hole and white hole opposite from each other. We demonstrate numerically that a pulse in one chiral phase is transmitted through the interface as thermal radiation to the other chiral phase, i.e., the pulse is thermalised when it passes through the event horizon. The temperature of the thermal radiation is approximated well by the Hawking temperature for a wide range of coupling profiles and initial conditions. Hence, the horizon physics can provide a high level modelling of the chiral interface that accurately predicts the evolution of the interacting chiral phases across a phase boundary. We envision that gravity at extreme curvatures can provide an elegant formalism that can efficiently model several strongly interacting systems and their interfaces in higher dimensions.

In this work we consider a chain of $N$ spin-$1/2$ particles with Hamiltonian
\begin{equation}
H = \sum_{n=1}^N \left[- \frac{u}{2}  \left(\sigma^x_n \sigma^x_{n+1} + \sigma^y_n \sigma^y_{n+1} \right) + \frac{v}{4} \chi_n \right] 
\label{eq:spin_ham},
\end{equation}
where $u,v \in \mathbb{R}$, $\{ \sigma^x_n , \sigma^y_n , \sigma^z_n \}$ are the Pauli matrices of the $n$th spin and $\chi_n$ is the spin chirality given by the three-spin interaction 
\begin{equation}
\chi_n = \boldsymbol{\sigma}_n \cdot \boldsymbol{\sigma}_{n+1} \times \boldsymbol{\sigma}_{n+2}, 
\end{equation}
where $\boldsymbol{\sigma}_n$ is the vector of Pauli matrices of the $n$th spin \cite{Pachos1,Pachos2}. The chirality operator is a measure of the solid angle spanned by three neighbouring spins. 

The model given in Eq. \eqref{eq:spin_ham} can be mapped to an interacting fermionic Hamiltonian via a Jordan-Wigner transformation. After application of self-consistent mean field (MF) theory to Eq.  \eqref{eq:spin_ham} (see Supplementary Material) we obtain 
\begin{equation}
H_\text{MF} = \sum_{n=1}^N \left( -u c^\dagger_n c_{n+1} - \frac{iv}{2} c_n^\dagger c_{n+2} \right) + \text{h.c.}.
\label{eq:bh_ham}
\end{equation}
This free Hamiltonian can be diagonalised exactly with a Fourier transform giving a gapless dispersion relation
\begin{equation}
E(p) = -2u \cos(p) + v\sin(2p), 
\label{eq:dispersion}
\end{equation}
which contains a single band.

\begin{figure}[t]
\begin{center}
\includegraphics[scale=1]{FIG1.pdf}
\end{center}
\caption{(a) The nearest- and next-to-nearest-neighbour tunnelling of fermionic Hamiltonian of Eq. \eqref{eq:bh_ham} with diatomic colouring, $A$ and $B$. The chirality operator, $\chi_n$, is defined for each triangular plaquette. (b) An example of interface profile for the couplings $u$ and $v$ across the chain with the corresponding ground state chirality, $\langle \chi_n\rangle$. The system changes from a chiral phase, $\langle \chi_n\rangle\neq 0$, when $|v| > |u|$, to non-chiral phase, $\langle \chi_n\rangle=0$, when $|v| < |u|$, across the interface, $n_\mathrm{h}$, corresponding to the event horizon in the continuum theory. The chirality, $\langle \chi_n\rangle$, is determined by mean field (MF) theory and by numerical modelling of the spin chain (DMRG with bond dimension $D=400$) for system size $N = 200$. (c) The corresponding mean field dispersion relations are obtained from a diatomic representation of Eq. (\ref{eq:dispersion}). The sign of $v$ determines the direction the cones tilt.} 
\label{fig_1} 
\end{figure}

Suppose we introduce two sublattices, $A$ and $B$, as shown in Fig. \ref{fig_1}(a). This has the effect of folding the Brillouin zone giving two bands, as shown in Fig. \ref{fig_1}(c). This interpretation reveales the relativistic dispersion exhibited by the model at low energy. There is a Dirac cone located at $p_0 = \pi$, which tilts as the magnitude of $v$ is increased. For $|v| < |u|$ we have a type-I Weyl fermion with a single Fermi point; for $|u| = |v|$ we have critical tilting giving rise to a type-III Weyl fermion, whereby one band is flat; and for $|v| > |u|$ we have an overtilted cone, giving rise to type-II Weyl fermions \cite{Volovik1,Volovik2,Volovik3}. In the latter case, additional Fermi points appear in the dispersion relation, due to the Nielsen-Ninomiya theorem~\cite{Nielsen1,Nielsen2}.

To find an effective geometric description of the chiral system we derive the corresponding field theoretic description. The continuum limit of $H_\text{MF}$ is found by Taylor expanding $H_\text{MF}$ in momentum space about $p_0$, the Fermi point where the energy becomes zero, to first order in $p$. This yields a Dirac Hamiltonian on a $(1+1)$D spacetime with metric 
\begin{equation}
\mathrm{d}s^2 =  \left( 1 -\frac{v^2}{u^2} \right) \mathrm{d}t^2 -  \frac{2v}{u^2} \mathrm{d}t \mathrm{d}x - \frac{1}{u^2}\mathrm{d}x^2, \label{eq:metric}
\end{equation}
which is the Schwarzschild metric written in Gullstrand-Painlev\'{e} coordinates \cite{Volovik_helium_droplet}. If the coupling $v$ is upgraded to a sufficiently slowly-varying function $v(x)$, such as in Fig. \ref{fig_1}(b), then this continuum description remains valid and an event horizon is located at the critical point $x_\mathrm{h}$, where $ |v(x_\mathrm{h})| = |u|$, which corresponds to the location of critical tilt in Fig. \ref{fig_1}(c). The corresponding Hawking temperature is given by~\cite{Volovik3} 
\begin{equation}
T_\text{H} = \frac{1}{2\pi} \left|\frac{\mathrm{d} v(x_\mathrm{h})}{ \mathrm{d}x} \right|. 
\label{eq:Hawking}
\end{equation}
Therefore, the original spin model of Eq. (\ref{eq:spin_ham}) with inhomogeneous couplings is effectively described by a $(1+1)$D black hole \cite{Volovik1,Volovik2,Volovik3,Volovik_helium_droplet}.


To investigate the nature of quantum phases supported by Eq. (\ref{eq:spin_ham}), and the transitions between them, we consider the case of homogeneous couplings $u$ and $v$. Beginning with the free fermion model derived using the mean field approximation of Eq. (\ref{eq:bh_ham}) reveals that for small $v$ the system is in a disordered, gapless, XY phase, while as $v$ increases it passes through a second-order phase transition into a gapless chiral phase, corresponding to a non-zero ground state chirality $\langle \chi_n \rangle$, as seen in Fig. \ref{fig_2}(a), so the chirality behaves as an order parameter. In the mean-field approximation the chiral phase transition is located at $|v|=|u|$, coinciding with the critical tilting of the Dirac cones and the appearance of additional Fermi points, as shown in Fig. \ref{fig_1}(c). Near the critical point, 
\begin{equation}
\langle \chi_n \rangle \sim (v-v_\mathrm{c})^\gamma,
\end{equation}
with critical point $v_c = u$ and critical exponent $\gamma = 1$. On the other hand, studying the full spin model of Eq. (\ref{eq:spin_ham}) using finite DMRG \cite{schollwock2011density}, we estimate that the phase transition of the full model is located at $v_\mathrm{c} \approx 1.12u$ with a critical exponent of  $\gamma \approx 0.39$. A comparison between the chirality of the full spin model at various system sizes and the mean-field model for inhomogeneous and homogenous couplings can be seen in Fig. \ref{fig_1}(b) and Fig. \ref{fig_2}(a) respectively, showing that the chirality behaves as an order parameter and the effectiveness of the mean field approximation.

To gain further insight into the nature of the chiral phase transition, we consider the behaviour of the bipartite entanglement entropy as $v$ is increased. As discussed above, the model is gapless for all $v$, i.e. it can be destined by a conformal field theory (CFT). In this case we expect the ground state entanglement entropy of a partition of $L \ll N$ spins to obey the Cardy formula
\begin{equation}
S_L = \frac{c}{3} \ln L + S_0, \label{eq:entropy}
\end{equation}
where $c$ is the central charge of the CFT and $S_0$ is a constant \cite{Calabrese}. Using this formula, we estimate the central charge as a function of $v$ for the full spin model and the mean field approximation. In Fig. \ref{fig_2}(b) we see that $ c \approx 1 $ in the XY phase which jumps to $ c \approx 2$ in the chiral phase, with good agreement between the spin and mean field results. We can clearly interpret this in the mean-field model: the additional Fermi points appearing when $|v| > |u|$ cause the model to transition from a $c=1$ CFT with a single Dirac fermion to a $c=2=1+1$ CFT with two Dirac fermions, as seen by the additional Fermi points of the dispersion in Fig. \ref{fig_1}(c). This can also be understood from the lattice structure of the MF model, as seen in Fig.~\ref{fig_1}(a), where for $|v|\ll |u|$ a single zig-zag fermionic chain dominates ($c=1$) while for $|v|\gg |u|$ two fermionic chains dominate, corresponding to the stands of the ladder, thus effectively doubling the degrees of freedom ($c=2$).

\begin{figure}
\begin{center}
\includegraphics[scale=1]{FIG2.pdf}
\end{center}
\caption{(a) The average ground state chirality, $\langle \chi \rangle = \sum_n \langle \chi_n \rangle/N $, for the mean-field (MF) model ($N=500$) and spin model found using DMRG ($N=200$, $D=300$) where $u =1$. (b) The central charge, $c$, obtained for the MF model ($N=500$) and the spin model found using DMRG ($v\leq u$: $N=200$, $D=300$ and $v > u$: $N=160$, $D=800$). (c)-(d) The Fermi velocities, $v_\mathrm{L}$ and $v_\mathrm{R}$ respectively ($u=1$) derived from the mean field (MF) and Luttinger liquid descriptions compared to the numerical results of the MPS excitation ansatz for the spin model at bond dimension $D=36$ in the thermodynamic limit.}
\label{fig_2}
\end{figure}

As we have seen numerically, the mean field approximation can faithfully reproduce many of the features of the fully interacting model, especially for $|v| < |u|$, suggesting that the interactions are not significant. We investigate this further by bosonising the spin Hamiltonian for $|v| < |u|$, mapping it to a Luttinger liquid \cite{Giamarchi,Miranda,Sreemayee}. After a Jordan-Wigner transformation, the spin Hamiltonian takes the form $H = H_\mathrm{MF} + H_\mathrm{int}$, where $H_\mathrm{MF}$ is the quadratic mean-field Hamiltonian of Eq. (\ref{eq:bh_ham}) and $H_\mathrm{int}$ is an interaction term containing quartic terms. For $|v| < |u|$, the single-band mean field dispersion relation of Eq. (\ref{eq:dispersion}) suggests that the spin model has two Fermi points located at $p_\mathrm{R,L} = \pm \pi/2$, with the Fermi velocities $v_\mathrm{R,L} = 2(\pm u - v)$. By expanding around these Fermi points and bosonising the interaction terms using the methods of \cite{Giamarchi,Miranda,Sreemayee}, we find the fully interacting Hamiltonian can be mapped to the free boson Hamiltonian
\begin{equation}
H = u \int \mathrm{d}x \left[ \Pi^2 + (\partial_x \Phi)^2 \right],
\end{equation}
where the fields obey the canonical commutation relations $[\Phi(x), \Pi(y)] = i\delta(x-y)$. The interactions of $H_\mathrm{int}$ simply rescaled the Fermi velocities $v_\mathrm{R,L} \rightarrow v'_\mathrm{R,L} = 2\left( \pm u - v\left(1-2/\pi \right) \right)$, but leaves the Luttinger parameter unchanged at $K = 1$ (see Supplementary Material), suggesting the model remains a non-interacting free fermion model.


The dispersion of the full spin model as a function of $v$ for $|v| < |u|$ can be calculated using the MPS excitation ansatz working in the thermodynamic limit \cite{haegeman2013post}. The full dispersion features unequal left- and right-moving Fermi velocities whose magnitudes change oppositely with $v$---the signature of tilting of the cones---and the appearance of additional Fermi points, similar to the mean field model. In Fig. \ref{fig_2}(c) and \ref{fig_2}(d), the Fermi velocities $v_\mathrm{L,R}$  obtained from the mean field approximation, the Luttinger liquid model and the spin Hamiltonian are compared. We see that the Luttinger liquid model is much more accurate than the mean field approximation. We expect that the quantitative disagreement and the observed asymmetry in the change to the left and right velocities to be lifted at higher order in perturbation theory.

\begin{figure}[t]
\begin{center}
\includegraphics[scale=1]{FIG3.pdf}
\end{center}
\caption{Thermalisation of the mean field (MF) model after a quench to the Hawking temperature $T_\mathrm{H}$. (a) the lattice wavefunction $\psi_n$ on the right half of the system ($n \in [n_\text{h},N]$)  transmitted through the horizon, for the couplings $u=1$ and $v$ given by Eq. \eqref{en:couplings} with $\alpha = 20$, $\beta = 0.1$, and the horizon at $n_\text{h} = N/2$ with $N=500$. The particle tunnels across at $t \approx 2$ and a small wavepacket escapes into the other half, which we interpret as Hawking radiation. (b) A snapshot of the overlap $-\ln P$ vs. the energy of the state $E$ at time $t = 4.5$. The system thermalises shortly after the particle passes through the interface, displaying a linear dependence on $E$, where the gradient is given by $1/T$. (c) The extracted temperature $T$ of the radiation vs. $\alpha$ extracted after a short time $t=4.5$. $T$ grows linearly with $\alpha$, very close to the predictions of the Hawking formula $T_\text{H} \approx \alpha \beta/2 \pi$.}
\label{fig_3}
\end{figure}

\begin{figure}[t]
\begin{center}
\includegraphics[scale=1]{FIG4.pdf}
\end{center}
\caption{(a) The measured temperature $T$ vs. the distance $|n_0 - n_\mathrm{h}|$ from the event horizon that the particle is released for the mean field (MF) system of size $N = 500$, $\alpha = 15$, $\beta = 0.1$ and $n_\mathrm{h} = N/2$. We see the temperature is insensitive to the initial position $n_0$ of the particle. This effect breaks down if $n_0$ is very close to the horizon or too far away, near the boundary of the system. (b) The measured temperature $T$ vs. the position of the horizon $n_\mathrm{h}$ for the same system, where $n_0 = n_\mathrm{h} - 25$.  The temperature is insensitive to the location of the horizon, except when it gets too close to the boundaries, in which case it breaks down.}
\label{fig_4}
\end{figure}

It has been shown in the literature that many analogue gravitational systems in condensed matter will exhibit a Hawking-like effect \cite{Volovik1,Volovik2,Volovik3,Yang,Sabsovich,Maertens,Huang,Hang,Guan,Retzker,Rodriguez,Roldan-Molina,Kosior,Steinhauer_2016}. Reversing the argument, we would like to see if the Hawking radiation can effectively describe quenched time evolutions across the chiral interface. To be able to simulate large system sizes and long evolution times, we resort to the mean field description of Eq. (\ref{eq:bh_ham}) rather than the full spin model. Consider an open, inhomogeneous system with couplings $u(x) = 1$ and
\begin{equation}
v(x) = \alpha \tanh[\beta(x-x_\mathrm{h})],
\label{en:couplings}
\end{equation}
where $\alpha ,\beta \in \mathbb{R}$ and $x_\text{h}$ is in the centre of the system. Here, we take $x$ as the unit cell coordinate in order to align with our continuum limit conventions (see Supplementary Material). This produces a positive and negative chiral region separated by a small zero-chirality region in the centre of the system. For large enough $\alpha$ and $\beta$, the zero-chirality region has arbitrarily small size, so the system effectively models an interface between two oppositely polarized chiral phases, corresponding to a black hole-white hole interface respectively in the continuum.

Following the method of Ref. \cite{Yang}, we initialise a single-particle state $|n_0 \rangle = c^\dagger_{n_0}|0\rangle$ on the $n_0$th lattice site inside the left half of the system, and let the wavefunction evolve freely across the boundary into the other half via the Hamiltonian $H_\mathrm{MF}$, as shown in Fig. \ref{fig_3}(a). We then measure the overlap of the wavefunction with localised energy modes that exist only on the other side of the boundary as
\begin{equation}
P(k,t) = |\langle k | e^{-iHt} |n_0\rangle|^2,
\end{equation}
where $|k \rangle$ are the single-particle eigenstates of the Hamiltonian $H_\mathrm{out}$, where $H_\mathrm{out}$ is the Hamiltonian of Eq. (\ref{eq:bh_ham}) truncated to the outside region only. This method utilises the result that Hawking radiation can be viewed as quantum tunnelling \cite{Wilczek}.

We find numerically that the interface between the two chiral phases thermalises the wavefunction: once the wavefunction evolves across the interface, the distribution takes the form $P(k,t) \propto e^{-E(k)/T}$, where $T$ is some effective temperature. Fig. \ref{fig_3}(b) shows the distribution $P(k,t)$ at time $t = 4.5$ for a system with parameters $N= 500$, $n_h=250$, $\alpha = 20$ and $\beta = 0.1$, where we prepared the particle at $n_0 = 230$. It is clear that $P(k,t)$ follows a Boltzmann distribution at some temperature $T$, where the gradient of the line is given by $1/T$. We observe the system strongly thermalises. The value $\beta = 0.1$ is taken to suppress the effects from having finite lattice spacing and finite system size. In Fig. \ref{fig_3}(c), we present the dependence of the measured temperature $T$ on the magnitude $\alpha$ for $\beta = 0.1$. We see it closely follows the predictions of the Hawking formula $T_\mathrm{H} = \alpha \beta/2 \pi$, obtained from Eq. \eqref{eq:Hawking} and Eq. \eqref{en:couplings}, for a wide range of couplings, $\alpha$, thus accurately modelling the physics of the chiral interface. The thermalisation to the Hawking $T$ breaks down when $\alpha < 4$ as the couplings will not be sharp enough to provide a sufficient interface, whereas for large $\alpha$ the couplings vary too fast for the continuum approximation to be valid, which is where the black hole physics should emerge.

The Hawking temperature $T_\text{H}= \alpha \beta/2 \pi$ of our system is a very simple analytic formula that describes a complex thermalisation process. In particular, it does not depend on the initial value $n_0$, i.e., the position where quench starts, nor the horizon location $n_\text{h}$ which effectively gives the size ratio between the two chiral phases in our system. To verify these properties, we numerically determine the dependence of $T$ on $n_0$ (Fig. \ref{fig_4}(a)) and $n_\text{h}$ (Fig. \ref{fig_4}(b)). We see that the measured temperature $T$ is largely insensitive to the initial position of the particle $n_0$. The black hole description only fails if $n_0$ is initially too close or too far away from the interface or when the interface $n_\mathrm{h}$ is too close the edges of the system. In all these cases boundary effects start to contribute and the exterior region which the overlap $P(k,t)$ is measured in becomes too small. These two observations show that the thermalisation across the interface is robust and will aid in any potential experimental realisation of the model.

In this work, we demonstrated that that low-energy behaviour of a chiral fermionic phase can be described by Dirac fermions on a black hole background. The interface of such a chiral phase with a non-chiral phase gives rise to an event horizon with quantum properties modelled by Hawking radiation. This analogy allowed us to predict the time evolution of a quenched system with a chiral interface which should thermalise to the Hawking temperature $T_\mathrm{H}$. We demonstrated that this is indeed the case for a wide variety of quenches, positions of the interface and coupling parameters, thus providing a faithful high level description of chiral interfaces. To demonstrate the faithfulness of our approach we employed mean field theory, matrix product state techniques and bosonisation to fully characterise the phase diagram of the model. We envision that this bridge between chiral systems and black holes can facilitate the quantum simulation of Hawking radiation, e.g. with cold atom technology \cite{Jaksch,Rodriguez,Kosior}. Moreover, our investigation opens the way for modelling certain strongly correlated systems by effective geometric theories with extreme curvature, thus providing an intuitive tool for their analytical investigation.

\textbf{Acknowledgements}: We thank Patricio Salgado-Rebolledo, Joe Barker and Diptiman Sen for insightful discussions. M.D.H., A.H. and J.K.P. acknowledge support by EPSRC (Grant No. EP/R020612/1).  


\section{Introduction}

Interfaces of quantum systems offer a fertile environment for rich and exotic physics to emerge that often cannot be met without the support of bulk systems. For example, domain walls between fractional quantum Hall states can give rise to parafermions, anyons with non-Abelian statistics \cite{Clarke}, while domain walls between 2D Heisenberg models can give rise to deconfined fractional excitations \cite{Batista}. Moreover, higher order topological phases support gapless edge states at boundary defects \cite{Benalcazar,Langbehn,Ezawa} and non-unitary conformal field theories can emerge at the boundary of interacting field theories \cite{Soderberg}. The complexity of interfaces, especially for interacting systems, is so high that simple and effective modelling is instrumental to obtaining a qualitative and quantitative understanding of their behaviour.

Here we investigate the interface between two spin-$1/2$ chains that are in different chiral phases. We consider the XY model Hamiltonian supplemented by a three-spin chiral operator, i.e., our system is intrinsically interacting. Such chiral systems are of much interest as they exhibit a rich spectrum of quantum correlations \cite{Pachos1} and can give rise to skyrmionic configurations \cite{Tikhonov}. We find that such a chiral interface can be effectively modelled by a black hole event horizon. Specifically, we demonstrate that the mean field theory of the system can be modelled as a Dirac fermion in the curved background geometry of a black hole, where the event horizon is positioned at the interface between the two phases. The chiral phase is identified with the inside of the black hole, whilst the non-chiral phase is identified with the asymptotically flat spacetime exterior to the black hole. 

To establish the validity of the mean field theory that gives rise to the back hole description and further understand the behaviour of the interacting chiral system, we model it both numerically with Matrix Product State (MPS) methods and analytically by bosonising the system and using the theory of Luttinger liquids. In this way, we show that the horizon description is quantitatively and quantitatively faithful to the one provided by mean field theory. We further demonstrate that the interface is between two conformal field theories with central charge $c=2$ inside the black hole and $c=1$ outside, thus identifying the change in the fermionic degrees of freedom across the interface. 

To test the faithfulness of the description of our system with black hole physics, we use the mean field description to investigate the time evolution following a quench that propagates through the horizon. We take an interface between two opposite chiralities that models a black hole and white hole opposite from each other. We demonstrate numerically that a pulse in one chiral phase is transmitted through the interface as thermal radiation to the other chiral phase, i.e., the pulse is thermalised when it passes through the event horizon. The temperature of the thermal radiation is approximated well by the Hawking temperature for a wide range of coupling profiles and initial conditions. Hence, the horizon physics can provide a high level modelling of the chiral interface that accurately predicts the evolution of the interacting chiral phases across a phase boundary. We envision that gravity at extreme curvatures can provide an elegant formalism that can efficiently model several strongly interacting systems and their interfaces in higher dimensions.

In this work we consider a chain of $N$ spin-$1/2$ particles with Hamiltonian
\begin{equation}
H = \sum_{n=1}^N \left[- \frac{u}{2}  \left(\sigma^x_n \sigma^x_{n+1} + \sigma^y_n \sigma^y_{n+1} \right) + \frac{v}{4} \chi_n \right] 
\label{eq:spin_ham},
\end{equation}
where $u,v \in \mathbb{R}$, $\{ \sigma^x_n , \sigma^y_n , \sigma^z_n \}$ are the Pauli matrices of the $n$th spin and $\chi_n$ is the spin chirality given by the three-spin interaction 
\begin{equation}
\chi_n = \boldsymbol{\sigma}_n \cdot \boldsymbol{\sigma}_{n+1} \times \boldsymbol{\sigma}_{n+2}, 
\end{equation}
where $\boldsymbol{\sigma}_n$ is the vector of Pauli matrices of the $n$th spin \cite{Pachos1,Pachos2}. The chirality operator is a measure of the solid angle spanned by three neighbouring spins. 

The model given in Eq. \eqref{eq:spin_ham} can be mapped to an interacting fermionic Hamiltonian via a Jordan-Wigner transformation. After application of self-consistent mean field (MF) theory to Eq.  \eqref{eq:spin_ham} (see Supplementary Material) we obtain 
\begin{equation}
H_\text{MF} = \sum_{n=1}^N \left( -u c^\dagger_n c_{n+1} - \frac{iv}{2} c_n^\dagger c_{n+2} \right) + \text{h.c.}.
\label{eq:bh_ham}
\end{equation}
This free Hamiltonian can be diagonalised exactly with a Fourier transform giving a gapless dispersion relation
\begin{equation}
E(p) = -2u \cos(p) + v\sin(2p), 
\label{eq:dispersion}
\end{equation}
which contains a single band.

\begin{figure}[t]
\begin{center}
\includegraphics[scale=1]{FIG1.pdf}
\end{center}
\caption{(a) The nearest- and next-to-nearest-neighbour tunnelling of fermionic Hamiltonian of Eq. \eqref{eq:bh_ham} with diatomic colouring, $A$ and $B$. The chirality operator, $\chi_n$, is defined for each triangular plaquette. (b) An example of interface profile for the couplings $u$ and $v$ across the chain with the corresponding ground state chirality, $\langle \chi_n\rangle$. The system changes from a chiral phase, $\langle \chi_n\rangle\neq 0$, when $|v| > |u|$, to non-chiral phase, $\langle \chi_n\rangle=0$, when $|v| < |u|$, across the interface, $n_\mathrm{h}$, corresponding to the event horizon in the continuum theory. The chirality, $\langle \chi_n\rangle$, is determined by mean field (MF) theory and by numerical modelling of the spin chain (DMRG with bond dimension $D=400$) for system size $N = 200$. (c) The corresponding mean field dispersion relations are obtained from a diatomic representation of Eq. (\ref{eq:dispersion}). The sign of $v$ determines the direction the cones tilt.} 
\label{fig_1} 
\end{figure}

Suppose we introduce two sublattices, $A$ and $B$, as shown in Fig. \ref{fig_1}(a). This has the effect of folding the Brillouin zone giving two bands, as shown in Fig. \ref{fig_1}(c). This interpretation reveales the relativistic dispersion exhibited by the model at low energy. There is a Dirac cone located at $p_0 = \pi$, which tilts as the magnitude of $v$ is increased. For $|v| < |u|$ we have a type-I Weyl fermion with a single Fermi point; for $|u| = |v|$ we have critical tilting giving rise to a type-III Weyl fermion, whereby one band is flat; and for $|v| > |u|$ we have an overtilted cone, giving rise to type-II Weyl fermions \cite{Volovik1,Volovik2,Volovik3}. In the latter case, additional Fermi points appear in the dispersion relation, due to the Nielsen-Ninomiya theorem~\cite{Nielsen1,Nielsen2}.

To find an effective geometric description of the chiral system we derive the corresponding field theoretic description. The continuum limit of $H_\text{MF}$ is found by Taylor expanding $H_\text{MF}$ in momentum space about $p_0$, the Fermi point where the energy becomes zero, to first order in $p$. This yields a Dirac Hamiltonian on a $(1+1)$D spacetime with metric 
\begin{equation}
\mathrm{d}s^2 =  \left( 1 -\frac{v^2}{u^2} \right) \mathrm{d}t^2 -  \frac{2v}{u^2} \mathrm{d}t \mathrm{d}x - \frac{1}{u^2}\mathrm{d}x^2, \label{eq:metric}
\end{equation}
which is the Schwarzschild metric written in Gullstrand-Painlev\'{e} coordinates \cite{Volovik_helium_droplet}. If the coupling $v$ is upgraded to a sufficiently slowly-varying function $v(x)$, such as in Fig. \ref{fig_1}(b), then this continuum description remains valid and an event horizon is located at the critical point $x_\mathrm{h}$, where $ |v(x_\mathrm{h})| = |u|$, which corresponds to the location of critical tilt in Fig. \ref{fig_1}(c). The corresponding Hawking temperature is given by~\cite{Volovik3} 
\begin{equation}
T_\text{H} = \frac{1}{2\pi} \left|\frac{\mathrm{d} v(x_\mathrm{h})}{ \mathrm{d}x} \right|. 
\label{eq:Hawking}
\end{equation}
Therefore, the original spin model of Eq. (\ref{eq:spin_ham}) with inhomogeneous couplings is effectively described by a $(1+1)$D black hole \cite{Volovik1,Volovik2,Volovik3,Volovik_helium_droplet}.


To investigate the nature of quantum phases supported by Eq. (\ref{eq:spin_ham}), and the transitions between them, we consider the case of homogeneous couplings $u$ and $v$. Beginning with the free fermion model derived using the mean field approximation of Eq. (\ref{eq:bh_ham}) reveals that for small $v$ the system is in a disordered, gapless, XY phase, while as $v$ increases it passes through a second-order phase transition into a gapless chiral phase, corresponding to a non-zero ground state chirality $\langle \chi_n \rangle$, as seen in Fig. \ref{fig_2}(a), so the chirality behaves as an order parameter. In the mean-field approximation the chiral phase transition is located at $|v|=|u|$, coinciding with the critical tilting of the Dirac cones and the appearance of additional Fermi points, as shown in Fig. \ref{fig_1}(c). Near the critical point, 
\begin{equation}
\langle \chi_n \rangle \sim (v-v_\mathrm{c})^\gamma,
\end{equation}
with critical point $v_c = u$ and critical exponent $\gamma = 1$. On the other hand, studying the full spin model of Eq. (\ref{eq:spin_ham}) using finite DMRG \cite{schollwock2011density}, we estimate that the phase transition of the full model is located at $v_\mathrm{c} \approx 1.12u$ with a critical exponent of  $\gamma \approx 0.39$. A comparison between the chirality of the full spin model at various system sizes and the mean-field model for inhomogeneous and homogenous couplings can be seen in Fig. \ref{fig_1}(b) and Fig. \ref{fig_2}(a) respectively, showing that the chirality behaves as an order parameter and the effectiveness of the mean field approximation.

To gain further insight into the nature of the chiral phase transition, we consider the behaviour of the bipartite entanglement entropy as $v$ is increased. As discussed above, the model is gapless for all $v$, i.e. it can be destined by a conformal field theory (CFT). In this case we expect the ground state entanglement entropy of a partition of $L \ll N$ spins to obey the Cardy formula
\begin{equation}
S_L = \frac{c}{3} \ln L + S_0, \label{eq:entropy}
\end{equation}
where $c$ is the central charge of the CFT and $S_0$ is a constant \cite{Calabrese}. Using this formula, we estimate the central charge as a function of $v$ for the full spin model and the mean field approximation. In Fig. \ref{fig_2}(b) we see that $ c \approx 1 $ in the XY phase which jumps to $ c \approx 2$ in the chiral phase, with good agreement between the spin and mean field results. We can clearly interpret this in the mean-field model: the additional Fermi points appearing when $|v| > |u|$ cause the model to transition from a $c=1$ CFT with a single Dirac fermion to a $c=2=1+1$ CFT with two Dirac fermions, as seen by the additional Fermi points of the dispersion in Fig. \ref{fig_1}(c). This can also be understood from the lattice structure of the MF model, as seen in Fig.~\ref{fig_1}(a), where for $|v|\ll |u|$ a single zig-zag fermionic chain dominates ($c=1$) while for $|v|\gg |u|$ two fermionic chains dominate, corresponding to the stands of the ladder, thus effectively doubling the degrees of freedom ($c=2$).

\begin{figure}
\begin{center}
\includegraphics[scale=1]{FIG2.pdf}
\end{center}
\caption{(a) The average ground state chirality, $\langle \chi \rangle = \sum_n \langle \chi_n \rangle/N $, for the mean-field (MF) model ($N=500$) and spin model found using DMRG ($N=200$, $D=300$) where $u =1$. (b) The central charge, $c$, obtained for the MF model ($N=500$) and the spin model found using DMRG ($v\leq u$: $N=200$, $D=300$ and $v > u$: $N=160$, $D=800$). (c)-(d) The Fermi velocities, $v_\mathrm{L}$ and $v_\mathrm{R}$ respectively ($u=1$) derived from the mean field (MF) and Luttinger liquid descriptions compared to the numerical results of the MPS excitation ansatz for the spin model at bond dimension $D=36$ in the thermodynamic limit.}
\label{fig_2}
\end{figure}

As we have seen numerically, the mean field approximation can faithfully reproduce many of the features of the fully interacting model, especially for $|v| < |u|$, suggesting that the interactions are not significant. We investigate this further by bosonising the spin Hamiltonian for $|v| < |u|$, mapping it to a Luttinger liquid \cite{Giamarchi,Miranda,Sreemayee}. After a Jordan-Wigner transformation, the spin Hamiltonian takes the form $H = H_\mathrm{MF} + H_\mathrm{int}$, where $H_\mathrm{MF}$ is the quadratic mean-field Hamiltonian of Eq. (\ref{eq:bh_ham}) and $H_\mathrm{int}$ is an interaction term containing quartic terms. For $|v| < |u|$, the single-band mean field dispersion relation of Eq. (\ref{eq:dispersion}) suggests that the spin model has two Fermi points located at $p_\mathrm{R,L} = \pm \pi/2$, with the Fermi velocities $v_\mathrm{R,L} = 2(\pm u - v)$. By expanding around these Fermi points and bosonising the interaction terms using the methods of \cite{Giamarchi,Miranda,Sreemayee}, we find the fully interacting Hamiltonian can be mapped to the free boson Hamiltonian
\begin{equation}
H = u \int \mathrm{d}x \left[ \Pi^2 + (\partial_x \Phi)^2 \right],
\end{equation}
where the fields obey the canonical commutation relations $[\Phi(x), \Pi(y)] = i\delta(x-y)$. The interactions of $H_\mathrm{int}$ simply rescaled the Fermi velocities $v_\mathrm{R,L} \rightarrow v'_\mathrm{R,L} = 2\left( \pm u - v\left(1-2/\pi \right) \right)$, but leaves the Luttinger parameter unchanged at $K = 1$ (see Supplementary Material), suggesting the model remains a non-interacting free fermion model.


The dispersion of the full spin model as a function of $v$ for $|v| < |u|$ can be calculated using the MPS excitation ansatz working in the thermodynamic limit \cite{haegeman2013post}. The full dispersion features unequal left- and right-moving Fermi velocities whose magnitudes change oppositely with $v$---the signature of tilting of the cones---and the appearance of additional Fermi points, similar to the mean field model. In Fig. \ref{fig_2}(c) and \ref{fig_2}(d), the Fermi velocities $v_\mathrm{L,R}$  obtained from the mean field approximation, the Luttinger liquid model and the spin Hamiltonian are compared. We see that the Luttinger liquid model is much more accurate than the mean field approximation. We expect that the quantitative disagreement and the observed asymmetry in the change to the left and right velocities to be lifted at higher order in perturbation theory.

\begin{figure}[t]
\begin{center}
\includegraphics[scale=1]{FIG3.pdf}
\end{center}
\caption{Thermalisation of the mean field (MF) model after a quench to the Hawking temperature $T_\mathrm{H}$. (a) the lattice wavefunction $\psi_n$ on the right half of the system ($n \in [n_\text{h},N]$)  transmitted through the horizon, for the couplings $u=1$ and $v$ given by Eq. \eqref{en:couplings} with $\alpha = 20$, $\beta = 0.1$, and the horizon at $n_\text{h} = N/2$ with $N=500$. The particle tunnels across at $t \approx 2$ and a small wavepacket escapes into the other half, which we interpret as Hawking radiation. (b) A snapshot of the overlap $-\ln P$ vs. the energy of the state $E$ at time $t = 4.5$. The system thermalises shortly after the particle passes through the interface, displaying a linear dependence on $E$, where the gradient is given by $1/T$. (c) The extracted temperature $T$ of the radiation vs. $\alpha$ extracted after a short time $t=4.5$. $T$ grows linearly with $\alpha$, very close to the predictions of the Hawking formula $T_\text{H} \approx \alpha \beta/2 \pi$.}
\label{fig_3}
\end{figure}

\begin{figure}[t]
\begin{center}
\includegraphics[scale=1]{FIG4.pdf}
\end{center}
\caption{(a) The measured temperature $T$ vs. the distance $|n_0 - n_\mathrm{h}|$ from the event horizon that the particle is released for the mean field (MF) system of size $N = 500$, $\alpha = 15$, $\beta = 0.1$ and $n_\mathrm{h} = N/2$. We see the temperature is insensitive to the initial position $n_0$ of the particle. This effect breaks down if $n_0$ is very close to the horizon or too far away, near the boundary of the system. (b) The measured temperature $T$ vs. the position of the horizon $n_\mathrm{h}$ for the same system, where $n_0 = n_\mathrm{h} - 25$.  The temperature is insensitive to the location of the horizon, except when it gets too close to the boundaries, in which case it breaks down.}
\label{fig_4}
\end{figure}

It has been shown in the literature that many analogue gravitational systems in condensed matter will exhibit a Hawking-like effect \cite{Volovik1,Volovik2,Volovik3,Yang,Sabsovich,Maertens,Huang,Hang,Guan,Retzker,Rodriguez,Roldan-Molina,Kosior,Steinhauer_2016}. Reversing the argument, we would like to see if the Hawking radiation can effectively describe quenched time evolutions across the chiral interface. To be able to simulate large system sizes and long evolution times, we resort to the mean field description of Eq. (\ref{eq:bh_ham}) rather than the full spin model. Consider an open, inhomogeneous system with couplings $u(x) = 1$ and
\begin{equation}
v(x) = \alpha \tanh[\beta(x-x_\mathrm{h})],
\label{en:couplings}
\end{equation}
where $\alpha ,\beta \in \mathbb{R}$ and $x_\text{h}$ is in the centre of the system. Here, we take $x$ as the unit cell coordinate in order to align with our continuum limit conventions (see Supplementary Material). This produces a positive and negative chiral region separated by a small zero-chirality region in the centre of the system. For large enough $\alpha$ and $\beta$, the zero-chirality region has arbitrarily small size, so the system effectively models an interface between two oppositely polarized chiral phases, corresponding to a black hole-white hole interface respectively in the continuum.

Following the method of Ref. \cite{Yang}, we initialise a single-particle state $|n_0 \rangle = c^\dagger_{n_0}|0\rangle$ on the $n_0$th lattice site inside the left half of the system, and let the wavefunction evolve freely across the boundary into the other half via the Hamiltonian $H_\mathrm{MF}$, as shown in Fig. \ref{fig_3}(a). We then measure the overlap of the wavefunction with localised energy modes that exist only on the other side of the boundary as
\begin{equation}
P(k,t) = |\langle k | e^{-iHt} |n_0\rangle|^2,
\end{equation}
where $|k \rangle$ are the single-particle eigenstates of the Hamiltonian $H_\mathrm{out}$, where $H_\mathrm{out}$ is the Hamiltonian of Eq. (\ref{eq:bh_ham}) truncated to the outside region only. This method utilises the result that Hawking radiation can be viewed as quantum tunnelling \cite{Wilczek}.

We find numerically that the interface between the two chiral phases thermalises the wavefunction: once the wavefunction evolves across the interface, the distribution takes the form $P(k,t) \propto e^{-E(k)/T}$, where $T$ is some effective temperature. Fig. \ref{fig_3}(b) shows the distribution $P(k,t)$ at time $t = 4.5$ for a system with parameters $N= 500$, $n_h=250$, $\alpha = 20$ and $\beta = 0.1$, where we prepared the particle at $n_0 = 230$. It is clear that $P(k,t)$ follows a Boltzmann distribution at some temperature $T$, where the gradient of the line is given by $1/T$. We observe the system strongly thermalises. The value $\beta = 0.1$ is taken to suppress the effects from having finite lattice spacing and finite system size. In Fig. \ref{fig_3}(c), we present the dependence of the measured temperature $T$ on the magnitude $\alpha$ for $\beta = 0.1$. We see it closely follows the predictions of the Hawking formula $T_\mathrm{H} = \alpha \beta/2 \pi$, obtained from Eq. \eqref{eq:Hawking} and Eq. \eqref{en:couplings}, for a wide range of couplings, $\alpha$, thus accurately modelling the physics of the chiral interface. The thermalisation to the Hawking $T$ breaks down when $\alpha < 4$ as the couplings will not be sharp enough to provide a sufficient interface, whereas for large $\alpha$ the couplings vary too fast for the continuum approximation to be valid, which is where the black hole physics should emerge.

The Hawking temperature $T_\text{H}= \alpha \beta/2 \pi$ of our system is a very simple analytic formula that describes a complex thermalisation process. In particular, it does not depend on the initial value $n_0$, i.e., the position where quench starts, nor the horizon location $n_\text{h}$ which effectively gives the size ratio between the two chiral phases in our system. To verify these properties, we numerically determine the dependence of $T$ on $n_0$ (Fig. \ref{fig_4}(a)) and $n_\text{h}$ (Fig. \ref{fig_4}(b)). We see that the measured temperature $T$ is largely insensitive to the initial position of the particle $n_0$. The black hole description only fails if $n_0$ is initially too close or too far away from the interface or when the interface $n_\mathrm{h}$ is too close the edges of the system. In all these cases boundary effects start to contribute and the exterior region which the overlap $P(k,t)$ is measured in becomes too small. These two observations show that the thermalisation across the interface is robust and will aid in any potential experimental realisation of the model.

In this work, we demonstrated that that low-energy behaviour of a chiral fermionic phase can be described by Dirac fermions on a black hole background. The interface of such a chiral phase with a non-chiral phase gives rise to an event horizon with quantum properties modelled by Hawking radiation. This analogy allowed us to predict the time evolution of a quenched system with a chiral interface which should thermalise to the Hawking temperature $T_\mathrm{H}$. We demonstrated that this is indeed the case for a wide variety of quenches, positions of the interface and coupling parameters, thus providing a faithful high level description of chiral interfaces. To demonstrate the faithfulness of our approach we employed mean field theory, matrix product state techniques and bosonisation to fully characterise the phase diagram of the model. We envision that this bridge between chiral systems and black holes can facilitate the quantum simulation of Hawking radiation, e.g. with cold atom technology \cite{Jaksch,Rodriguez,Kosior}. Moreover, our investigation opens the way for modelling certain strongly correlated systems by effective geometric theories with extreme curvature, thus providing an intuitive tool for their analytical investigation.

\textbf{Acknowledgements}: We thank Patricio Salgado-Rebolledo, Joe Barker and Diptiman Sen for insightful discussions. M.D.H., A.H. and J.K.P. acknowledge support by EPSRC (Grant No. EP/R020612/1).  

