\section{Introduction}
\label{sec:intro}

The last ten years or so have seen a revolution in the ways we probe strong-field relativistic gravity. The main breakthrough came in 2015  
with the first observation of gravitational waves (GWs) from merging black holes by the LIGO-Virgo Collaboration~\cite{Abbott:2016blz} . 
Since then these detectors have chalked up many more merging compact binary systems, allowing for new precision tests of General Relativity (GR)
and new astrophysical information on compact objects (for a review, see~\cite{isi2021}). 
The second most important milestone, and the one most relevant to this paper, was reached by the spectacular high-resolution 
image of the supermassive black hole in the M87 galactic center (usually dubbed M87*), obtained by the Event Horizon Telescope (EHT) 
Collaboration~\cite{EHT2019ApJ.V}. During the preparation of this paper a second sister image was released, picturing the  supermassive black hole SgrA*
in our galactic center~\cite{EHT22}. These images serve as direct evidence for the existence of black holes and can also be used
as probes of GR  (for a review, see, e.g., Ref.~\cite{Broderick2014ApJ}). Motivated by this exciting possibility,  a significant amount of work over the last decade or so, 
has focused on the calculation of shadows of black holes beyond GR 
(e.g., Refs.~\cite{Johannsen2010ApJ, Johannsen2013ApJ, Cunha2015PhRvL, Cunha2016IJMPD, Cunha2016PhRvD, 
Cunha2017PhysRevD, Cunha2017PhysRevLett, Cunha2018GReGr, Medeiros:2019cde, wielgus21, junior21, lara21})  
as well as on improving our understanding of the image produced by general relativistic Kerr black holes (e.g., Refs.~\cite{Gralla_etal2019, Gralla2020, 
Gralla_Lupsasca2020, bronzwaer21, paugnat22}).


The centerpiece in these images is the shadow cast by the black hole as silhouetted against its luminous accretion flow. The shape and overall scale of the shadow, 
as projected onto the `optical plane' of a distant observer, is formed by photons freely moving along geodesics of the black hole's spacetime.  
In principle then, a shadow image like that of M87* is a \emph{geodetic} `experiment' that could enable tests of GR via the so-called `Kerr hypothesis', that is, the 
theoretically predicted uniqueness of the Kerr metric as the correct description of astrophysical  black holes. 
As far as Kerr black holes are concerned, and assuming a source of illumination of angular size much greater than the hole itself,  it is known that light rays that are 
asymptotically captured at the location of the unstable photon orbit give rise to a shadow that is nearly circular-shaped provided the black hole spin is not close to the 
maximum allowed limit~\cite{bardeen73, takahashi04}. This `Bardeen shadow', which also coincides with the black hole's capture cross section for light rays incoming 
from infinity and moving parallel to the equatorial plane, can de defined in a mathematically invariant way by expressing its equatorial and polar radii in terms of the 
corresponding impact parameters which themselves are combinations of the geodesic constants of motion~\cite{bardeen73}. 


These considerations were exploited in a recent EHT paper~\cite{Psaltis_etal2020} where the shape of the M87* shadow was used to set limits on the
non-GR parameters of Johannsen's deformed Kerr metric~\cite{Johannsen:2013pca} 
(see, however, Refs.~\cite{kggp21,voelkel21} for a discussion on the limitations of black hole shadows as probes of GR gravity, when the deviations 
arise purely from non-GR gravitational degrees of freedom and are completely unrelated to accretion physics~\cite{ozel22, younsi23}).

A closely related notion has to do with the connection between the shape of a black hole shadow and the hole's mass and spin \emph{multipole moments} 
$\{ M_\ell, S_\ell \}$ (where the index $\ell$ is a non-negative integer).  According to GR's no-hair theorem, Kerr black holes are characterised by finely tuned set of multipole 
moments that are fully determined algebraically by the first two, the mass $M=M_0$ and spin angular momentum $J=S_1$. The theorem is encapsulated in the
formula~\cite{hansen74},
\begin{equation}
M_\ell + i S_\ell = M ( i J/M )^\ell.
\label{nohair}
\end{equation}
With the help of a Kerr-like metric with an adjustable quadrupole moment $M_2$, previous work has shown that a deviation from the Kerr quadrupole moment $M_2 = -J^2/M$ 
manifests itself as an oblate or prolate deformation of the Kerr shadow~\cite{Johannsen2010ApJ, Broderick2014ApJ}. These results have led to the commonly stated claim that 
the no-hair theorem itself is testable by the shape of a black hole shadow.


A nearly circular shadow is not an exclusive characteristic of Kerr black holes; in fact there are known non-Kerr black hole spacetimes that enjoy the same property. 
Examples include the aforementioned Johannsen metric~\cite{Johannsen:2013pca, Johannsen2013ApJ} and the more general metric of Carson \& Yagi~\cite{Carson_Yagi2020}. 
These spacetimes, however, are special in the sense that they are separable, i.e. they admit a third integral of motion for geodesic motion (in Kerr this is the 
well-known Carter constant~\cite{MTW1973}). This observation could be taken as evidence for  a plausible intrinsic relation between the near circularity of the shadow 
and the spacetime's separability.
 

In this paper, we make contact with the above considerations and assess to what extent a shadow image like that of M87*
could provide a reliable test of the no-hair theorem relation~\eqref{nohair}. To this end, we revisit the dependence of the shadow shape 
on the multipole moments by making use of stationary-axisymmetric spacetime metrics with an \emph{arbitrary} structure in the quadrupole and higher 
multipole moments. 

Our analysis is performed within the framework of GR and the employed metrics are vacuum solutions of the theory, taking the form of an expansion 
in the spin $J$ or in the inverse radial distance $1/r$ and the multipole order. It should be emphasised that, as a consequence of GR's 
uniqueness theorems for black holes~\cite{MTW1973}, these metrics do not represent true black holes in the sense that they are infested with 
horizon-piercing curvature singularities. Nevertheless, as long as the spacetime's circular photon orbit does not approach the event horizon,
this pathology has little bearing on the properties of the shadow. Alternatively, the central body could have a material surface (with negligible emission) 
instead of a horizon, i.e. something akin to a `gravastar'~\cite{mazur04}. In both cases the aforementioned spacetimes can describe the exterior 
spacetime of the putative non-Kerr compact object or even serve as proxies for the spacetime of non-GR black holes. These are precisely the 
systems likely to violate the no-hair theorem of canonical Kerr black holes. As a secondary topic of our paper we provide an analysis of the multipolar 
structure of the Johannsen metric in order to explore a connection (if any) to the spacetime's Kerr-like shadow. 

A key assumption underpinning our analysis is that the deformation away from Kerr, and the ensuing `decircularisation' of the shadow, is caused by the rotation 
of the  `black hole'. As a consequence, the effect is \emph{maximised} at the equatorial plane where the shadow radius of a backlit system is identical to the impact 
parameter $b_{\rm ph}$ associated with the radius $r_{\rm ph}$ of the unstable equatorial photon orbit (the so-called light ring). This implies that as long as we are limited 
to modest deviations from the Kerr spacetime, the photon orbit impact parameter provides an accurate measure of the shadow's non-circularity (as, in fact, it does for 
Kerr black holes). 

 
The remainder of the paper is organised as follows. In Section~\ref{sec:axistationary} we describe the formalism for the calculation of the 
radius and associated impact parameter of the equatorial light ring of a general stationary-axisymmetric spacetime. The following two sections
comprise the paper's main calculation and results. In Section~\ref{sec:HTshadow} we obtain relations for the shadow equatorial radius 
as a function of the first few multipole moments in the Hartle-Thorne vacuum spacetime. A similar calculation is repeated in Section~ \ref{sec:axisym} 
in the context of another general stationary-axisymmetric spacetime. Section~\ref{sec:Jmetric} is dedicated to the multipole
moment analysis of the Johanssen metric. Our concluding remarks can be found in Section~\ref{sec:conclusions}. 
Throughout the paper we adopt relativistic units $G=c=1$ and use a prime to denote a radial derivative. 



\section{Photon ring and impact parameter: general formalism}
\label{sec:axistationary}

The formation of a shadow is the manifestation of extreme light bending in the spacetime of a massive `dark' body,
caused by the presence of an unstable photon orbit, the so-called `light ring'. Therefore, the first step of our analysis 
is the calculation of the light ring radius $r_{\rm ph}$ and the associated impact parameter $b_{\rm ph}$ which describes photons that approach 
the black hole from infinite distance and get trapped at the light ring.

Here we consider an arbitrary axisymmetric and stationary metric of the form
\begin{equation}
ds^2 = g_{tt} dt^2 + g_{rr} dr^2 + 2g_{t\varphi} dt d\varphi + g_{\theta\theta} d\theta^2 + g_{\varphi\varphi} d\varphi^2,
\label{metric}
\end{equation}
with $ g_{\alpha\beta} = g_{\alpha\beta} (r,\theta)$, assuming a spherical-like coordinate system. 
 %
