\section{Introduction}

Einstein's theory of general relativity (GR) predicts the existence of compact dark objects such as black holes (BHs). Under the assumption of strong energy condition, collapsing matter in GR leads to spacetime singularities \cite{hawking_ellis_1973, Wald:1984rg, PhysRevLett.14.57}. These set of results are known as Hawking-Penrose singularity theorems that suggest a generic existence of singularities in BH solutions of GR. At the location of a spacetime singularity, various curvature scalars diverge, indicating breakdown of the underlying gravity theory and its predictive power. It is widely believed that at small enough length scales, general relativity would be replaced by a quantum theory of gravity, which will ultimately resolve this issue. However, in the absence of a fully consistent quantum gravity framework, one may consider a phenomenological approach in which the known BH solutions are modified to regularize the central singularity with a non-singular core. In general, such approaches are highly non-unique and theory agnostic in nature. Several interesting regularized BH solutions (both static and rotating) have been proposed over years \cite{Dymnikova:1992ux, Modesto:2008im, Bonanno:2000ep, Hayward:2005gi, Platania:2019kyx, Bambi:2013ufa, Azreg-Ainou:2014pra, PhysRevD.89.104017, Eichhorn:2021etc, Delaporte:2022acp}. Such solutions contain one or more additional parameters over the usual hairs, namely the mass and spin of a Kerr BH. Usually these regularized BHs possess Cauchy horizons with nonzero surface gravity causing an exponential growth of perturbations at its vicinity known as the mass inflation singularity \cite{PhysRevD.101.084047, Carballo-Rubio:2019nel, PhysRevD.41.1796, PhysRevLett.67.789, Bhattacharjee:2016zof}. However, the recently proposed regularized BH solutions in Ref.\cite{Carballo-Rubio:2022kad} (spherically symmetric) and in Ref.\cite{Franzin:2022wai} (rotating) circumvent both the issues of singular core and mass inflation instability. Due to its astrophysical relevance, we focus on this regularized rotating BH geometry which we refer as regularized stable Kerr BH and study various  properties in the subsequent sections.\\

\noindent
As discussed previously, the approach of regularizing the singularity is highly non-unique and as a result, it leads to a plethora of regularized BHs of different types. In the absence of a unique approach, an effective way to check the viability of these models could be to perform various consistency tests. In this paper, we explore the regularized BH solution presented in Ref.\cite{Franzin:2022wai} and attempt to constrain the deviations from Kerr BH by theoretical and observational studies. The central theme of this article is of two folds. Firstly, inspired by seminal works of Regge, Wheerle and Vishweswara, we check the stability of such solutions by analysing the scalar quasi-normal mode (QNM) spectrum \cite{PhysRev.108.1063, PhysRevLett.24.737, Vishveshwara:1970zz, PhysRevLett.29.1114, Leaver:1985ax, PhysRevD.34.384}. Ensuring the stability of a BH solution is essential not only from a mathematical perspective but also due to its astrophysical relevance. For a comprehensive review on the subject of BH QNMs and stability, we refer the reader to Refs. \cite{Berti:2009kk, Kokkotas:1999bd, Konoplya:2011qq}.
If the regularized BH under consideration is found to be unstable in certain region of the parameter space, one can rule out those regions. This is a very powerful method to check the viability of any classical BH solution and constrain any additional parameter present in the theory. Conversely, any well behaved BH solution must also be stable against small perturbation in the entire range of the parameter space. One of the aims of this work is to look for any instabilities present in this BH spacetime. Our study suggests that the regularized BH under consideration is a stable configuration against scalar perturbation. As a result, to constrain possible deviation from Kerr solution, we resort to observational techniques, which corresponds to the second part of our work.
\\

\noindent
The bending of light by a strong gravitational field is one of the fascinating features of the spacetime curvature, which gives rise to the existence of spherical/circular null orbits around BHs. As far as observational implications are concerned, these null orbits further lead to BH shadow, a dark patch around BHs seen by a distant observer. The recent observations of $M87^*$ and $Sgr\, A^*$ BH shadows by the Event Horizon Telescope (EHT) collaboration~\cite{EventHorizonTelescope:2019dse, EventHorizonTelescope:2019ths, EventHorizonTelescope:2019pgp, EventHorizonTelescope:2019ggy, EventHorizonTelescope:2022xnr, EventHorizonTelescope:2022xqj} are the first direct detection of isolated BHs, providing us with an unprecedented platform to test the viability of modified Kerr solutions. However, since the EHT observations are in strong agreement with the Kerr paradigm, any deviation from Kerr present in the theory are highly constrained~\cite{Jusufi:2022loj, Hendi:2022qgi, Banerjee:2022iok, KumarWalia:2022aop, Banerjee:2022bxg, Banerjee:2022jog, Shaikh:2022ivr, Li:2022eue}. To that end, we compute the shadow cast by the regularized stable Kerr BH and compare the angular shadow diameter with the EHT observations for both $M87^*$ and $Sgr\, A^*$ shadows and obtain upper bounds on the additional parameters.\\

\noindent
The rest of the article is organized as follows: In \ref{Sec_2} we briefly review the regularized stable Kerr solution and discuss some of its properties. In \ref{Sec_3} we discuss the consistency conditions on the mass profile suggested in Ref. \cite{Franzin:2022wai}. This in some sense represents a study similar to the weak cosmic censorship~\cite{WALD1974548, Hubeny:1998ga, Sorce:2017dst}. The computation of scalar quasi-normal modes and stability analysis is presented in \ref{Sec_4}. In \ref{Sec_5}, we study the shadow cast by this BH and constrain various model parameters requiring consistency with EHT observations. Finally, we conclude with a brief discussion of main results and possible future extensions of our work.\\

\noindent
Notations and conventions: In this work we follow the mostly positive signature convention $(-, +, +, +)$ for the metric. Indices referring the four dimensional spacetime are represented by Greek letters. Unless specified otherwise, we also set the fundamental constants to unity, i.e., $c = 1 = G$.

\section{Constructing the metric} \label{Sec_2}
An astrophysical BH with Arnowitt-Deser-Misner (ADM) mass $M$ and spin $a<M$ is well described by the Kerr metric, which is an exact solution of Einstein field equations. However, it is infected with three major issues: (i) the central singularity where various curvature scalars blow up and GR looses its predictive power, (ii) the instability originating due to the mass inflation at the Cauchy horizon, and (iii) the existence of closed timelike curves in the spacetime. Over years, several attempts have been made to get rid of one such issue at a time by introducing some modified phenomenological metrics. However, in a recent work \cite{Franzin:2022wai}, the authors proposed the following ``inner-degenerate'' regularized stable Kerr metric that is free from all these issues :

\begin{align} \label{g}
ds^2 = \mathcal{C}(r,\theta) \Bigg[-\left(1-\frac{2\, m(r)\, r}{\Sigma(r,\theta)}\right) \mathrm{d} t^{2}-\frac{4\, a\, m(r)\, r\, \sin ^{2} \theta}{\Sigma(r,\theta)}\, \mathrm{d} t\, \mathrm{d} \phi\, +\, &\frac{\Sigma(r,\theta)}{\Delta(r)}\, \mathrm{d} r^{2}\,  + \\ \nonumber
&\Sigma(r,\theta)\, \mathrm{d} \theta^{2}+\frac{A(r,\theta) \sin ^{2} \theta}{\Sigma(r,\theta)}\, \mathrm{d} \phi^{2}\Bigg] ,
\end{align}

\noindent
with $\Sigma(r,\theta)=r^{2}+a^{2} \cos ^{2} \theta,\, \Delta(r)=r^{2}-2\, m(r)\, r+a^{2},\, A(r,\theta)=\left(r^{2}+a^{2}\right)^{2}-\Delta\,  a^{2} \sin ^{2} \theta$, and $m(r)$ is the mass profile introduced to remove the mass inflation instability at the Cauchy horizon $r_{-}$. Note that the locations of the horizons $r_{+}$ and $r_{-} < r_{+}$ are given by the real positive roots of $\Delta(r) = 0$. On the other hand, the conformal factor $\mathcal{C}(r,\theta)$ is there to tame the central curvature singularity of the Kerr metric.\\

