\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).


\section{Introduction}

The study of the structure of null and timelike circular orbits around a given generic ultracompact object is important to access several of its phenomenology properties. These properties can, in principal, be very different from the ones found for the paradigmatic General Relativity black hole (BH) -- the Kerr BH \cite{Kerr:1963ud}.
Such differences are opportune to know given the recent measurements of gravitational waves (GWs) by the LIGO/Virgo/KAGRA collaboration~\cite{LIGOScientific:2016aoc,LIGOScientific:2018mvr,LIGOScientific:2020ibl,LIGOScientific:2021djp} and the shadow and lensing of light around the M87* and Sgr A* supermassive BHs provided by the EHT collaboration~\cite{EventHorizonTelescope:2019dse,EventHorizonTelescope:2022xnr}.
Furthermore, with the advances towards the construction and implementation of the future Laser Interferometer Space Antenna (LISA), designed to be able to detect GW signals with much lower frequencies than the LIGO/Virgo/KAGRA collaboration, such as the ones produced by extreme mass ratio inspirals (EMRIs)~\cite{Barausse:2020rsu}, it will be possible to probe the structure of timelike circular orbits around these ultracompact objects, that can be BHs or other ultracompact objects that could mimic a BH -- \textit{eg.} scalar and vector boson stars \cite{PhysRev.172.1331,PhysRev.187.1767,Schunck_2003,Herdeiro:2017fhv,Brito:2015pxa,Herdeiro:2019mbz,Bustillo:2020syj,Herdeiro:2021lwl}.

To get a sense of what structures of circular orbits one can find, we start by recalling two remarkable works developed recently in~\cite{Cunha:2017qtt,Cunha:2020azh}. In the first work~\cite{Cunha:2017qtt}, the authors showed, using a topological argument, that a generic stationary, axisymmetric and asymptotically flat ultracompact horizonless object must have at least two null circular orbits -- otherwise known as light-rings (LRs) --, where one of them is stable. In the second work~\cite{Cunha:2020azh}, the authors used a similar topological argument as in the previous one to prove that a stationary, axisymmetric, asymptotically flat and non-extremal BH must have at least one unstable LR outside of its event horizon.
Both these results define precisely the structure of null circular orbits for generic (within the assumptions mentioned) ultracompact objects, with or without an horizon.

A natural question to ask now is the following: does the existence of a LR establish any structure for timelike circular orbits (TCOs) around it?
An answer to this question was given in \cite{Delgado:2021jxd}, where it was shown that, for a generic ultracompact object with the same assumptions as mentioned previously, together with a $\mathbb{Z}_2$ symmetry fixing an equatorial plane, the radial stability of the existing LRs determines the localisation and radial stability of TCOs in their vicinity. However, this work only mentions and studies the radial stability of the circular orbits, and it leaves open an opportunity to perform a more complete study where the vertical stability of the circular orbits is also included.

In this paper, we shall take the aforementioned opportunity.
For that, we will study the radial and vertical epicyclic frequencies of null and timelike particles following circular orbits on a spacetime generated by a generic ultracompact object with the same assumptions mentioned in the previous paragraph.
These epicyclic frequencies are computed by perturbing the circular orbits in the radial and vertical directions, and are tightly connected with the stability of the orbits.
With them, we shall arrive at the same result present in \cite{Delgado:2021jxd}, together with its generalisation, where we also look at the vertical stability of the orbits.

Throughout the first part of this paper, we do not impose any further assumption upon the ultracompact object in question. The generic ultracompact object may or may not be a BH, and it is a solution of an undefined set of equations of motion obtained from an undefined theory of gravity. However, it is expected that the matter which composes the object obeys the energy conditions. Therefore, in the second part of this paper, we establish connections between the epicyclic frequencies, and consequently the structure of circular orbits, with the energy conditions.

The first study of its kind was done in~\cite{10.1093/mnras/stt1185}, where the authors found that the sum of the squared epicyclic frequencies of TCOs around a Maclaurin spheroid~\cite{maclaurin1742treatise,chandrasekhar1969ellipsoidal} in the Newtonian regime, could be written only in terms of the angular velocity of the timelike particles, $\Omega$, and the density of the spheroid, $\rho$, as~\cite{10.1093/mnras/stt1185},
\begin{equation}\label{Eq:SumEpicyclicFreqNewtonian}
	\omega_r^2 + \omega_\theta^2 = 2 \Omega^2 + 4\pi G \rho~,
\end{equation}
where $\omega_r$ and $\omega_\theta$ are the radial and vertical epicyclic frequencies and $G$ is Newton's constant.
In the same work, the authors also found that, for the case of a Kerr BH, the sum vanishes when computed at the LR.

Some years later, a proof developed in~\cite{Vieira:2017zau} showed that, for a very generic static, spherically symmetric and asymptotically flat object, the sum of the squared epicyclic frequencies of TCOs could be written at the expense of a linear combination of the components of the Ricci tensor~\cite{Vieira:2017zau},
\begin{equation}\label{Eq:VieiraResult}
	\omega_r^2 + \omega_\theta^2 = \left[ R_{tt} + \Omega^2 R_{\varphi\varphi} \right] + \frac{\Omega^2}{2g_{rr}}  \frac{d g_{\varphi\varphi}}{d r} \frac{1}{\tilde{r}^2} \frac{d\tilde{r}^2}{dr^2}~,
\end{equation}
where, $\tilde{r}^2 = -g_{\varphi\varphi}/g_{tt}$. They also found, through the above equation, that if the object obeys the strong energy condition (SEC) then $\omega_r^2 + \omega_\theta^2 > 0$.

In the second part of this paper, we shall generalise the work done in~\cite{Vieira:2017zau} by considering a more generic stationary, axisymmetric ultracompact object. We shall obtain a similar result to the one presented in Eq. \eqref{Eq:VieiraResult} and we shall also consider the SEC as well as the null energy condition (NEC). If the matter which compose the object in question obeys the SEC and NEC, then the sum of the epicyclic frequencies of LRs and TCOs will always be non negative.

This work is organised as follows.
In Section 2, we start by introduction and defining all quantities of interest for the null and timelike circular orbits, namely, their radial and vertical epicyclic frequencies.
Then, in Section 3, we show, together with the work done in \cite{Delgado:2021jxd}, that the epicyclic frequencies of the timelike particle coincide with the epicyclic frequencies of the null particles when both of them are computed at the LRs. The four possible structure of TCOs surrounding a LR are presented.
Section 4 provides the relation between the sum of the squared epicyclic frequencies and the Ricci tensor. Connections with the null and strong energy condition are also provided and analysed.
Finally, in Section 5, we close the present paper with the conclusions and final remarks.

\section{Epicyclic frequencies of equatorial circular orbits}
\label{Sec:Definitions}

To study the epicyclic frequencies, we shall follow a similar setup as the one described in \cite{Delgado:2021jxd}. Let us consider a stationary, axisymmetric, asymptotically flat, 1+3 dimensional spacetime, $(\mathcal{M}, g)$, that describes a generic ultracompact object with a $\mathbb{Z}_2$ symmetry that may or may not have an event horizon. No assumption is made on the field equations $(\mathcal{M},g)$ solves.

