\section{Introduction}

Scattering between excitons and free carriers has been observed experimentally since the 1960s in highly excited bulk semiconductors\cite{PhysRev.177.567,10.1063/1.1653814,10.1143/JPSJ.34.693,10.1002/pssb.2220580237,10.1002/pssa.2210220102,10.1002/pssb.2220710216,10.1002/pssb.2220780219,Levy1975ExperimentalIO}. In these bulk semiconductor systems, the scattering cross sections have been studied previously\cite{ELKOMOSS1977557,ELKOMOSS1979431}, with the distinction between elastic and inelastic scattering fundamental for the interpretation of the experimental data.

Recently, the quality of monolayer semiconductor samples has drastically increased\cite{Reynolds1986,Museur2008,Cadiz2017,Elias2019,Li2021a}, with advances in techniques such as molecular beam epitaxy\cite{Vuong2017}, chemical vapor deposition\cite{Kim2015} or solution methods\cite{Yano2000}. The improved quality of these samples allowed for a much more detailed study of excitons in mono- and few--layer materials, where well--defined resonance peaks are experimentally identifiable even at room temperature\cite{PhysRevLett.33.827,10.1063/1.93648,Poellmann2015,Merkl2019}. This then creates the necessity of calculating and measuring the excitonic linewidth.

Several mechanisms contribute to the excitonic linewidth, such as acoustic\cite{PhysRevB.33.5512} and optical\cite{PhysRevB.34.2554} phonon scattering \cite{PhysRevB.103.235402}, radiative recombination\cite{10.1080/00018739100101522,10.1021/nl503799t}, as well as scattering in semiconductor alloys\cite{10.1063/1.94649,10.1063/1.102583}. Besides these mechanisms, others can play a part in determining the exciton linewidth. In this paper, we focus our attention on one of these mechanisms, namely scattering with free carriers. This scattering mechanism can play an important part in determining the excitonic linewidth, especially in systems with high density of free carriers and excitons\cite{10.1063/1.360578}, high pump fluences\cite{Djurii2005}, in tunneling experiments\cite{Young1990,Teran2009} or when an electric field is applied to the semiconductor\cite{10.1063/1.363972,10.1002/adom.202102132}.  Additionally, electron exciton scattering  processes also play an important role when studying exciton - polaron systems\cite{PhysRevB.103.075417}.

This paper is structured as follows. In Sec. \ref{sec:sec2}, we quickly review the approach outlined by Feng and Spector in Ref. \cite{FENG1987593} to the scattering between free carriers and excitons in semiconducting quantum wells. This derivation was performed in the central field \cite{Massey1951} and Born approximations \cite{Bates1953}.
In Sec. \ref{sec:sec3}, we turn our discussion to the total scattering cross section. We discuss the distinction between elastic and inelastic scattering, briefly reviewing the variational exciton wave functions for various states. We then explicitly compute the total cross section for a few select transitions, discussing the thresholds present in inelastic scattering processes. Finally, in Sec. \ref{sec:sec4}, we compute the contribution to excitonic linewidth from the scattering cross section with free carriers, analyzing its dependence on both the temperature of the system and the free carrier density in the monolayer.
\section{Free Carrier - Exciton Scattering\label{sec:sec2}}

In this section, we will follow the expressions derived by Feng and Spector in Ref. \cite{FENG1987593} for free carrier - exciton scattering, based on assuming two--dimensional (2D) gases of free carriers (electrons or holes) and excitons interacting with one another. The obtained expressions are derived following the same approach as those discussed by Mott and Massey for the general theory for tridimensional (3D) two--body collisions \cite{Mott1965}. The cross sections due to collisions between the carriers and excitons are then calculated using the central field \cite{Massey1951} and Born approximations \cite{Bates1953}.

\subsection{Differential Scattering Cross Section}

Let us begin by considering a two body system consisting of an exciton and a free carrier (electron or hole). The reduced mass of such a system is given by
\begin{align}
	M &=\frac{m_c\left(m_e + m_h\right)}{m_c + m_e + m_h},
\end{align}
where $ m_{c/e/h} $ is the mass of the free carrier/ electron/hole.

Following the derivations by Feng and Spector\cite{FENG1987593,FENG7096} in--depth in Appendix \ref{app:derivation}, the differential scattering cross section for our free carrier - exciton system is written as
\begin{equation}
	I_{fi}\left(\theta\right)=\frac{M^{2}}{2\pi\hbar^{4}k_{i}}\left|\int d^{2}\mathbf{r}\,d^{2}\mathbf{R}\,e^{i\mathbf{q}\cdot\mathbf{R}} V\left(\mathbf{r},\mathbf{R}\right)\chi_{f}^{\dagger}\left(\mathbf{r}\right)\chi_{i}\left(\mathbf{r}\right)\right|^{2},\label{eq:generic_crossection}
\end{equation}
where $ \mathbf{q}=\mathbf{k}_i-\mathbf{k}_f $ is the difference between the initial and final relative momentum of the system, $ \chi_{i/f} $ represent the initial/final exciton wave functions, and $ V\left(\mathbf{r},\mathbf{R}\right) $ the interaction potential between the free carrier and the exciton. In Eq. (\ref{eq:generic_crossection}), $ \mathbf{r} $ is the relative position vector of the electron and the hole in the exciton, and $ \mathbf{R} $ is the relative position vector from the free carrier to the center of mass of the exciton.

The interaction potential between the free carrier and the exciton in the central field approximation\cite{Massey1951} will be modeled by the Rytova--Keldysh potential\cite{rytova1967,keldysh1979coulomb}, usually employed to describe excitonic phenomena in mono- and few--layer materials and obtained by solving the Poisson equation for a charge embedded in a thin film of vanishing thickness. In real space, the Rytova--Keldysh potential is given by
\begin{equation}
	V_{RK}(\mathbf{r})=\frac{\hbar c \alpha}{\epsilon}\frac{\pi}{2 r_{0}}\left[H_{0}\left(\epsilon\frac{r}{r_{0}}\right)-Y_{0}\left(\epsilon\frac{r}{r_{0}}\right)\right], \label{eq:RK}
\end{equation}
where $\alpha=1/137$ is the fine--structure constant, $\epsilon$ the mean dielectric constant of the medium above/below the layered material, $H_{0}\left(x\right)$ is the zeroth-order Struve function and $Y_{0}\left(x\right)$ is the zeroth-order Bessel function of the second kind. The parameter $r_{0}$ corresponds to an in--plane screening length related to the 2D polarizability of the material and can be calculated from the single particle Hamiltonian of the system \cite{PhysRevB.99.035429}. In the limit of zero screening length, the Rytova--Keldysh potential becomes the Coulomb potential. Considering, again, the interaction between the free carrier and the exciton, the interaction potential will be
\begin{equation}\label{eq:pot_old}
	V\left(\mathbf{r},\mathbf{R}\right)=\pm \left[V_{RK}(\mathbf{r_{ch}})-V_{RK}(\mathbf{r_{ce}})\right],
\end{equation}
where $ \pm $ distinguishes between the free carrier being a hole ($ + $) or an electron ($ - $), $ r_{ch} $ is the distance between the free carrier and the hole of the exciton, and $ r_{ce} $ is the distance between the free carrier and the electron of the exciton. These two vectors can be written from $ \mathbf{r} $ and $ \mathbf{R} $ as \cite{FENG1987593}
\begin{align}\label{eq:vectors_from_masses}
	\mathbf{r}_{ch} &= \mathbf{R}-\frac{\sigma}{1+\sigma}\mathbf{r},
	&\mathbf{r}_{ce} &= \mathbf{R}+\frac{1}{1+\sigma}\mathbf{r},
\end{align}
where $ \sigma=m_e / m_h $ is the ratio between effective electron and hole masses.
Returning to the discussion on the scattering cross section from Eq. (\ref{eq:generic_crossection}), the integration over $ R $ reads
\begin{equation}\label{eq:R_integration}
	\int d^{2}\mathbf{R}\,e^{i\mathbf{q}\cdot \mathbf{R}} V\left(\mathbf{r},\mathbf{R}\right)=\pm \left[\int d^{2}\mathbf{R}\,e^{i\mathbf{q}\cdot \mathbf{R}}V_{RK}(\mathbf{r_{ch}})-\int d^{2}\mathbf{R}\,e^{i\mathbf{q}\cdot \mathbf{R}}V_{RK}(\mathbf{r_{ce}})\right],
\end{equation}
which can be computed directly\cite{FENG7096} by performing a change of integration variables back to $\mathbf{r}_{ch},\,\mathbf{r}_{ce} $ and reads
\begin{equation}\label{eq:change_variables_int}
	\pm\left[e^{iq\frac{\sigma}{1+\sigma}r\cos\left(\phi_{r}\right)}-e^{-iq\frac{1}{1+\sigma}r\cos\left(\phi_{r}\right)}\right]2\pi\frac{\hbar c\alpha}{\epsilon}\frac{1}{q\left(1+r_{0}q\right)},
\end{equation}
with $ \phi_{r} $ the angle between $ \mathbf{r} $ and $ \mathbf{q}$, and
\begin{equation}
	V_{RK}\left(\mathbf{q}\right)=2\pi\frac{\hbar c\alpha}{\epsilon}\frac{1}{q\left(1+r_{0}q\right)}
\end{equation}
the Fourier transform of the Rytova--Keldysh potential. Finally, the differential scattering cross section can be written as

\begin{align}
	I_{fi}\left(\theta\right)=\frac{M^{2}}{2\pi\hbar^{4}k_{i}}\left|2\pi \frac{\hbar c\alpha}{\epsilon}\frac{1}{q\left(1+r_{0}q\right)}J\left(i\rightarrow f\right)\right|^{2}.\label{eq:cross-section_integrated_RK}

\end{align}
with the dependence on the initial and final exciton states included in $ J\left(i\rightarrow f\right) $, defined as
\begin{equation}\label{eq:J}
	J\left(i\rightarrow f\right)=\int_{0}^{+\infty}r\,dr\,\int_{0}^{2\pi} d\phi_{r}\,\left[e^{i\frac{\sigma}{1+\sigma}qr\cos \phi_{r}}-e^{-i\frac{1}{1+\sigma}qr\cos \phi_{r}}\right]\chi_{f}^{\dagger}\left(\mathbf{r}\right)\chi_{i}\left(\mathbf{r}\right).
\end{equation}

\section{Total Cross Section\label{sec:sec3}}

Knowing the differential cross section given by Eq. (\ref{eq:cross-section_integrated_RK}), we can now compute the full scattering cross section. To this effect, we must simply perform an angular integration in $ \theta $ as
\begin{align}\label{eq:full_scatter}
	Q_{\mathrm{i\rightarrow f}}&=\int_{-\pi}^{\pi}d\theta\,I_{\mathrm{i\rightarrow f}}\left(\theta\right).