\noindent
We are interested in investigating the QNM stability and the shadow observables of this modified Kerr BH and put constraints on various parameters of the metric. To keep our analysis as general as possible, we shall work with the mass profile $m(r)$ and the conformal factor $\mathcal{C}(r,\theta)$ unspecified, and only fix them when necessary. However, in the same spirit of Ref.\cite{Franzin:2022wai}, we impose two minimal criteria on the conformal factor in order to avoid any additional singularity: $\mathcal{C}(r,\theta) > 0$ everywhere, and to keep the spacetime indistinguishable from Kerr metric to a distant observer: $\mathcal{C}(r,\theta) \to 1 + \mathcal{O}(r^{-n})$ with $n \geq 2$. The second condition implies that the ADM mass and the specific angular momentum of the black hole are $M = \text{lim}_{r \to \infty} m(r)$ and $a$, respectively. \\

\noindent
Also, we follow the same efficient way of Ref.\cite{Franzin:2022wai} for parameterizing $m(r) = F(r; r_{\pm})$ in terms of the horizons $r_{+} = M+\sqrt{M^2-a^2}$ and $r_{-}=a^2\, \left[ M+(1-e)\sqrt{M^2-a^2}\right]^{-1}$. Then, for a given value of $e$, the extremal limit $r_{-} \to r_{+}$ remains same as the Kerr case (namely $a=M$), which we shall refer as ``a-extremality". Moreover, to ensure $0< r_{-} < r_{+}$, we need to choose the parameter $e$ such that $0 \neq e < 2$. The metric in Eq.(\ref{g}) reduces to a conformal Kerr BH by setting $e=0$. On the other hand, $e \to 2$ is the ``e-extremal" limit for any fixed $a \leq M$. In the next section, we discuss the consistency of this parametrization as the BH absorbs test particles.

\section{Consistency criterion on the mass profile} \label{Sec_3}
As discussed in the previous section, we want to impose a specific parameterization on the mass profile $m(r) = F(r; r_{\pm})$. This is only possible upto the a-extremal ($a=M$) configuration of the BH given by Eq.(\ref{g}). Now consider the process in which such a BH with $ a \leq M$ absorbs a test particle with energy $E << M$ and angular momentum $L << M^2$. This may change both the ADM mass and the angular momentum of the initial black hole in such a way that the final configuration has $(a\, M + L) > (M+E)^2$ and $r_{\pm}$ becomes complex. As a result, occurrence of such a process puts the viability of the aforesaid parameterization in question. Thus, we want to put some constraints on $m(r)$ that can prevent such difficulties.\\

\noindent
The above inequality holds if the particle's energy satisfy the upper bound $E < E_{max} = aL/(M^2+a^2)$. On the other hand, such a particle will be absorbed by the BH if it does not come across any turning points before crossing the event horizon at $r_{+}$. Neglecting back reactions and using geodesic equations, the entering condition can be translated to a lower bound on the test particle's energy: $E \geq E_{min} = \left(a\, \sigma\, L \right)/\left[2 M r_+ + a^2(\sigma-1)\, sin^2\theta \right]$, where $\sigma = m(r_+)/M$. Thus, to prevent such a process we must require, $E_{min} \geq E_{max}$. It implies the consistency condition: $\sigma \geq r_+^2/M^2$, that is strictest for all allowed values of $\theta \in [0,\pi)$. Therefore, if the mass profile satisfy the inequality $m(r_+) \geq r_+^2/M$, the BH will not absorb those particle that could have made the parametrization inconsistent in the final configuration.\\

\noindent
For the mass profile given in Ref.\cite{Franzin:2022wai}, this consistency criteria is not satisfied if the initial BH is sub-extremal with $a<M$. Therefore, one may choose a different mass profile which satisfies the consistency criteria under test particle approximation. However, in the same spirit of \cite{Sorce:2017dst}, it is reasonable to expect that the situation improves and the apparent inconsistency in the parametrization of Ref.\cite{Franzin:2022wai} goes away as one properly takes into account the back reaction effects. In contrast, if we start with an a-extremal BH, the consistency condition becomes $m(r_+) \geq M$, which is satisfied by the choice of $m(r)$ given in Ref.\cite{Franzin:2022wai} even under test particle approximation.\\

\noindent
Though it looks similar to the calculation done to check the validity of the weak cosmic censorship (WCC) conjecture for a Kerr BH~\cite{WALD1974548, Hubeny:1998ga, Sorce:2017dst}, there is a subtle difference. WCC protects the predictive power the underlying gravity theory, it is necessary to censor out those BH solutions in which there is a possibility of destroying the horizon and exposing the central singularity by test particle absorption. However, for the metric in \ref{g}, there is no such curvature singularity inside the event horizon. Thus, the inequality derived above is not a strict necessity, it is just to make the aforesaid parameterization consistent. If we use different parameterization, we would have got a different criteria on $m(r)$. However, for the purpose of this paper, we shall stick to the same parameterization for the mass profile as given in Ref.\cite{Franzin:2022wai}.

\section{Scalar QNM spectrum} \label{Sec_4}
Note that, so far we have not specified any particular form for the conformal factor $\mathcal{C}(r,\theta)$ or the mass profile $m(r)$ except putting some conditions on them. However, to calculate the QNM spectrum and shadow of this BH in later sections, we have to specify their functional forms. For the mass profile, we consider the same form as given in Eqs.(28-32) of Ref.\cite{Franzin:2022wai}.\\

\noindent
However, we choose a different conformal factor $\mathcal{C}(r,\theta) = 1+b/r^{2z}$ than what are proposed in Ref.\cite{Franzin:2022wai}. This choice for $\mathcal{C}(r,\theta)$ is motivated by two facts. Firstly, it is much simpler to work with and performs the job that it is designed for, namely it regularizes the central curvature singularity of the Kerr metric if $z \geq 2$. In addition, such a choice of the conformal factor makes the scalar wave equation separable, which we shall exploit in finding the QNM spectrum. However, for the purpose of this paper, we shall consider the simplest case $z=2$.\\

\noindent
 Consider that a massless scalar field $\Phi(x^\mu)$ perturbs the the stationary and axisymmetric BH geometry given in \ref{g} with $\mathcal{C}(r,\theta) =\psi(r)= 1+b/r^{4}$, and $m(r)$ as prescribed in Ref.\cite{Franzin:2022wai}. Owing to the symmetries of the background spacetime, we can decompose the field as $\Phi(x^\mu)=  R_{lm}(r)\, S_{lm}(\theta) e^{-i\,( \omega\, t-\, m\, \phi)}$, where $S_{lm}(\theta)$ is the spheroidal harmonics which satisfies the angular Teukolsky equation \cite{1973ApJ...185..635T},

\begin{equation}\label{ang_Teuk}
\left[\frac{d}{dx}\left(1-x^2\right)\frac{d}{dx}+\left(E_{lm}+a^2\, \omega^2\, x^2-\frac{m^2}{1-x^2}\right)\right]S_{lm}(x)=0,
\end{equation}

\noindent
with $x=\cos\theta$ and $E_{lm}$ is the separation constant. Note that even in the presence of conformal factor $\psi(r)$, we can cast the radial equation into the standard Teukolsky form \cite{1973ApJ...185..635T}:

\begin{equation} \label{rad_teuk}
	\left[\frac{d}{dr}\, \Big(\Delta _1\, \frac{d}{dr} \Big)+	 \left(\frac{  K^2}{\Delta _1}-\lambda_{lm}\, \psi \right)\right] R_{lm}(r)=0 ,

\end{equation}

\noindent
by introducing a function $\Delta_1 (r) =\psi(r)\,  \Delta(r)$. Here, $\lambda_{lm}=E_{lm}+a^2\, \omega ^2-2\, a\, m\, \omega $ and $K=\psi(r)\left[\omega  \left(a^2+r^2\right)-a\, m\right]$. The QNMs are the solution of the above equations with ingoing boundary conditions at the event horizon and outgoing boundary conditions at infinity,
\begin{align}
R(r\rightarrow r_h)\sim e^{-i\omega r_*}\qquad \text{and} \qquad R(r\rightarrow \infty)\sim e^{i\omega r_*}\ .
\end{align}