Stationarity and axial symmetry imply the existence of two Killing vectors field, $\{\eta_1, \eta_2\}$, that commune, $[\eta_1, \eta_2] = 0$, thanks to a theorem develop by Carter (together with asymptotic flatness) \cite{Carter:1970ea}.
Such result open the possibility to choose a coordinate system, $(t,r,\theta,\varphi)$, that can be adapted to the Killing vector fields such that $\eta_1 = \partial_t$ and $\eta_2 = \partial_\varphi$.
We also assume that the metric associated to the spacetime is, at least, $C^2$-smooth on and outside the possible horizon, and circular. This implies, for asymptotically flat spacetimes, that the geometry possesses a 2-space orthogonal to the Killing vector fields (\textit{c.f.} theorem 7.11 in \cite{Wald:1984rg}). Therefore, the discrete symmetry $(t,\varphi) \rightarrow (-t,-\varphi)$ is present on the spacetime.

By a gauge choice, we can defined spherical-like coordinates $(r,\theta)$ in the orthogonal 2-space that are orthogonal to each other. An extra gauge choice can be used to fix the localisation of the horizon at a constant positive radial coordinate, $r = r_H$, for the case of ultracompact objects with a horizon. With these choices, $g_{r\theta} = 0$, $g_{rr} > 0$ and $g_{\theta\theta} > 0$ (outside the possible horizon). One can also impose that $(r,\theta)$ must reduce to the standard spherical coordinates as one approaches spatial infinity, $r \rightarrow \infty$. The range of the coordinates are, $t \in (-\infty, +\infty)$; $r \in (r_H,+\infty)$, if a horizon exists, or $r \in [0, +\infty)$, if no horizon exists; $\theta \in [0, \pi]$; and $\varphi \in [0, 2\pi)$. The rotating axis is located at $\theta = \{0,\pi\}$, and the equatorial plane can be found at $\theta = \pi/2$. Outside of a possible horizon, causality implies that $g_{\varphi\varphi} \geq 0$.

In the end, due to all the assumptions, gauge choices, and using a Lorentzian signature $(-,+,+,+)$, we can write the following line element,
\begin{equation}
	ds^2 = g_{tt}(r,\theta) dt^2 + 2 g_{t\varphi}(r,\theta) dt d\varphi + g_{\varphi\varphi}(r,\theta) d\varphi^2 + g_{rr}(r,\theta) dr^2 + g_{\theta\theta}(r,\theta) d\theta^2 ~.
\end{equation}
Note that we shall consider that the radial coordinate is a faithful measurement of the distance to the ultracompact object.

The motion of test particles in the above geometry, can be described through the following effective Lagrangian,
\begin{equation}
	2 \mathcal{L} = g_{\mu\nu} \dot{x}^\mu \dot{x}^\nu = \epsilon ~,
\end{equation}
where the dot denotes the derivation with respect to an affine parameter, and $\epsilon = \{-1, 0\}$ for timelike and null particles, respectively.
The effective Lagrangian will depend on the spherical-like coordinate, and can be written as follows
\begin{equation}
	2 \mathcal{L} = 	g_{tt}(r,\theta) \dot{t}^2 + 2 g_{t\varphi}(r,\theta) \dot{t} \dot{\varphi} + g_{\varphi\varphi}(r,\theta) \dot{\varphi}^2 + g_{rr}(r,\theta) \dot{r}^2 + g_{\theta\theta}(r,\theta) \dot{\theta}^2 = \epsilon ~.
\end{equation}

The existence of Killing vector fields give raise to constants of motion that we can be introduce into the Lagrangian. Those are the energy, $E$, and angular momentum, $L$, of the test particle,
\begin{equation}
	-E \equiv g_{t\mu} \dot{x}^\mu = g_{tt} \dot{t} + g_{t\varphi} \dot{\varphi} ~, \hspace{10pt} L \equiv g_{\varphi \mu} = g_{\varphi t} \dot{t} + g_{\varphi\varphi} \dot{\varphi} ~.
\end{equation}
The effective Lagrangian reads now,
\begin{equation}\label{Eq:LagrangianAB}
	2 \mathcal{L} = - \frac{A(r,\theta,E,L)}{B(r,\theta)} + g_{rr}(r,\theta) \dot{r}^2 + g_{\theta\theta}(r,\theta) \dot{\theta}^2= \epsilon ~,
\end{equation}
where
\begin{equation}
	A(r,\theta,E,L) = g_{\varphi\varphi}(r,\theta) E^2 + 2 g_{t\varphi}(r,\theta) E L + g_{tt}(r,\theta) L^2 ~, \hspace{10pt} \text{and} \hspace{10pt} B(r,\theta) = g_{t\varphi}(r,\theta)^2 - g_{tt}(r,\theta) g_{\varphi\varphi}(r,\theta) ~.
\end{equation}
Eq. \eqref{Eq:LagrangianAB} suggests the introduction of an effective potential, $V_\epsilon(r,\theta)$, defined in the following way,
\begin{equation}
	V_\epsilon(r,\theta) \equiv g_{rr}(r,\theta) \dot{r}^2 + g_{\theta\theta}(r,\theta) \dot{\theta}^2= \epsilon + \frac{A(r,\theta,E,L)}{B(r,\theta)} ~.
\end{equation}

In this work, we are interested in objects with a $\mathbb{Z}_2$ symmetry and in circular orbits in the equatorial plane, $\theta = \pi/2$, meaning that $\dot{\theta} = 0$. With these additional assumptions, the effective potential will only have a radial dependency, and can be written as,
\begin{equation}
	V^r_\epsilon(r) \equiv g_{rr}(r,\pi/2) \dot{r}^2 = \epsilon + \frac{A(r,\pi/2,E,L)}{B(r,\pi/2)} ~.
\end{equation}

An equatorial circular orbit can be easily computed through this new effective potential, by imposing that the following two equations must be true simultaneously,
\begin{equation}\label{Eq:FirstEquationPotential}
	V^r_\epsilon(r_\text{cir}) = 0 \hspace{20pt} \Leftrightarrow \hspace{20pt} A(r_\text{cir},\pi/2,E,L) = - \epsilon B(r_\text{cir},\pi/2) ~,
\end{equation}
and
\begin{equation}\label{Eq:SecondEquationPotential}
	\partial_r V^r_\epsilon(r_\text{cir}) = 0 \hspace{20pt} \Leftrightarrow \hspace{20pt} \partial_r A(r_\text{cir},\pi/2,E,L) = -\epsilon \partial_r B(r_\text{cir},\pi/2) ~.
\end{equation}
Note that we have used Eq. \eqref{Eq:FirstEquationPotential}, to obtain the Eq. \eqref{Eq:SecondEquationPotential}.

Beside knowing where and how to find circular orbits, it is equally important to know their stability. This can be analysed by computing both their radial and vertical epicyclic frequencies.

\subsection{Radial and vertical epicyclic frequencies}

The epicyclic frequencies can be computed by perturbing a circular orbit in either the radial or vertical direction. Through the former, we compute the radial epicyclic frequency, whereas through the latter, we obtain the vertical epicyclic frequency.