The assumed symmetries allow us to write the following equations for the $u^t, u^\varphi$  four-velocity components,
\begin{equation}
u^t = \frac{1}{{\cal D}} \left (\,  g_{t\varphi} b + g_{\varphi\varphi}  \, \right ), \quad
u^\varphi = -\frac{1}{{\cal D}} \left (\,  g_{t\varphi}  + g_{tt} b \, \right ),
\end{equation}
where $  {\cal D} = g_{t\varphi}^2  - g_{tt} g_{\varphi\varphi} $. The orbital constants $E$ (energy per unit mass) and 
$L$ (angular momentum per unit mass) enter through the impact parameter $b = L/E$. Upon inserting these in the 
normalisation condition $u^\mu u_\mu =0$ and setting $u^\theta=0$ we end up with a radial motion equation,
\begin{equation}
 g_{rr} (u^r )^2  = \frac{1}{{\cal D}} \left (\, g_{tt} b^2 + 2 g_{t\varphi} b + g_{\varphi\varphi} \, \right ) \equiv   V_{\rm eff} (r ,b).
\end{equation}
A circular orbit obeys 
\begin{equation}
 V_{\rm eff} ( r_{\rm ph} ,b_{\rm ph}) = 0, \qquad   V_{\rm eff}^\prime (r_{\rm ph} ,b_{\rm ph}) = 0.
\end{equation}
These two conditions lead to the light ring equation,
\begin{align}
& \Big \{ 4 \left (  g_{t\varphi}  g_{tt}^\prime   -  g_{tt}  g_{t\varphi}^\prime  \right )   
\left (  g_{t\varphi} g_{\varphi\varphi}^\prime    - g_{\varphi\varphi} g_{t\varphi}^\prime +  \right ) 
\nonumber \\
&  + \left ( g_{\varphi\varphi} g_{tt}^\prime - g_{tt} g_{\varphi\varphi} ^\prime \right )^2 \Big \}_{r=r_{\rm ph}} =0,
\label{photonEq}
\end{align}
and the associated impact parameter,
\begin{equation}
b_{\rm ph}= \frac{1}{2} \frac{ g_{\varphi\varphi} g_{tt}^\prime - g_{tt} g_{\varphi\varphi} ^\prime}{g_{tt}  g_{t\varphi}^\prime -g_{t\varphi}  g_{tt}^\prime }.
\label{bgeneral}
\end{equation} 
These formulae describe both prograde and retrograde photon motion, the two cases being distinguished by the sign of the spin parameter
(which does not appear explicitly here). As a benchmark example we may consider the standard Kerr metric in Boyer-Lindquist parameters. We find,
\begin{align}
& r_{\rm ph} ( r_{\rm ph}  -3M)^2 = 4 a^2 M^3, 
\\
& b_{\rm ph} = M \left [ a + \left ( \frac{r_{\rm ph} }{M} \right )^{3/2} \right ],
\label{bphKexact}
\end{align}
where $a = J/M^2$ is the dimensionless spin parameter. Prograde (retrograde) motion would correspond to $a>0$  $(a<0)$.

From its definition (and its association with the light ring), $b_{\rm ph}$ coincides with the equatorial radius of the shadow cast
by a black hole that is backlit from a source of large angular size (e.g., a distant luminous plane) and viewed from the opposite direction 
by an equatorial observer. As discussed, for example, in Ref.~\cite{bardeen73}, this is a  coordinate \emph{invariant} identification. 
Rotation causes the bifurcation of the light ring into prograde (co-rotating) and retrograde (counter-rotating) branches and as a consequence
the shadow develops a left-right asymmetry.
The two orbits vary asymmetrically with the black hole spin, resulting in a relative displacement between the centre of the black hole (at $r=0$) and 
that of the non-circular shadow.  Nevertheless, knowledge of the prograde and retrograde impact parameters $ b_{\rm ph} = \{ b_{\rm pro}, b_{\rm retro} \}$ 
allows us to calculate an invariant equatorial shadow radius (as viewed by an equatorial observer) defined by the geometric average,
\begin{equation}
\bar{b}_{\rm ph} = \frac{1}{2} \left ( b_{\rm pro} + b_{\rm retro} \right ).
\label{bphaver}
\end{equation}
The equatorial plane (or, more generally, the latitude slice  $\theta \sim \pi/2$) is where the shadow's shape is expected to show
the most pronounced deviation from circularity because of the maximum differential dependence of $r_{\rm ph}$ and $ b_{\rm pro},  b_{\rm retro} $
on the spin.  The opposite arrangement is expected to happen near the symmetry axis $\theta = \{ 0, \pi \}$ where the impact
of rotation is at its minimum.  
This situation is exemplified by the shadow shape found in a number of deformed Kerr spacetimes 
(see e.g.~\cite{Johannsen2013ApJ, Carson_Yagi2020}) and of course by the Kerr shadow itself, where as the spin parameter $a \to 1$ it 
acquires a characteristic `D' shape as a result of the rapid  shrinkage of $b_{\rm pro}$ with respect to the much slower changing $b_{\rm retro}$, 
see, e.g.,  Ref.~\cite{Johannsen2010ApJ}. The example most relevant to this paper is the shadow obtained in Ref.~\cite{kostaros22}
\footnote{In that work the Hartle-Thorne metric is used `as it is', without further expansion of the geodesic equations with respect
to the spin. An interesting consequence of this approach is the emergence of non-equatorial light rings, a property not present in our model
where all equations are consistently spin expanded.}  
using the Hartle-Thorne spacetime (the subject of the next section); it is found that, in spite of the presence of moderate rotation, 
the shadow's polar radius remains very close to the radius of a Schwarzschild black hole.   





\section{Shadow radius in the Hartle-Thorne spacetime} 
\label{sec:HTshadow}

\subsection{The ${\cal O} (J^4)$ Hartle-Thorne spacetime}

The purpose of this and the following section is to put on a quantitative basis the relation between the black hole shadow
size and a given spacetime's multipole moment structure. Perhaps the best available tool for undertaking this task in
GR is the celebrated \emph{exterior} Hartle-Thorne (HT) metric~\cite{HT68} in its modern incarnation of ${\cal O} (J^4)$ 
precision~\cite{Yagi_etal:2014}, where $J$ represents the central body's angular momentum.

The HT metric has the advantage of allowing arbitrary values for the first five mass and spin multipole moments, $\{M, S_1, M_2, S_3, M_4 \}$;
these are Geroch-Hansen moments, see Ref.~\cite{quevedo90} for a review. In addition, it is generically non-separable with respect to the Hamilton-Jacobi 
equation for geodesic motion, that is, it does not admit a Carter constant.  As already mentioned in the Introduction, we employ the HT metric to describe the
exterior spacetime of a rotating ultracompact body with a surface or as an effective black hole spacetime, ignoring the fact that it is not well behaved in the 
vicinity of the event horizon. This is not a problem as long as the light ring stays well clear of the horizon radius. It should be pointed out that the HT metric 
includes the Kerr spacetime as a special limit (and as an expansion in the spin) with the added property of being singularity-free (except at the center) and separable. 

The HT metric is an expansion in the spin and takes the following functional form~\cite{Yagi_etal:2014}
(here $\epsilon$ is the spin-order bookkeeping parameter, $\mu =\cos\theta$ and $\nu = \{t, r,\mu,\varphi \}$):
\begin{align}
g_{\nu\nu} (r,\mu) &= g_{\nu\nu}^{(0)} (r) + \left[ g_{\nu\nu}^{(20)} (r) + g_{\nu\nu}^{(22)} (r) P_2 (\mu) \right] \epsilon^2 
 \nonumber \\
&+ \left[ g_{\nu\nu}^{(40)} (r)  + g_{\nu\nu}^{(42)} (r) P_2 (\mu) + g_{\nu\nu}^{(44)} (r) P_4 (\mu) \right]\epsilon ^4 
\nonumber \\
& + {\cal O} (\epsilon^6 ),
\\
\nonumber \\
g_{t \varphi} (r,\mu) & = (1-\mu^2 ) \left \{ g_{t\varphi}^{(1)} (r)   \epsilon  \right.
\nonumber \\
& \left. +\left[ g_{t \varphi}^{(31)} (r) + g_{t\varphi}^{(33)} (r) \frac{d P_3}{d\mu} (\mu) \right]\epsilon^3 \right \} 
+ {\cal O} (\epsilon^5 ),
\end{align}
where $P_n (\mu)$ is the standard Legendre polynomial. The coordinate change $\mu \to \theta$ in the line element results in 
$ g_{\theta \theta} (r,\theta) =  g_{\mu \mu } \sin ^2 \theta$. The radial metric functions $g_{\mu\nu}^{(nm)}(r)$ are given by rather lengthy expressions 
and therefore are not presented here. The interested reader can find them in Ref.~\cite{Yagi_etal:2014}.