\noindent
The radial function $R_{lm}$, subject to the above QNM boundary conditions, can be expressed as follows \cite{Rahman:2019uwf, Rahman:2020guv},
\begin{equation}\label{frobenious}
    R_{lm}=\frac{u_{lm}(z)}{r}\, e^{i\, B_1(r)}\, z^{-2\, i\,  B_2(r_+)}\ .
\end{equation}

\noindent
Here, $z=(r-r_+)/(r-r_-)$, $B_1(r)=K(r)/\Delta_1(r)$, and  $B_2(r_+)=K(r_+)/\Delta'_1(r_+)$. Replacing the above expression in \ref{rad_teuk}, we obtain an equation for $u_{lm}(z)$. Note that both the radial and angular equation depends on $\omega$ and $\lambda_{lm}$. Thus, to obtain the desired QNMs, we need to solve these two equations simultaneously. We use the \texttt{Black Hole Perturbation Toolkit} \cite{BHPToolkit} to solve the angular Teukolsky equation \ref{ang_Teuk}, and  \texttt{QNMSpectral} package \cite{Jansen:2017oag} to solve the equation for $u_{lm}(z)$ numerically. For given values of $(a, b, e, \ell)$, the output of this analysis is the set of QNM frequencies and separation constants. The real and imaginary part of the scalar QNMs are shown in \ref{fig_stability_Scalar} for various choice of parameters $(e, b, a)$ for the case of $\ell = 0$ fundamental mode. As it can be seen from \ref{fig_stability_Scalar}, the imaginary part of QNM remains negative for the entire range of the additional parameters $e$ and $b$, indicating the regularized black hole to be stable in the entire parameter space.
Furthermore, we interestingly find that with increasing values of the conformal parameter ($b$), the imaginary part of the scalar QNM frequency decreases for the entire range of $e$. As the conformal factor is introduced to tame the central curvature singularity, we can conclude that BHs with regular cores are more stable as compared to the singular ones. As represented in \ref{fig_stability_Scalar}, this behaviour is seen for both rotating ($a\neq 0$) as well as non-rotating case ($a=0$). Also,  \ref{fig_stability_Scalar} clearly demonstrates the dependence of QNMs on the parameter $e$. For a given value of $b$, the decay time of the perturbation increases (or, stays the same for non-rotating case) with increasing values of $e$. \\

\begin{figure}

	\centering
	\minipage{0.48\textwidth}
	\includegraphics[width=\linewidth]{ReQNMa0.pdf}
	\endminipage\hfill

	\minipage{0.48\textwidth}
	\includegraphics[width=\linewidth]{ImQNMa0.pdf}
	\endminipage
	\hfill

	\minipage{0.48\textwidth}
	\includegraphics[width=\linewidth]{ReQNMa04.pdf}
	\endminipage\hfill

	\minipage{0.48\textwidth}
	\includegraphics[width=\linewidth]{ImQNMa04.pdf}
	\endminipage
	\hfill

	\minipage{0.48\textwidth}
	\includegraphics[width=\linewidth]{ReQNMa08.pdf}
	\endminipage\hfill

	\minipage{0.48\textwidth}
	\includegraphics[width=\linewidth]{ImQNMa08.pdf}
	\endminipage
	\caption{In this figure, we illustrate the variation of both real (left panel) and imaginary parts (right panel) of scalar QNMs for $\ell = 0$ mode. Refer to text for more discussion.}\label{fig_stability_Scalar}
\end{figure}

\noindent
Because of the complicated structure of the perturbation equations in a general curved spacetime, it is often difficult to obtain any analytical results. In such a scenario, one resorts to numerical techniques, as we did previously. However, in some limited cases, it is indeed possible to solve the perturbation equation analytically. One such approximation is the eikonal limit or geometric optics limit. We have cross-verified our results for eikonal modes (presented in \ref{fig_Lyapunov})  with analytical Lyapunov methods \cite{Cardoso:2008bp}, which are in excellent agreement with each other. In the eikonal limit ($\ell >> 1$ keeping $ m = \ell)$, the QNMs can be interpreted in terms of various properties of the equatorial null geodesics trapped at the photon sphere at $r = r_c$. The real parts of the QNMs are proportional to the angular frequency $\Omega_c$ of these null geodesics at the photon sphere. In contrast, the imaginary part is related to the Lyapunov exponent $\Lambda_{\textrm{Lya}}$, which dictates the instability timescale of the equatorial null geodesics at the photon sphere. Thus, the eikonal QNMs take the form \cite{Cardoso:2008bp}:

\begin{equation}\label{lya}
    \omega_{n\ell}=\ell\, \Omega_c-i\,  (n+\frac{1}{2})\, |\Lambda_{\textrm{Lya}}|\, ,
\end{equation}

\noindent
where $n$ is the overtone number. The expression for the $\Omega_c$ and $\Lambda_{\textrm{Lya}}$ can be written in terms of the metric coefficients as \cite{Rahman:2018oso},
\begin{equation}
    \begin{aligned}
    & \Omega_c =\frac{g_{t\phi}+g_{tt}\, y}{g_{\phi\phi}+g_{t\phi}\, y}\Bigg\rvert_{r = r_c}\, ,\\
    & \Lambda_{\textrm{Lya}}^2 =\frac{g_{t\phi}^2-g_{tt}\, g_{\phi\phi}}{g_{rr}}\left\{\frac{g_{tt}\, g''_{\phi\phi}-g''_{tt}\, g_{\phi\phi}-2\, y \left(g''_{tt}\, g_{t\phi}-g_{tt}\, g''_{t\phi}\right)}{2\, g_{tt}\left(g_{\phi\phi}+y\,  g_{t\phi}\right)^2}\right\}\Bigg\rvert_{r = r_c}\, ,
    \end{aligned}