Let us start by consider that $x$ is one of the spheroidal coordinates, $x = \{r,\theta\}$, and $y$ is the remaining spheroidal coordinate. Moreover, let us also consider that we are on a circular orbit, such that $x = x_\text{c}$ and $y = y_\text{c}$. We shall fix $y$ and introduce a perturbation in $x$ around $x_\text{c}$, such that $y = y_\text{c} \Leftrightarrow \dot{y} = 0$ and $x = x_\text{c} + \delta x  \Leftrightarrow \dot{x} = \dot{\delta x}$. Under these assumptions, the perturbed effective potential reads,
\begin{equation}\label{Eq:PerturbedPotential}
	V^x_\epsilon(x_\text{c} + \delta x,y_\text{c}) = g_{xx}(x_\text{c} + \delta x,y_\text{c}) \dot{\delta x}^2 ~.
\end{equation}
The left-hand side of the above equation can be expanded to,
\begin{equation}
	V^x_\epsilon(x_\text{c} + \delta x,y_\text{c}) = V^x_\epsilon(x_\text{c},y_\text{c}) + \partial_x V^x_\epsilon(x_\text{c},y_\text{c}) \delta x + \frac{1}{2} \partial_x^2 V^x_\epsilon(x_\text{c},y_\text{c}) \delta x^2 + \mathcal{O}(\delta x^3) ~.
\end{equation}
The first two terms will always vanish because we are both on a circular orbit (if $x=r: V^r_\epsilon(r_\text{cir}) = \partial_r V^r_\epsilon(r_\text{cir}) = 0$) and on the equatorial plane of a $\mathbb{Z}_2$ symmetric object (if $x=\theta: V^\theta_\epsilon(\pi/2) = \partial_\theta V^\theta_\epsilon(\pi/2) = 0$). Thus, equation \eqref{Eq:PerturbedPotential} becomes\footnote{Note that, for simplicity and notation ease, we dropped the error term, $\mathcal{O}(\delta x^3)$.},
\begin{equation}
	\frac{1}{2} \partial^2_x V^x_\epsilon(x_\text{c},y_\text{c}) \delta x^2 = g_{xx}(x_\text{c} + \delta x,y_\text{c}) \dot{\delta x}^2 ~.
\end{equation}
Expanding now the right-hand size of the above equation, and considering only terms up to second order in $\dot{\delta x}$, we arrive at,
\begin{equation}
	\frac{1}{2} \partial^2_x V^x_\epsilon(x_\text{c},y_\text{c}) \delta x^2 = g_{xx}(x_\text{c},y_\text{c}) \dot{\delta x}^2 ~.
\end{equation}
By performing a dot derivative of both sizes, we can write a more instructive result,
\begin{equation}
	\ddot{\delta x} + (\omega^x_\epsilon)^2 \delta x = 0  ~,\hspace{20pt} \text{where} \hspace{20pt} (\omega^x_\epsilon)^2 \equiv -\frac{1}{2} \frac{\partial_x^2 V^x_\epsilon(x_\text{c},y_\text{c})}{g_{xx}(x_\text{c},y_\text{c})}
\end{equation}
This is an harmonic oscillator, which means that if the frequency $\omega_x^2$ is negative (positive), the perturbation grows exponentially (remains a small perturbation) leading to unstable (stable) circular orbits. This frequency is known as the epicyclic frequency.
Going back to the spheroidal coordinates, we have the radial and vertical epicyclic frequencies, respectively,
\begin{equation}\label{Eq:EpicyclicFrequenciesComovingObs}
	(\omega^r_\epsilon)^2 \equiv -\frac{1}{2} \frac{\partial_r^2 V^r_\epsilon(r_\text{cir},\pi/2)}{g_{rr}(r_\text{cir},\pi/2)} ~, \hspace{20pt} (\omega^\theta_\epsilon)^2 \equiv -\frac{1}{2} \frac{\partial_\theta^2 V^\theta_\epsilon(r_\text{cir},\pi/2)}{g_{\theta\theta}(r_\text{cir},\pi/2)} ~.
\end{equation}
These epicyclic frequencies are measured with respect to the proper time of a comoving observer. However, in the following computations, it proved useful to use the epicyclic frequencies measured by an observer at spatial infinity, rather than by a comoving observer.
To obtain the frequencies measured by such observer, one only has to dividing the frequencies in Eq. \eqref{Eq:EpicyclicFrequenciesComovingObs} by the squared redshift factor, $\dot{t}^2$,
\begin{equation}\label{Eq:EpicyclicFrequencyRadial}
	(\nu^r_\epsilon)^2 \equiv \frac{(\omega^r_\epsilon)^2}{\dot{t}^2} = -\frac{1}{2} \frac{\partial_r^2 V^r_\epsilon(r_\text{cir},\pi/2)}{g_{rr}(r_\text{cir},\pi/2)} \left[ \frac{B(r_\text{cir},\pi/2)}{E g_{\varphi \varphi}(r_\text{cir},\pi/2) + L g_{t\varphi}(r_\text{cir},\pi/2) }\right]^2 ~,~
\end{equation}
\begin{equation}\label{Eq:EpicyclicFrequencyAngular}
	(\nu^\theta_\epsilon)^2 \equiv \frac{(\omega^\theta_\epsilon)^2}{\dot{t}^2} = -\frac{1}{2} \frac{\partial_\theta^2 V^\theta_\epsilon(r_\text{cir},\pi/2)}{g_{\theta\theta}(r_\text{cir},\pi/2)} \left[ \frac{B(r_\text{cir},\pi/2)}{E g_{\varphi \varphi}(r_\text{cir},\pi/2) + L g_{t\varphi}(r_\text{cir},\pi/2) }\right]^2 ~.
\end{equation}
Henceforth, we shall drop the explicit dependency of the several functions, and it shall be understood that all quantities are computed at $r=r_\text{cir}$ and $\theta = \pi/2$, unless stated otherwise.

We shall now compute the epicyclic frequencies for null and timelike circular orbits.

\subsection{Null particles}

For null particles, $\epsilon = 0$, circular orbits are known as LRs.
In order to obtain them, we will use Eqs. \eqref{Eq:FirstEquationPotential} and \eqref{Eq:SecondEquationPotential}.
Both equations can be rewritten in terms of the inverse impact parameter, $\sigma_\pm = E_\pm/L_\pm$, where $\pm$ represents the two possible solutions due to the rotation of the ultracompact object and they are associated with prograde ($+$) and retrograde ($-$) orbits,
\begin{eqnarray}
	A = 0 \hspace{10pt} &\Leftrightarrow& \hspace{10pt} \left[ g_{\varphi\varphi} \sigma_\pm^2 + 2 g_{t\varphi} \sigma_\pm + g_{tt} \right]_\text{LR} = 0~,\label{Eq:LRFirstEquationPotential} \\
	\partial_r A = 0 \hspace{10pt} &\Leftrightarrow& \hspace{10pt} \left[\partial_r g_{\varphi\varphi} \sigma_\pm^2 + 2 \partial_r g_{t\varphi} \sigma_\pm + \partial_r g_{tt} \right]_\text{LR} = 0 ~.\label{Eq:LRSecondEquationPotential}