\end{align}
Explicitly substituting Eq. (\ref{eq:cross-section_integrated_RK}), the full cross section is given by
\begin{align}\label{eq:scatter_theta_RK_tot}
	Q_{\mathrm{i\rightarrow f}}\left(k_{i}\right)&=\frac{2\pi M^{2}\left(\hbar c\alpha\right)^{2}}{\hbar^{4}\epsilon^{2}k_{i}}\int_{-\pi}^{\pi}d\theta\,\left|\frac{J\left(i\rightarrow f\right)}{q(1+r_{0}q)}\right|^{2}.
\end{align}
To compute this integral, however, we must first define the exciton wave functions which we will consider when computing Eq. (\ref{eq:J}). We must also define the type of scattering in question, as it will introduce both the specific $ \theta $ dependence in $ \mathbf{q} $ as well as specific thresholds for the relative momentum of the free carrier - exciton system from conservation of energy.

\subsection{Elastic Scattering}

In an elastic scattering process, the exciton remains in its ground state after the collision, meaning $ \left|\mathbf{k}_i\right|=\left|\mathbf{k}_f\right| $. As such, we can write
\begin{equation}
	q= 2 \left|\sin\left(\frac{\theta}{2}\right)\right| k_i.\label{eq:elast_q}
\end{equation}
Additionally, we also have $ \chi_{f,i}\left(\mathbf{r}\right)=\chi_{1s}\left(r\right) $, where we consider a simple variational \emph{ansatz} \cite{Quintela,Gomes2021} based on the eigenfunctions of the two--dimensional Hydrogen atom\cite{Lee1979,PhysRevA.43.1186} and given by
\begin{equation}\label{eq:ansatz}
	\chi_{1s}\left(r\right)=\mathcal{N}_{1s}e^{-r\gamma_{1s}/2},
\end{equation}
with $ \mathcal{N}_{1s} $ a normalization constant given by
\begin{equation}\label{eq:normalization}
	\mathcal{N}_{1s}=\left[\int rdrd\theta\,\left(e^{-r\gamma_{1s}/2}\right)^{2}\right]^{-1/2}=\frac{\gamma_{1s}}{\sqrt{2\pi}}
\end{equation}
and $ \gamma_{1s} $ a variational parameter. This variational parameter is computed by minimization of the energy expectation value of the Wannier equation\cite{Wannier1937}
\begin{equation}\label{eq:wannier}
	H=-\frac{\hbar^{2} c^{2}}{2\mu}\nabla^2+V_{RK}(\mathbf{r}),
\end{equation}
with $ 	V_{RK}(\mathbf{r}) $ the Rytova--Keldysh potential.

With this \emph{ansatz}, we can directly substitute the wave function into $ J\left(i\rightarrow f\right) $, given by Eq. (\ref{eq:J}), and obtain
\begin{align}
	J_{\mathrm{elast}}\left(q\right)=J\left(\mathrm{1s}\rightarrow \mathrm{1s}\right)
	&=\gamma_{1s}^{3}\left[\frac{1}{\left[\left(\frac{\sigma}{1+\sigma}q\right)^{2}+\gamma_{1s}^{2}\right]^{3/2}}-\frac{1}{\left[\left(\frac{1}{1+\sigma}q\right)^{2}+\gamma_{1s}^{2}\right]^{3/2}}\right]\label{eq:J_transition_elastic_exponential}
\end{align}
after integration. This is then substituted into Eq. (\ref{eq:scatter_theta_RK_tot}), reading
\begin{align}\label{eq:scatter_theta_RK_tot_elast}
	Q_{\mathrm{elast}}\left(k_{i}\right)&=\frac{2\pi M^{2}\left(\hbar c\alpha\right)^{2}}{\hbar^{4}\epsilon^{2}k_{i}}\int_{-\pi}^{\pi}d\theta\,\left|\frac{J_{\mathrm{elast}}\left(q\right)}{q(1+r_{0}q)}\right|^{2},
\end{align}
where $ q $ is given by Eq. (\ref{eq:elast_q}). This integral has no analytical solution and must be computed numerically.

\begin{figure}
	\centering{}\includegraphics{Fig1}
	\caption{Total elastic cross section for electron--exciton (black) and hole--exciton (red) scattering in $ \mathrm{hBN} $ encapsulated in fused quartz as a function of the initial relative wave vector. \label{fig:cross-section_coulomb_GaAs}}
\end{figure}

To finalize the computation of the elastic cross section, we must choose a set of material specific parameters. We consider those corresponding to monolayer hexagonal Boron--Nitride ($ \mathrm{hBN} $) encapsulated in fused quartz, with dielectric constant $ \epsilon=3.8 $ \cite{Serway2003}. The electron and hole masses in this material are $ m_e = 0.83\, m_0$, $ m_h = 0.63\, m_0$ \cite{Ferreira19}, with $ m_0 $ the electron rest mass, and $ r_0 = 10\,\text{\AA}$ \cite{henriques_optical_2020}. The obtained variational energy for the $ 1s $ excitonic state is $ E_{1s}=-58.9\,\mathrm{meV} $.

Varying the initial relative wave vector, we obtain the plot of the total elastic cross section from Eq. (\ref{eq:scatter_theta_RK_tot}) in Fig. (\ref{fig:cross-section_coulomb_GaAs}). We can see that the cross section for electron scattering is always larger than that for hole scattering, as expected from the fact that the reduced mass of the system is larger when the free carrier considered is an electron. A very quick increase from zero relative momentum up to a global maximum is also observed, consistent with the results of Feng and Spector \cite{FENG7096} for elastic scattering.

\subsection{Inelastic Scattering}

Let us now consider inelastic scattering between free carriers and the exciton. The relative momenta is now given by
\begin{equation}\label{eq:momenta_1}
	q^{2} = k_{f}^{2} + k_{i}^{2}-2k_{i}k_{f}\cos(\theta),
\end{equation}
with $ k_{f} $ obtained from conservation of energy as
\begin{equation}\label{eq:conservation_energy}
	k_{f}^{2}=k_{i}^{2}-\frac{2M}{\hbar^{2}}(E_{f}-E_{i}).
\end{equation}
Here, $ E_{f/i} $ is the energy of the final/initial state of the exciton, respectively. Substituting this relation into Eq. (\ref{eq:momenta_1}), we obtain
\begin{equation}\label{eq:momenta_2}
	q^{2} = 2k_{i}^{2}-\frac{2M}{\hbar^{2}}\Delta_{f,i}-2k_{i}^{2}\sqrt{1-\frac{2M}{\hbar^{2}}\frac{\Delta_{f,i}}{k_{i}^{2}}}\cos(\theta),
\end{equation}
with $ \Delta_{f,i}=E_{f}-E_{i} $. A threshold in $ k_i $ below which no scattering is allowed is immediately evident, obtained from Eq. (\ref{eq:momenta_2}) as
\begin{equation}\label{eq:threshold}
	k_{\mathrm{min}}=\sqrt{\frac{2M}{\hbar^{2}}(E_{f}-E_{i})}.
\end{equation}
Below this threshold, there is not enough energy in the scattering process to allow the jump between excitonic states $ i\rightarrow f $.

Besides knowing the energies of the final states, we must also know their wave functions. These are obtained\cite{Lee1979,Quintela,Gomes2021} in a similar form to Eq. (\ref{eq:ansatz}) and are given by
\begin{align}\label{eq:2s_2p_ansatz}
	\chi_{2s}\left(r\right)&=\mathcal{N}_{2s}\left(1-\frac{r}{d}\right)e^{-r\gamma_{2s}/2}, &
	\chi_{2p_{\pm}}\left(r\right)&=\mathcal{N}_{2p}r e^{\pm i \theta}e^{-r\gamma_{2p}/2},
\end{align}
where $ \gamma_{2s},\gamma_{2p} $ are variational parameters, $ \mathcal{N}_{2s},\mathcal{N}_{2p} $ are normalization constant given by
\begin{align}
	\mathcal{N}_{2s}&=\frac{2\gamma_{2s}^{2}}{\sqrt{\pi}\sqrt{3\gamma_{1s}^{2}-2\gamma_{1s}\gamma_{2s}+3\gamma_{2s}^{2}}},&
	\mathcal{N}_{2p}&=\frac{\gamma_{2p}^{2}}{2\sqrt{3\pi}},
\end{align}
and $ d $ is a parameter obtained by imposing orthogonality between $ \chi_{1s} $ and $ \chi_{2s} $, given by
\begin{equation}
	d=\frac{4}{\gamma_{1s}+\gamma_{2s}}.
\end{equation}
These wave functions, together with the $ 1s $ wave function, are plotted in Fig. (\ref{fig:wave_functions}) for monolayer $ \mathrm{hBN} $ encapsulated in fused quartz.
\begin{figure}
	\centering{}\includegraphics{Fig2}
	\caption{Radial probability density for the $ 1s $ (blue), $ 2s $ (orange) and $ 2p $ (green) states as described by the variational wave functions of Eq. (\ref{eq:ansatz}) and Eq. (\ref{eq:2s_2p_ansatz}) in a $ \mathrm{hBN} $ monolayer encapsulated in fused quartz. \label{fig:wave_functions}}
\end{figure}

\subsubsection{$ 1s \rightarrow 2s $ Transitions}
We will first consider $ 1s \rightarrow 2s $ transitions. To compute $ J_{2s}=J\left(1s\rightarrow 2s\right) $, we recall Eq. (\ref{eq:J}) and, after integration, obtain
\begin{align}\label{eq:J_1s2s}
	J_{2s}
	&= \frac{3\gamma_{1s}\left(\gamma_{1s}+\gamma_{2s}\right)\gamma_{2s}^{2}}{\sqrt{6\gamma_{1s}^{2}-4\gamma_{1s}\gamma_{2s}+6\gamma_{2s}^{2}}} \left[\frac{q^{2}\left(\frac{\sigma}{1+\sigma}\right)^2}{\left[q^2\left(\frac{\sigma}{1+\sigma}\right)^2+\left(\frac{\gamma_{1s}+\gamma_{2s}}{2}\right)^{2}\right]^{5/2}}-\frac{q^{2}\left(\frac{1}{1+\sigma}\right)^2}{\left[q^2\left(\frac{1}{1+\sigma}\right)^2+\left(\frac{\gamma_{1s}+\gamma_{2s}}{2}\right)^{2}\right]^{5/2}}\right],\nonumber
\end{align}
where $ q $ is obtained from Eq. (\ref{eq:momenta_2}) as
\begin{equation}\label{eq:q_1s2s}
	q^{2}=2k_{i}^{2}-\frac{2M}{\hbar^{2}}\Delta_{2s,1s}-2k_{i}^{2}\sqrt{1-\frac{2M}{\hbar^{2}}\frac{\Delta_{2s,1s}}{k_{i}^{2}}}\cos(\theta).