\end{equation}
where $g_{\mu\nu}$ denotes the metric coefficient given in \ref{g}, prime ($'$) denotes derivative with respect to $r$ and $y$ is the impact factor, the expression of which can be written as follows \cite{Rahman:2018oso},

\begin{equation}
    y=-\frac{g_{t\phi}}{g_{tt}}+\sqrt{\left(\frac{g_{t\phi}}{g_{tt}}\right)^2-\frac{g_{\phi\phi}}{g_{tt}}}\Bigg\rvert_{r = r_c}\, .
\end{equation}

\begin{figure}[ht]

	\centering
	\minipage{0.48\textwidth}
	\includegraphics[width=\linewidth]{reLypunov.pdf}
	\endminipage\hfill

	\minipage{0.48\textwidth}
	\includegraphics[width=\linewidth]{imLypunov.pdf}
	\endminipage
	\caption{In this figure we represent the real and imaginary part of the QNM computed analytically in the eikonal limit ($\ell = 10$) for the parameters indicated in the plot. }\label{fig_Lyapunov}
\end{figure}

\noindent
To summarize, this section considers the scalar perturbation of the regularized stable Kerr metric and computed the QNM spectrum both analytically and numerically. We found that both for small and large values of $\ell$, the imaginary part of the QNM modes are always negative, implying the stability of the underlying BH in the entire parameter space of $b$ and $e$. Thus, it is not possible to constrain these extra parameters from the stability analysis. As a result, we resort to shadow observations in order to obtain any possible bounds on the additional parameters, which is discussed in the next section.

\section{Bounds from Shadow Observations} \label{Sec_5}
BH shadow provides us with a powerful tool to constrain various parameters of the metric under consideration. In this section, we shall study the shadow cast by the regularized stable Kerr BH and compare it with the EHT observations for $M87^*$ and $Sgr\, A^*$ shadows \cite{EventHorizonTelescope:2019dse, EventHorizonTelescope:2019pgp, EventHorizonTelescope:2022xnr}. The goal is to restrict the parameter space of these regularized BHs such that their shadow sizes match with EHT observations.\\

\noindent
Although we have specified a form of the conformal factor in the previous section, we shall see that for the study of the BH shadow, no particular form for $\mathcal{C}(r,\theta)$ needs to be assumed. It is because for any choice of the conformal factor $\mathcal{C}(r,\theta)$, null geodesic equations are separable and take the following form,

\begin{align} \label{geod}
\mathcal{C}^2(r,\theta)\, \Sigma^2(r,\theta)\, \dot{r}^2 = -Q\, \Delta(r) + \Big[ (r^2+a^2)\, E - a\, L_z \Big]^2 := R(r), \\
\text{and}\quad \mathcal{C}^2(r,\theta)\, \Sigma^2(r,\theta)\, \dot{\theta}^2 = Q - \Big[ \frac{L_z}{\text{sin}\theta} - a\, E\, \text{sin}\theta \Big]^2 := \Theta(\theta)\ .
\end{align}

\noindent
Here, $E$ and $L_z$ are two constants of motion associated to photon's energy, and the $z$-components of the angular momentum. The separability implies the existence of another constant of motion $Q$, namely the Carter constant. Now, to construct the BH shadow, we require the presence of a photon region filled with unstable spherical null geodesics satisfying $R(r_p) = R'(r_p) = 0$ with $R''(r_p) < 0$. Solving these two equations simultaneously, we obtain two useful quantities,

\begin{equation}
\xi := \frac{L_z}{E} = \frac{1}{a} \Big[ r_p^2 +a^2 -\frac{2\, r_p\, \Delta(r)}{r_p-M(r_p)} \Big],\, \, \text{and} \quad \eta := \frac{Q}{E^2} =  \frac{4\, r_p\, \Delta(r_p)}{\left[ r_p - M(r_p) \right]^2}\, ,
\end{equation}

\noindent
where $M(r) = m(r) + r\, m'(r)$ with the mass profile $m(r)$ given in Ref.\cite{Franzin:2022wai}. However, one can not put any arbitrary values of $r_p$ in the above set of equations. It is because the spherical null geodesics are confined only in a finite portion of the spacetime outside the BH known as the photon region. This region is represented by the inequality $\Theta(\theta, \xi, \eta) \geq 0$. At any fixed value of $\theta$, this inequality is satisfied when $r_p$ varies between two end points $\left[r^p_{-}, r^p_{+}\right]$ where equality holds true. \\

\noindent
One can use $(\xi, \eta)$ to span the observer's sky with suitably defined celestial coordinates $(\alpha, \beta)$. For a distant observer, these coordinates are given in terms of the inclination angle $\theta_i$:

\begin{equation} \label{celes}
\alpha = - \xi\, \text{cosec}\theta_{i},\, \, \text{and}\, \, \, \, \beta = \pm \sqrt{\eta - \Big(\xi\, \text{cosec}\theta_{i} - a\, \text{sin}\theta_i\Big)^2}\ .
\end{equation}

\begin{figure}[ht]
\begin{center}
\includegraphics[width=0.4\textwidth]{Shadow_img.pdf}
\caption {This figure shows $M87^*$ BH shadow ($M=1$) for $a_* = 0.5$ and different values of $e$. Beyond a critical value $e > e_c \approx 1.75$, the shadow size becomes bigger than $45\, \mu as$ and becomes inconsistent with EHT observations\cite{EventHorizonTelescope:2019dse}.}
\label{sh}
\end{center}
\end{figure}

\noindent
For a Schwarzschild BH, the shadow boundary is a circle of radius $R_{sh} = 3\sqrt{3}\, M$. In contrast, rotation makes the size of a Kerr shadow varies away from $R_{sh}$. Usually, such variation is negligible if the BH is not rotating very rapidly. However, in our case, the shadow of a regularized stable Kerr BH depends on an extra parameter $e$. For a fixed values of $(M, a)$, the BH shadow will have different angular sizes as $e$ varies in the range $( -3 - 3/ \sqrt{1-a_*^2},\, 2]$ (with $a_* := a/M$) as prescribed in Ref.\cite{Franzin:2022wai}. For example, see \ref{sh} in which $M87^*$ shadow sizes (in $M=1$ unit) are shown for $a_* = 0.5$ and different choices of $e$ values.\\

\begin{center}
\begin{table}[ht!]
\begin{tabular}{|c|c|c|c|c|c|c|}
\Xhline{1\arrayrulewidth}
Central BHs  &    $M\, ( \times 10^6\, M_\odot)$  &  $a_*$   &  $D$                            &  $\theta_i$                        &  $\theta_d$                        &   References
\\\Xhline{1\arrayrulewidth}
$M87^*$	       &    $6.5 \times 10^3 $   &	 $0.5-0.94$	   &	$16.8\, \text{Mpc}$	&	$ 17 ^{\circ}$               & 	$42 \pm 3\, \mu$as        &   \cite{EventHorizonTelescope:2019dse, EventHorizonTelescope:2019pgp}
\\\Xhline{1\arrayrulewidth}
$Sgr\, A^*$      &    $4.0$                       &	 $0.44$	   &	$8\, \text{kpc}$	        &	$< 50 ^{\circ}$             & 	$51.8 \pm 2.3\, \mu$as  &   \cite{EventHorizonTelescope:2022xnr, Kato_2010}
\\\Xhline{1\arrayrulewidth}
\end{tabular}
\caption{EHT observational data for $M87^*$ and $Sgr\, A^*$ shadows.}
\label{T1}
\end{table}
\end{center}

\noindent
Thus, for fixed values of $(M, a_*, \theta_i, D)$, there can be a range of $e$-values for which the shadow will be inconsistent with the angular size $\theta_d$ \cite{Kumar:2020owy} as observed by EHT, check \ref{T1}. For this purpose, let us first consider the $M87^*$ shadow. In this case, since the angular momentum parameter consistent with observation can take a range of values $a_* \in [0.5,0.94]$, we have shown below the variation in shadow size using contour plots \ref{M87plt}. From these plots, it is very easy to identify the disallowed regions which are inconsistent with EHT observations. For example, if by some independent measurement $a_* =0.6$ is fixed for $M87^*$ BH, then all values of $e$ lies in the range $( -3 - 3/ \sqrt{1-a_*^2}, 1.1799] \cup (1.7105,2]$ are disallowed by EHT observations.\\

\begin{figure}[ht]
  \subfloat[For $-5 \leq e \leq 0$]{
	\begin{minipage}[c][1\width]{
	   0.3\textwidth}
	   \centering
	   \includegraphics[width=1.1\textwidth]{M87_1.pdf}
	\end{minipage}}
 \hfill
  \subfloat[For $0 \leq e \leq 1.68$]{
	\begin{minipage}[c][1\width]{
	   0.3\textwidth}
	   \centering
	   \includegraphics[width=1.1\textwidth]{M87_2.pdf}
	\end{minipage}}
 \hfill
  \subfloat[For $1.68 \leq e \leq 2$]{
	\begin{minipage}[c][1\width]{
	   0.3\textwidth}
	   \centering
	   \includegraphics[width=1.1\textwidth]{M87_3.pdf}
	\end{minipage}}
\caption{Contour plots for angular size of $M87^*$ shadow. Values of $(a_*,\, e)$ for which the shadow size is less than $39\, \mu$as or greater than $45\, \mu$as are disallowed.}
\label{M87plt}
\end{figure}

\noindent
For the case of $Sgr\, A^*$, a unique spin parameter $a_* \approx 0.44$ can be fixed by studying the multiple quasi-periodic oscillations of radio emissions \cite{Miyoshi, Kato_2010}. However, EHT shadow observation is consistent with all inclination angles that are less than $50 ^{\circ}$. Therefore, we can show the variation in shadow size using contour plots for different choices of $(\theta_i, e)$, see \ref{SgrAplt}. It is interesting to note that the $Sgr\, A^*$ shadow size is consistent with EHT observations for all values of $e \in \Big(-3 - 3/ \sqrt{1-(0.44)^2} \approx -6.341, 1.7\Big]$ with $\theta_i < 50^{\circ}$ and $a_* \approx 0.44$. Whereas all values of $e > 1.7$ are disallowed by EHT observations \cite{EventHorizonTelescope:2022xnr}.\\

\begin{figure}[ht]
  \subfloat[For $-5 \leq e \leq 0$]{
	\begin{minipage}[c][1\width]{
	   0.3\textwidth}
	   \centering
	   \includegraphics[width=1.1\textwidth]{SgrA_1.pdf}

	\end{minipage}}
 \hfill
  \subfloat[For $0 \leq e \leq 1.68$]{
	\begin{minipage}[c][1\width]{
	   0.3\textwidth}
	   \centering
	   \includegraphics[width=1.1\textwidth]{SgrA_2.pdf}
	\end{minipage}}
 \hfill
  \subfloat[For $1.68 \leq e \leq 2$]{
	\begin{minipage}[c][1\width]{
	   0.3\textwidth}
	   \centering
	   \includegraphics[width=1.1\textwidth]{SgrA_3.pdf}
	\end{minipage}}
\caption{Contour plots for angular size of $Sgr\, A^*$ shadow. Values of $(\theta_i,\, e)$ for which the shadow size is less than $49.5\, \mu$as or greater than $54.1\, \mu$as are disallowed.}
\label{SgrAplt}
\end{figure}

\section{Discussion and Conclusion}
BHs are inevitable consequence of Einstein's field equations. The recent observation of shadow by the EHT along with the detection of gravitational wave by the LIGO-Virgo-KAGRA (LVK) collaboration provide direct and compelling evidence for the existence of such objects in nature. These astrophysical observations are in consistent with the Kerr paradigm, namely the central compact object being a Kerr BH completely described by the mass (M) and spin (a). However, as discussed in the introduction, the Kerr solution suffers from various pathologies such as the existence of spacetime singularities and mass inflation singularity at the Cauchy horizon. It is widely believed that such issues are consequence of the limitations of classical GR and will be resolved by a full quantum gravity formalism. However, in the absence of a consistent quantum theory of gravity, one may try to tackle the these problems using phenomenological models. One such model of regularized stable Kerr BH has been recently developed in Refs. \cite{Carballo-Rubio:2022kad, Franzin:2022wai}. In addition to mass and spin, these BH solutions depend on two additional parameters: a conformal factor $b$, and a deviation parameter $e$ that regularizes and stabilizes the Kerr solution. In this work, we have studied some observational signatures of these parameters on BH QNMs and shadow observations. \\

\noindent
Astrophysical BHs are expected to be stable under perturbations. As a result, stability turns out to be an important criteria for any well behaved BH spacetime. One of the main goal of this work is to put bounds on $(b,e)$ requiring stability of regularized stable Kerr BH when perturbed. For this purpose, we compute the scalar QNM modes and investigate whether the imaginary part of  these modes are negative needed for the stability. Our analysis shows that the regularized Kerr BHs are stable under scalar perturbation for all allowed values of $(b,e)$. Interestingly, for a fixed value of $e$, the imaginary parts $\omega_I$ of QNM frequencies decrease as we increase the values of $b$, which makes the BHs ``more" stable as it will take less relaxation time to shake off the perturbation. Also, for a fixed $a < M$, in the ``e-extremal" limit $(e \to 2)$, both imaginary and real parts of QNM frequencies tend towards zero. This gives us a way to identify ``e-extremality". Since small values $\omega_I$ lead to longer relaxation time, near-extremal values of $e$ make the BH ``less" stable.\\

\noindent
We have also studied the shadow structure of these regularized BH and its dependence on the additional parameters $(e,b)$. Since the parameter $b$ appears in the conformal factor, it doesn't affect the null geodesics and the shadow radius. In contrast, choice of $e$ values has a sharp impact on the angular shadow size, which we compare with $M87^*$ and $Sgr\, A^*$ images. Consistency with EHT shadow observations put stringent restrictions on the parameter space of $e$, but keep $b$ unconstrained. For the allowed values of the $e$, these regularized BHs mimic the Kerr shadow of corresponding mass and spin, but can still be distinguished from Kerr using QNM structure. To the best of our knowledge, this is the first observational bounds on the parameters of regularized BHs in consideration.\\

\noindent
In future, one may try to extend the consistency criteria discussed in \ref{Sec_3} on the mass profile $m(r)$ by taking back reaction into account. Also, we want to study other important properties of these BHs, such as their response to an external tidal field and check whether these BHs have zero Love number like Kerr. Further constraints in the parameter space can come from the study of gravitational lensing as well.\\

\section*{Acknowledgement}
The research of R.G. is supported by the Prime Minister Research Fellowship (PMRF Reference Number: 192002-120 and ID: 1700531), Government of India. Research of M.R. is supported by SERB, Government of India through the National Post Doctoral Fellowship grant (Reg. No. PDF/2021/001234). Research of AKM is supported by SERB, Government of India through the National Post Doctoral Fellowship grant (Reg. No. PDF/2021/003081). R.G. wants to thank Sudipta Sarkar for many useful discussions.


\section{Introduction \& Related Work}

Robots have a love-hate relationship with contact. On the one hand, substantial research effort has been invested over the years to develop algorithms that allow robots to \textit{avoid} colliding with nearby humans and obstacles~\cite{lavallePlanningAlgorithms2006,dawson22perception,singletary_cbf_arm}. On the other hand, robots must eventually make contact with their environment in order to be useful. Robots cannot help with household chores or manufacturing tasks without the ability to touch their environment.

That said, we cannot embrace the idea of robots intentionally making contact with their environment without answering important safety questions. For instance, a certain amount of force is required to complete most useful tasks, but too much force can cause damage to the humans or objects on the receiving end of the interaction. Over the years, roboticists have approached these questions in a variety of ways~\cite{manipulation}, but perhaps the most successful is the family of strategies known as ``stiffness''~\cite{singletary_cbf_arm} or ``impedance'' control~\cite{hogan_impedance_control,hogan_survey}. These controllers are conceptually simple but have seen great success; the core idea is to replace (via feedback control) the complex dynamics of a robotic arm with those of a much simpler mechanical system, such as a spring-mass-damper~\cite{manipulation}. As a result, when a human pushes on the robot it will respond in a predictable and physically intuitive way, and the force the robot exerts in response starts small and only increases gradually as it is pushed away from some target pose.

\begin{figure}[tb]
    \centering
    \includegraphics[width=\linewidth]{figs/task.jpg}
    \caption{A human-robot collaboration task that motivates our work, in which the human and robot must work together to move a large object. The challenge of this task lies in making the robot responsive to feedback from the human, which comes in the form of forces and torques exerted on the object.}
    \label{fig:task}
\end{figure}

Stiffness controllers have a number of appealing theoretical and practical properties~\cite{hogan_survey}, but they can be difficult to tune for practical use. For example, what is the appropriate stiffness to use when carrying an object? If the object is particularly heavy, then a large stiffness would be required, but this large stiffness could lead to large and potentially-unsafe motions if the object is dropped accidentally, where a large stiffness would cause the robot's arm to violently spring back to its equilibrium position. Tuning the stiffness and damping parameters of these controllers is non-trivial (and becomes more difficult if the user decides to incorporate nonlinear stiffness and damping laws), particularly if one set of parameters must be chosen to perform safely in a range of interactions~\cite{manipulation,Villani2016}. This tuning problem becomes even more complicated when humans enter the picture; some works attempt to estimate (at run-time) approximate mechanical models for the human, but it is difficult to construct an approximate model of human manipulation~\cite{cos_adaptive_cooperative_control}.

In this paper, we propose an alternate approach to programming robots for safe, compliant contact. Instead of specifying the desired behavior \textit{by analogy}, as stiffness control does by specifying the parameters of a spring-mass-damper system that the robot attempts to imitate, we allow the user to specify the desired behavior \textit{directly} by constraining the maximum force the robot is allowed to exert on its environment. To impose these constraints, we develop a force feedback controller based on control barrier functions (CBFs), which allow us to directly specify a set of safe and unsafe behaviors and derive a controller that is guaranteed to respect these constraints. This controller is simple to implement and tune, and we demonstrate its performance on an example task where a human and robot must collaboratively lift and maneuver a large object.

\subsection{Contributions}

Our main contribution in this paper is a novel architecture for force-feedback control based on control barrier functions. Our proposed controller is simple to implement and easy to tune, allowing the user to proceed directly from a specification of the robot's behavior (in terms of desired pose and maximum allowable force) to a compliant controller. Our controller requires no sensor feedback other than joint positions and force-torque measurements at the end effector, it is compatible with commonly-available robotic hardware, and it does not require any system identification for either the robot or its environment. We compare the performance of our controller with that of stiffness control to highlight the advantages and drawbacks of our approach, and we showcase the effectiveness of our controller in programming a robot to carry a bulky object collaboratively with a human. We hope that this novel, easy-to-use force feedback controller will be a useful building block for future robots to collaborate more safely and efficiently with humans.

\section{Problem Statement}\label{problem}

In this work, we develop a force-feedback controller that enables responsive, intuitive, and safe human-robot interaction. We motivate this work using the collaborative carrying task shown in Fig.~\ref{fig:task}, inspired by construction tasks, where a human and a robot must work together to manipulate a large object. The robot does not receive explicit information about the desired pose of the object, but it receives feedback from the human in the form of forces and torques exerted on the object (measured using a force-torque sensor in the robot's wrist).

This task highlights three requirements that drive the design of our force-feedback controller. First, the controller must be \textit{safe}: it must not cause the robot to exert too much force on the object, since excessive force could damage to the object, cause the robot's or human's grasp to fail, or push the human off balance. Second, the controller must be \textit{responsive}; since the robot has no explicit knowledge of where the human intends to move the object, it must be able to react when the human signals (by applying a force to the object) their intent to move in a certain direction. Finally, the controller must be \textit{stable}: in the absence of any force feedback from the human, the robot should hold the object level and avoid drifting or oscillation.

In traditional approaches to force-feedback control, e.g. stiffness or impedance control, these requirements are considered implicitly by replacing the dynamics of the robot's end effector with those of a spring-mass-damper surrogate and tuning the parameters of that surrogate to achieve empirically good performance. Our goal in this paper is to eliminate tedious tuning of stiffness parameters and instead derive a controller directly from the safety, responsiveness, and stability requirements. Our key insight is that satisfying the safety requirement (i.e. avoiding applying excess force on the carried object) will naturally lead to responsive behavior. When the human pushes on the object, signaling an intent to move in some direction, they cause the robot to approach its force limits, prompting it to respond appropriately; a \textit{safe} collaborator is thus necessarily a \textit{responsive} one. The rest of this section will provide a formal problem statement of the safety and stability properties that give rise to an easy-to-implement responsive force-feedback controller described in Section~\ref{approach}.

To simplify the presentation, we will derive our controller in task space (considering the Cartesian position, velocity, and wrench at the end effector), but a similar controller could be developed in joint space with suitable assumptions on the rank of the Jacobian. Let us denote by $\mathbf{w}_R = [\mathbf{f}, \mathbf{\tau}]$ the wrench measured at the end effector and expressed in the robot's base frame, comprising the concatenation of the force $\mathbf{f}$ and torque $\mathbf{\tau}$ vectors. Let $\mathbf{p}_R^G$ and $\mathbf{V}_R^G$ denote the pose and spatial velocity of the gripper in the robot's base frame, respectively, and let us denote by $\hat{\mathbf{p}}_R^G$ the desired end effector pose, which the robot will attempt to track in the absence of any force input from the human (e.g. to keep the carried object level). We assume the ability to command $\mathbf{V}_R^G$ directly, which in practice requires a) access to an underlying differential inverse kinematics~\cite{manipulation} or whole-body controller~\cite{wbc_paper} and b) that the robot is not operating near its joint limits or any kinematic singularity. Appropriate controllers are available for most commercially-available robots (e.g. the PAL Robotics TIAGo~\cite{Pags2016TIAGoTM}, which we use for our experiments, ships with an suitable whole-body controller), and the stability requirement encourages the robot to avoid joint limits and singularities.