\end{eqnarray}
The first equation gives an algebraical equation for the inverse impact parameter,
\begin{equation}
	\sigma_\pm = \left[ \frac{-g_{t\varphi} + \sqrt{B}}{g_{\varphi\varphi}} \right]_\text{LR}~,
\end{equation}
whereas the second equation gives the radial coordinate of the LR, $r = r_\text{LR}$.

The epicyclic frequencies of the LRs, can be easily computed through Eqs. \eqref{Eq:EpicyclicFrequencyRadial} and \eqref{Eq:EpicyclicFrequencyAngular},
\begin{equation}
	(\nu^r_{0})^2 = - \frac{1}{2} \left[ \frac{\partial_r^2 g_{\varphi\varphi} \sigma_\pm^2 + 2 \partial_r^2 g_{t\varphi} \sigma_\pm + \partial_r^2 g_{tt}}{g_{rr}} \right]_\text{LR}
\end{equation}
\begin{equation}
	(\nu^\theta_{0})^2 = - \frac{1}{2} \left[ \frac{\partial_\theta^2 g_{\varphi\varphi} \sigma_\pm^2 + 2 \partial_\theta^2 g_{t\varphi} \sigma_\pm + \partial_\theta^2 g_{tt}}{g_{\theta\theta}} \right]_\text{LR}
\end{equation}
Since $g_{rr}$ and $g_{\theta\theta}$ are always positive outside of a possible horizon, the epicyclic frequencies are real (complex) if the their numerator are negative (positive).

\subsection{Timelike particles}

For timelike particles, $\epsilon = -1$, the same analysis can be done. In order to simplify the following computations, it is convenient to introduce the angular velocity of timelike particles along circular orbits (measured with respect to an observer at infinity),
\begin{equation}
	\Omega = \frac{d\varphi}{dt} = - \frac{E g_{t\varphi} + L g_{tt}}{E g_{\varphi\varphi} + L g_{t\varphi}} ~.
\end{equation}
With this result and Eq. \eqref{Eq:FirstEquationPotential}, one can obtain an analytical expression for the energy and angular momentum of the timelike particle written in terms of the metric functions and the angular velocity,
\begin{equation}\label{Eq:EnergyAngMomTCOs}
	E_\pm = - \left[ \frac{g_{tt} + g_{t\varphi}\Omega_\pm}{\sqrt{\beta_\pm}} \right]_{r_\text{cir}}~,\hspace{10pt} L = \left[ \frac{g_{t\varphi} + g_{\varphi\varphi}\Omega_\pm}{\sqrt{\beta_\pm}} \right]_{r_\text{cir}} ~,
\end{equation}
where $\beta_\pm \equiv -g_{tt} - 2 g_{t\varphi} \Omega_\pm - g_{\varphi\varphi} \Omega_\pm^2$.
An expression for the angular velocity can be obtained by solving equation \eqref{Eq:SecondEquationPotential},
\begin{equation}\label{Eq:AngularVelTCOs}
	\Omega_\pm = \left[ \frac{-\partial_r g_{t\varphi} \pm \sqrt{C}}{\partial_r g_{\varphi\varphi}} \right]_{r_\text{cir}} ~,
\end{equation}
where $C \equiv (\partial_r g_{t\varphi})^2 - \partial_r g_{tt} \partial_r g_{\varphi\varphi}$.

With Eqs. \eqref{Eq:EnergyAngMomTCOs} and \eqref{Eq:AngularVelTCOs}, we have defined all possible timelike circular orbits (TCOs) in a given spacetime (with the assumptions presented in the beginning of this section). Their stability can be analysed through their epicyclic frequencies,
\begin{equation}
	(\nu^r_{-1})^2 = - \frac{1}{2} \left[ \frac{(\partial_r^2 g_{\varphi\varphi} \Omega_\pm^2 + 2 \partial_r^2 g_{t\varphi} \Omega_\pm + \partial_r^2 g_{tt})B - 2 C \beta_\pm}{B g_{rr}} \right]_{r_\text{cir}} ~,
\end{equation}
\begin{equation}
	(\nu^\theta_{-1})^2 = - \frac{1}{2} \left[ \frac{\partial_\theta^2 g_{\varphi\varphi} \Omega_\pm^2 + 2 \partial_\theta^2 g_{t\varphi} \Omega_\pm + \partial_\theta^2 g_{tt}}{g_{\theta\theta}} \right]_{r_\text{cir}} ~.
\end{equation}
Similar as to the case of LRs, the epicyclic frequencies of TCOs are real (complex) is the numerator of the above expressions are negative (positive), since $B$, $g_{rr}$ and $g_{\theta\theta}$ are always positive outside a possible horizon.

\section{Connection between null and timelike particles}

We shall show now that the epicyclic frequencies of TCOs coincide with the epicyclic frequencies of LRs when we compute the former on a LR.

We start by recalling the results derived in \cite{Delgado:2021jxd}. Firstly, the authors showed that the function $\beta_\pm$ always vanishes on a LR, by proving that Eqs. \eqref{Eq:LRFirstEquationPotential} and \eqref{Eq:LRSecondEquationPotential} identically vanish when $\beta_\pm = 0$. A consequence of this result is that the angular velocity of a timelike particle on such a TCO is the same as the inverse impact parameter of a null particle on the LR, $[\Omega_\pm = \sigma_\pm]_\text{LR}$.
Secondly, they showed that the localisation and \textit{radial} stability of the regions that can harbour TCOs around the LR, depend exclusively on the \textit{radial} stability of said LR. Such statement is verified by looking at the sign of $\beta_\pm$ to find the regions where TCOs are allowed, and by studying the stability of those allowed regions. In short, a radially stable (unstable) LR accommodates radially stable (unstable) TCOs in the regions radially below (above) the LR.
A schematic representation of these results is presented in Fig. \ref{Fig:CompIllustStableLR}.

\begin{figure}[h!]
	\centering
	\begin{tikzpicture}
		\draw[ultra thick,->] (-0.05,-0.05) -- (10.45,-0.05) node[anchor=north west]{$r$};
		\filldraw[red!10, draw=red!85!black,ultra thick] (0,0) rectangle (5,1.5) node[pos=0.5,red!85!black,font=\large\bfseries]{TCOs forbidden};
		\draw[dash pattern=on 3pt off 4pt,ultra thick] (5.1,0) -- (5.1,1.5);
		\filldraw[yellow!10, draw=yellow!20!orange,ultra thick] (5.2,0) rectangle (10.2,1.5) node[pos=0.5,yellow!20!orange,align=center,font=\large\bfseries]{\textit{Radially}\\ Unstable TCOs};
		\node[font=\bfseries,align=center] at (5.1,1.95){\textit{Radially}\\ Unstable LR};
		\node at (5.1,-0.35){$r_\pm^{\text{LR}}$};
	\end{tikzpicture}
	\begin{tikzpicture}
		\draw[ultra thick,->] (-0.05,-0.05) -- (10.45,-0.05) node[anchor=north west]{$r$};
		\filldraw[green!10, draw=green!70!black,ultra thick] (0,0) rectangle (5,1.5) node[pos=0.5,green!70!black,align=center,font=\large\bfseries]{\textit{Radially}\\ Stable TCOs};;
		\draw[dash pattern=on 3pt off 4pt,ultra thick] (5.1,0) -- (5.1,1.5);
		\filldraw[red!10, draw=red!85!black,ultra thick] (5.2,0) rectangle (10.2,1.5) node[pos=0.5,red!85!black,font=\large\bfseries]{TCOs forbidden};
		\node[font=\bfseries,align=center] at (5.1,1.95){\textit{Radially}\\ Stable LR};
		\node at (5.1,-0.35){$r_\pm^{\text{LR}}$};
	\end{tikzpicture}
	\caption{Structure of the equatorial TCOs in the immediate vicinity of an unstable (top panel) and stable (bottom panel) LR. Adapted from \cite{Delgado:2021jxd}.}
	\label{Fig:CompIllustStableLR}