Apart from their dependence on the coordinates $\{r,\mu\}$, the metric components depend on a number of parameters 
associated with the central body. These are
\begin{equation}
\{ M_{\rm s}, J, C_{20}, C_{40}, C_{22}, C_{42}, C_{44}, C_{31}, C_{33} \},
\end{equation}
where $M_{\rm s}$ is the mass of the spherical body in the limit of zero rotation. A given $C_{nm}$ parameter is of order $\epsilon^n$
in the spin and appears in tandem with the $P_m (\mu)$ polynomial. As we are about to see, these are related to the source's first five
multipole moments, $\{M, S_1, M_2, S_3, M_4 \}$. The simplest way for extracting the HT moments is to expand $g_{tt}, g_{t\varphi}$ in powers of $1/r$,
\begin{widetext}
\begin{align}
g_{tt} & = -1 + \frac{2 M_{\rm s}}{r} ( 1 + C_{20} \epsilon^2 + C_{40} \epsilon^4) - \frac{2}{r^3} \left \{ C_{42} M_{\rm s}^3 \epsilon^4
+ \left ( \frac{J^2}{M_{\rm s}} + \frac{8}{5} C_{22} M^3_{\rm s} \right ) \epsilon^2  \right \} P_2 
\nonumber \\
&  + \frac{2 P_4}{r^5} \left [ \frac{107}{105} \frac{J^4}{M_{\rm s}^3} + \frac{9428}{735} C_{22} J^2 M_{\rm s} - \frac{4}{7} C_{33} J M^3_{\rm s}
+ \left (  C_{44} + \frac{3144}{245} C_{22}^2 \right ) M_{\rm s}^5  \right ] \epsilon^4  + {\cal O} \left (\frac{\epsilon^2}{r^4}, \frac{1}{r^6}, \epsilon^6 \right ),
\label{gexpand1}
\\
\nonumber \\
g_{t\varphi} & = \sin^2\theta \left \{  -\frac{2}{r} \left ( J \epsilon +  C_{31} M_{\rm s}^2 \epsilon^3 \right ) 
- \frac{2}{3} \frac{C_{33} M_{\rm s}^4}{r^3} \frac{d P_3}{d\mu}   \epsilon^3 \right \}  + {\cal O} \left ( \frac{1}{r^5}, \epsilon^5 \right ).
\label{gexpand2}
\end{align}
\end{widetext}
Following the same practice as in Newtonian gravity, we can read off the multipoles from the coefficients of  the $1/r$ powers (see discussion in Ref.~\cite{Pappas:2012}).
 We find,
\begin{align}
M &= M_0 =  M_{\rm s} \left ( 1 + C_{20} \epsilon^2 + C_{40} \epsilon^4 \right )  + {\cal O} \left (\epsilon^6 \right ), 
\\
\nonumber \\
S_1 &=   J \epsilon +  C_{31} M_{\rm s}^2 \epsilon^3  + {\cal O} \left (\epsilon^5 \right ),  
\\
\nonumber \\
M_2 & =  -\left ( \frac{J^2}{M_{\rm s}} + \frac{8}{5} C_{22} M^3_{\rm s} \right ) \epsilon^2 -  C_{42} M_{\rm s}^3 \epsilon^4 + {\cal O} \left (\epsilon^6 \right ), 
\\
\nonumber \\
S_3 & = C_{33} M_{\rm s}^4 \epsilon^3  + {\cal O} \left (\epsilon^5 \right ), 
\\
\nonumber \\
M_4 & =  \left [ \frac{107}{105} \frac{J^4}{M_{\rm s}^3} + \frac{9428}{735} C_{22} J^2 M_{\rm s} - \frac{4}{7} C_{33} J M^3_{\rm s}
\right.
\nonumber \\
&\left. \quad + \left (  C_{44} + \frac{3144}{245} C_{22}^2 \right ) M_{\rm s}^5 \right ] \epsilon^4 + {\cal O} \left (\epsilon^6 \right ).
\end{align}
The mass correction terms can be set to zero, $C_{20} = C_{40} = 0$, so that we can replace  $M_{\rm s} \to M$ in all expressions henceforth. 
This amounts to a simple shift in the mass scale of the system without any further physical importance. 


The Kerr limit of the HT spacetime is retrieved for
\begin{align}
C_{22}  &= C_{42} = C_{31} = 0, 
\\
C_{33} &= - \frac{J^3}{M_{\rm s}^6}, ~C_{44}= - \frac{62}{105} \frac{J^4}{M^8_{\rm s}},
\end{align}
and corresponds to the following set of multipole moments [see Eq.~\eqref{nohair}]: 
\begin{equation}
 S_1^{\rm K} = J, ~M_2^{\rm K} = - \frac{J^2}{M}, ~S_3^{\rm K} = - \frac{J^3}{M^2}, ~M_4^{\rm K} = \frac{J^4}{M^3}.
\end{equation}
(In this and other formulae below, the index `K'  labels a Kerr metric quantity.)

Instead of the original $C_{nm}$ parameters, it is more convenient to work with a new set of dimensionless  
parameters $c_{nm}$, defined as
\begin{equation}
C_{nm} = c_{nm} \frac{J^n}{M^{2n}}.
\end{equation}
We then define the Kerr \emph{deviation} parameters $\varepsilon_{nm}$ as
\begin{equation}
c_{nm} = c_{nm}^{\rm K} + \varepsilon_{nm}, 
\label{Kerrdeviation_params}
\end{equation}
with $\left (  c_{22}^{\rm K},  c_{42}^{\rm K} ,  c_{31}^{\rm K} ,  c_{33}^{\rm K},  c_{44}^{\rm K} \right ) = ( 0, 0, 0, -1, - 62/105) $.


\subsection{The equatorial shadow radius}

It is straightforward to use the general expressions~\eqref{photonEq} and \eqref{bgeneral} in the HT spacetime and obtain the light ring radius $r_{\rm ph}$ and 
the associated impact parameter  $b_{\rm ph} $ as expansions in the dimensionless spin parameter $a = J/M^2$ (or in $\epsilon$). For example, for the latter parameter 
we find
\begin{widetext}
\begin{align}
\frac{b_{\rm ph}}{M} &= 3\sqrt{3} - 2 a + \frac{1}{2\sqrt{3}} \left [ -1 + 18 c_{20} + 9 c_{22} (-16 + 15\log3) \right ] a^2  
 \nonumber \\
 & + \left  [ \frac{425957471}{15552} + \frac{605}{8} c_{20} + \frac{27}{2} c_{40} + \frac{585}{8} c_{42} + \frac{1534869}{128} c_{44} 
 - \frac{6380937}{256} \log 3 - \frac{2025}{32} \log 3\, c_{20}  - \frac{2025}{32} \log 3 \, c_{42} \right.
 \nonumber \\
& \left.   - \frac{5587785}{512} \log 3\, c_{44}  + \frac{9}{2}  ( 45 \log 3 -52) c_{20} c_{22} + \frac{c_{31}}{48} ( 6075 \log 3 -7132) 
- \frac{c_{33}}{384}  ( 7096815 \log 3 -7798684 ) \right.
\nonumber \\
& \left. - \frac{13}{40320}  ( 510678081\log 3 -561143860  ) c_{22} + \frac{3}{1120}  \left (41913428 - 39495033 \log 3 + 1225854 \log 3^2 \right ) c_{22}^2
\right ] a^3 
\nonumber \\
& + \frac{a^4}{19595520 \sqrt{3}} \Big [ -176359680 c_{20}^2 - 11664 \left ( 2520553852 - 3237502203 \log 3 + 858538926 \log 3^2 \right )   c_{22}^2
\nonumber \\
& + 20412 c_{20} \left \{ 8417925 \log 3 -9265540 + 8640 c_{31}   +   ( 8037225 \log 3 - 8848980 )  c_{33} + 48 c_{22} ( 500175 \log 3 -550652) \right \} 
\nonumber \\
& - 7 \left \{ 588719779672 + 50388480 c_{40} + 1214222400 c_{42} + 272213921700 c_{44} - 535888505070 \log 3 
\right. 
\nonumber \\
&  \left. - 1110121200  \log 3\, c_{42} - 247788794925 \log 3\, c_{44}  + 38880  ( 57105 \log 3 -62576) c_{31} - 16200 ( 24047847 \log 3 -26418796)  c_{33} \right \}
\nonumber \\
&  + 108 c_{22} \left \{ -292460091548 + 296252520363 \log 3 - 27342639450 \log3^2 + 181440 c_{31} ( 297 \log 3 -310 ) 
\right.
\nonumber \\
&  \left. - 39690 c_{33} \left (842680 - 1523826 \log 3 + 688905 \log3^2 \right ) \right \} \Big ] + {\cal O} (a^5).
\label{bphHT}
\end{align}
\end{widetext}
With the help of the Kerr limit of this expression we can benchmark the precision of the HT metric if we compare it against the full Kerr result~\eqref{bphKexact}. 
Although these formulae assume different coordinate systems, the comparison is nevertheless meaningful because $b_{\rm ph}$ is itself a gauge-invariant
quantity. The outcome of this exercise is shown in Fig.~\ref{fig:KerrTest} where we plot $b_{\rm ph} (a)$.
\begin{figure}[htb!]
\begin{center}
\includegraphics[width=\columnwidth]{KerrHTvsExact.pdf}
\end{center}
\caption{Comparison of the light ring impact parameter $b_{\rm ph}$, calculated using the Kerr limit of the HT spacetime, Eq.~\eqref{bphHT},
and the exact Kerr metric, Eq.~\eqref{bphKexact}, across the allowed spin range $-1 \leq  a \leq 1$. The prograde (retrograde) $b_{\rm ph}$ 
corresponds to $a >0$ ($a<0$).}
\label{fig:KerrTest}
\end{figure}
It is clear that, at least as far as $b_{\rm ph}$ is concerned, the ${\cal O} (J^4 )$ HT metric performs extremely well for
 $ -0.8 \lesssim a \lesssim 0.8 $. We expect the same level of accuracy to be representative of the 
 HT spacetime in general. 


Moving away from Kerr, in Fig.~\ref{fig:bphHT} we show $b_{\rm ph}$ as a function of $a$ for fixed $ \varepsilon_{nm}$
(top panel) and as a function of the deviation parameters for $a = 0.5$  (bottom panel). 
According to the displayed results, a moderate deviation $ | \varepsilon_{nm} | \sim 1 $ from Kerr is more than
enough to cause a notable change in $b_{\rm ph}$, especially for the prograde case, even for a moderate value of the spin. 
As expected from the expansions~\eqref{gexpand1} and \eqref{gexpand2}, the quadrupole parameter $\varepsilon_{22}$ is the one 
typically associated with the largest deviation. 
\begin{figure}[htb!]
\begin{center}
\includegraphics[width=\columnwidth]{bphHT_spin.pdf}
\\ \vspace{0.5cm}
\includegraphics[width=\columnwidth]{bphHT_deform.pdf}
\end{center}
\caption{The light ring impact parameter $b_{\rm ph}$ of the HT spacetime as a function of the dimensionless spin $a=J/M^2$
and the Kerr deviation parameters $\varepsilon_{nm}$ defined in Eq.~\eqref{Kerrdeviation_params}. Top panel: 
individually fixed $\varepsilon_{nm} <0$ and varying $a$ (both prograde $a>0$ and retrograde $a<0$ cases). 
The plot also shows the Schwarzschild impact parameter $b_{\rm Schw}$ (thin horizontal line) as well as the Kerr limit of the HT $b_{\rm ph}$ (dashed curve).
Bottom panel: fixed spin $a=0.5$ and individually varying deviation parameters $ x=\varepsilon_{nm}$.}
\label{fig:bphHT}
\end{figure}

The key parameter in relation with the shape of the shadow is the averaged prograde-retrograde impact parameter $\bar{b}_{\rm ph}$
as defined in Eq.~\eqref{bphaver}. This parameter is obtained most easily by simply removing the odd-order spin terms in an expansion like 
the one in Eq.~\eqref{bphHT}; these are the terms that cancel out when we sum $b_{\rm pro}$ and $b_{\rm retro}$. 
 As already discussed, $\bar{b}_{\rm ph}$ coincides with the equatorial radius of the shadow cast by a backlit black hole, as viewed by an equatorial observer.  
This is displayed in Fig.~\ref{fig:meanbHT1} in the form  $\bar{b}_{\rm ph}/b_{\rm Schw}$ (where $ b_{\rm Schw} = 3 \sqrt{3} M$ is the radius of the Schwarzschild 
black hole's circular shadow) as a function of $a$ 
and for different values of the deviation parameters $ \varepsilon_{nm} $ (each one `switched on' individually). The top (bottom) panel shows 
results for $\varepsilon_{nm} >0$ ( $\varepsilon_{nm} < 0$); in general, flipping sign in a given $\varepsilon_{nm}$ causes $\bar{b}_{\rm ph}$ 
to move to the opposite side of the circular $b_{\rm Schw}$ radius. All $\varepsilon_{nm}$ parameters appear to be correlated in this respect, 
with the exception of $\varepsilon_{31}$ which is anticorrelated. 
\begin{figure}[htb!]
\begin{center}
\includegraphics[width=\columnwidth]{meanRadiusHT1.pdf}
\\ \vspace{0.2cm}
\includegraphics[width=\columnwidth]{meanRadiusHT2.pdf}
\end{center}
\caption{The HT spacetime's equatorial shadow radius $\bar{b}_{\rm ph} = (b_{\rm pro} + b_{\rm retro} )/2$ (normalised to the Schwarzschild 
radius $b_{\rm Schw}$) as a function of the spin $a$ for $\varepsilon_{nm} <0$ and  $\varepsilon_{nm} > 0$ Kerr
deviation parameters (top and bottom, respectively). In all cases only a specific $\varepsilon_{nm}$ is set to a non-zero value. 
For the purpose of comparison the figure includes the radius  $\bar{b}_{\rm ph}$ for the exact Kerr metric (solid black curves) and for
the Kerr limit of the HT metric (dashed curves).}
\label{fig:meanbHT1}
\end{figure}