\end{equation}

As discussed above, the energy of the $ 2s $ state is obtained by minimization of the Wannier equation with the variational wave functions and its value is $ E_{2s}=-8.83\,\mathrm{meV} $ for our system. Explicitly computing the thresholds from Eq. (\ref{eq:threshold}), we obtain $ k_{\mathrm{min}}= 0.0833\,\text{\AA}^{-1}$ for electron--exciton scattering and $ k_{\mathrm{min}}= 0.0760\,\text{\AA}^{-1}$ for hole--exciton scattering.

\subsubsection{$ 1s \rightarrow 2p $ Transitions}

For computing $ J\left(1s\rightarrow 2p_{\pm}\right) $ following Eq. (\ref{eq:J}), we must quickly consider the distinction between $ p_{\pm} $ states. This is, however, not important, as the two states are degenerate and the two integrals $ J\left(1s\rightarrow 2p_{+}\right) $, $ J\left(1s\rightarrow 2p_{-}\right) $ are, in fact, equal. As such, we take into account the two $ 2p $ states by multiplying the total cross section by an angular momentum degeneracy factor $ g_{\ell}=2 $.

Explicitly, $ J_{2p}=J\left(1s\rightarrow 2p\right) $ is given by
\begin{align}\label{eq:J_1s2p}
	J_{2p}
	&=\frac{3\gamma_{1s}\gamma_{2p}^{2}\left(\gamma_{1s}+\gamma_{2p}\right)}{2\sqrt{6}}\left[\frac{iq\frac{\sigma}{1+\sigma}}{\left[q^2\left(\frac{\sigma}{1+\sigma}\right)^2+\left(\frac{\gamma_{1s}+\gamma_{2p}}{2}\right)^{2}\right]^{5/2}}+\frac{iq\frac{1}{1+\sigma}}{\left[q^2\left(\frac{1}{1+\sigma}\right)^2+\left(\frac{\gamma_{1s}+\gamma_{2p}}{2}\right)^{2}\right]^{5/2}}\right].\nonumber
\end{align}
The relative momentum $ q $ is defined analogously to Eq. (\ref{eq:q_1s2s}), with the only factor missing being the energy of the $ 2p_{\pm} $ states.

For our system, this energy is $ E_{2p}=-10.2\,\mathrm{meV} $. As such, the thresholds from Eq. (\ref{eq:threshold}) are now $ k_{\mathrm{min}}= 0.0822\,\text{\AA}^{-1}$ for electron--exciton scattering and $ k_{\mathrm{min}}=0.0749\,\text{\AA}^{-1}$ for hole--exciton scattering.

\subsection{Joint Elastic and Inelastic Scattering}

Finally, we consider the joint contribution to the scattering cross section from both elastic and inelastic scattering. This cross section will, therefore, involve a sum over final states, where only the $ 1s\rightarrow 1s $ contribution originates from elastic scattering processes. Explicitly, and for an arbitrary set of final exciton states $ f $, the total scattering cross section is given by
\begin{equation}\label{eq:scatt_total_1}
	Q_{\mathrm{Total}}=\sum_{f}Q_{1s\rightarrow f}.
\end{equation}
For the three transitions discussed previously, the sum in Eq. (\ref{eq:scatt_total_1}) is restricted and is explicitly written as
\begin{equation}\label{eq:scatt_total}
	Q_{\mathrm{Total}}=Q_{1s\rightarrow 1s}+Q_{1s\rightarrow 2s}+Q_{1s\rightarrow 2p}.
\end{equation}
This total scattering cross section is plotted in Fig. (\ref{fig:cross-section_total}) for both types of free carriers, together with the dashed lines representing the thresholds for the inelastic scattering processes considered.

Analyzing Fig. (\ref{fig:cross-section_total}), we can see that, when the same scattering process is allowed for both types of free carriers, electron--exciton scattering has a cross section roughly $ 1.5\times $ larger. This trend is, however, inverted between the threshold momentum for hole--exciton $ 1s \rightarrow 2p $ scattering and the threshold momentum for electron--exciton $ 1s \rightarrow 2p $, \emph{i.e.}, between $ k_i = 0.0749\,\text{\AA}^{-1} $ and $ k_i = 0.0822\,\text{\AA}^{-1}$. In this momentum range, the dominant $ 1s \rightarrow 2p $ process is already allowed for hole--exciton scattering, leading to a vastly superior cross section relative to that for electron--exciton scattering.

\begin{figure}
	\centering{}\includegraphics{Fig3}
	\caption{Total cross section for electron--exciton (black) and hole--exciton (red) scattering in $ \mathrm{hBN} $ encapsulated in fused quartz as a function of the initial relative wave vector. Vertical dashed lines represent the electron--exciton and hole--exciton scattering thresholds for the $ 1s \rightarrow 2p $ transition, located at $ k_{\mathrm{min},e^{-}}= 0.0822\,\text{\AA}^{-1}$ and $ k_{\mathrm{min},h}= 0.0749\,\text{\AA}^{-1}$.

		\label{fig:cross-section_total}}
\end{figure}

The final threshold included, visible in Fig. (\ref{fig:cross-section_total}) as the dashed purple lines, originates from $ 1s \rightarrow 3d $ scattering, as $ E_{3d}=-3.74\,\mathrm{meV} $ is the lowest energy state after $ 2s $. These take place at $ k_{\mathrm{min}}=0.0875\,\text{\AA}^{-1} $ for electron--exciton scattering and $ k_{\mathrm{min}}=0.0798\,\text{\AA}^{-1} $ for hole--exciton scattering. These are, however, much smaller than the peaks in Fig. (\ref{fig:cross-section_total}), similarly to what is observed Ref. \cite{FENG1987593}, and are invisible in Fig. (\ref{fig:cross-section_total}).

\section{Scattering Contribution to Exciton Linewidth\label{sec:sec4}}

To conclude our study of the scattering of excitons with free carriers in layered materials, we will now discuss the contribution from these scattering processes to the excitonic linewidth. This will provide a point of comparison against experimental studies\cite{Selig2016}. Although other phenomena will also contribute to this linewidth, such as radiative lifetimes\cite{Palummo2015} and phonon scattering\cite{Henriques2021}, the dependence of the free carrier scattering on both temperature and carrier density should provide a good distinction of the various contributing processes.

The contribution to the excitonic linewidth from free carrier scattering is given by\cite{PhysRevB.34.2554,FENG1987459,feng1989excitons,10.1063/1.364274}
\begin{equation}\label{eq:linewidth}
	\Gamma_{\mathrm{Total}}=\sum_{f}\frac{2\hbar^{2}}{\pi M}\int_{0}^{\infty}dk\,k^{2}f\left(\frac{m_c+m_e+m_h}{m_e+m_h}k\right)Q_{1s\rightarrow f},
\end{equation}
where $ Q_{1s\rightarrow f} $ is, as described earlier, the scattering cross section associated with a specific transition from the excitonic ground state to a final state $ f $. As we are summing over final states $ f $ and only $ Q_{1s\rightarrow f} $ depends on the final state, this is equivalent to switching the sum over final states and the integral and writing
\begin{equation}\label{eq:linewidth_tot}
	\Gamma_{\mathrm{Total}}=\frac{2\hbar^{2}}{\pi M}\int_{0}^{\infty}dk\,k^{2}n_{F}\left(\frac{m_c+m_e+m_h}{m_e+m_h}k\right)Q_{\mathrm{Total}},
\end{equation}
where $ Q_{\mathrm{Total}} $ is the total scattering cross section as plotted in Fig. (\ref{fig:cross-section_total}).
Here, $ n_{F}(k) $ is the Fermi--Dirac distribution for free carriers, given by
\begin{equation}\label{eq:fermi_dirac-distr}
	n_{F}(k)=\frac{1}{e^{\frac{E_k-E_F}{k_B T}}+1}
\end{equation}
where the dispersion relation is given by
\begin{equation}\label{eq:free-carrier_dispersion}
	E_k=\frac{\hbar^{2}}{2 m_{c}}k^{2}
\end{equation}
and the Fermi energy is
\begin{equation}\label{eq:fermi_energy_free}
	E_F=2\pi \frac{\hbar^{2}}{2 m_{c}}n,
\end{equation}
with $ n $ the area density of free carriers. We consider carrier densities up to a maximum of $ 10^{12}\,\mathrm{cm}^{-2} $, where the average separation of free carriers $ d_{c}= 2 (\pi n)^{-1/2} $ is still larger but already of the order of the root mean square (RMS) exciton radius, given by \begin{equation}\label{eq_rms}
	r_{\mathrm{RMS};n}=\int_{0}^{\infty}\int_{0}^{2\pi} rdrd\theta\,\psi_{n}(r,\theta)^{\dagger} r^{2} \psi_{n}(r,\theta).
\end{equation}
For the two excitonic states most relevant to the scattering cross section, $ 1s $ and $ 2p $, the RMS exciton radius is $ r_{\mathrm{RMS};1s}=19.5\,\text{\AA} $ and $ r_{\mathrm{RMS};2p}=73.2\,\text{\AA} $, respectively, while the average separation between free carriers is $ d_{c}=113\,\text{\AA} $.

\begin{figure}
	\includegraphics[scale=0.78]{Fig4_a}\includegraphics[scale=0.78]{Fig4_b}
	\caption{Contribution to exciton linewidth from scattering with free carriers as a function of the temperature at constant free carrier area density $n=10^{9}\,\text{cm}^{-2}$ (left) and $n=10^{12}\,\text{cm}^{-2}$ (right). Insets show the low--temperature region of the plot for clarity.\label{fig:linewidth_const_density}}
\end{figure}

\begin{figure}
	\centering{}\includegraphics[scale=0.9]{Fig5_a}\includegraphics[scale=0.9]{Fig5_b}
	\includegraphics[scale=0.9]{Fig5_c}\includegraphics[scale=0.9]{Fig5_d}
	\caption{Contribution to exciton linewidth from scattering with free carriers as a function of the free carrier area density at constant temperature $ T=10\,\mathrm{K}$ (top-left), $  T=100\,\mathrm{K}$ (top-right), $  T=200\,\mathrm{K}$ (bottom-left), and $  T=300\,\mathrm{K}$ (bottom-right). \label{fig:linewidth_const_temp}}
\end{figure}