\end{figure}

These results can now pave the way for us to generalise them, by also including the vertical stability of the allowed TCOs.  For that, let us consider a spacetime with a LR. If we compute the epicyclic frequencies of a TCO infinitely close to the LR, such that we can use $\beta_\pm|_\text{LR} = 0$ and $[\Omega_\pm = \sigma_\pm]_\text{LR}$, we obtain,

\begin{equation}\label{Eq:EpicyclicFreqRadialLR}
	(\nu^r_{-1})^2 = - \frac{1}{2} \left[ \frac{\partial_r^2 g_{\varphi\varphi} \sigma_\pm^2 + 2 \partial_r^2 g_{t\varphi} \sigma_\pm + \partial_r^2 g_{tt}}{g_{rr}} \right]_\text{LR} = (\nu^r_{0})^2 ~,
\end{equation}
\begin{equation}\label{Eq:EpicyclicFreqAngularLR}
	(\nu^\theta_{-1})^2 = - \frac{1}{2} \left[ \frac{\partial_\theta^2 g_{\varphi\varphi} \sigma_\pm^2 + 2 \partial_\theta^2 g_{t\varphi} \sigma_\pm + \partial_\theta^2 g_{tt}}{g_{\theta\theta}} \right]_\text{LR} = (\nu^\theta_{0})^2 ~.
\end{equation}
Here we see that the epicyclic frequencies, and hence the stability, of a TCO infinitely close to a LR are precisely the same as the ones for the same LR.
Therefore, the full stability (both radial and vertical) of a LR determines the full stability of a TCO infinitely close to that LR.

This result, however, it is only true infinitely close to the LR, where $\beta_\pm|_\text{LR} = 0$ and $[\Omega_\pm = \sigma_\pm]_\text{LR}$ are valid. To analyse the stability of the allowed TCOs in the vicinity of the LR, where $\beta_\pm|_\text{LR} = 0$ and $[\Omega_\pm = \sigma_\pm]_\text{LR}$ are no longer valid, we need to investigate what happens to the epicyclic frequencies on those regions.
Fortunately, this analyse is rather simple to do. One only has to use the continuity properties of the epicyclic frequencies.

Imagine a spacetime with a radial unstable LR. By Eq. \eqref{Eq:EpicyclicFreqRadialLR} we know that a timelike particle on a TCO infinitely close to the LR is radially unstable, hence, $\nu^r_{-1}(r_\text{LR})^2 < 0$. If we analyse the epicyclic frequency of an allowed TCO adjacent to the LR, one can argue that the epicyclic frequency will continue to be negative (however one does not know if it increases or decreases). Such argument can be done by assuming that the epicyclic frequencies are continuous on the allowed regions for TCOs, which is a fair assumption to make since we assumed that all metric functions are, at least, $C^2$-smooth and all remaining quantities are well-behaved. Hence, the allowed region of TCOs adjacent to the radially unstable LR will harbour radially unstable TCOs. Such result is consistent with the results provided in \cite{Delgado:2021jxd}.

The same argument can be done regarding a radially stable LR, recovering the same results found in \cite{Delgado:2021jxd}, or regarding a vertically stable and unstable LR. In the last two cases, the allowed region of TCOs continues to the determined by the \textit{radial} stability of the LR, but their vertical stability is entirely determined by the vertical stability of the LR. We can summarise all the main results in the following list:
\begin{itemize}
	\item Spacetime with a \textit{radially} unstable LR.

	For such spacetime, the region immediately inwards of the LR can not accommodate TCOs, whereas the outward region can. To know the full stability of the allowed TCOs we need to know the full stability of the LR, therefore, let us assume that the LR is,

	\begin{itemize}
		\item \textit{Vertically} unstable.

			In this case, the allowed TCOs will be both \textit{radially} and \textit{vertically} unstable.

			\begin{figure}[H]
				\centering
				\begin{tikzpicture}
					\draw[ultra thick,->] (-0.05,-0.05) -- (11.85,-0.05) node[anchor=north west]{$r$};
					\filldraw[red!10, draw=red!85!black,ultra thick] (0,0) rectangle (5.7,1.5) node[pos=0.5,red!85!black,font=\large\bfseries]{TCOs forbidden};
					\draw[dash pattern=on 3pt off 4pt,ultra thick] (5.8,0) -- (5.8,1.5);
					\filldraw[yellow!10, draw=yellow!20!orange,ultra thick] (5.9,0) rectangle (11.6,1.5) node[pos=0.5,yellow!20!orange,align=center,font=\large\bfseries]{\textit{Radially} and \textit{Vertically}\\ Unstable TCOs};

					\node at (5.8,-0.35){$r_\pm^{\text{LR}}$};
				\end{tikzpicture}
				\caption{Structure of the equatorial TCOs in the immediately vicinity of a \textit{radially} and \textit{vertically} unstable LR.}
				\label{Fig:StructureFullyUnstable}
			\end{figure}

		\item \textit{Vertically} stable.

			Here, the allowed TCOs will be \textit{radially} unstable but \textit{vertically} stable.

			\begin{figure}[h!]
				\centering
				\begin{tikzpicture}
					\draw[ultra thick,->] (-0.05,-0.05) -- (11.85,-0.05) node[anchor=north west]{$r$};
					\filldraw[red!10, draw=red!85!black,ultra thick] (0,0) rectangle (5.7,1.5) node[pos=0.5,red!85!black,font=\large\bfseries]{TCOs forbidden};
					\draw[dash pattern=on 3pt off 4pt,ultra thick] (5.8,0) -- (5.8,1.5);
					\filldraw[yellow!30!blue!10, draw=yellow!30!blue!70!blue,ultra thick] (5.9,0) rectangle (11.6,1.5) node[pos=0.5,yellow!30!blue!70!blue,align=center,font=\large\bfseries]{\textit{Radially} Unstable and \\ \textit{Vertically} Stable TCOs};

					\node at (5.8,-0.35){$r_\pm^{\text{LR}}$};
				\end{tikzpicture}
				\caption{Structure of the equatorial TCOs in the immediately vicinity of a \textit{radially} unstable and \textit{vertically} stable LR.}
				\label{Fig:StructureRadiallyUnstableAngularlyStable}
			\end{figure}
	\end{itemize}

	\item Spacetime with a \textit{radially} stable LR.

	For such spacetime, the region immediately outwards of the LR can not accommodate TCOs, whereas the inward region can. Similar as the previous case, the full stability of the allowed region is determined by the full stability of the LR. Let us assume that the LR is,

	\begin{itemize}
		\item \textit{Vertically} unstable.

			For this case, the allowed TCOs with be \textit{radially} stable but \textit{vertically} unstable.
			\begin{figure}[h!]
				\centering
				\begin{tikzpicture}
					\draw[ultra thick,->] (-0.05,-0.05) -- (11.85,-0.05) node[anchor=north west]{$r$};
					\filldraw[green!30!blue!10, draw=green!30!blue!70,ultra thick] (0,0) rectangle (5.7,1.5) node[pos=0.5,green!30!blue!70,align=center,font=\large\bfseries]{\textit{Radially} Stable and\\ \textit{Vertically} Unstable TCOs};;
					\draw[dash pattern=on 3pt off 4pt,ultra thick] (5.8,0) -- (5.8,1.5);
					\filldraw[red!10, draw=red!85!black,ultra thick] (5.9,0) rectangle (11.6,1.5) node[pos=0.5,red!85!black,font=\large\bfseries]{TCOs forbidden};

					\node at (5.8,-0.35){$r_\pm^{\text{LR}}$};
				\end{tikzpicture}
				\caption{Structure of the equatorial TCOs in the immediately vicinity of a \textit{radially} stable and \textit{vertically} unstable LR.}
				\label{Fig:StructureRadiallyStableAngularlyUnstable}
			\end{figure}

		\item \textit{Vertically} stable.

			In this particular case, the allowed TCOs will be both \textit{radially} and \textit{vertically} stable.

			\begin{figure}[h!]
				\centering
				\begin{tikzpicture}
					\draw[ultra thick,->] (-0.05,-0.05) -- (11.85,-0.05) node[anchor=north west]{$r$};
					\filldraw[green!10, draw=green!70!black,ultra thick] (0,0) rectangle (5.7,1.5) node[pos=0.5,green!70!black,align=center,font=\large\bfseries]{\textit{Radially} and \textit{Vertically}\\ Stable TCOs};;
					\draw[dash pattern=on 3pt off 4pt,ultra thick] (5.8,0) -- (5.8,1.5);
					\filldraw[red!10, draw=red!85!black,ultra thick] (5.9,0) rectangle (11.6,1.5) node[pos=0.5,red!85!black,font=\large\bfseries]{TCOs forbidden};

					\node at (5.8,-0.35){$r_\pm^{\text{LR}}$};
				\end{tikzpicture}
				\caption{Structure of the equatorial TCOs in the immediately vicinity of a \textit{radially} and \textit{vertically} stable LR.}
				\label{Fig:StructureFullyStable}
			\end{figure}
	\end{itemize}