The results of Fig.~\ref{fig:meanbHT1} suggest that a moderate $ | \varepsilon_{nm} | \sim 0.5$ deviation could produce a minimum $ 10\%-20 \%$ 
variation with respect to a circular shadow for $ 0.6 \lesssim a \lesssim 0.8$. At the same time the Kerr `yardstick' discussed earlier suggests 
that the accuracy of the HT shadow radius begins to deteriorate at  $a \approx 0.8$. 
Once again, the quadrupole deformation $\varepsilon_{22}$ stands out as the shadow's dominant decircularisation factor. 
The moderate degree of variation of the equatorial radius with $a$ implies that the corresponding change in the shadow's \emph{polar} radius 
should be even less pronounced, thus essentially retaining its zero-rotation value $ b_{\rm Schw}$. 


A perhaps surprising situation could arise when several of the deformation parameters are present at the same time. 
For instance, it is fairly easy to produce a very Kerr-like shadow that remains nearly circular for the entire spin range  
despite having a strongly non-Kerr multipolar structure. An example is shown in Fig.~\ref{fig:meanbHT2} where a 
dominant $\varepsilon_{22} < 0$ is counterbalanced by the combined presence of  $\varepsilon_{42}, \varepsilon_{44}, \varepsilon_{33} >0$.
What this really means is that the connection between the shadow shape and the multipole moments of the spacetime (or their proxies 
$\varepsilon_{nm}$) is not as direct as suggested in previous work \cite{Johannsen2010ApJ, Broderick2014ApJ}.  This important issue is further 
explored in the following section where we express $b_{\rm ph}$  and the shadow in terms of the multipole moments themselves. 
\begin{figure}[htb!]
\begin{center}
\includegraphics[width=\columnwidth]{meanRadiusHT3.pdf}
\end{center}
\caption{An example of a HT spacetime with a nearly circular shadow despite a strongly non-Kerr structure.  We show the
normalised equatorial shadow radius $\bar{b}_{\rm ph}/b_{\rm Schw}$ as a function of the spin $a$ for Kerr deviation parameters
$ ( \varepsilon_{22},  \varepsilon_{33}, \varepsilon_{42}, \varepsilon_{44} ) = (-0.1, 0.5, 0.5, 0.5) $.  The HT result (designated as `non-Kerr') 
is compared against the shadow radius of the full Kerr spacetime (solid curve) and the Kerr limit of the HT spacetime (dashed curve).
}
\label{fig:meanbHT2}
\end{figure}



\subsection{The shadow radius as a function of the multipole moments  $M_\ell, S_\ell$.}
\label{sec:bphMM}

The previous section provided us with some understanding of the dependence of the equatorial 
shadow radius as a function of the original $C_{nm}$ parameters of the HT metric, expressed as deviations from 
the Kerr spacetime. Here we complete this analysis by expressing the shadow radius directly in terms of the Geroch-Hansen 
multipole moments $\{ M_\ell, S_\ell \}$. 


As a first step, we define a new set of dimensionless moments, $ \{ {\cal M}_\ell, {\cal S}_\ell \}$,
by factoring out the spin dependence together with the suitable powers of $M$,
\begin{equation}
S_\ell = {\cal S}_\ell \frac{J^\ell}{M^{\ell-1}}, \quad M_\ell = {\cal M}_\ell \frac{J^\ell}{M^{\ell-1}}.
\label{MMnondim}
\end{equation}
The next step is to express $J, C_{nm} $ in terms of these new parameters and the spin $a$. 
In particular, $J$ is replaced with the help of  $J \to S_1 - C_{31} M^2 \epsilon^2 $ which means that the parameter 
$c_{31}$ is expected to be present in $b_{\rm ph}$. The same is true for $c_{42}$ after replacing $ C_{22} \to C_{22} ( {\cal M}_2, J, c_{42} )$. 