As before, we must compute the integral of Eq. (\ref{eq:linewidth_tot}) numerically. The specific methodology for the discretization of Eq. (\ref{eq:linewidth_tot}) is discussed in Sec. \ref{app:numeric_integral}. Choosing a Gauss--Legendre quadrature\cite{Kythe2002} of size $ N=450 $, the results for the contribution of scattering with free carriers to the excitonic linewidth are presented in Fig. (\ref{fig:linewidth_const_density}) as a function of temperature for free carrier area densities of $ n = 10^{9}\,\,\text{cm}^{-2} $ and $ n = 10^{12}\,\,\text{cm}^{-2} $. In Fig. (\ref{fig:linewidth_const_temp}), we present the excitonic linewidth as a function of the free carrier area density for four distinct values of the temperature $ T $ between $ 10\,\mathrm{K} $ and $ 300\,\mathrm{K}  $.
\section{Conclusions}

In this paper, we studied the effects of scattering between free carriers and excitons in monolayer materials and its contribution to the excitonic linewidth. To this end, we began by reviewing the general form of the expressions for the differential cross section between free carriers and excitons in two dimensions\cite{FENG1987593}, as well as the inclusion of screening in the differential cross section. Quickly reviewing the computation of variational functions for the exciton wave functions, we discussed both elastic and inelastic scattering processes.

With the differential cross section known for both elastic and inelastic scattering, we proceeded to the computation of the total scattering cross section.
While elastic scattering is the sole contributor for low relative momentum, inelastic $ 1s\rightarrow 2p $ scattering dominates the total cross section after it becomes allowed. This dominant scattering has a maximum cross section roughly $ 40 $ times larger than that of elastic scattering, and at least $ 400 $ times larger than both $ 1s\rightarrow 2s $ and $ 1s\rightarrow 3d $ inelastic scatterings. This dominant behavior of $ 1s\rightarrow 2p $ transitions is very similar to what is is presented in Figs. (2), (9) of Ref. \cite{FENG1987593}. This behavior is, of course, dependent on the ratio between electron and hole masses, although we only focus our attention on the masses characteristic of single layer hBN.

After discussing the relation between the total scattering cross section and the excitonic linewidth, we considered both its dependence on both temperature and free carrier density. First looking at the exciton linewidth for fixed values of the carrier density, we observe that the difference of linewidth between the two densities studied ($ 10^{9}\,\,\text{cm}^{-2} $ and $ 10^{12}\,\,\text{cm}^{-2} $) is noticeable in both the finite value of the linewidth as $ T$ approaches $0\,\mathrm{K} $ and in the crossing of the two curves happening slightly earlier at higher carrier density (at $ T\gtrsim150\,\mathrm{K} $ for $ n=10^{9}\,\,\text{cm}^{-2} $ and at $ T\approx140\,\mathrm{K} $ for $ n=10^{12}\,\,\text{cm}^{-2} $). An even greater noticeable change in behavior would be present at densities above $ 10^{13}\,\,\text{cm}^{-2} $, reasonably larger than those usually obtained in experimental works. These higher carrier densities would already imply an average free carrier separation $ d_{c} $ much smaller than the RMS exciton radius for the $ 2p $ excitonic state, which would make more complex excitonic phenomena, such as biexcitons and trions\cite{Lampert1958,Kheng1993,Chen2002,Li2003}, increasingly significant.

Regarding the excitonic linewidth at fixed temperatures, the free carrier scattering contribution remained in the $ 0.5-20\,\mathrm{meV} $ range for temperatures between $ 100-300\,\mathrm{K} $, although a much faster increase for higher densities is observed at $ T=10\,\mathrm{K} $. For temperatures in the range $ 100-300\,\mathrm{K} $, the computed linewidth remained essentially constant as carrier density increases until roughly $ n\approx10^{11}\,\,\text{cm}^{-2} $, although the actual value of the linewidth is strongly dependent on the temperature. Past $ n\approx10^{11}\,\,\text{cm}^{-2} $, the computed linewidth rapidly increases for all values of the temperature, although the growth occurs sooner for lower temperatures. It should, therefore, be feasible to measure this contribution to the excitonic linewidth, as this is both an attainable carrier density \cite{Du2016,Kuiri2021,Li2021}, and good quality samples of both hBN and transition metal dichalcogenides have been grown in the past which presented excitonic linewidths in this order of magnitude\cite{Reynolds1986,Museur2008,Cadiz2017,Elias2019,Li2021a}.

\funding{M. F. C. M. Q. acknowledges the International Nanotechnology Laboratory (INL) and the Portuguese Foundation for Science and Technology (FCT) for the Quantum Portugal Initiative (QPI) grant SFRH/BD/151114/2021.
N. M. R. P. acknowledges support by the Portuguese Foundation for Science and Technology (FCT) in the framework of the Strategic Funding UIDB/04650/2020, COMPETE 2020, PORTUGAL 2020, FEDER, and  FCT through projects PTDC/FIS-MAC/2045/2021, EXPL/FIS-MAC/0953/ 2021, and from the European Commission through the project Graphene Driven Revolutions in ICT and Beyond (Ref. No. 881603, CORE 3).}

\dataavailability{Not applicable.}

\conflictsofinterest{The authors declare no conflict of interest.}

\abbreviations{Abbreviations}{
The following abbreviations are used in this manuscript:\\

\noindent
\begin{tabular}{@{}ll}
hBN & hexagonal Boron Nitride\\
RMS & root mean squared
\end{tabular}
}


\section{Introduction}

There is nowadays a compelling direct evidence of the existence of dark compact objects.
While general-relativistic black holes (BHs) --- characterized solely by their mass~$M$ and specific angular momentum~$a$ --- are favorite to describe these astrophysical objects, beneath their horizon they house a classically inevitable singularity.
It is widely believed, though, that quantum-gravity effects change the internal structure of BHs at some scale $\ell$ and cure the singularity issue.
Without specifying the actual theory responsible for the singularity regularization, there is only a handful of viable alternatives to classical BHs, i.e.~regular black holes (RBHs), wormholes, black bounces, and quasi black holes~\cite{Carballo-Rubio:2018jzw,carballo-rubio_opening_2020,carballo-rubio_geodesically_2020}.

Here we focus on RBHs, i.e.~geometries that are almost indistinguishable from a classical BH for a distant observer, yet with a smooth regular core. (For a classification of possible asymptotic geometries for the core in the static case see ref.~\cite{bargueno_global_2020}.)
For this reason, it is often useful to analyze these spacetimes through the lens of general relativity (GR), as we expect quantum-gravity-induced deviations only close to the core, and there being described by some effective-matter stress--energy tensor distribution making the regular geometry a solution of Einstein's field equations. As a quite general feature, such effective stress--energy tensor sports almost always a violation of the weak energy condition (and possibly of other energy conditions). This is not surprising given the expected quantum-gravitational origin of this feature. In several cases, one can even interpret such matter content as a consequence of coupling gravity to some non-linear electrodynamics theory~\cite{AyonBeato:2000,AyonBeato:2004,Balart:2014,Rodrigues:2018,Bronnikov:2000,Bolokhov:2012,Bronnikov:2017,Bokulic:2022cyk}, albeit one has to take into account that such an interpretation is not strictly speaking necessary as it is not obvious that quantum-gravitational effects should sum up to some classical model.

For spacetimes which are everywhere well-defined and analytical, RBHs can have a non-singular core by sporting inside the outer horizon or a wormhole throat or an inner horizon~\cite{carballo-rubio_opening_2020,carballo-rubio_geodesically_2020}. Here we shall focus on the second class of regular solutions.\footnote{See however~\cite{simpson_black-bounce_2019,Franzin:2021vnj,mazza_novel_2021} for static or rotating examples of the first class and~\cite{Franzin:2022iai, Riaz:2022rlx} for some of their phenomenological implications.}

The viability of such class of solutions has been recently questioned by the realization that their inner horizon(s) are generically associated to an exponential instability \cite{Carballo-Rubio:2018pmi,carballo-rubio_inner_2021,Carballo-Rubio:2022kad} known as mass inflation~\cite{Poisson:1989zz} regulated by their surface gravity, which is generically non-zero.\footnote{While it has been claimed that such an exponential instability is replaced by a milder polynomial growth in some RBH solutions~\cite{Bonanno:2020fgp}, it was later shown that the onset of such behavior is generically preceded by an exponential, and hence fatal, phase~\cite{carballo-rubio_inner_2021,DiFilippo:2022qkl}.} Remarkably, this instability has been proven to be avoidable in a recent spherically symmetric RBH proposal where a solution endowed with a degenerate inner horizon, characterized by a vanishing surface gravity, was constructed and shown to be stable~\cite{Carballo-Rubio:2022kad}. In what follows we shall generalize this idea and extend it to the phenomenologically more relevant rotating RBH spacetimes.

Indeed, the first RBH models have been introduced assuming spherical symmetry~\cite{ansoldi_spherical_2008}, which of course led to solutions lacking a key feature for being astrophysical GR-black-hole mimickers: rotation. Over the years, there have been several proposals for rotating RBHs, often by means of the Janis--Newman procedure. Even in this case, the required matter content to sustain the solution can be interpreted as associated to some non-linear electrodynamics~\cite{dymnikova_regular_2015,turimov_generic_2021,toshmatov_rotating_2014} and suitable criteria for reasonable RBHs were proposed so to rule out some models, or support others~\cite{maeda_quest_2021}.

Apart from the above-mentioned mass-inflation instability, an unpleasant common feature of the proposed rotating spacetimes is that the curvature invariants are found to be finite but not continuous at the would-be ring singularity \cite{bambi_rotating_2013}.
In addition, as for Kerr BHs, rotating RBHs generically still allow for closed timelike curves (CTCs).
Nonetheless, these models represented the simplest alternatives to Kerr BHs and for this reason their phenomenological features have been extensively studied~\cite{abdujabbarov_shadow_2016, Tsukamoto:2017fxq, ghosh_ergosphere_2020,Kumar:2020yem, akiyama_first_2022,kumar_walia_testing_2022, banerjee_signatures_2022}. In order to overcome the above-mentioned shortcomings of the so-far proposed RBH solutions, it is worth analyzing in more detail how these are normally constructed, i.e.~how singularities are regularized.

The most common approach to model RBH geometries is by considering the Schwarzschild or Kerr metric and promoting its mass $M$ to be a function of the radial coordinate, $M\mapsto m(r)$. This is done in most of the above-mentioned proposals and generally amounts to requiring a $m(r)$ which goes to zero suitably fast close to the singularity and that after a transient region goes to the prescribed black-hole Arnowitt--Deser--Misner (ADM) mass $M$ at infinity.

Alternatively, it is also possible to remove the singularity by means of a suitable Weyl rescaling~\cite{narlikar_space-time_1977} of the metric, i.e.~by multiplying the Schwarzschild or Kerr metric by a conformal factor~\cite{bambi_spacetime_2017,jusufi_rotating_2020}, or simply of some metric components~\cite{Canate:2022zst}.