\end{itemize}

\section{Epicyclic frequencies and energy conditions}

We have showed, so far, that the structures of TCOs around a LR for a generic stationary, axisymmetric and asymptotically flat object with a $\mathbb{Z}_2$ symmetry can be characterised in 4 different ways. This, however, does not assume any type of constrain regarding the matter composing the object. If one wants to consider a object that is compose of real physical matter, then one assumes that the matter obeys the energy conditions.

The energy conditions that we will take special attention to are the NEC and SEC.
Regarding the former, a proof given in \cite{Cunha:2017qtt} showed that an ultracompact horizonless object, that obeys the NEC, can never have a fully unstable LR. Likewise, following \cite{Cunha:2020azh}, a BH can not have fully unstable LRs if its matter obeys the NEC. Here, we shall arrive to the same result using a different approach than the one used in \cite{Cunha:2017qtt,Cunha:2020azh}.
Regarding the latter, we shall show a similar conclusion as the one stated previously but related with TCOs. In particular, we shall present a proof that the existence of a region with fully unstable TCOs implies the violation of the SEC.
We will also look into the more simple case of a Ricci-flat object, where we will show that such object can never have neither fully unstable nor stable LRs, and, consequently, can never have neither fully stable or unstable TCOs close to the LRs. However, fully stable TCOs can still exist in regions far away from the LRs, such as, \textit{e.g.}, in the asymptotically flat region.

All previous statement can be demonstrated by relating the sum of the squared of the epicyclic frequencies and the Ricci tensor.

Following a similar work developed in~\cite{Vieira:2017zau}, where the authors constructed a relation between the epicyclic frequencies and the Ricci tensor, for a generic static, spherically symmetric and asymptotically flat spacetime, we need to find a linear combination of the Ricci tensor components such that we can write it in terms of the epicyclic frequencies.
From the connection with the Newtonian limit \cite{Abramowicz:2004tm,binney2011galactic}, we expect that $R_{tt}$ must be connected to the epicyclic frequencies, since, in this limit\cite{Vieira:2017zau,Abramowicz:2004tm,binney2011galactic}, $(\nu_{-1}^r)^2 + (\nu_{-1}^\theta)^2 = 2\Omega_\pm + R_{tt}$, where $R_{tt} = 4\pi G \rho$ -- \textit{cf.} Eq. \eqref{Eq:SumEpicyclicFreqNewtonian}. Furthermore, since the epicyclic frequencies do not depend on the derivatives of $g_{rr}$ and $g_{\theta\theta}$, the linear combination must account that fact.

The linear combination of Ricci tensor components that satisfy the requirements presented above must be of the form of $R_{tt} + 2 f R_{t\varphi} + f^2 R_{\varphi\varphi}$, where $f$ is either the inverse impact parameter of null particles, $\sigma_\pm$, or the angular velocity of timelike particles, $\Omega_\pm$.
This way, the correct terms to write the epicyclic frequencies appear, and all the terms with derivatives of $g_{rr}$ and $g_{\theta\theta}$ disappear.

Let us focus in the particular case of TCOs. In this case, $f=\Omega_\pm$, and we start by expanding the linear combination,
\begin{eqnarray}\label{Eq:RicciCombination}
	R_{tt} + 2\Omega_\pm R_{t\varphi} + \Omega_\pm^2 R_{\varphi\varphi} &=& \left(\partial_\mu \Gamma^\mu_{tt} + 2 \Omega_\pm \partial_\mu \Gamma^\mu_{t\varphi} + \Omega_\pm^2 \Gamma^\mu_{\varphi\varphi} \right) + \nonumber \\
	&& \Gamma^\mu_{\mu\nu} \left( \Gamma^\nu_{tt} + 2 \Omega_\pm \Gamma^\nu_{t\varphi} + \Omega_\pm^2 \Gamma^\nu_{\varphi\varphi} \right) - \nonumber \\
	&& \left( \Gamma^\mu_{\nu t} \Gamma^\nu_{\mu t} + 2 \Omega_\pm \Gamma^\mu_{\nu t} \Gamma^\nu_{\mu \varphi} + \Omega_\pm^2 \Gamma^\mu_{\nu \varphi} \Gamma^\nu_{\mu \varphi} \right) ~,