Once $b_{\rm ph} (a,  {\cal M}_\ell, {\cal S}_\ell) $ is calculated, we can easily obtain the shadow radius $\bar{b}_{\rm ph}$. The 
final ${\cal O} (a^4)$ expression is
\begin{widetext}
\begin{align}
\frac{\bar{b}_{\rm ph}}{M} &= 3\sqrt{3} + \frac{a^2}{16 \sqrt{3}} \left [\, {\cal S}_1^2 ( 712 - 675 \log 3 )
+ 45 {\cal M}_2 (16 - 15 \log3 )\, \right ] + \frac{a^4}{\sqrt{3}} \left [\,  \frac{5}{8}  c_{42} (-622 + 567\log 3 )
\right.
\nonumber \\
\nonumber \\
& \left. + \frac{315}{256} {\cal M}_4 (-79028 + 71937 \log 3)   + {\cal S}_1^4 \frac{ \left (403363330924 + 2226769271865 \log 3 - 2361245507100 \log 3^2 \right )}{22394880}
\right.
\nonumber \\
\nonumber \\
& \left. + {\cal M}_2 {\cal S}_1^2 \frac{\left ( 2538443668 + 25485296247 \log 3 - 25302890160 \log 3^2 \right )}{82944}
+ {\cal S}_1 c_{31} \frac{(70118 - 64395 \log 3)}{72}
\right.
\nonumber \\
\nonumber \\
& \left. + \frac{35}{6912} {\cal S}_1 {\cal S}_3 \left (-18414152 - 3671082 \log 3 + 18600435 \log 3^2 \right)
+ \frac{35}{256} \frac{{\cal M}_2 {\cal S}_3}{{\cal S}_1} \left (842680 - 1523826 \log 3 + 688905 \log 3^2 \right ) 
\right.
\nonumber \\
\nonumber \\
& \left. + \frac{5}{6144} {\cal M}_2^2 \left (121182452 - 379769553 \log 3 + 245296836 \log 3^2 \right)
+ \frac{5}{8} \frac{ {\cal M}_2 c_{31}}{{\cal S}_1} (310 - 297 \log 3)
\, \right ].
\label{bphMM1}
\end{align}
The numerical value of this result is given by the much shorter expression:
\begin{align}
\frac{\bar{b}_{\rm ph}}{M} &\approx 5.19615  - a^2 \left ( 0.778098\, {\cal M}_2   + 1.06677\, {\cal S}_1^2 \right ) 
+ a^4 \Big (  0.329511\, c_{42} - 11.1427 \,{\cal M}_2^2 + 2.04045\, {\cal M}_4 
\nonumber \\
& - 16.6617\, {\cal M}_2 {\cal S}_1^2    - 4.72798 \, {\cal S}_1^4  -5.87737\, c_{31}  \frac{{\cal M}_2}{{\cal S}_1}
+  4.67328 \,  \frac{{\cal M}_2 {\cal S}_3 }{{\cal S}_1} 
-5.02887 \, c_{31}  {\cal S}_1 + 7.39031\,   {\cal S}_1 {\cal S}_3   \Big ).
\label{bphMM2}
\end{align}
\end{widetext}
Let us first focus on the ${\cal O} (a^2)$ portion of this formula; we can see that 
 a ${\cal M}_2 <0$ quadrupole counteracts the shadow shaping action of the `frame-dragging' multipole ${\cal S}_1$.
This is indeed what happens in Kerr, which has ${\cal S}_1^{\rm K} =1$ and ${\cal M}_2^{\rm K} =-1$ and is in agreement with
the analysis of Ref.~\cite{Broderick2014ApJ}.  A near perfect cancellation of the two effects takes place  for
\begin{equation}
{\cal M}_2 \approx - 1.371 {\cal S}_1^2,
\label{M2circHT}
\end{equation}
which means that non-Kerr bodies too can produce quasi-circular shadows. At the other end of the spectrum,
prolate bodies have ${\cal M}_2 > 0$ and are likely to cast a markedly non-circular shadow.  

The inclusion of the ${\cal O} (a^4 )$ term in the shadow radius opens the door to more possibilities. 
An interesting exercise in this respect is to set some of the parameters to their Kerr values,  
$S_1 =1, {\cal M}_2 = -1, c_{42} = c_{31}=0$. The resulting shadow radius now depends on the higher multipoles ${\cal S}_3, {\cal M}_4$,
\begin{align}
\frac{\bar{b}_{\rm ph}}{M} & \approx 5.19615 - 0.288675 a^2  +  \left (0.79105  \right.
\nonumber \\
& \left. \quad + 2.04045 {\cal M}_4  + \,2.71703 {\cal S}_3 \right )  a^4.
\label{baverM4S3}
\end{align}
The extent to which these two multipoles can decircularise the shadow can be understood by looking at the top panel of 
Fig.~\ref{fig:bHTcontour} where we show the contour plot of the fractional difference $ \Delta_b =  \bar{b}_{\rm ph}/b_{\rm Schw} -1$ (for $a=0.7$). 
For the $ (-1,3) \times (-3,1)$ `box' shown in the figure, the shadow radius can vary up to $\sim 40\%-50\%$ with respect to the circular 
Schwarzschild value. For the Kerr spacetime, ${\cal M}_4^{\rm K} =1, {\cal S}_3^{\rm K} = -1$, the deviation from a circular shadow is small, $| \Delta_b | \lesssim 0.1 $
(this is indicated in~Fig.~\ref{fig:bHTcontour} by a black dot). 
However, it is also clear that there is a high degree of degeneracy in the sense that  we can have $| \Delta_b | \ll 1 $ even for a markedly non-Kerr 
multipolar structure; an example is the point $({\cal M}_4, {\cal S}_3) = (-1, 0.5)$. 

We can repeat the same exercise by varying the mass moments ${\cal M}_2, {\cal M}_4$ while setting the rest of the parameters 
equal to their Kerr values (i.e., ${\cal S}_1= -{\cal S}_3 =1, c_{31}=c_{42}=0$). The associated shadow radius is
\begin{align}
\frac{\bar{b}_{\rm ph}}{M} & \approx 5.19615 -  ( 1.06677 +  0.778098 {\cal M}_2) a^2
\nonumber \\
&  + \left (\,  2.04045 {\cal M}_4 -12.1183 - 21.335 {\cal M}_2  \right.
\nonumber \\
& \left. - 11.1427 {\cal M}_2^2 \, \right  ) a^4.
\label{baverM2M4}
\end{align}
Compared to~\eqref{baverM4S3}, this radius displays a markedly larger ${\cal O} (a^4)$ piece. This effectively limits how much we can 
vary ${\cal M}_2$ while being consistent with the spin expansion character of $b_{\rm ph}$. For example, ${\cal M}_2 \gtrsim 1$ could easily result in a 
negative radius, which is clearly unphysical. In order to account for this limitation, the corresponding contour plot of $\Delta_b$ is calculated 
for a somewhat lower spin, $a=0.5$, see middle panel of Fig.~\ref{fig:bHTcontour}.
The dominant role of ${\cal M}_2$ as a shadow decircularisation factor is clearly visible in this figure: the contour lines are almost vertical and a variation 
across the $-2 \leq {\cal M}_2 \leq 0.5$ range causes a reduction in the equatorial radius up to $\approx 40\%$ with respect to the Kerr radius. 
At the same time, however, we can tweak ${\cal M}_4$ so that $\Delta_b$ remains close to its Kerr value even if ${\cal M}_2$ deviates from Kerr. 
An example of this is the point $({\cal M}_2,{\cal M}_4 ) = (-1.5,1.5)$ in Fig.~\ref{fig:bHTcontour}. 

Finally, we consider the shadow radius as a function of  ${\cal M}_2, {\cal S}_3$ with the rest of the parameters set equal to their Kerr values
(${\cal S}_1^{\rm K}= {\cal M}_4^{\rm K} =1, c_{31}=c_{42}=0$). The associated equatorial shadow radius is,
\begin{align}
\frac{\bar{b}_{\rm ph}}{M} & \approx 5.19615 - (1.06677 +  0.778098 {\cal M}_2 ) a^2 + 
\nonumber \\
& + ( -2.68753 - 11.1427 {\cal M}_2^2 + 7.39031 {\cal S}_3 
\nonumber \\
& -16.6617  {\cal M}_2  + 4.67328  {\cal M}_2  {\cal S}_3) a^4.
\end{align}
The quadrupole ${\cal M}_2$ is still the dominant factor (leading to a similar deviation from a spherical shadow as the previous case) 
but its decircularising influence  is much more easily counteracted by a variation in ${\cal S}_3$  than in ${\cal M}_4$. 

The upshot of this analysis, and in combination with the results of the previous section, is that a quasi-circular shadow
does \emph{not} provide a reliable test of the no hair-theorem [Eq.~\eqref{nohair}] as different combinations of multipole moments 
beyond the quadrupole can produce the same, more or less, deviation from circularity as the Kerr spacetime. 
In terms of the dimensionless parameters ${\cal M}_\ell, {\cal S}_\ell$, this statement entails comparable magnitude shifts
away from the Kerr moments. 

\begin{figure}[htb!]
\begin{center}
\includegraphics[width=0.73\columnwidth]{RshadowHTM4S3_fractional.pdf}
\\ \vspace{0.1cm}
\includegraphics[width=0.73\columnwidth]{RshadowHTM2M4_fractional.pdf}
\\ \vspace{0.1cm}
\includegraphics[width=0.73\columnwidth]{RshadowHTM2S3_fractional.pdf}
\end{center}
\caption{Contour plot of the equatorial HT shadow radius expressed as a fractional difference 
$ \Delta_b =  \bar{b}_{\rm ph}/b_{\rm Schw} -1 $ with respect to the circular Schwarzschild radius. 
Top panel: $\Delta_b$ is shown as a function of the ${\cal M}_4, {\cal S}_3$ multipole moments for spin  $a=0.7$,
see Eq.~\eqref{baverM4S3}. The remaining moments have been set equal to their Kerr values, 
$ {\cal M}_2 = -1, {\cal S}_1 = 1$.
Middle panel: $\Delta_b$ is shown as a function of the ${\cal M}_2, {\cal M}_4$ mass moments for $a=0.5$,
see Eq.~\eqref{baverM2M4}. The remaining moments have been set equal to their Kerr values,  
${\cal S}_1 =1, {\cal S}_3= -1$.
Bottom panel: $\Delta_b$ is shown as a function of the ${\cal M}_2, {\cal M}_4$ mass moments for $a=0.5$,
see Eq.~\eqref{baverM2M4}. The remaining moments  ${\cal S}_1, {\cal S}_3$  have been set equal to their Kerr values. 
In all cases, we have set the remaining metric constants $c_{42}, c_{31}$ equal to zero.
The black dot marks the HT Kerr  limit of $\Delta_b$.}
\label{fig:bHTcontour}
\end{figure}