In this paper we propose to combine these two approaches to obtain a rotating RBH which is free of all the issues discussed above.
In particular the conformal factor will be used to regularize the singularity, while the mass function will be chosen so to avoid the mass inflation along the idea proposed in~\cite{Carballo-Rubio:2022kad} (which cannot be as easily implemented in the rotating case by using only a mass function).
In summary, we shall show that our candidate spacetime has a number of desirable features:
the curvature invariants are everywhere finite and continuous;
the region of the spacetime at $r<0$, responsible for CTCs, is inaccessible to timelike and null geodesics;
the inner horizon has zero surface gravity and as such is stable;
the outer/trapping horizon is located as in the Kerr spacetime;
the ergoregion, light ring and innermost stable circular orbit are similar, in shape and location, to their Kerr counterparts, yet allowing for, in principle, detectable differences.

\medskip

The paper is structured as follows. \Cref{sec:RBHs} introduces the basic ingredients for regularizing BHs, discusses some common pitfalls of traditional RBHs and introduces our strategy for avoiding them. \Cref{s:conformal} explains how a conformal factor can cure singularities and presents our own choice. \Cref{sec:stabilizing} shows how a carefully picked mass function can stabilize the inner horizon; we then recap the metric we are putting forward with \cref{eq:finalmetric}. In \cref{sec:mimicker} we investigate the extent to which the object we propose deviates from a Kerr BH by exploring some simple phenomenological applications: first we describe the ergosurfaces, then compute the coordinate location of the light ring and ISCO. Finally, \cref{sec:conclusions} reports our conclusions.

\section[Basic ingredients of regular black holes]{Basic ingredients of regular black holes\label{sec:RBHs}}

Most static and spherically symmetric RBHs can be described by the line element
\be\label{eq:seed}
ds^2 = -\left(1-\frac{2m(r)}{r} \right) \dd{t}^2 + \frac{\dd{r}^2}{\left(1-\frac{2m(r)}{r} \right)} + r^2 \left[ \dd{\theta}^2 + \sin^2 \theta \dd{\phi}^2 \right]\, .
\ee

For large $r$, $m(r)\to M$, i.e.~to the ADM mass, and the spacetime is asymptotically flat.
Regularity is ensured if the mass function $m(r) = \order{r^3}$ as $r\to0$~\cite{frolov_notes_2016} and $\lim_{r\to 0^+} m'(r)/r^2 \gtreqqless 0$ corresponds respectively to a de Sitter, Minkowski and anti-de Sitter core.

Some widely studied choices of $m(r)$ are summarized in \cref{tab:RBHs}. For all of these RBH models, thanks to the introduction of some characteristic (supposedly quantum-gravity-induced) length scale $\ell$, one finds that $r=0$ is a (regular) ``point'' --- it is the degenerate sphere obtained by shrinking constant-$r$ surfaces.
An interesting alternative, which we will not explore in this article, consists in endowing the center with a non-zero size: popular examples are black bounces \cite{simpson_black-bounce_2019,mazza_novel_2021,Franzin:2021vnj} and Bronnikov's black universes \cite{Bronnikov:2005, Bronnikov:2006}.

\begin{table}[ht]
\centering
\begin{tabular}{ll}
\toprule
Model & $m(r)$\\
\midrule
Bardeen \cite{bardeen_non-singular_1968} & $M \frac{r^3}{(r^2+\ell^2)^{3/2}}$\\
Hayward \cite{hayward_formation_2006} & $M \frac{r^3}{r^3 + 2M\ell^2}$ \\
Dymnikova \cite{dymnikova_vacuum_1992} & $M \left[1-\exp(\frac{r^3}{\ell^3}) \right]$\\
Fan--Wang \cite{Fan:2016hvf} & $M \frac{r^3}{(r+\ell)^3}$\\
\bottomrule
\end{tabular}
\caption{Some of the most popular RBH models. More examples e.g.\ in refs.~\cite{Ghosh:2014pba,Tsukamoto:2017fxq,Kumar:2018ple,Kumar:2020yem,simpson_eye_2022}}\label{tab:RBHs}
\end{table}

Perhaps the most straight-forward way to add rotation to these spacetimes is to repeat the same ``regularization'' on the Kerr metric, replacing its mass parameter $M$ with the function $m(r)$. The resulting line element is known as G\"urses--G\"ursey metric~\cite{Gurses:1975vu}:
\be\label{metricGG}
ds_\text{GG}^2 = -\left(1-\frac{2m(r)r}{\Sigma}\right) \dd{t}^2 - \frac{4a \, m(r) r \sin^2 \theta }{\Sigma} \dd{t} \dd{\phi} + \frac{\Sigma}{\Delta} \dd{r}^2 + \Sigma \dd{\theta}^2 + \frac{A \sin^2 \theta}{\Sigma} \dd{\phi}^2
\ee
with
\be
\Sigma = r^2 + a^2\cos^2\theta, \quad \Delta=r^2-2m(r)r +a^2, \quad A= (r^2+a^2)^2 -\Delta a^2 \sin^2 \theta\, .
\ee

Usually, the same $m(r)$ is used in the rotating and non-rotating case, although there might be reasons to introduce a dependence on the angle \cite{eichhorn_locality-principle_2021, eichhorn_image_2021-1}.
In the G\"urses--G\"ursey spacetime, a necessary and sufficient condition to avoid scalar curvature singularities is $m(0)=m'(0)=m''(0)=0$~\cite{torres_regular_2017}.
Basic properties of the line element \eqref{metricGG} as a rotating RBH are discussed in ref.~\cite{dymnikova_basic_2017}.

Killing horizons are located at values of $r$ that solve $\Delta = 0$; since this equation must have more than one root, these spacetimes generically possess inner horizons. As in the Kerr spacetime, such inner Killing horizons coincide with Cauchy horizons, hence these spacetimes are not globally hyperbolic.

Metrics of this form have been extensively used for phenomenological applications. From a theoretical point of view, however, these examples are not completely satisfying, for several reasons.

First of all, their inner horizons have non-zero surface gravity and for this reason they are prone to a phenomenon known as mass inflation. As a consequence of this phenomenon, any perturbation --- no matter how small --- will lead to an accumulation of energy close to the inner horizon and, eventually, to a large backreaction on the geometry. Since this phenomenon proceeds exponentially fast, it is regarded as a source of instability. Hence, RBHs as those described above can at most be metastable.

Secondly, the interpretation of the region close to the center is far from clear. As we said, in the non-rotating case the metric approaches that of an (anti) de Sitter or Minkowski spacetime close to $r=0$. As rotation is included, however, the would-be singularity becomes a ring located at $r=0, \ \theta=\pi/2$ and the simple interpretation in terms of maximally symmetric spaces is lost.
A particularly puzzling issue concerns the behavior of the curvature scalars close to the would-be singularity: the limits $r\to 0$ and $\theta \to \pi/2$ often do not commute (typically, if the $r$ limit is taken first the result is zero, but if the $\theta$ limit is taken first the result is a number proportional to inverse powers of $\ell$).
Hence, strictly speaking, the limit does not exist.\footnote{The $\Sigma\to0$ limit of the curvature invariants is well defined, and zero, if the mass function goes to zero faster than $r^3$. This is not the case for the most popular RBH models, but it is so for e.g.\ a Bardeen-like RBH with $m(r)=r^4/(r^2+\ell^2)^2$, as remarked also in ref.~\cite{Toshmatov:2017zpr}, and for the metric introduced in ref.~\cite{simpson_eye_2022}.}

Moreover, these spacetimes usually admit an analytical extension to a region of negative radial coordinate, where other pathologies might still arise. In particular, the $r<0$ region is responsible for the existence of CTCs and of a chronological horizon in the Kerr spacetime: since there seems to be no general mechanism to prevent the formation of CTCs in these RBHs, their existence becomes a model-dependent question.

\medskip

The line element \eqref{metricGG} can also be obtained via the Newman--Janis procedure with appropriate choices of complexification.
To avoid the intrinsic ambiguities of the standard implementation, Azreg-A\"inou~\cite{azreg-ainou_generating_2014, azreg-ainou_regular_2014, azreg-ainou_static_2014} introduced a modified version of the procedure that leads to a metric whose components are the same as in \cref{metricGG}, up to a free multiplicative function $\Psi(r, \theta)$,
\be\label{eq:PsiGG}
ds^2 = \frac{\Psi}{\Sigma}\,ds_\text{GG}^2\,,
\ee
where the free function $\Psi$ can then be constrained by some physical arguments. In particular, if the seed metric can be interpreted as a solution of some non-linear electrodynamics, the field equations imply $\Psi=\Sigma$.
Moreover, within the Newman--Janis framework, in the $a\to0$ limit, the line element \eqref{eq:PsiGG} must reduce to the seed metric, i.e.~$\lim_{a\to 0} \Psi(r, \theta) = r^2$.
The choice $\Psi = \Sigma$ thus seems the most natural, but it entails all the issues we just discussed.

On the other hand, a large body of work in the context of conformal gravity has shown that appropriate choices for the conformal factor can lead to geodesically complete spacetimes. In light of this fact, \cref{eq:PsiGG} is particularly suggestive, as it encodes a significantly larger amount of freedom than \cref{metricGG}. In particular, when $\Psi/\Sigma \neq 1$ it becomes possible to disentangle the regularization of the central singularity from the choice of $m(r)$, which instead determines the location and properties of the horizons.

The goal of this paper is therefore to show how one can exploit such freedom to build a RBH that is free of the issues presented above: in our proposal, $\Psi$ will be used to improve the appearances of the spacetime close to $r=0$, while $m(r)$ will be chosen so as to trim the properties of the inner horizon.