Given this notation, we can formalize the \textit{safety} and \textit{stability} requirements as follows.
\begin{itemize}
    \item \textit{Safety}: the robot should respect user-specified force and torque limits; i.e., $\forall t$, $|f_i| \leq f_{i,\ max}$ and $|\tau_i| \leq \tau_{i,\ max}$ for $i \in \{x, y, z\}$. For many tasks, including the collaborative carrying task shown in Fig.~\ref{fig:task}, the limit for each axis can be estimated in a straightforward manner from the task specification (we provide an example of this in our experiments in Section~\ref{experiments}).
    \item \textit{Stability}: if it is possible to do so without violating the safety constraint, the robot should attempt to track the desired end effector pose; i.e., $\norm{\mathbf{p}_R^G - \hat{\mathbf{p}}_R^G} \to 0$ so long as the safety constraint is not active.
\end{itemize}

In the next section, we introduce a novel force-feedback controller based on control barrier functions (CBFs, for respecting the safety constraint) and control Lyapunov functions (CLFs, to encode the stability constraint); this controller is easy to implement and does not require us to estimate the impedance of the gripper or its coupling to the environment.

\section{Approach}\label{approach}

A common approach to constructing force-feedback controllers, stiffness and impedance control, involves replacing the dynamics of the robot with those of a simple damped spring and then carefully tuning the stiffness and damping parameters by hand until the safety and stability requirements are met~\cite{manipulation,Villani2016}. These methods approach the safety and stability requirements indirectly. In contrast, in this paper we take a \textit{requirements-first} approach to constructing a force-feedback controller that explicitly models the safety and stability requirements and acts accordingly.