\section{Shadow radius in a general stationary-axisymmetric spacetime}
\label{sec:axisym}

As a second case study of the shadow produced by a non-Kerr spacetime, in this section we consider a general
stationary-axisymmetric spacetime in GR with an arbitrary multipolar structure. This model is hinged on a double expansion; 
first with respect to the multipole moment order (which in practice is a slow rotation approximation if an $\ell$-pole is assumed to scale 
as $\sim J^\ell$) and subsequently with respect to $M/r$ for a given multipole moment order. 
The moments themselves are written with the Kerr and non-Kerr parts separated, 
\begin{equation}
M_\ell = M_\ell^{\rm K} + \delta M_{\ell},  \quad S_\ell = S_\ell^{\rm K} + \delta S_\ell, \nonumber
\end{equation}
and the former part is resummed so that the full Kerr part of the metric is recovered. This step ensures that the metric has the correct 
Schwarzschild and Kerr limits when, respectively,  rotation is turned off or the moments are set to their Kerr values.
In essence, this metric is what comes out when we superimpose Ryan's post-Newtonian metric of arbitrary multipole moments~\cite{Ryan:1995}
with the Kerr spacetime, making sure not to count terms twice (details on the construction of this spacetime will appear elsewhere). 
This makes it an ideal complement to the HT metric which is fully accurate at each spin order but is truncated to a lower multipole order than 
the metric discussed here. 
This spacetime is a direct descendant of the metric first introduced in Ref.~\cite{Maselli:2019qbf}, which was designed with an arbitrary set of 
multipole moments up to $M_4$ and the correct Schwarzschild limit. The present model improves on that earlier construction by pushing the 
expansion beyond $M_4$ and having the correct Kerr limit. Therefore, the new metric can be used as a parametrised non-Kerr solution within GR, 
with arbitrary deformations in the multipole moments.
  

A general stationary-axisymmetric vacuum spacetime can be constructed algorithmically via the Ernst potential formalism~\cite{ernst68}. 
For the case at hand, we employ the cylindrical-like Weyl-Papapetrou coordinates $\{t,\rho,z, \phi\}$ and the resulting line element takes the 
form~\cite{papapetrou53}
\begin{align}
ds^2  &=  -f \left(dt-\omega d\phi\right)^2  \nonumber
\\ 
& \quad + f^{-1} \left[ e^{2\gamma} (d\rho^2+dz^2)+\rho^2d\phi^2 \right].
\label{ds2axisym}
\end{align}
The three metric potentials $\{f, \omega, \gamma\}$ are functions of $\{ \rho, z \}$; they can be written in a separable form, comprising Kerr and 
non-Kerr parts (their full functional forms can be found in Appendix~\ref{sec:Ryanmetric}),
\begin{align}
 f &= f_{\rm K}  +\delta f,
 \label{eq:axisym1}
 \\
 \omega &=  \omega_{\rm K} + \delta\omega,
\label{eq:axisym2}
\\
e^{2  \gamma} &=  e^{2 \gamma_{\rm K}} \left ( 1+ \delta \gamma \right ).
\label{eq:axisym3}
\end{align}
The functions $\{ f_{\rm K}, \omega_{\rm K}, \gamma_{\rm K} \}$ represent the baseline Kerr solution in Weyl-Papapetrou coordinates and depend on $M$ 
and the dimensionless  Kerr spin parameter  $a=S_1/M^2$ (not to be confused with the Kerr parameter $a = J/M^2$ used elsewhere in the paper).  

The arbitrary deformations $\{ \delta M_\ell, \delta S_\ell \}$ in the multipole moments enter through the non-Kerr corrections 
$\{ \delta f, \delta \omega, \delta \gamma \}$. Once the metric is known, we can use the general formulae~\eqref{photonEq} and \eqref{bgeneral} 
to obtain $r_{\rm ph}$ and $b_{\rm ph}$ as expansions in the spin/multipolar order parameter $a$. The results take a somewhat simpler form when written 
in terms of the dimensionless multipolar deformation (this definition assumes that any deformation away from Kerr is rotation induced),
\begin{equation}
\delta {\cal M}_\ell \equiv \frac{\delta M_\ell}{a^\ell M^{\ell+1}}, \qquad  \delta {\cal S}_\ell \equiv \frac{\delta S_\ell}{a^\ell M^{\ell+1}}.
\end{equation}
After a numerical evaluation, we arrive at the following expressions:
\begin{widetext}
\begin{align}
\frac{b_{\rm ph}}{M} &=3 \sqrt{3}-2 a  -1.19772 \Big ( 0.646842 \, \delta {\cal M}_2 +0.24102  \Big )  a^2 
-0.0800267  \Big( 16.1044  \, \delta {\cal M}_2  -4.47113 \, \delta {\cal S}_3 
\nonumber \\
& + \, 1.85123 \Big )  a^3 
+ 6.16711  \Big ( 0.0321565 \, \delta {\cal M}_4 -0.0838228\, ( \delta {\cal M}_2) ^2 -0.342846 \, \delta {\cal M}_2   
+ 0.126939\, \delta {\cal S}_3
\nonumber\\ 
& -\, 0.0151695 \Big ) a^4  -31.0214   \Big( 0.0951536\, ( \delta {\cal M}_2)^2  -0.0199417 \, \delta {\cal M}_2\, \delta {\cal S}_3 +0.163008\, \delta {\cal M}_2  
\nonumber \\
& - \, 0.0189\, \delta {\cal M}_4 -0.0357977 \, \delta {\cal S}_3 + 0.00330828\, \delta {\cal S}_5 + 0.00212252 \Big ) a^5  +O\left(a^6\right),
 \\
 \nonumber \\

\frac{\bar{b}_{\rm ph}}{M} &=3 \sqrt{3}  -1.19772 \Big ( 0.646842 \,\delta {\cal M}_2 + 0.24102  \Big ) a^2  
 + 6.16711  \Big (-0.0838228\, ( \delta {\cal M}_2 )^2 -0.342846\, \delta {\cal M}_2 
 \nonumber \\
&  \quad +\, 0.0321565 \,\delta {\cal M}_4 + 0.126939 \, \delta {\cal S}_3 - 0.0151695 \Big ) a^4  +O\left(a^6 \right).
\end{align}
\end{widetext} 
As pointed out earlier, $\bar{b}_{\rm ph}$ is given by the same expression as $b_{\rm ph}$ after removing the odd powers of $a$.


These results look very similar to those obtained in the HT spacetime [e.g., Eq.~\eqref{bphMM2}] but, as expected, are not identical. 
By construction, the general metric of this section is a function of the multipole moments only, while the HT metric contains the structure constants $C_{nm}$.
As we have seen, it is not possible to replace all of these parameters with $M_\ell, S_\ell$. 
The proximity of the two $\bar{b}_{\rm ph}$ results can be gauged if we calculate the ${\cal M}_2$ that makes the shadow circular at ${\cal O} (a^2)$. 
We find 
\begin{equation}
\delta {\cal M}_2 \approx - 0.373 ~\Rightarrow ~ {\cal M}_2 \approx -1.373 {\cal S}_1^2,
\end{equation}
which lies very close to our earlier HT result~\eqref{M2circHT}.


Without further ado we calculate the degree of deviation of this new $\bar{b}_{\rm ph}$ from the circular radius $b_{\rm Schw}$ and 
the Kerr shadow radius. In Fig.~\eqref{fig:shadowG_frac} we show contour plots of the fractional difference  $ \Delta_b =  \bar{b}_{\rm ph}/b_{\rm Schw} -1$ 
in the shadow radius for a dimensionless spin $a=0.7$ and three complementary choices for the multipole moments: 
(i) We fix $M_2 = M_2^{\rm K}$ and vary $\delta {\cal S}_3,  \delta {\cal M}_4$ (top panel). In both cases, the Kerr limit 
is indicated by a black dot and is clearly seen to represent a nearly circular shadow for the chosen spin value.
(ii) We fix $S_3 = S_3^{\rm K}$ and vary $\delta {\cal M}_2, \delta {\cal M}_4$ (middle panel).
(iii) We fix $M_4 = M_4^{\rm K}$ and vary $\delta {\cal M}_2, \delta {\cal S}_3$ (bottom panel).