\end{eqnarray}
where $\Gamma^\alpha_{\mu\nu} = \frac{1}{2} g^{\alpha \rho}\left(\partial_\mu g_{\rho \nu} + \partial_\nu g_{\mu \rho} - \partial_\rho g_{\mu\nu} \right)$ are the Christoffel symbols.
Since we are computing this quantities in a circular orbit at the equatorial plane, we can use the relations discussed in previous sections. This way, the first term (first line) in Eq.  \eqref{Eq:RicciCombination} reduces to,
\begin{equation}
	\partial_\mu \Gamma^\mu_{tt} + 2 \Omega_\pm \partial_\mu \Gamma^\mu_{t\varphi} + \Omega_\pm^2 \Gamma^\mu_{\varphi\varphi} = (\nu_{-1}^r)^2 + (\nu_{-1}^\theta)^2 - \frac{1}{g_{rr}} \frac{C \beta_\pm}{B} ~.
\end{equation}
Likewise, the second term (second line) vanishes, whereas the third and final term (third line) reduces to,
\begin{equation}
	\Gamma^\mu_{\nu t} \Gamma^\nu_{\mu t} + 2 \Omega_\pm \Gamma^\mu_{\nu t} \Gamma^\nu_{\mu \varphi} + \Omega_\pm^2 \Gamma^\mu_{\nu \varphi} \Gamma^\nu_{\mu \varphi} = -\frac{1}{2 g_{rr}} \frac{C \beta_\pm}{B} ~.
\end{equation}
Therefore, we can write the sum of the squared epicyclic frequencies through the linear combination of Ricci tensor components in the following way,
\begin{equation}\label{Eq:SumEpicyclicFrequenciesTCOs}
	(\nu_{-1}^r)^2 + (\nu_{-1}^\theta)^2 = R_{tt} + 2\Omega_\pm R_{t\varphi} + \Omega_\pm^2 R_{\varphi\varphi} + \frac{1}{2g_{rr}} \frac{C \beta_\pm}{B}~.
\end{equation}
We have successfully generalised the result presented in \cite{Vieira:2017zau} for the more general case of a stationary and axisymmetric spacetime. In the limit of static and spherically symmetric spacetimes our result converges to the one obtained in \cite{Vieira:2017zau}.

The same procedure can be done for the case of LRs. In this case, $f = \sigma_\pm$ and, after expanding and simplifying the linear combination, one can obtain a similar expression as seen in Eq. \eqref{Eq:SumEpicyclicFrequenciesTCOs},
\begin{equation}\label{Eq:SumEpicyclicFrequenciesLRs}
	(\nu_0^r)^2 + (\nu_0^\theta)^2 = R_{tt} + 2\sigma_\pm R_{t\varphi} + \sigma_\pm^2 R_{\varphi\varphi}~.
\end{equation}
Thus, for null particles, the sum of the squared epicyclic frequencies can be entirely written in terms of the Ricci tensor components and the inverse impact parameter.
Note that the same result can be obtained by simply using some results from the previous section, in particular, $\beta_\pm = 0$, $\sigma_\pm = \Omega_\pm$, $\nu_{-1}^r = \nu_0^r$ and $\nu_{-1}^\theta = \nu_0^\theta$, and using them in Eq. \eqref{Eq:SumEpicyclicFrequenciesTCOs}.

We can now investigate what happen to Eqs. \eqref{Eq:SumEpicyclicFrequenciesTCOs} and \eqref{Eq:SumEpicyclicFrequenciesLRs} when we impose some conditions to the Ricci tensor.

\subsection{Ricci-flat object}

For a Ricci-flat object, all components of the Ricci tensor vanish, $R_{\mu\nu} = 0$. For such objects, the sum of the squared epicyclic frequencies of null and timelike particles simplifies to,
\begin{eqnarray}
	\text{LR:}&~ (\nu_0^r)^2 + (\nu_0^\theta)^2 = 0~, \label{Eq:RicciFlatLRs} \\
	\text{TCOs:}&~ (\nu_{-1}^r)^2 + (\nu_{-1}^\theta)^2 = \dfrac{1}{2g_{rr}} \dfrac{C \beta_\pm}{B}~. \label{Eq:RicciFlatTCOs}
\end{eqnarray}

From the first equation, one can conclude that a LR can never be neither fully stable, $(\nu_0^r)^2 > 0 ~ \land ~ (\nu_0^\theta)^2 > 0$, or unstable, $(\nu_0^r)^2 < 0 ~ \land ~ (\nu_0^\theta)^2 < 0$, since its sum must vanish. This implies that the LR must have opposite radial and vertical stabilities.
This is consistent with what is already known for Ricci-flat solutions that are stationary, axisymmetric, asymptotically flat and in 4 dimensions, such as the Kerr solution~\cite{abramowicz1998theory,abramowicz2013foundations,Torok:2005ct}. Both the prograde and retrograde LRs present in the Kerr spacetime are always radially unstable and vertically stable.

For the case of TCOs, one can prove that the right-hand side of Eq. \eqref{Eq:RicciFlatTCOs} is always positive. The proof rely on the fact that each individual quantity is positive. The $g_{rr}$ component of the metric and the $B$ function are always positive when we are outside of a possible horizon -- \textit{cf.} Sec. \ref{Sec:Definitions}. The function $C$ is always positive on the regions where TCOs are allowed, otherwise the angular velocity of the timelike particle is complex, and we no longer have TCOs -- see \cite{Delgado:2021jxd} for a more detailed analysis. The $\beta_\pm$ function is also always positive on the regions where TCOs are allowed, as we discussed in the previous section. Therefore, we have proved that the right-hand size of Eq. \eqref{Eq:RicciFlatTCOs} is always positive on the regions where we can find TCOs.

A direct consequence of this proof is that TCOs around a Ricci-flat solution can never be fully unstable. However, their stability can be any of the three remaining possibilities combinations. For the well-known case of Kerr BHs, it is known that sufficiently far away from the BH, TCOs are fully stable, whereas, if we approach the BH, we start to find radially unstable but vertically stable TCOs~\cite{abramowicz1998theory,abramowicz2013foundations,Torok:2005ct}. Thus, our result is consistent with the well-known results from the Kerr case.

\subsection{NEC and SEC obeying object}