To develop this controller, we turn to two powerful methods from control theory: control barrier functions (CBFs, for safety) and control Lyapunov functions (CLFs, for stability). CBFs and CLFs are certificate functions that allow a controller to not only find safe and stable control actions but also prove that those actions are safe and stable. We will begin by briefly introducing the theory of CBFs and CLFs and then demonstrate how they can be combined into a single force-feedback controller.

Consider a generic continuous-time dynamical system $\dot{x} = f(x, u)$, and suppose we are interested in imposing a safety constraint defined by the sub-level set of a function $h(x) \leq 0$ (i.e. all safe states result in a negative value of $h$ and all unsafe states result in a positive value of $h$) and a stability constraint about $x = 0$. We say that $h$ is a CBF if, for all $x$ and some $\alpha > 0$ there exists a $u$ such that
\begin{equation}
    \nabla_x h(x) \dot f(x, u) \leq -\alpha h(x) \label{eq:generic_cbf_condition}
\end{equation}

If $h$ is in fact a CBF (and a $u$ satisfying~\eqref{eq:generic_cbf_condition} always exists), then any $u$ that satisfies~\eqref{eq:generic_cbf_condition} is guaranteed to respect the safety constraint; i.e. if $h < 0$, then $h$ will never be allowed to become positive (safe states will never be driven into unsafe states), and if $h > 0$ then it must become negative exponentially quickly (unsafe states must be driven into safe states~\cite{ames_cbf,ferraguti_cbf_survey}). Just as a CBF lets us enforce a safety constraint on our controller, CLFs allow us to enforce stability constraints: if we can find some function $V(x)$ such that $V(\bar{x}) = 0$, $V(x \neq \bar{x}) > 0$, and for all $x \neq \bar{x}$ and some $\lambda > 0$ there exists a $u$ such that
\begin{equation}
    \nabla_x V(x) \dot f(x, u) \leq -\lambda V(x) \label{eq:generic_clf_condition}
\end{equation}
then we say that $V$ is a CLF and any $u$ that satisfies~\eqref{eq:generic_clf_condition} will cause the system to exponentially stabilize around the desired state $x \to \bar{x}$~\cite{ames_cbf,ferraguti_cbf_survey}. We can combine these two certificates into a single controller that respects both safety and stability requirements by solving an optimization problem for $u$:
\begin{subequations}
\begin{align}
    \min_{u, \gamma} \quad& \norm{u}^2 + k\gamma \label{eq:cbf_clf_qp}\\
  \text{s.t.} \quad& \nabla_x h(x) f(x, u) \leq -\alpha h(x) \label{eq:cbf_constraint}\\
              & \nabla_x V(x) f(x, u) \leq -\lambda V(x) + \gamma \label{eq:clf_constraint}\\
              & \gamma \geq 0 \label{eq:gamma_constraint}