Note, incidentally, that adding a conformal factor to \cref{metricGG} is by no means disruptive: the spacetime described by the metric \eqref{eq:PsiGG} is of Petrov type D, exactly as \cref{metricGG}, with double null directions given by
\be
l^\mu = \sqrt{\frac{\Sigma}{\Psi}} \frac{1}{\Delta} \left\{r^2+a^2, \Delta, 0, a \right\}, \quad
n^\mu = \sqrt{\frac{\Sigma}{\Psi}} \frac{1}{\Delta} \left\{r^2+a^2, -\Delta, 0, a \right\};
\ee
these two null vectors can be complemented with
\be
m^\mu = \sqrt{\frac{\Sigma}{\Psi}} \frac{1}{\sqrt{2}\left(r + \mathrm{i} a \cos \theta\right)} \left\{\mathrm{i} a \sin \theta, 0, 1, \mathrm{i} \csc\theta \right\}
\ee
and its complex conjugate $\Bar{m}^\mu$ to form a Kinnersley-like tetrad. When  $\Psi(r, \theta)=\psi_r(r) + \psi_\theta (\theta)$ (i.e.\ it is ``separable''), the geometry admits the non-trivial Killing tensor
\be\label{eq:KT}
K_{\mu \nu} = \Psi(r, \theta) \left[l_\mu n_\nu + l_\nu n_\mu \right] + \psi_r (r) g_{\mu \nu}\, .
\ee
In this case, the equations of motion for a test particle of Killing energy per unit mass $E$ and Killing angular momentum along the axis of rotation per unit mass $L$ are
\be
\Dot{t} &= \frac{1}{\Psi\Delta}\left[A E - 2m(r) a r L \right]\,, \\
\Dot{\phi} &= \frac{1}{\Psi\Delta} \left[ \frac{L}{\sin^2\theta}\left(\Sigma-2m(r) r\right) +2m(r) a r E\right] \,, \\
\Psi^2 \Dot{r}^2 &= \left[(r^2 +a^2)E -aL\right]^2 - \Delta \left(\delta \psi_r +K \right)  \label{rgeod} \,, \\
\Psi^2 \Dot{\theta}^2 &= K-\delta \psi_\theta - \left(aE\sin\theta -\frac{L}{\sin\theta} \right)^2 = Q + \cos^2\theta \left(E^2a^2 - \frac{L^2}{\sin^2\theta} \right) - \delta \psi_\theta \,,
\ee
where $\delta=0$ for massless particles and $\delta=1$ for massive ones, while $K$ is the conserved quantity associated to the Killing tensor \eqref{eq:KT}, $K=K_{\mu \nu}\Dot{x}^\mu \Dot{x}^\nu$, and $Q = K -(Ea-L)^2$. Clearly, planar equatorial orbits are possible only if $\psi_\theta(\pi/2)=0$.

In the more general case in which $\Psi$ is not separable, the equations of motion are more involved and not separable. Motion with $\Ddot{\theta} = \Dot{\theta}=0$ can take place on the equator and on the axis of symmetry if $\partial_\theta \Psi= 0$ there. Note that if $\Psi$ is a function of $\Sigma$ only this is always the case, since $\partial_\theta \Psi = \Psi' \partial_\theta \Sigma = 2a^2 \cos\theta \sin\theta \Psi'$.

\section{Regularizing the singularity with $\Psi$}\label{s:conformal}

In this section we discuss how the function $\Psi$ can regularize the spacetime, regardless of the specific choice of $m(r)$. We assume such $\Psi$ will satisfy a very minimal set of requirements, namely: $\Psi(r,\theta)>0$ everywhere, in order to ensure that no additional singularities are introduced; and
\be
\frac{\Psi}{\Sigma} = 1+\order{\frac{1}{r^2}} \qq{as} r\to \infty\, ,
\ee
so that the spacetime's ADM mass and specific angular momentum are still given by the parameters $M =\lim_{r\to \infty}m(r)$ and $a$, respectively --- this is tantamount to a slightly stricter version of the usual asymptotic-flatness condition. (Note, in particular, that we do not follow the physical interpretation of \cite{azreg-ainou_generating_2014, azreg-ainou_regular_2014,azreg-ainou_static_2014} and hence we do not impose the partial differential equations that descend from that reasoning.)

Let us now look for a $\Psi$ that regularizes the singularity of the Kerr BH, i.e.~one for which
\be
ds^2 = \frac{\Psi}{\Sigma}\,ds^2_\text{Kerr}
\ee
is the line element of a spacetime free of scalar polynomial curvature singularities. The same $\Psi$ will also regularize more general metrics characterized by a generic (analytic) $m(r)$. As will become clear momentarily, the function $\Psi$ can also ``remove'' regions of the spacetime with undesirable features.

A simple example of such $\Psi$ is
\be\label{eq:Psi}
\Psi = \Sigma + \frac{b}{r^{2 z}}\, ,
\ee
with $z$ a real number, which we will further constrain in a moment, and $b$ a positive constant with dimensions $[M]^{2z+2}$. Note that if $b \to 0$ as $M\to 0$ or as $a \to 0$ one may recover respectively the Minkowski or the Schwarzschild metric. The Ricci scalar has the form
\be
R = -\frac{6b\, r^{2 z}}{r^2\Sigma^2 \left(b+r^{2z} \Sigma \right)^3}P_z\left(r, \cos\theta\right)\, ,
\ee
where $P_z\left(r,\cos\theta\right)$ is an expression (a polynomial in $r$ and $\cos\theta$ when $z$ is an integer) that goes to zero at least as fast as $\Sigma^2$ in the limit $r\to 0,\ \theta\to \pi/2$.
Hence, the Ricci scalar never blows up for $z \geq 1$. However, $z=1$ still does not yield a well-defined limit, while for $z>1$ the limit exists and is zero, irrespective of the path taken to reach the would-be singularity in the ($r$-$\theta$) space. Similar remarks hold for the Ricci tensor squared $R^{\mu \nu}R_{\mu \nu}$ and the Kretschmann scalar $R^\mu_{\ \nu \rho \sigma}R_\mu ^{\ \nu \rho \sigma}$.

The ansatz in \cref{eq:Psi} can be written as $\psi_r(r) + \psi_\theta (\theta)$, i.e.\ it is ``separable'' in the terminology introduced at the end of \cref{sec:RBHs}, and hence has the advantage of leading to separable equations of motion.

Note that, with $z>1$, $\Psi$ is divergent on the whole disk $r=0$ --- which will have consequences for CTCs. The fact that this divergence can in fact cancel the divergences in the curvature scalars is quite remarkable. For these reasons, \cref{eq:Psi} is the choice we will mostly explore in the remainder of the paper: in particular, we will often consider the ``minimal'' choice $z=3/2$, corresponding to the smallest integer exponent of $r$ that yields a well defined limit.

It is also worth mentioning that, if one focuses on the non-spinning case only, lower values of the exponent $z$ are required.
Indeed, in order to regularize the metric
\be
ds^2 = \frac{\Psi}{r^2}\,ds^2_\text{Schw}\,,
\ee
with the $a\to 0$ limit of \cref{eq:Psi}
\be
\Psi = r^2 + \frac{b}{r^{2z}}\,,
\ee
one must have $z \geq 1/2$.

Finally, before moving on, let us add that an interesting alternative to \cref{eq:Psi} can be represented by the ansatz
\be\label{eq:Psi2}
\Psi = \Sigma + \frac{b}{\Sigma^z} \,.
\ee
In this case it is easy to check that the Ricci scalar tends to zero for $r\to 0, \ \theta \to \pi/2$ for any $z>1$. The same holds true for the Ricci tensor squared and the Kretschmann scalar, hence \cref{eq:Psi2} seems equivalent to \cref{eq:Psi}. Notably, however, in this second case $z=1$ too yields a well-defined, and finite, limit
\be
\lim_{\Sigma \to 0} R = - \frac{24 a^2}{b}
\ee
and similar results can be found for $R^{\mu \nu}R_{\mu \nu}$ and the Kretschmann.\footnote{To our knowledge, those built with a $\Psi$ are the only examples of rotating RBHs whose curvature scalars are continuous and non-zero at the would-be singularity.}
With this choice, $\Psi$ only diverges on the ring $r=0, \ \theta=\pi/2$, but not on the disk $r=0, \ \theta\neq \pi/2$.
\Cref{eq:Psi2} will be juxtaposed to \cref{eq:Psi} in \cref{subsec:r0} to highlight the properties that make us prefer the latter.

\subsection{The spacetime close to $r=0$ \label{subsec:r0}}

Although the scalar curvatures we computed are everywhere finite, the components of the metric still diverge for $\Sigma = 0$.
Previous works~\cite{bambi_spacetime_2017} have argued that the resulting spacetime is in fact geodesically complete, since the would-be singularity is reached in infinite proper time. Since our choice of conformal factor is slightly different from that discussed in ref.~\cite{bambi_spacetime_2017}, we sketch the relevant computations below.

Consider first a particle moving on the equatorial plane $\theta=\pi/2$ and falling radially towards $r=0$. With $E$ and $L$ being the particle's energy and angular momentum per unit mass, its radial velocity satisfies
\be
\Psi^2 \Dot{r}^2 = \mathcal{R}_\text{Kerr} - \delta \Delta (\Psi-r^2)\, .
\ee
Here $\delta = 0$ or $1$ for massless or massive particles respectively, and $\mathcal{R}_\text{Kerr}$ is the right-hand side of \cref{rgeod} with $\Psi=\Sigma$. The proper time it takes for the particle to fall from $r_0$ to $r$ is
\be
\Delta \tau = - \int^r_{r_0} \frac{\Psi}{\left[\mathcal{R}_\text{Kerr} - \delta \Delta (\Psi-\tilde r^2) \right]^{1/2}} \dd{\tilde r}\, .
\ee
For both our ansatz (\cref{eq:Psi} or \cref{eq:Psi2}), one finds that on the equatorial plane $\Psi - r^2 = b\, r^{-2z} > 0$ ; therefore the infall time for massless particles ($\delta=0$) is shorter than that for massive particles ($\delta =1$).
(Obviously, this is true as long as $\mathcal{R}_\text{Kerr} - \delta \Delta (\Psi- r^2)>0$, i.e.\ only where the trajectory is classically allowed: where the condition is not met, such motion could not take place.)

Let us then focus on massless particles.
At $r=0$, $\mathcal{R}_\text{Kerr}=a^2(Ea-L)^2$, while $\Psi$ diverges at least as fast as $1/r^2$.
We conclude that massless particles reach the would-be singularity in an infinite amount of proper time. Given the inequality above, the conclusion remains true for massive particles.

Next, consider a particle that falls along the axis of symmetry $\theta=0$.
Such particle could reach the disk $r=0$ without encountering the would-be singularity, and potentially cross it through its center. On-axis motion requires $L=0$, so the radial velocity now satisfies
\be
\Psi^2 \Dot{r}^2= -\delta \Psi \Delta + E^2(r^2+a^2)^2
\ee
where $\Psi$ is now evaluated at $\theta=0$. The infall proper time becomes in this case
\be
\Delta \tau = -\int_{r_0}^r \frac{\Psi}{\sqrt{(\tilde r^2 + a^2)^2E^2-\delta \Delta \Psi}} \dd{\tilde r}\, .
\ee
First of all, we can see that it is still true that massless particles fall in a shorter time than massive ones, therefore we again focus on the former. We have
\be
E \Delta \tau_\text{light} = -\int \frac{\Psi}{r^2+a^2} \dd{r}\,
\ee
and with our ansatz we have
\be
\frac{\Psi}{r^2+a^2} = 1+b \begin{cases} (r^2+a^2)^{-1}r^{-2z} &\qq{for \cref{eq:Psi}}\\
(r^2+a^2)^{-(z+1)} &\qq{for \cref{eq:Psi2}.}
\end{cases}
\ee
In the first case, the integrand diverges faster than $r^{-2}$ as $r\to 0$, hence the particle will reach the would-be singularity in an infinite time. In the second case, instead, the integrand is everywhere finite. For massive particles, one can show that the infall time remains finite in the second case but, according to the inequality above, it is infinite in the first.