If an object possesses matter such that its energy-momentum tensor (and, by Einstein's equations, the Ricci tensor) is non-vanishing, a way to inform ourselves about its exotic properties is by evaluating possible violation of the energy conditions.
One of them is the SEC. In short, this condition is defined as $R_{\mu\nu} t^\mu t^\nu \geq 0$, for any timelike vector field, $t^\mu$, and encapsulated the condition that matter must gravitate towards matter.
To study this energy condition, let us consider a timelike vector tangent to a TCO,
\begin{equation}
	t^\mu = a \left( \eta_1^\mu + \Omega_\pm \eta_2^\mu \right)~,
\end{equation}
where $a = a(r_\text{cir},\pi/2) > 0$ is a constant. By computing the SEC, we arrive to,
\begin{equation}
	R_{\mu\nu} t^\mu t^\nu \geq 0 ~~\Leftrightarrow~~ a^2 \left( R_{tt} + 2 \Omega_\pm R_{t\varphi} + \Omega_\pm^2 R_{\varphi\varphi} \right) \geq 0~.
\end{equation}
Therefore, we can written the sum of the squared epicyclic frequencies of a TCO as,
\begin{equation}\label{Eq:SECObeyingTCOs}
	(\nu_{-1}^r)^2 + (\nu_{-1}^\theta)^2 = \frac{1}{a^2} R_{\mu\nu} t^\mu t^\nu + \dfrac{1}{2g_{rr}} \dfrac{C \beta_\pm}{B}~.
\end{equation}
Thus, for a object whose matter obeys the SEC, the first term of the right-hand size of Eq. \eqref{Eq:SECObeyingTCOs} is always positive. Furthermore, using the same argument used in the previous subsection where we prove that the right-hand size of Eq. \eqref{Eq:RicciFlatTCOs} is always positive, we can conclude that the sum of the squared epicyclic frequencies of TCOs around an object obeying the SEC is always positive.
This implies that no fully unstable TCOs exist on the spacetime generated by such object.
Or, in different words, if a generic stationary, axisymmetric, asymptotically flat spacetime has a region with fully unstable TCOs, then we can guarantee that the SEC is violated.

A second energy condition that we can also study is the NEC. This condition is defined as $T_{\mu\nu} k^\mu k^\nu \geq 0$, for any null vector field $k^\mu$, and it is intertwined with the SEC, in the sense that one can imply the other.
In particular, if the SEC is obeyed then the NEC is also obeyed, or, if the NEC is violated, the SEC is also violated.
We will now show that the interlacement between both energy conditions is consistent with all results discussed so far.

Let us consider the following null vector field,
\begin{equation}
	k^\mu = b \left( \eta_1^\mu + \sigma_\pm \eta_2^\mu \right)~,
\end{equation}
where $b = b(r_\text{cir},\pi/2) > 0$. With this null vector, the NEC can be written as\footnote{We are assuming that, within a given theory of gravity, we can defined a effective energy-momentum tensor such that we can use the Einstein's equations, $R_{\mu\nu} - \frac{1}{2}g_{\mu\nu} R = T_{\mu\nu}^\text{eff}$.},
\begin{equation}
	T_{\mu\nu} k^\mu k^\nu \geq 0 ~~\Leftrightarrow~~ b^2 (R_{tt} + 2\sigma_\pm R_{t\varphi} + \sigma_\pm^2 R_{\varphi\varphi}) \geq 0 ~.
\end{equation}
Therefore, the sum of the squared epicyclic frequencies of the LRs take a very simple form,
\begin{equation}\label{Eq:SECObeyingLRs}
	(\nu_0^r)^2 + (\nu_0^\theta)^2 = \frac{1}{b^2} T_{\mu\nu} k^\mu k^\nu~.
\end{equation}
From this result we can conclude that fully unstable LRs never exist for objects whose matter obeys the NEC.
Or, in different words, if a generic stationary, axisymmetric, asymptotically flat spacetime has fully unstable LRs, then we can guarantee that the NEC is violated.
This is precisely the same result that was first obtained in \cite{Cunha:2017qtt} for an ultracompact horizonless object and \cite{Cunha:2020azh} for a BH.

Finally, let us comment on the interlacement of the energy conditions and all results discussed so far. Consider now that we have a hypothetical spacetime such that it possesses a fully unstable LR, hence it violated the NEC. From the result obtained in the previous section, we know that in the immediately vicinity of the LR we have a region which harbours fully unstable TCOs -- \textit{cf.} Fig. \ref{Fig:StructureFullyUnstable}. Thus, by Eq. \eqref{Eq:SECObeyingTCOs} and the following argument, the SEC is also violated.
This is consistent with the implication tree of the NEC and SEC. Namely, that the violation of the NEC also implies the violation of the SEC.

\section{Conclusions and final remarks}

In this work, we have analysed the radial and vertical epicyclic frequencies of null and timelike particles following circular orbits in a spacetime generated by a generic stationary, axisymmetric, asymptotically flat ultracompact object with a $\mathbb{Z}_2$ symmetry. We have shown that the epicyclic frequencies of LRs coincide with the epicyclic frequencies of TCOs when the latter approach the former. This implies, by continuity, that the full stability of a LR determines the full stability of TCOs in the immediately vicinity of the LR.
The results present in the first part of this paper generalise the ones obtain in \cite{Delgado:2021jxd}, where the author only analysed the radial stability of the LRs and TCOs.
Together with the results of \cite{Delgado:2021jxd}, we determined the 4 possible structures that one can found for a generic spacetime with a LR -- \textit{cf.} Figs. \ref{Fig:StructureFullyUnstable} -- \ref{Fig:StructureFullyStable}.

The results of the first part of this paper can be summarised as follows: the \textit{radial} stability of a LR determines the localisation and \textit{radial} stability of TCOs in its vicinity, whereas, the \textit{vertical} stability of a LR only determines the \textit{vertical} stability of TCOs in its vicinity.

In the second part of this work, we established relations between the sum of the squared epicyclic frequencies with linear combinations of the Ricci tensor components, that generalise the results obtained in \cite{Vieira:2017zau}. Such relations opened the possibility of introduce some energy conditions regarding the matter that the generic ultracompact object may be compose of. In particular, we showed that if the object obeys the SEC, which implies that the NEC is also obeyed, LRs and TCOs can never be fully unstable. Reciprocally, if a generic ultracompact object possesses fully unstable LRs, it violated the NEC, and consequently the SEC. This is consistent with the works presented in \cite{Cunha:2017qtt,Cunha:2020azh}. Furthermore, if the same object only possesses a region with fully unstable TCOs, then the object only violates the SEC.

We also study the more simpler case of a Ricci-flat ultracompact object. In this case, we showed that the LRs can never be neither fully stable or fully unstable, since the sum of the squared epicyclic frequencies of each LR always vanishes -- \textit{cf.} Eq. \eqref{Eq:RicciFlatLRs}. Regarding the TCOs, the sum of its epicyclic frequencies is always positive, hence, a Ricci-flat object never has regions of fully unstable TCOs.

As future work, it would be interesting to generalise further the results presented here to more generic ultracompact objects that do not have a $\mathbb{Z}_2$ symmetry. For such objects, circular orbits will not always be found at the equatorial plane, but will be found at an angular coordinate $\theta_\text{cir}$ that will be a function of the radial coordinate, $\theta_\text{cir} = f(r_\text{cir})$. Due to the missing $\mathbb{Z}_2$ symmetry, the relation between the full stability of a LR and the full stability of TCOs in its vicinity will, most likely, be more convoluted and less obvious. Nevertheless, one could obtain all the possible structures of TCOs around LRs and also study which are connected to violations of the energy conditions.

\section{Acknowledgements}

We would like to thank C. Herdeiro and E. Radu for comments on a draft of this paper.
This work is supported by the Center for Research and Development in Mathematics and Applications (CIDMA) and Center for Astrophysics and Gravitation (CENTRA) through the Portuguese Foundation for Science and Technology (FCT - Fundação para a Ciência e a Tecnologia), references UIDB/04106/2020, UIDP/04106/2020 and UIDB/00099/2020.
We acknowledge support from the projects PTDC/FIS-OUT/28407/2017, CERN/FISPAR /0027/2019,
PTDC/FIS-AST/3041/2020. This work has further been supported by the European Union’s Horizon 2020 research and innovation (RISE) programme H2020-MSCA-RISE-2017 Grant No. FunFiCO-777740.

\bibliographystyle{ieeetr}