\end{align}
\end{subequations}
As discussed in Section~\ref{problem}, the robot may temporarily ignore the stability requirement if necessary to preserve safety; the additional decision variable $\gamma$ allows the controller to relax the CLF condition~\eqref{eq:clf_constraint} governing stability as needed to satisfy the safety-critical CBF condition~\eqref{eq:cbf_constraint}. When the system has control-affine dynamics (i.e. $f(x, u)$ is affine in $u$), then this optimization problem takes the form of a quadratic program (QP) and can be solved efficiently online~\cite{dawson2021safe,ames_cbf}. QP-based CBF/CLF controllers like this have emerged in recent years as powerful tools for deriving feedback controllers directly from safety and stability constraints~\cite{ames_cbf,ferraguti_cbf_survey,singletary_cbf_arm,dawson2021safe,dawson22perception,dawson2022survey}.

Next, we must adapt this CBF/CLF framework to design a controller that responds safely to forces applied to the robot's end-effector. Recall that we are concerned with controlling the velocity $\mathbf{V}_R^G$ of the robot's end-effector. This leaves us with a simple continuous-time control-affine dynamical system involving the end-effector pose $\mathbf{p}_R^G$ (for notational convenience, both linear and angular velocities are concatenated into a single vector $\mathbf{V}_R^G$): $\der{}{t} \mathbf{p}_R^G = \mathbf{V}_R^G$.

Deriving a CLF for this system is straightforward; it is easy to show that $V(x) = \frac{1}{2}\norm{\mathbf{p}_R^G - \hat{\mathbf{p}}_R^G}^2$ is a valid control Lyapunov function for these dynamics. Deriving a CBF is slightly more involved and requires some assumptions about the nature of the robot's contact interactions with its environment. Let us model the contact wrench between the robot and its environment using a simple linear stiffness model $\mathbf{w}_R = -\mathbf{K} (\mathbf{p}_R^G - \mathbf{p}_R^{G_0})$, where $\mathbf{p}_R^{G_0}$ is the pose at which contact is made and $\mathbf{K}$ is an \textit{unknown} stiffness matrix which we assume to have no coupling between the various degrees of freedom (i.e. $\mathbf{K}$ is a matrix with unknown positive diagonal entries and zero off-diagonal entries, as in~\cite{Villani2016}). In our experiments, we will show that despite the simplicity of this contact model, our system is nevertheless able to remain safe; moreover, we are able to achieve this safe performance \textit{without explicitly estimating the parameters of $\mathbf{K}$}, as we show in the following theorem.

\begin{theorem}\label{thm:cbf}
Recall that our safety requirement involved bounding the magnitude of the wrench along each degree of freedom (i.e. $|w_i| \leq w_{i, max}$ for $i \in \set{fx, fy, fz, \tau x, \tau y, \tau z}$). We can exploit the fact that our control over each degree of freedom (DoF) is decoupled and define a corresponding set of six CBFs, one for each DoF:
\begin{align}
    h_i(x) = \frac{1}{2}\pn{w_{i}^2 - w_{i, max}^2};\quad i \in \set{fx, fy, fz, \tau x, \tau y, \tau z} \label{eq:contact_cbf}
\end{align}

Each individual $h_i(x)$ is a valid CBF; moreover, we can find a single control input $\mathbf{V}_R^G$ that satisfies all of these CBFs \textit{without} any knowledge of the stiffness matrix $\mathbf{K}$.
\end{theorem}
\begin{proof}
First, let us prove that each $h_i$ is a valid CBF by showing that Eq.~\eqref{eq:generic_cbf_condition} is satisfiable for each DoF. Let $k_i$ denote the $i$-th diagonal entry of $\mathbf{K}$, and let $v_i$ denote the $i$-th component of $\mathbf{V}_R^G$.
\begin{align}
    \nabla_x h_i(x) \dot f(x, u) &\leq -\alpha h_i(x) \nonumber\\
    -w_i k_i v_i &\leq -\alpha h_i(x) \nonumber \\
    -w_i v_i &\leq -\frac{\alpha}{k_i} h_i(x) \label{eq:contact_cbf_condition}
\end{align}

Thus, any $v_i$ that satisfies $-w_i v_i \leq -\alpha' h(x)$ for some $\alpha' \geq 0$ will also satisfy the CBF condition \eqref{eq:contact_cbf_condition} with $\alpha = k_i\alpha'$. Since $k_i$ is unknown but assumed to be positive, $k_i \alpha' > 0$, so the safety properties of the CBF will still hold. As a result, it is possible to choose the control input for each DoF $v_i$ without any knowledge of the stiffness parameters. Moreover, since the control inputs for each DoF are uncoupled, it is sufficient to satisfy~\eqref{eq:contact_cbf_condition} individually for each DoF and then concatenate the resulting $v_i$ into a single control input $\mathbf{V}_R^G$.
\end{proof}

By exploiting the decoupling between the various degrees of freedom at the end effector, we can thus create a simple QP-based controller that explicitly respects both the safety and stability requirements of our problem:
\begin{subequations}
\begin{align}
        \min_{\gamma\geq 0, \mathbf{V}_R^G} \quad& \norm{\mathbf{V}_R^G}^2 + k\gamma \label{eq:contact_cbf_clf_qp}\\
  \text{s.t.} \quad& -w_i k_i v_i \leq -\frac{\alpha}{2}\pn{w_{i}^2 - w_{i, max}^2}; \quad \text{(CBF \textit{i})} \label{eq:contact_cbf_constraint}\\
             & \quad\qquad i \in \set{fx, fy, fz, \tau x, \tau y, \tau z} \nonumber \\
              & \pn{\mathbf{p}_R^G - \hat{\mathbf{p}}_R^G}\mathbf{V}_R^G \leq -\frac{\lambda}{2} \norm{\mathbf{p}_R^G - \hat{\mathbf{p}}_R^G}^2 + \gamma \quad \text{(CLF)} \label{eq:contact_clf_constraint}
\end{align}
\end{subequations}

This controller is parameterized only by the user-specified limits on the wrench in each axis, the desired end effector pose, and the CBF and CLF parameters $\alpha$ and $\lambda$. We find that tuning $\alpha$ and $\lambda$ is typically not necessary ($\alpha = \lambda = 1$ is a sensible default); these parameters do not affect the safety of the system and govern only how fast the system may allow the wrench to increase and how fast it should converge towards the desired pose (in contrast to the parameters of a traditional stiffness controller, which can be difficult to tune and have a direct impact on the safety and stability of the resulting controller).

\section{Experiments}\label{experiments}

To demonstrate the effectiveness of our force-limiting CBF controller in hardware, we first validate the performance and safety of the controller in simple experiments before highlighting the use of our controller to solve the collaborative carrying task from Fig.~\ref{fig:task}. In all hardware experiments, we implement our controller using Python and Casadi~\cite{Andersson2019} to solve the QP in~\eqref{eq:contact_cbf_clf_qp}. We use a 2-armed TIAGo mobile robot developed by PAL Robotics~\cite{Pags2016TIAGoTM} for our experiments, which runs ROS middleware and provides 6-DoF wrench measurements at the wrist of each arm and control of the end-effector pose via an onboard whole body kinematic controller. Before running our controller, we measure the gravitational wrench on the end-effector and subtract this from the measured wrench to estimate $\mathbf{w}_R$. We deploy our controller on the robot's onboard computer at \SI{30}{Hz}.

\subsection{Validation of the CBF force controller}

To validate the performance of our CBF-based force controller, we a simple experiment where the robot is commanded to maintain a pose $\hat{\mathbf{p}}_R^G = [0.5, -0.3, 1, 0, 0, 0]$ (for readability, we write poses as $[x, y, z, roll, pitch, yaw]$, but our implementation relies on quaternions to represent rotations). We then hang a weighted bag from the robot's end effector to induce a wrench and observe the robot's response.