Therefore, the two choices of $\Psi$ lead to a very different structure of the region close to the would-be singularity: in the first case, the whole disk $r=0$ is (regularized and) ``sent to infinity''; in the second case, only the ring $r=0,\ \theta=\pi/2$ is pushed away, so that particles can still cross the disk inside the ring. This is a non-negligible difference as in the case of ansatz \cref{eq:Psi} we end up precluding, to light or matter, access to that region of the Kerr geometry ($r<0$) characterized by the presence of CTCs.

Let us stress that while usually such region is taken to be nonphysical in the Kerr geometry, due to the fact that it is shielded by a Cauchy horizon which is widely (albeit non-unanimously) considered unstable, the same region would represent a problem for us once we shall have proceeded to stabilize the RBH inner horizon by making it degenerate. It is henceforth even more pressing for stable RBH to chose an ansatz such as \cref{eq:Psi} over that of \cref{eq:Psi2}.

\section{Stabilizing the inner horizon with $m(r)$}\label{sec:stabilizing}

In ref.~\cite{Carballo-Rubio:2022kad} it was shown that the mass-inflation instability can be turned off if the surface gravity of the inner horizon $\kappa_{-}$ vanishes thanks to a wise choice of the mass function. The problem to extend this idea to the rotating case consists in the fact that adding this condition (i.e.~having a degenerate inner horizon), in addition to the conditions necessary to remove the ring singularity, avoid CTCs as well as have a well defined limit to the regularized region, makes the task daunting if not impossible. While we do not have a no-go theorem in this sense, it is rather clear to a first investigation that, even if viable, such regular metrics would be too cumbersome for any application towards phenomenology.

We shall then pursue a different path here, starting from the realization that if we regularize the singularity with the conformal factor as above, the functional form of the mass function is left with very few constraints (namely it must be everywhere finite and it must reduce to the ADM mass $M$ at infinity), and can be easily shaped so to stabilize the inner horizon.

Since the surface gravity of the inner horizon $r_-$ depends on $m(r)$ as
\be
\kappa_{-} \propto \partial_r \Delta |_{r=r_-}\, ,
\ee
if we assume a rational-function form for the mass function, to have a vanishing $\kappa_-$ the inner horizon must be a degenerate root of $\Delta$
\be
\Delta \equiv r^2-2 m(r) r+ a^2 = 0 \implies (r-r_+)(r-r_-)^d = 0\, ,
\label{delta}
\ee
for some $d \in\mathbb{N}_{\geq 2}$. $d=2$ is not viable, since it implies that $m(r)$ has a pole at some positive $r$. Thus the minimal choice ends up being $d=3$ which implies (given also the required asymptotic behavior) a mass function of the form
\be\label{m}
m(r) = M\,\frac{r^2+\alpha r + \beta}{r^2 + \gamma r+ \mu}.
\ee

From \cref{delta}, it can been shown that $\beta$ cannot be zero and thus the limit of $m(r)$ for $r \to 0$ is not zero but the finite value $M \beta/\mu$.
In this form, $m(r)$ is parametrized by four coefficients, two of dimension $[M]$ ($\alpha$ and $\gamma$) and two of dimension $[M]^2$ ($\beta$ and $\mu$).
However, through \cref{delta}, they can all be expressed as functions of the position of the two horizons
\be
\alpha &=\frac{a^4+r_-^3 r_+ - 3 a^2 r_- (r_- + r_+)}{2 a^2 M}\,,\\
\beta &=\frac{a^2 (2 M - 3 r_- - r_+) + r_-^2 (r_- + 3 r_+) }{2 M}\,,\\
\gamma &= 2M - 3r_- - r_+\,,\\
\mu &=\frac{r_-^3 r_+}{a^2}\,.
\ee
If we choose $r_+=M+\sqrt{M^2-a^2}$, i.e.~the outer horizon to coincide with its Kerr analog, our family of solutions can be parametrized in terms of $r_-$ only.
It is quite remarkable, and very relevant for phenomenological studies, that in spite of being located beyond a trapping horizon, the position of the inner horizon can matter for observables in the outside geometry. An example of this can be exposed by looking at the large-$r$ behavior of the mass function:
\be\label{masympt}
m(r) \sim M +\frac{M (\alpha -\gamma)}{r} + \order{1/r^2}\,,
\quad (r\to\infty)\,.
\ee
The second term in the above expansion could be interpreted as an electric charge, and could lead to a different quadrupole moment with respect to a Kerr BH.

The choice $\alpha=\gamma$ must be discarded as it forces the inner horizon to coincide with the Kerr one, and in turn implies a non-zero inner horizon surface gravity (actually the usual one for the Kerr geometry) making the conformal Kerr metric still unstable to mass inflation.

Nonetheless, we can introduce a parameter controlling the difference between the inner-horizon position in our geometry and in Kerr. This parameter will in turn control the difference $\alpha-\gamma$. Let us write then
\be
r_-\equiv\frac{a^2}{M+(1-e)\sqrt{M^2-a^2}}\,,
\label{eq:epar}
\ee
with $e\neq0$ and $e<2$ in order to ensure $0<r_-<r_+$. Further requiring the mass function to have no poles implies
\be\label{ebounds}
-3 - \frac{3M}{\sqrt{M^2 - a^2}} < e < 2\,,
\ee
where in the positive (negative) part of the interval $r_-$ is larger (smaller) than the Kerr inner horizon.

With the above choice, it follows that $\alpha-\gamma = \order{e^3}$ --- the same holds true for all the other coefficients in the large-$r$ expansion. This suggests that sizable deviations of $r_-$ from its Kerr value could translate into measurable differences in the value of the quadrupole moment, or in the periastron precession and the orbital frequency in a binary system \cite{Ryan:1995wh}. Such differences would all be $\order{e^3}$, which entails that values of $\abs{e}$ close to one or smaller might be phenomenologically favored; but the possible impact of $e$ on astrophysical observables certainly deserves further scrutiny, which we leave for the future.

Let us also note that, with the parametrization \eqref{eq:epar}, the mass function becomes $m(r) = M + \order{e^3}$ and in particular $m(r_+)=M$. This entails, among other things, that the outer-horizon angular velocity is the same as in Kerr, while its surface gravity is
\be
\kappa_+ = \frac{\partial_r \Delta(r_+)}{2(r_+^2+a^2)} = \kappa_+^\text{Kerr} + \order{e^3}\, .
\ee
Moreover, $e\to2$ is an extremal limit similar to $a\to M$, since in this limit $r_-\to r_+$ and $\kappa_+\to 0$.

Of course, different choices from Eq.~\eqref{eq:epar} for $r_-$ are in principle possible but they are strongly limited by a series of sanity requirements:
the inner horizon must lie within the outer horizon for all values of $a$;
$m(r)$ must go to $M$ asymptotically;
the denominator of $m(r)$ must have no zeros (for all $r>0$), that is $\gamma^2<4\mu$;
all the coefficients of $m(r)$ must be finite for all values of $a$;
the extremal limit $a \to M$ should remain thermodynamically unattainable and thus also the surface gravity of $r_+$ should become zero in this limit --- indeed this is possible only if $r_- \to r_+$ for $a \to M$.

\medskip

In conclusion, the complete form of our rotating ``inner-degenerate'' metric is
\be\label{eq:finalmetric}
ds^2=\frac{\Psi}{\Sigma} \left[-\left(1-\frac{2m(r)r}{\Sigma}\right) \dd{t}^2 - \frac{4a \, m(r) r \sin^2 \theta }{\Sigma} \dd{t} \dd{\phi} + \frac{\Sigma}{\Delta} \dd{r}^2 + \Sigma \dd{\theta}^2 + \frac{A \sin^2 \theta}{\Sigma} \dd{\phi}^2 \right],
\ee
with $m(r)$ given in \cref{m} and
\be
\Psi=\Sigma+\frac{b}{r^3}, \quad \Sigma = r^2 + a^2\cos^2\theta, \quad \Delta=r^2-2m(r)r +a^2, \quad A= (r^2+a^2)^2 -\Delta a^2 \sin^2 \theta\,,
\ee
where for the power law of $\Psi$ we have chosen the lowest integer that makes the curvature scalars continuous and finite (see \cref{s:conformal}).

Fixing $r_+={r_{+}}^{\rm Kerr}$ and choosing $r_-$ as in Eq.~\eqref{eq:epar}, this metric represents a family of stable, rotating, CTCs-free, regular solutions with two free parameters (beyond the usual spin one): the ``Kerr-deviation parameter'' $e$  and the ``regularization parameter" $b$.
Notice that for $a \to M$ the metric becomes conformal to the extremal Kerr, while for $a \to 0$ the metric becomes conformal to Schwarzschild.

\subsection{The rotating ``inner-degenerate'' RBH as a Kerr black hole mimicker\label{sec:mimicker}}

In this section we investigate the extent to which our metric \cref{eq:finalmetric} can mimic a Kerr BH: we will briefly describe the effective matter content; then the position of ergosurfaces; and finally the location of the light ring and the ISCOs.

\subsubsection{Effective matter content}