Considering first $\bar{b}_{\rm ph}$ as a function of the higher multipoles $M_4, S_3$, we can observe a much smaller
variation in $\Delta_b$ ($\lesssim 15\% $) with respect to the previous HT case, for the same$(-2,2)\times (-2,2)$ `box' and spin. 
In other words, for the spacetime discussed in this section the shadow remains nearly circular for a wide range of deviation
around the $S_3, M_4$ moments. The degree of decircularisation of the shadow becomes much higher 
($\lesssim 50 \%$ for the range shown in Fig.~\ref{fig:shadowG_frac}), and similar to the one found in the HT spacetime, 
when the mass quadrupole $M_2$ is varied. In all cases, however, we can draw the same conclusion as before: we can shift two 
or more multipoles away from Kerr while maintaining a Kerr-like shadow radius.  
\begin{figure}[htb!]
\begin{center}
\includegraphics[width=0.36\textwidth]{RshadowG_S3M4_fractional.pdf}
\\ \vspace{0.45cm}
\includegraphics[width=0.36\textwidth]{RshadowG_M2M4_fractional.pdf}
\\ \vspace{0.45cm}
\includegraphics[width=0.36\textwidth]{RshadowG_M2S3_fractional.pdf}
\end{center}
\caption{Contour plot of the equatorial shadow radius expressed as a fractional difference $ \Delta_b =  \bar{b}_{\rm ph}/b_{\rm Schw} -1 $
with respect to the circular Schwarzschild radius. Top panel:  $\Delta_b$ as a function of the $\delta {\cal S}_3, \delta {\cal M}_4$ deviation 
parameters, assuming $\delta {\cal M}_2= 0$ (i.e. ${\cal M}_2 = {\cal M}_2^{\rm K} = -1$). 
Middle panel: $\Delta_b$ as a function of the $\delta {\cal M}_2, \delta {\cal M}_4$ deviation parameters, assuming $\delta {\cal S}_3 = 0$ 
(i.e., ${\cal S}_3 = {\cal S}_3^{\rm K} =-1$). 
Bottom panel: $\Delta_b$ as a function of the $\delta {\cal M}_2, \delta {\cal S}_3$ deviation parameters, assuming $\delta {\cal M}_4 = 0$ 
(i.e., ${\cal M}_4 = {\cal M}_4^{\rm K} =1$). 
In all cases, the spin is set to $a=0.7$. The black dot marks the Kerr limit of $\Delta_b$.}
\label{fig:shadowG_frac}
\end{figure}


\section{The Johannsen metric: multipole moments and shadow radius}
\label{sec:Jmetric}


\subsection{Extracting the multipole moments}
\label{sec:Jmoments}

The widely used Johannsen metric~\cite{Johannsen2013} (hereafter `J metric', see Appendix~\ref{sec:Jexpansions}) is an example of a deformed Kerr spacetime 
with the characteristic property of separability, that is, it admits a Carter-like constant. The purpose of this section is to study the multipolar structure of the J spacetime 
in relation to its black hole shadow.

In reality, this is an ill-defined objective; the formal calculation of a spacetime's multipole moments in a given theory of gravity requires the use of 
field equations and appropriate asymptotic conditions~\cite{Pappas:2012, pappas15a}. No such field equations are available for the J spacetime for the simple reason that
it is not a solution of GR or any other known theory of gravity. 

In the absence of field equations that could determine the multipole moments, we follow a more practical approach by 
deriving post-Newtonian expansions for the orbital frequencies of a test body in a small eccentricity/inclination orbit.
These expressions can be compared against the ones obtained by Ryan~\cite{Ryan:1995} for a stationary-axisymmetric spacetime 
of arbitrary multipolar structure in GR with the aim of extracting some information about the multipolar structure of the J spacetime 
from the  coefficients of the post-Newtonian expansion.
Two of the body's equations of motion are
\begin{equation}
u^t = \frac{1}{{\cal D}} \left (  g_{t\varphi} L + g_{\varphi\varphi} E \right ), ~
u^\varphi = -\frac{1}{{\cal D}} \left (  g_{t\varphi} E  + g_{tt} L  \right ).
\end{equation}
In addition, we can define an effective potential (not to be confused with the one used earlier for null geodesics)
\begin{align}
 & g_{rr} (u^r )^2  + g_{\theta\theta} \left ( u^\theta \right )^2 
  = \frac{1}{{\cal D}} \left (\, g_{tt} L^2 + 2 g_{t\varphi} L E  + g_{\varphi\varphi}  E^2\, \right )
 \nonumber  \\
 & -1  \equiv   V_{\rm eff} (r,\theta,E,L).