We attach the bag so that it induces an approximately \SI{20}{N} force in the $-z$ direction. To isolate the controller's response in linear degrees of freedom, we set the maximum allowable torques $\tau_{x, max} = \tau_{y, max} = \tau_{z, max} = \SI{10}{Nm}$, larger than the torque exerted by the test mass, and set the maximum allowable force in all linear directions to be \SI{25}{N}. We also set parameters $\lambda = 10$, $\alpha = 1$, and $k = 1$. We then start the controller, add the \SI{20}{N} load, wait several seconds, and then add an additional \SI{5}{N} load to the bag. The test setup is shown in Fig.~\ref{fig:linear_test_setup}.

\begin{figure}[tb]
    \centering
    \includegraphics[width=\linewidth]{figs/test_linear_setup.jpg}
    \caption{The test setup used to assess the CBF controller's response to a load exceeding its linear load limit. First, a \SI{20}{N} test load is attached to the arm, inducing a force primarily along the robot's $-z$ axis; then, a \SI{5}{N} mass is added to the test load. The additional load approaches the force limit of the CBF controller, prompting it to respond.}
    \label{fig:linear_test_setup}
\end{figure}

The results of this test are shown in Figs.~\ref{fig:linear_test_results} and~\ref{fig:linear_test_comparison} (and in our supplementary video): as intended, the robot does not move when the \SI{20}{N} load is applied, but when the \SI{25}{N} limit is exceeded the robot moves to lower the test mass. Since the test mass is hanging from the robot, the robot is not able to immediately reduce the force (as we assumed with the linear stiffness model used in Theorem~\ref{thm:cbf}), but it still reacts appropriately. We compare our CBF controller with a traditional stiffness controller, as shown in Fig.~\ref{fig:linear_test_comparison} and in the supplementary video.

Our controller is able to support the \SI{20}{N} load but moves in response to a load that exceeds the safety threshold (i.e. it respects the safety constraints but is able to track the desired position with a reduced load). In contrast, even though we tuned the parameters of the stiffness controller to be as stiff as possible (without causing instability due to measurement noise and the bandwidth of the controller), it was unable to support even the base load of \SI{20}{N}. This comparison highlights a significant advantage of our approach over stiffness control: it is in general difficult to balance safety (i.e. not applying excessive force) with performance (in this case, goal tracking under a lighter load) without complicated system identification modules, and even then the process requires significant manual tuning. In contrast, we are able to proceed nearly directly from a specification (e.g. ``track this position but avoid applying any force in excess of \SI{25}{N}'') to a controller that satisfies that specification, simplifying the process of developing a safe controller for human-robot interaction applications.

\begin{figure}[tb]
    \centering
    \includegraphics[width=\linewidth]{figs/test_linear_results_2.jpg}
    \caption{The measured net force along the robot's $z$ axis, showing our controller's response to applied \SI{20}{N} (applied at \SI{10}{s}) and \SI{25}{N} (at \SI{17}{s}) loads. Since we set the safety limit in the $z$ direction at \SI{25}{N}, the robot does not move until the \SI{25}{N} load is applied; once the additional load is applied it smoothly drops the load. The spikes occurring at $>$\SI{26}{s} are caused by the test mass touching the ground, de-loading the robot and causing it to try to move back towards its setpoint.}
    \label{fig:linear_test_results}
\end{figure}

\begin{figure*}[tb]
    \centering
    \includegraphics[width=\linewidth]{figs/linear_cbf_vs_impedance.jpg}
    \caption{A comparison of our CBF controller and a traditional stiffness controller in the test described in Fig.~\ref{fig:linear_test_setup}. Our controller maintains its position when the \SI{20}{N} force is applied and only moves when the \SI{25}{N} force is applied to exceed its prescribed safety limit. In contrast, the stiffness controller fails to support the \SI{20}{N} load, despite the stiffness being tuned as high as possible without becoming unstable (\SI{40}{N/m}).}
    \label{fig:linear_test_comparison}
\end{figure*}

\subsection{Collaborative carrying task}

We can now apply our CBF force-feedback controller to the motivating human-robot collaboration example from Fig.~\ref{fig:task}. In this task, the human and the robot must collaboratively lift a bulky object; however, to avoid pushing the human off balance or causing their grasp to fail, we wish to impose safety limits on the forces and torques that the robot can apply to the carried object. We set the safety limits as follows: the force in the $x$ and $y$ directions (in the robot's frame) must not exceed \SI{10}{N}, the force in the $z$ direction (less 50\% of the carried object's weight) must not exceed \SI{10}{N}, the torques about the $y$ and $z$ axes must not exceed \SI{3}{Nm}, and the torque around the $x$ axis must not exceed \SI{0.5}{Nm}. We use CBF parameters $\lambda = 10$, $k = 1$, $\alpha = 1$ for forces, and $\alpha = 10$ for torques. To add an additional layer of responsiveness, we add a simple inverse-dynamics controller to control the movement of the robot's base to that of a free-floating rigid body. This allows the human to not only guide the robot's end-effector but also move the robot around the room. This inverse-dynamics controller only affects the motion of the base; the arm and torso motions are driven by the CBF controller, and we include a video of our controller operating without any base motion in our supplementary video.

Using this capability, we simulate a basic construction task where the human and the robot must collaboratively position a panel of foam insulation. We use an ArUco marker~\cite{aruco} to localize the panel and establish a grasp, but we cover this marker once the robot has grasped the panel (it is not used while carrying; only force and joint position measurements are used). We also program the robot to recognize when the human has rotated the panel by 90 degrees and lowered it to the ground; this is the signal to release the board and end the test.

\begin{figure*}[tb]
    \centering
    \includegraphics[width=\linewidth]{figs/carrying_snapshots.jpg}
    \caption{Using our CBF-based force-feedback controller to control the right arm and torso of the TIAGo robot (the base is controlled separately with a velocity controller that makes the base behave like a free-floating, damped rigid body). Our supplementary video shows the robot responding to push/pull, side-to-side, and twisting motions from the human; our CBF controller allows the robot to respond without applying too much force to the carried object, which could cause either the human's or the robot's grasp to fail.}
    \label{fig:carrying_snapshots}
\end{figure*}

As shown in Fig.~\ref{fig:carrying_snapshots} and in our supplementary video, our controller enables the robot to successfully perform this collaborative carrying task, allowing the human to guide the robot through back-and-forth, side-to-side, and twisting motions of the carried panel. We found that for the majority of the task, the robot's end effector was pushed away from the desired target pose, indicating that the CBF constraints were active when the human was guiding the motion of the carried object.

\section{Limitations \& Conclusion}

In our experiments, we found that the human's qualitative experience of working with our CBF-based controller was sensitive to the bandwidth of the controller and underlying hardware. Our experiments on the TIAGo robot were limited by the bandwidth of the force-torque sensor (which publishes measurements at \SI{30}{Hz}) and the whole-body controller that we used to implement velocity control of the robot's end effector. For safety reasons, the speed of the TIAGo's end effector is limited when using the provided whole-body controller, which means that human testers reported that the robot sometimes felt ``slow'' or ``sluggish'' to work with. We anticipate that the user experience could be improved by using more performant hardware for both sensing and actuation; we do not believe that the CBF-based controller will cause issues in this domain, as previous work has shown that control frequencies of \SI{>100}{Hz} are possible with CBF-based controllers (even with un-optimized Python implementations~\cite{dawson2021safe}).

Despite these limitations, we present an easy-to-implement and simple-to-tune force feedback controller that enables safe, responsive human-robot interaction. Our controller uses control barrier functions to ensure that the robot does not exert too much force on its environment while still allowing useful contact to occur. Moreover, the use of our controller does not require any system identification of either the robot or its environment, and our controller does not require any sensor feedback beyond joint positions and wrench measurements at the end effector (both of which are commonly available for commercial robotic arms). We use our controller to accomplish a simulated human-robot collaborative construction task, cooperatively maneuvering a large bulky object without subjecting either the human or the carried object to excessive force. We hope that our controller will help enable solutions to a wide range of similar problems where safe, responsive force-feedback can help robots collaborate more effectively with humans.

\section*{Acknowledgments}

We would like to sincerely thank the MIT-IBM Watson AI Lab team for providing access to and support for the TIAGo robot, particularly Jason Phillips, Luke Inglis, Mark Weber, and David Cox.

\bibliographystyle{IEEEtran}