Clearly, the metric we are considering is not a vacuum solution of general relativity. The Einstein's equations can still be used to characterize the spacetime by interpreting the Einstein tensor $G^\mu_{\ \nu}= R^\mu_{\ \nu}-\frac{1}{2}R\, \delta^\mu_{\ \nu}$ as an effective stress--energy tensor. To properly characterize the effective matter content, one first needs to project the Einstein's tensor onto an orthonormal tetrad, e.g.~the one of Refs.~\cite{azreg-ainou_generating_2014, azreg-ainou_regular_2014, azreg-ainou_static_2014}.
The behavior of the orthonormal components close to spatial infinity is particularly relevant: since the spacetime is asymptotically flat, they must all tend to zero as $r\to \infty$, but they do so at different rates. In particular, the slowest decaying (non-zero) components are those on the diagonal, all the others being of higher order in powers of $1/r$. Such components, at infinity, are the effective energy density and pressures:\footnote{Technically, the energy density and pressures are defined in terms of the eigenvalues of the orthonormalized Einstein's tensor, when these are real. In asymptotically flat spacetimes, this procedure and the one presented in the text agree at leading order.}
\be\label{eq:rhop}
\rho = -p_r = p_\theta = p_\phi = -\frac{2M(\alpha-\gamma)}{r^4} +\order{1/r^5}\, .
\ee

Note that these quantities fall off quickly as $r\to \infty$, meaning that deviations from a vacuum solution (or alternatively from GR, given the supposed non-classical origin of the effective stress--energy tensor) are sizable only in a region close to the object. Moreover, they are $\order{e^3}$ and do not depend on $b$; the next-to-leading order $\order{1/r^5}$ also does not depend on $b$.
\Cref{eq:rhop} can lead to violations of the null energy condition (NEC), which requires $\rho + p_i \geq 0$, if $\alpha-\gamma >0$. When the NEC is violated, all the other classical energy conditions are violated too. When instead $\gamma > \alpha$, not only the null but also the weak (NEC + $\rho \geq 0$) and dominant ($\rho\geq \abs{p_i}$) energy conditions are met; the strong energy condition (NEC + $\rho+3p_i\geq 0$) instead is always violated.

Moving closer to $r=0$, the simple interpretation in terms of energy density and pressures is not always viable, since there are regions in which the Einstein's tensor cannot be diagonalized over the real numbers: in these regions, the effective matter content is of type IV in the Hawking--Ellis classification \cite{hawking_large_1973}. (The existence of these regions is entirely due to the presence of the conformal factor: when $\Psi=\Sigma$, the effective stress--energy tensor is of Hawking--Ellis type I for any $m(r)$.)

\subsubsection{Ergosurfaces}

The ergosurfaces are defined by the roots of $g_{tt}=0$, or equivalently of $r^2 - 2m(r)r + a^2\cos^2\theta = 0$, whose solution can be given in closed form.
Since the result is cumbersome, in \cref{fig:embed} we show the embedding in Euclidean space of the horizons and ergosufaces for some illustrative choice of the parameters.
The main difference with respect to a Kerr BH is the shape of the inner ergosurface around the poles: values of $e$ closer to the upper and lower bounds in \cref{ebounds} correspond to a more pronounced cuspid around the poles; for values of $e$ closer to the lower bound, the inner horizon and ergosurface move close and eventually touch also at the equator; for values of $e$ closer to the upper bound the horizons move closer as previously said. The conformal factor does not affect the ergosurfaces at all.

\begin{figure}[t]
\centering
\includegraphics[width=0.5\textwidth]{structure.pdf}
\caption{Embedding in Euclidean space of the horizons (green and blue surfaces) and ergosurfaces (red and purple surfaces) for $a/M = 0.95$ and $e=1$.}
\label{fig:embed}
\end{figure}

Finally, let us notice that, since with our choice $m(r_+)=M$, the textbook expression for the maximal efficiency of the Penrose process \cite{penrose_extraction_1971, bardeen_rotating_1972, wald_energy_1974, kovetz_efficiency_1975} seems to yield the same result as in Kerr:
\be
\eta_\text{max} = 1-\frac{2m(r_+)}{r_+} = 1-\frac{2M}{M+\sqrt{M^2-a^2}}\, .
\ee
Checking whether this is actually the case would require a more careful analysis of the motion of test particles in our spacetime --- an interesting question which however lies outside the scope of this work.

\subsubsection{Notable equatorial orbits}

In order to characterize the spacetime and its deviations away from Kerr from a phenomenological point of view, we compute the coordinate location of the light ring (LR) and the innermost stable circular orbit (ISCO).
We focus on the equatorial plane, where the radial motion is governed by the function (cf.\ \cref{rgeod})
\be
\mathcal{R} = E^2r^2(r^2+a^2) - r^2 L^2 + 2m(r) r (aE+L)^2 - \delta\Psi \Delta\,,
\ee
with $\delta=0$ or $1$ for null and timelike geodesics, respectively.
Circular orbits correspond to $\mathcal{R}=\mathcal{R}'=0$ and are stable if $\mathcal{R}''\leq 0$. Since the analytical expressions are not particularly illuminating, the values of $r_\text{LR}$ and $r_\text{ISCO}$ are computed numerically.

The location of the light ring, which is a null geodesic, does not depend on $\Psi$. Its fractional deviation from its Kerr analog is shown in \cref{lightring}, as a function of the spin, for some choices of the parameter $e$. The extrema and sign changes displayed by the curves of \cref{lightring} are ultimately determined by the behavior of the function $m(r)$ (and its derivative), which is not monotonic.

\begin{figure*}[h]
\centering
\includegraphics[width=0.7\textwidth]{Lightring.pdf}
\caption{Relative difference in the position of the light ring as a function of the spin between Kerr and our spacetime with $e=0.9$ (dotted lines), $e=1$ (dashed lines) and $e=1.1$ (continuous lines). We do not need to specify a value for the regularization parameter $b$ as null geodesics are insensitive to the conformal factor $\Psi$. While we display only values of $e $ near 1 (see previous discussion), corrections to the light-ring position actually grow very fast with $ e $ and they can be up to order $60\% $ for $e \rightarrow 2$.  Note also that the extrema and sign changes displayed by the curves are ultimately determined by the behavior of the function $m(r)$ (and its derivative), which is not monotonic. }
\label{lightring}
\end{figure*}

The analogous plot for the ISCO is reported in \cref{ISCO}.
\begin{figure*}[ht!]
\centering
\includegraphics[width=1.\textwidth]{ISCO.pdf}
\caption{Difference in the position of the ISCO between Kerr and our spacetime with $e=1$ and $b/M^5=0.8$ (continuous lines), $b/M^5=1$ (dashed lines) and $b/M^5=1.2$ (dotted lines). The prograde orbit, being in the more internal region of the spacetime where the conformal factor is greater, presents larger deviations, particularly for high spin.}
\label{ISCO}
\end{figure*}
Contrary to the previous case, $r_\text{ISCO}$ depends on $\Psi$, hence the curves in the figure correspond to specific choices of $b$. In fact, varying the parameter $b$ substantially affects the location of the ISCO, particularly for high spin.
The peculiar spike associated to prograde orbits and high spin, in particular, can be entirely explained in terms of the behavior of $\Psi$: since, as the spin increases, the prograde ISCO shrinks, $r_\text{ISCO}$ enters deeper into the region where $\Psi$ is markedly different from $r^2$. In order to further explore the parameter space in the high-spin regime, we set $a=0.998M$ (roughly the Thorne limit) and let the parameters vary in the ranges $b\in [0,1]$ and $e \in [-3-3M/\sqrt{M^2-a^2},2]$, thereby producing the contour plots of \cref{fig:ISCOspan}.
\begin{figure}[h!]
    \centering
    \subfloat[Prograde orbits.\label{fig:ISCOspanPR}]{\includegraphics[width=0.48\textwidth]{ISCO998PlotPr.pdf}}
    \hfill
    \subfloat[Retrograde orbits.\label{fig:ISCOspanRE}]{\includegraphics[width=0.48\textwidth]{ISCO998PlotRe.pdf}}\\
    \subfloat[Retrograde orbits, zoom to the region $\abs{e}<1$ of \cref{fig:ISCOspanPR}.\label{fig:ISCOspanZOOM}]{\includegraphics[width=0.48\textwidth]{ISCO998PlotReZoom.pdf}}
    \caption{Fractional deviation of the ISCO from its Kerr analog, computed as $r_\text{ISCO}/r_\text{ISCO}^\text{Kerr}-1\, .$ Spin $a/M=0.998$.}
    \label{fig:ISCOspan}
\end{figure}
Despite the much larger interval spanned by $e$, the gradient of the deviation is dominated by the $b$ component: this is clear for prograde orbits (\cref{fig:ISCOspanPR}), but is also true for retrograde orbits (\cref{fig:ISCOspanRE}) if $e$ is restricted to take reasonably small values as in \cref{fig:ISCOspanZOOM}. Note, however, that even for spins as high as $a=0.998M$, except for rather extreme values of the parameters, the ISCO moves less than a few percent in the prograde case and less than a few per mil in the retrograde case.

\section[Conclusions]{Conclusions\label{sec:conclusions}}

In this work we built and studied a new regular alternative to Kerr BHs that is stable under mass inflation. To construct it, we combined the two most common tools for regularization, in a novel way: we used a mass function to construct a degenerate (zero surface gravity) and thus stable inner horizon, and a conformal factor to regularize the singularity. In general, this procedure leads to a family of metrics, depending on the precise choice of the conformal factor and of $m(r)$.

We decided to focus on a particular form of the conformal factor that accomplishes the regularization in a minimal way and at the same time ensures the non-existence of CTCs and the separability of the equations of motion for test particles. With this choice, the curvature scalars are continuous and tend to zero on the would-be singularity\footnote{With a slightly different choice, however, the limit can also be made non-zero.} thereby solving a long-standing issue that affects many rotating RBHs. The regularization is controlled by a scale that we parametrize in terms of the quantity $b$, with dimensions $[M]^5$.

We further took an ansatz for $m(r)$ that is again minimal, in a suitable sense, and fixed the coordinate location of the outer horizon so that it coincides with its Kerr analog.
The resulting mass function can be expressed entirely in terms of the coordinate location of the inner horizon, whose deviation from that of Kerr is encoded by the dimensionless quantity $e$.
In the limit $e \rightarrow 0$ we obtain the conformal Kerr metric that, though regular, is characterized by the usual surface gravity at the inner horizon (as it should be given that the surface gravity is conformal invariant) and hence is again unstable under mass inflation. However, it is important to notice that our metric cannot indefinitely deviate from the conformal Kerr one since the deviation parameter $e$ must lie in a specific interval in order for the mass function to be everywhere finite and for the horizons to be well ordered ($0<r_-<r_+$).

Our metric thus depends on a total of four real parameters: the ADM mass $M$, the spin $a$, the regularization parameter $b$ and the deviation parameter $e$. The two additional parameters $b$ and $e$ can be constrained by observations, at least in principle. In particular, $e$ enters at low order in the PPN expansion of this object's gravitational field and thus influences its multipolar structure; moreover, it affects the orbits of massless test particles and therefore shifts the position of the light ring. Finally $e$ and $b$ both affect the motion of massive test particles, with $b$ having the dominant effect on the location of the ISCO (at least when $e$ is taken to vary in reasonably small ranges) especially at high spins.

Let us stress the remarkable relevance of this fact: given that the regularization and Kerr-deviation parameters might be directly related to quantum-gravitational effects, the possibility to constrain them via observations on the exterior geometry of the black hole is further evidence that a new window for quantum-gravity phenomenology might be opening via astrophysical observations.

In conclusion, the rotating RBH geometry proposed here is characterized by a wealth of physically desirable features that make it a plausible candidate for the end point of a gravitational collapse to be contrasted with the Kerr geometry to which it can be made arbitrarily close. Constraining the values of the parameters $e$ and $b$ with astrophysical data now becomes an interesting challenge and might represent a new window over testing possible quantum-gravitational scenarios. We leave this for future work.

\section*{Acknowledgments}
The authors acknowledge funding from the Italian Ministry of Education and Scientific Research (MIUR) under the grant PRIN MIUR 2017-MB8AEZ.

\bibliographystyle{JHEP.bst}