\label{norm_gen}
\end{align}
Circular equatorial orbits ($\theta=\pi/2$) are required to solve
\begin{equation}
V_{\rm eff} (r_0, E, L) = V_{\rm eff}^\prime (r_0, E, L) = 0,
\end{equation}
where $r_0$ denotes the orbital radius. These two conditions can be solved for any pair of parameters; in
the present case, we are interested in the angular frequency $\Omega = u^\varphi/u^t$. Making the substitutions 
\begin{equation}
E= - u^t \left  ( g_{tt} + \Omega g_{t\varphi} \right ), \qquad  L=  u^t \left  ( g_{t\varphi} + \Omega g_{\varphi\varphi} \right ).
\end{equation}
we can solve the above system for $\{ \Omega, u^t \}$. The $V^\prime_{\rm eff} =0$ equation becomes a binomial,
\begin{equation}
g_{\varphi\varphi}^\prime \Omega^2 + 2 g_{t\varphi}^\prime \Omega + g_{tt}^\prime = 0, 
\end{equation}
with solutions
\begin{equation}
\Omega_{\pm} = \frac{1}{g_{\varphi\varphi}^\prime} \left ( -  g_{t\varphi}^\prime \pm \sqrt{  g_{t\varphi}^{\prime 2} -  g_{tt}^\prime  g_{\varphi\varphi}^\prime }\right ).
\end{equation}
The condition $ V_{\rm eff} = 0$ yields the `redshift' formula,
\begin{equation}
u^t = \left (- g_{tt}  -g_{\varphi\varphi} \Omega^2 - 2 g_{t\varphi} \Omega  \right )^{-1/2}.
\end{equation}
A perturbed circular orbit is characterised by a small inclination/eccentricity. In addition to $\Omega$, we now have to consider
the epicyclic orbital frequencies $\Omega_r, \Omega_\theta$. To leading order in the perturbation, these are given by,
\begin{equation}
\Omega_r = \frac{1}{u^t} \sqrt{-\frac{V_{\rm eff}^{\prime\prime} }{2 g_{rr}}}, \qquad  
\Omega_\theta = \frac{1}{u^t} \sqrt{-\frac{ \partial_\theta^2 V_{\rm eff}}{2 g_{\theta\theta}}},
\label{epicyclicOm}
\end{equation}
where the right-hand-side terms are to be evaluated for the unperturbed circular equatorial orbit.  
The orbital precession frequencies are defined as
\begin{equation}
\omega_i \equiv \Omega -\Omega_i, \qquad i = \{r,\theta\}.
\label{precessOm}
\end{equation}
These general formulae can now be applied to the J metric (the detailed form of this metric is given in Appendix~\ref{sec:Jexpansions}). 
Using $\epsilon$ as a collective bookkeeping parameter for the metric's deformation parameters $ \{ \alpha_{13}, \alpha_{22}, \alpha_{52}, \varepsilon_3 \}$, 
we can write $\Omega$ as an expansion in $\epsilon$ and $M/r_0$. The result is Eq.~\eqref{Omexpan} in Appendix~\ref{sec:Jexpansions}.
The $\Omega$ expansion can be inverted and furnishes $M/r_0$ as an expansion in $\epsilon$ and the orbital velocity $v^3 = M \Omega$, see Eq.~\eqref{vexpan}.

This latter result can be subsequently used in Eqs.~\eqref{epicyclicOm}-\eqref{precessOm} to obtain expansions for the precession frequencies.
For the normalised frequencies  $\tilde{\omega}_i = \omega_i /\Omega$ we find (here $a = J/M^2$ is the Kerr spin parameter)
\begin{widetext}
\begin{align}
\tilde{\omega}_r & = 3 v^2 - 4a  v^3  + \frac{1}{2}  \left [ 3 \left ( 3 + a^2 \right ) + \epsilon ( 6\alpha_{13} -\alpha_{52} - 3\varepsilon_3 ) \right ] v^4
+ \frac{1}{2}  \left [ 27 + 17 a^2 + 3 \epsilon ( \alpha_{52} -4\alpha_{13} + 5 \varepsilon_3 ) \right ] v^6 
\nonumber \\
&  -5a  \left ( 2 + \epsilon \alpha_{22} \right ) v^5  + a \left [ -4 \left ( 12 + a^2 \right ) + \frac{1}{3} \epsilon ( 3 \alpha_{13} -21\alpha_{22} \
- 8\alpha_{52} -15 \varepsilon_3 )  \right ] v^7  
+ {\cal O} \left (\epsilon^3, v^8 \right ), 
\\
\nonumber \\
\tilde{\omega}_\theta & =  2a v^3 - \frac{3}{2} a^2 v^4 + \epsilon  a \alpha_{22}  v^5
+ 4a^2 v^6 + a \left [ - 5a^2 + \epsilon ( 3\varepsilon_3 -5 \alpha_{13} ) \right ] v^7 
+ {\cal O} \left (\epsilon^3, v^8 \right ).
\end{align}
\end{widetext}
Note that no ${\cal O} (\epsilon^2)$ terms appear in these results. These expressions can be compared against 
the multipole moment expansions of Ref.~\cite{Ryan:1995} for $\tilde{\omega}_i$, derived for orbits in an 
arbitrary axisymmetric-stationary spacetime in GR:
\begin{widetext}
\begin{align}
( \tilde{\omega}_r)_{\rm GR}  & = 3v^2 - \frac{4 S_1}{M^2} v^3 + \frac{3}{2} \left ( 3-\frac{M_2}{M^3} \right ) v^4
- 10 \frac{S_1}{M^2} v^5 + \left ( \frac{27}{2} - 2  \frac{S_1^2}{M^4} - \frac{21}{2}  \frac{M_2}{M^3}  \right ) v^6 
\nonumber \\
& \quad + \left ( -48  \frac{S_1}{M^2}  - 5   \frac{S_1}{M^2} \frac{M_2}{M^3} + 9 \frac{S_3}{M^4}   \right ) v^7 
 + {\cal O} \left ( v^8 \right ), 
\\
\nonumber \\
( \tilde{\omega}_\theta )_{\rm GR} & = \frac{2 S_1}{M^2} v^3 + \frac{3}{2} \frac{M_2}{M^3} v^4 + 
\left ( 7 \frac{S_1^2}{M^4} + 3 \frac{M_2}{M^3} \right ) v^6 
+ \left (  11 \frac{S_1}{M^2} \frac{M_2}{M^3} - 6 \frac{S_3}{M^4} \right ) v^7 + {\cal O} \left ( v^8 \right ).
\end{align}
\end{widetext}
Inspection of the first four $v^n$ pairs leads to the following identifications [the first (second) entry in each line corresponds to 
 $\tilde{\omega}_r$ ($\tilde{\omega}_\theta$)]:
\begin{align}
{\cal O} (v^3): & ~ S_1 = a M^2,
\\
{\cal O} (v^4): &~ \frac{M_2}{M^3} =  a^2  -\frac{1}{3} \epsilon( 6 \alpha_{13} - \alpha_{52} - 3 \varepsilon_3 ), 
\nonumber \\
& ~ \frac{M_2}{M^3} = - a^2
\\
{\cal O} (v^5): & ~ \mbox{unbalanced terms}~ \left \{ 5 a \alpha_{22}, -a  \alpha_{22} \right \},
\\
{\cal O} (v^6): &~ \frac{M_2}{M^3} = - a^2 - \frac{1}{7}  \epsilon ( \alpha_{52} -4 \alpha_{13} + 5 \varepsilon_3 ),
\nonumber \\
& ~  \frac{M_2}{M^3} = - a^2.
\end{align}
The predicted quadrupole moment from the ${\cal O} (v^4)$ terms is unique and identical to the Kerr value, $M_2 = -a^2 M^3$,  
provided we set $\alpha_{52} = 6 \alpha_{13}  - 3 \varepsilon_3$.
Similarly, in order to avoid inconsistency in the ${\cal O} (v^5) $ terms we must set $\alpha_{22} = 0$. Using this information in the ${\cal O} (v^6)  ~\mbox{terms}$
fixes one more parameter, $\alpha_{13} = - \varepsilon_3$. Finally, the highest-order terms lead to
\begin{equation}
{\cal O} (v^7): ~ \frac{S_3}{M^4} = - a^3  +  2a  \varepsilon_3,  \quad \frac{S_3}{M^4} = - a^3  - \frac{4}{3} a \varepsilon_3.
\end{equation}
We have thus hit an inconsistency wall: setting $\varepsilon_3 =0$ makes the rest of the parameters vanish and the metric
reduces to Kerr. Based on our earlier comment this result is hardly surprising; the J metric is not a solution of the vacuum GR equations
and therefore we should not have expected a fully consistent correspondence with  Ryan's post-Newtonian multipolar expansion. 



\subsection{Shadow radius in the Johannsen metric}
\label{sec:Jshadow}

The second (and last) part of our discussion of the J spacetime is concerned with the black hole's shadow (previous work 
on the subject can be found in~\cite{Johannsen2013ApJ}). 

The light ring radius and impact parameter are calculated with the help of Eqs.~\eqref{photonEq} and \eqref{bgeneral}, 
without expanding in the deformation parameters. As it turns out both quantities depend only on two parameters, 
$\{\alpha_{13}, \alpha_{22}\}$. 

As already discussed, the equatorial shadow radius is given by the prograde-retrograde averaged impact parameter,
$ \bar{b}_{\rm ph} =  \left ( b_{\rm pro} + b_{\rm retro} \right )/2$. The ratio $\bar{b}_{\rm ph}/b_{\rm Schw}$ is shown in 
Fig.~\ref{fig:Jshadow} as a function of the dimensionless spin $a$ (top panel) and of the deformation parameters (bottom panel). 
\begin{figure}[htb!]
\begin{center}
\includegraphics[width=\columnwidth]{meanRadiusJ1.pdf}
\\ \vspace{0.5cm}
\includegraphics[width=\columnwidth]{meanRadiusJ2.pdf}
\end{center}
\caption{The normalised equatorial shadow radius $\bar{b}_{\rm ph}/b_{\rm Schw}$ of a black hole 
in the J spacetime (solid curves). Top panel: as a function of the dimensionless spin $a=J/M^2$ for fixed deformation parameters 
$\alpha_{13} = \alpha_{22} = 1.5$.
Bottom panel: as a function of  $\alpha_{13},  \alpha_{22}$ for $a=0.9$. For comparison, we also show the Kerr equatorial shadow 
(dashed curves). }
\label{fig:Jshadow}
\end{figure}

These results suggest that not all deformation parameters lead to the same degree of deviation from a circular-shaped 
shadow. In the example shown in Fig.~\ref{fig:Jshadow},  the departure from the $a=0$ shadow radius is negligible when
$\alpha_{13} > 0$; in contrast, a positive $\alpha_{22}$ makes the shadow markedly less circular than its Kerr counterpart 
for any spin $a \gtrsim 0.5$. 




\section{Concluding remarks}
\label{sec:conclusions}

The purpose of this paper was to build a bridge between the no-hair theorem's multipolar relation~\eqref{nohair} and the equatorial radius of the 
shadow cast by a Kerr black hole mimicker (this term includes ultracompact objects with a material surface but negligible surface emission as well as 
non-Kerr black holes) within the framework of GR.  

With the help of two stationary-axisymmetric vacuum metrics endowed with rotation and an arbitrary set of multipole moments, we have
been able to produce analytic formulae, in the form of slow rotation expansions, for the shadow radius as an explicit function 
of the moments themselves (or, equivalently, their deviation from the Kerr moments). These formulae have subsequently allowed us to explore the 
variation in the equatorial shadow radius (which serves as a measure of the shadow's degree of circularity or lack thereof) in the multipole moment 
parameter space. 
Our model represents a significant improvement on previous work on the subject~\cite{Johannsen2010ApJ, Broderick2014ApJ} by 
making use of GR spacetimes with rotation several non-Kerr moments beyond quadrupole order and by directly expressing the shadow radius
as a function of these multipoles.

The main result of our analysis is rather clear: a quasi-circular shadow which typically characterises a Kerr black hole could also be the result of 
light propagating in a spacetime with a significant degree of  deviation from the Kerr multipolar structure.  Therefore, addressing the question posed in the title
of this paper, we can say that a black hole shadow \emph{ may not} necessarily be a reliable test of the no-hair theorem. 
This conclusion appears to be at odds with Refs.~\cite{Johannsen2010ApJ, Broderick2014ApJ}, but there is no
real contradiction here because those earlier papers assumed a spacetime with a single non-Kerr multipole moment (the quadrupole). 
This would correspond to a variation along a horizontal line passing through the `Kerr point' in some of the panels of Figs.~\ref{fig:bHTcontour} and
\ref{fig:shadowG_frac}. 
It should be emphasised that our analysis does not imply that a black hole shadow image cannot be used as
a test of the Kerr spacetime. In fact, the results shown in aforementioned figures do allow for an appreciable deviation from a 
Kerr-like shadow if one `moves' along a suitable direction in the multipole moment plane.  The real impact of our results is to
weaken the link between the moments and the shadow shape in the sense that an observation of a Kerr-like shadow 
does not necessarily imply a small deviation from the Kerr moments.



A secondary implication of our results is that a nearly circular shadow does \emph{not} appear to be an exclusive property of separable stationary-axisymmetric 
metrics with a Carter-like geodesic constant of motion. Indeed, the two spacetimes explored in this paper (and in contrast to the Kerr and J metrics) 
are not separable with respect to the Hamilton-Jacobi equation for point particle motion.   

It is conceivable that the analysis presented in this paper could be extended to non-GR theories of gravity, thus enabling 
a connection between the black hole shadow shape and the multipole moments of genuine non-Kerr black holes. 
For example, Refs.~\cite{pappas15a,pappas15b} have developed an Ernst formalism-based framework for the definition and calculation of 
multipole moments of asymptotically flat, stationary-axisymmetric spacetimes in scalar-tensor gravity. The resulting metric could be
used in the same way as the two metrics of the present paper for the calculation of a black hole shadow radius in a slow rotation
approximation. It is unclear, however, if the same method could work equally well for other theories of gravity whose field equations cannot
be reduced to an Ernst potential formalism. 



\acknowledgments

K. G. acknowledges support from research Grant No. PID2020-1149GB-I00 of the Spanish Ministerio de Ciencia e Innovaci{\'o}n.  

