# Approximate solutions for Dorodnitzyn’s gaseous boundary layer limit formula

Carla V. Valencia-Negrete Department of Physics and Mathematics
Universidad Iberoamericana Ciudad de México
Prolongación Paseo de la Reforma 880, Mexico City —01219, MEXICO
<EMAIL_ADDRESS>
###### Abstract

Oleinik’s _no back-flow_ condition ensures the existence and uniqueness of
solutions for the Prandtl equations in a rectangular domain
$R\subset\mathbb{R}^{2}$. It also allowed us to find a limit formula for
Dorodnitzyn’s stationary compressible boundary layer with constant total
energy on a bounded convex domain in the plane $\mathbb{R}^{2}$. Under the
same assumption, we can give an approximate solution $u$ for the limit formula
if $|u|<\\!\\!<\\!\\!<1$ such that:

$u(z)\cong\delta*c*\left[z+\frac{6}{25}\cdot\frac{1}{2i_{0}}\cdot\frac{4U^{2}}{3}z^{4}\right]+o(z^{5}),$

that corresponds to an approximate horizontal velocity component when a small
parameter $\epsilon$ given by the quotient of the maximum height of the domain
divided by its length tends to zero. Here, $c>0$, $\delta$ is the boundary
layer’s height in Dorodnitzyn’s coordinates, $U$ is the _free-stream_ velocity
at the upper boundary of the domain, and $T_{0}$ is the absolute surface
temperature.

_K_ eywords Boundary layer theory $\cdot$ theory Gas dynamics

##  1 Introduction

First, the limit formula is rewritten as a non-linear ordinary differential
equation of order 1 in Lemma 1,

$\frac{\partial u}{\partial
s}=c\left[1-\left(u^{2}/2i_{0}\right)\right]^{-6/25},$

where $i_{0}=c_{p}\ T_{0}$ the constant total energy value, $c_{p}$ is
_specific heat at constant_ _pressure_ , $T_{0}$ is the absolute surface
temperature, $c=\left.\partial u/\partial s\right|_{s=0}>0$ is a strictly
positive constant that represents the _no back-flow_ condition continuous
extension to the lower boundary, and the fractional exponent $-6/25=19/25-1$
comes from the empirical _Power-Law_
$\mu/\mu_{h}=\left(T/T_{h}\right)^{\frac{19}{25}}$ for the absolute
temperature $T$ and the dynamic viscosity $\mu$ with the _free-stream dynamic
viscosity_ $\mu_{h}>0$, and the _free-stream temperature_ $T_{h}>0$ at the
upper boundary [3, p. 46].

Then, it is shown in Lemma 2 that if $|u|<\\!\\!<\\!\\!<1$, the non-linear
term can be expressed as a power series
$\left[1-\left(u^{2}/2i_{0}\right)\right]^{-6/25}=1+(6/25)(1/2i_{0})\
u^{2}+o(u^{4})$. In Theorem 1, we seek a solution of the form of a fourth
degree polynomial, $Az+Bz^{2}+Cz^{3}+Dz^{4}$, in the same way it is done in
[1, 2]. The existence and unicity of the solutions comes from the Darboux sums
limit of the Riemann integral taken to solve Eq. (2).

##  2 Problem Statement

The limit formula (Eq. (1) in the next page) is presented as a non-linear
ordinary differential equation of order $1$.

######  Definition 1.

Let $\tilde{D}\subset\mathbb{R}^{3}$ be a bounded open subset of the three
dimensional Euclidean space such that its projection $\tilde{D}\ \cap\
\mathbb{R}^{2}$ is a planar region $D$,

$D\colon=\left\\{(x,y)\in\mathbb{R}^{2}\ |\ 0<x<L\hskip 2.0pt\&\hskip
4.0pt0<y<h(x)\right\\},$

where $L>\\!\\!>\\!\\!>h>0$, and $h\colon[0,L]\to(0,\infty)\in C^{2}([0,L])$
has derivatives up to order two in $(0,L)$ with extension to the corresponding
extremes, $\\{0,L\\}$, of this interval, and $h(0)=h(L)=\delta$. Moreover, the
upper boundary $h$ has only one critical point $c$ which is a global maximum,
$h(c)=H\geq h(x)$ $\forall x\in[0,L]$. We denote the topological boundary of
$D$ by $\partial D$.

######  Definition 2.

Consider: $(i)$ the first coordinate $\ell$ of Dorodnitzyn’s diffeormorphism
such that for all $\left(x,y\right)\in D$, $\ell\left(x,y\right)=c_{1}\hskip
2.0pt\left[1-\left(U^{2}/2i_{0}\right)\right]\hskip 2.0ptx$, where the _free-
stream_ velocity $u|_{y=h(x)}=U>0$, $c_{1}=p_{0}\hskip 2.0ptT_{0}^{2b/b-1}$,
$p|_{y=0}=p_{0}$ is the pressure at the lower boundary, $T_{0}$ is absolute
surface temperature, and $b=1.405$ is the exponent of the _adiabatic_ and
_polytropic_ atmosphere constant, $p\ V^{b}=cte.$, for the volume
$V=\int\\!\\!\int\\!\\!\int\hat{D}d\mathbf{x}$ [4, p. 35]; $(ii)$ the second
coordinate $s$ of Dorodnitzyn’s diffeomorphism

$\begin{split}s\left(x,\hat{y}\right)&=\int_{0}^{\hat{y}}\rho(x,y)dy,\end{split}$

where $\rho$ is the density; and $(iii)$ the upper boundary in the
Dorodnitzyn’s domain, $\delta(x,h(x))=\int_{0}^{h(x)}\rho(x,y)dy$. Then,
$\mathbf{s}=(\ell,s)$ is a diffeomorphism defined in the domain $D$ [5].

######  Lemma 1.

Let $D$ be a convex domain and $\mathbf{s}:D\to\mathbf{s}(D)$ as described in
Definition 1 and 2. Assume that $(i)$ $u\in C^{2}(D)$, $(ii)$ $\partial
u/\partial s>0$ in $D$ with a continuous extension to $\partial u/\partial
s|_{s=0}=c>0$ to a positive constant $c$, $(iii)$ $u(\ell,0)=0$ for all
$\ell\in[0,1]$, and $(iv)$ the following relation:

$\displaystyle f\hskip 2.0pt\frac{\partial^{2}u}{\partial s^{2}}$
$\displaystyle=$ $\displaystyle\frac{\partial f}{\partial s}\hskip
2.0pt\frac{\partial u}{\partial s},$ (1)

where $f=\left[1-\left(u^{2}/2i_{0}\right)\right]^{-6/25}>0$, and $i_{0}$ is
an strictly positive constant. Then, $u$ is a solution of Eq. (1) if and only
if it verifies Eq. (2):

$\displaystyle\frac{\partial u}{\partial s}$ $\displaystyle=$ $\displaystyle
c\left[1-\left(u^{2}/2i_{0}\right)\right]^{-6/25}.$ (2)

###### Proof.

If $g=\partial u/\partial s$, then Eq. (1) becomes $(\partial g/\partial
s)/g=(\partial f/\partial s)/f$. This is,
$\int_{0}^{y}\frac{\partial}{\partial
s}\left[ln\left(g(s)\right)\right]ds=\int_{0}^{y}\frac{\partial}{\partial
s}\left[ln\left(f(s)\right)\right]ds$. Thus,
$exp\left[ln\left(g(y)/g(0)\right)\right]=exp\left[ln\left(f(y)/f(0)\right)\right]$.
Therefore, $g(y)=\left(g\left(0\right)/f\left(0\right)\right)f(y)$. The _no-
slip_ condition $u(x,0)=0$ for all $x\in[0,L]$ implies that $f(0)=1$. Finally,
we substitute $g=\partial u/\partial s$, and we obtain Eq. (2). ∎

##  3 Analytic solutions for the limit formula

In Lemma 2, we show that there is an expression of the non-linear term
$\left[1-\left(u^{2}/2i_{0}\right)\right]^{-6/25}$ as a convergent power
series if $|u|<\\!\\!<\\!\\!<1$. Finally, following [1, 2], in Theorem 1, we
may seek a solution of the form of a fourth degree polynomial with
coefficients chosen so that the boundary conditions of the problem are
satisfied.

######  Lemma 2.

Let $|u|<\\!\\!<\\!\\!<1$. Then,

$\left[1-\left(u^{2}/2i_{0}\right)\right]^{-6/25}=1+\frac{6}{25}\cdot\frac{u^{2}}{2i_{0}}+\frac{6}{25}\cdot\frac{u^{4}}{2(2i_{0})^{2}}+\frac{1}{2}\cdot\frac{6}{25}\cdot\frac{u^{4}}{(2i_{0})^{2}}+o(u^{6}).$
(3)

###### Proof.

If $\sigma=1-(u^{2}/2i_{0})\cong 1$, then
$\theta=(-6/25)\ln\left(1-(u^{2}/2i_{0})\right)<\\!\\!<\\!\\!<1$. We apply the
_Maclaurin formula_ for the exponential,
$e^{\theta}=1+\theta+(1/2)\theta^{2}+(1/3!)\theta^{3}+\cdot\cdot\cdot$, and
the Taylor series expansion for the logarithmic function if $|\sigma|<1$,
$log(1+\sigma)=\sigma-(1/2)\sigma^{2}+(1/3!)\sigma^{3}+\cdot\cdot\cdot$, so
that:

$\begin{split}\left[1-\left(u^{2}/2i_{0}\right)\right]^{-6/25}&=\exp\left[\frac{-6}{25}\
\ln\left(1-\frac{u^{2}}{2i_{0}}\right)\right],\\\ &=\exp\left[\frac{-6}{25}\
\left\\{\ln\left(1+\frac{u}{\sqrt{2i_{0}}}\right)+\ln\left(1-\frac{u}{\sqrt{2i_{0}}}\right)\right\\}\right],\\\
&=1+\frac{6}{25}\cdot\frac{u^{2}}{2i_{0}}+\frac{6}{25}\cdot\frac{u^{4}}{2(2i_{0})^{2}}+\frac{1}{2}\cdot\frac{6}{25}\cdot\frac{u^{4}}{(2i_{0})^{2}}+\\\
&\quad+\frac{6}{25}\cdot\frac{u^{6}}{3(2i_{0})^{3}}+\frac{6}{25}\cdot\frac{u^{6}}{2(2i_{0})^{3}}+o(u^{8}).\end{split}$

∎

######  Remark 1.

In order to identify the coefficients of a polynomial expression for $u$, we
observe Dorodnityzn’s steps [1, 2]. If we take the upper boundary conditions
for the second order derivative from the problem stated by Dorodnitzyn on the
domain $\mathbf{s}(D)$, then the boundary conditions establish a system of
four equations for the coefficients $A$, $B$, $C$, $D$ so that $u$ can be
expressed as a fourth degree polynomial $Az+Bz^{2}+Cz^{3}+Dz^{4}$ for
$z=s/\delta$ if $A=U\left[2+\left(\lambda/6\right)\right]$,
$B=-U\left(\lambda/2\right)$, $C=U\left[\left(\lambda/2\right)-2\right]$, and
$D=U\left[1-\left(\lambda/6\right)\right]$, where $\lambda$ is the Pohlhausen
coefficient $\lambda=\left[-\delta/\left(1-U^{2}\right)\right]\partial
U/\partial\ell$, and $\lambda=0$ if $U$ is constant.

######  Remark 2.

For each fixed value of $x\in[0,L]$, the height’s normalization
$z(x)=s(x)/\delta(x)\in[0,1]$. Additionally, $\partial u/\partial
z=\delta*\partial u/\partial s$ because $\partial u/\partial z=(\partial
u/\partial s)(\partial s/\partial z)+(\partial
u/\partial\ell)(\partial\ell/\partial z)$ and $\ell$ is independent of $z$.

######  Theorem 1.

Under the same conditions of Lemma 1, Eq. (2) has a unique solution

$u(z)=\delta*c*\left[z+\frac{6}{25}\cdot\frac{1}{2i_{0}}\cdot\frac{4U^{2}}{3}z^{4}\right]+o(z^{5}),$

for the normalized height $z=s/\delta\in[0,1]$.

###### Proof.

Given Eq. (2), we can substitute Eq. (3), so that:

$\begin{split}\frac{\partial u}{\partial s}&=\frac{1}{\delta}\frac{\partial
u}{\partial z},\\\ &=c\left[1-\left(u^{2}/2i_{0}\right)\right]^{-6/25},\\\
&=c\left[1+\frac{6}{25}\cdot\frac{u^{2}}{2i_{0}}+\frac{6}{25}\cdot\frac{u^{4}}{2(2i_{0})^{2}}+\frac{1}{2}\cdot\frac{6}{25}\cdot\frac{u^{4}}{(2i_{0})^{2}}+o(u^{6})\right].\\\
\end{split}$ (4)

Because the _stream-flow velocity_ $U$ is constant, then the Pohlhausen
coefficient $\lambda=0$. Thus, we can seek a solution with a fourth degree
polynomial form $u(z)=2Uz-2Uz^{3}+Uz^{4}$. If we replace this expression in
the Eq. (4):

$\begin{split}u(z)&=\delta*c*\int_{0}^{z}\left[1+\frac{6}{25}\cdot\frac{u(\tau)^{2}}{2i_{0}}+\frac{6}{25}\cdot\frac{u(\tau)^{4}}{2(2i_{0})^{2}}+\frac{1}{2}\cdot\frac{6}{25}\cdot\frac{u(\tau)^{4}}{(2i_{0})^{2}}+o(u^{6})\right]d\tau,\\\
&=\delta*c*\left[z+\frac{6}{25}\cdot\frac{1}{2i_{0}}\cdot\frac{4U^{2}}{3}z^{4}\right]+o(z^{5}).\\\
\end{split}$

∎

##  4 Conclusion

Olga Oleinik’s _no back-flow_ condition is sufficient to show the existence
and unicity of solutions for Dorodnitzyn’s gaseous boundary layer limit when
the horizontal velocity component is small. The fact that one can find
analytic approximations for in a compressible boundary layer when the
horizonal component of the velocity is small allows us not only to empirically
see that the velocity profile is stable at low velocities, but also to
analytically understand that, in this case, one can arrive at flow expressions
that linearly depend on the height.

##  5 Appendix

Funding: This research did not receive any specific grant from funding
agencies in the public, commercial, or not-for-profit sectors.

Declarations of interest: none

## References

  * [1] Dorodnitzyn, A. A.  Laminar Boundary Layer in Compressible Fluid.  In Comptes Rendus de l’Académie des Sciences de l’URSS (Doklady), 34:8, pages 213–219, 1942.
  * [2] Lees, L.  Stability of the Supersonic Laminar Boundary Layer With a Pressure Gradient.  In Princeton University NJ James Forrestal Research Center, 1950-11-20.
  * [3] Smits, A. J. Dussauge, J.-P.  Turbulent Shear Layers in Supersonic Flow. Springer-Verlag, 2006.
  * [4] Tietjens, O. G.  Fundamentals of Hydro- and Aerodynamics. Based on Lectures by L. Prandtl, Ph. D. Dover Publications Inc., 1934.
  * [5] Valencia-Negrete, C. V., et al.  Reynolds’ Limit Formula for Dorodnitzyn’s Atmospheric Boundary Layer In Convective Conditions.  In International Journal of Applied Mathematics, Vol. 31, No 4 (2018), pp. 673–695, DOI: 10.12732/ijam.v31i4.12


# Constraints on the Inner Regions of Lensing Galaxies from Central Images
using a Recent AGN Offset Distribution

D. Perera1, L. L. R. Williams1, C. Scarlata1
1School of Physics and Astronomy, University of Minnesota, Minneapolis, MN,
55455, USA.
E-mail<EMAIL_ADDRESS>
(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

In gravitational lensing, central images in quads can serve as a powerful
probe of the inner regions of lens galaxies. The presence of an offset central
supermassive black hole (SMBH) has the potential to distort the time-delay
surface in a way such that 3 central images form: a strongly demagnified image
near the SMBH, and two less demagnified (and potentially observable) images at
a central maximum and saddle point. Using a quad-lens macro-model, we simulate
the constraints that could be placed on various lens galaxy parameters based
on their central images’ probability of detection or non-detection. Informed
by a recent low-redshift distribution of off-nucleus AGN, we utilize Bayesian
inference to constrain the mean SMBH off-nucleus distance and galactic core
radius for a sample of 6 quads. In general, we find that a detection of the
central image in any quad would favor larger SMBH off-nucleus distances and
galaxy core sizes. Assuming a linear relationship between core radii and
velocity dispersion $r_{c}=b\sigma$, these results similarly imply strong
constraints on $b$, where the likely case of a central image non-detection in
each quad constraining $b$ to $3.11^{+2.72}_{-2.26}\times 10^{-4}$ kpc km-1 s.
Our results show that tight constraints on lens galaxy parameters can be made
regardless of a detection or non-detection of a central image. Therefore, we
recommend observational searches for the central image, possibly using our
suggested novel detection technique in UV filters, to formalize stronger
constraints on lens galaxy parameters.

###### keywords:

gravitational lensing: strong – galaxies: general – quasars: supermassive
black holes

††pubyear: 2022††pagerange: Constraints on the Inner Regions of Lensing
Galaxies from Central Images using a Recent AGN Offset Distribution–A

##  1 Introduction

Gravitational lensing theory predicts that the number of multiple images
formed by non-singular mass distributions must always be odd. Following
Fermat’s Principle, the locations of these images correspond to stationary
points on the time-delay surface. Three image systems (known as "doubles"),
form visible images at a minimum and saddle point, where the saddle point
image is usually demagnified relative to its minimum counterpart. Five image
systems (known as "quads") form 2 images at minima and 2 images at saddle
points. In both cases, a demagnified central image is formed at a maximum near
the center of the lens. The central image is usually demagnified beyond
visibility due to steep central lens density profiles causing sharp peaked
behavior of the time-delay surface. Almost all existing searches for the
central image rely on the optical, or radio wavelengths. In the optical, the
central image is drowned by the light of the lensing galaxy, while detections
in the radio are limited because most QSO sources are radio quiet.

Out of about 200 known doubles, only 2 have observed central images: PMN
J1632-0033, demagnified to 0.004 (by 6 mag.), compared to the brightest image
(Winn et al., 2004), and PKS 1830-211, demagnified to 0.007 (by 5.4 mag.
Muller et al., 2020). Of the $\sim$50 known quads lensed by an isolated
galaxy, no detections exist in the optical, radio, or at ALMA wavelengths
(Wong et al., 2015; Tamura et al., 2015; Wong et al., 2017b). The more
recently discovered systems from  Gaia (Lemon et al., 2022) do not appear to
contain central images, but that is yet to be confirmed with further
observations.

Even though the central image is seldom detected, its study is of great
interest because it can serve as a useful probe of the central regions of lens
galaxies and Supermassive Black Hole (SMBH) properties. Most, and probably all
galaxies have SMBHs at their centers (Kormendy & Richstone, 1995; Ferrarese &
Merritt, 2000). The black holes are tied to two important aspects of galaxies:
SMBH growth is closely linked to galaxy formation, including galaxy mergers
(Di Matteo et al., 2005; Koss et al., 2018) and formation of density cores
(Nasim et al., 2021), and inspiralling and mergers of binary SMBH (Begelman et
al., 1980) leading to the emission of gravitational waves.

An SMBH at kpc distances from the galaxy center will slowly spiral inwards
through dynamical friction from collisionless particles, namely stars and dark
matter, as well as gas (Chen et al., 2022). Observationally determined
distribution of galaxy host-SMBH separations will constrain these dynamical
friction timescales, and test hydrodynamical cosmological simulations
(Volonteri et al., 2020; Katz et al., 2020). In massive elliptical galaxies,
this process will carve out a density core. To get a better physical
understanding of the dynamical friction within $\sim$1 kpc of the galaxy
center it is crucial to constrain the sizes of galaxies’ density cores. Nearby
($\lesssim$ 100 Mpc) massive galaxies with large ($\sim$700 pc) density cores
can be resolved (e.g. Rantala et al., 2018; Thomas et al., 2016), but this is
more difficult for more distant galaxies and smaller cores. Since the central
lensed image is affected by the SMBH and the galaxy mass density near the
center (Mao et al., 2001; Rusin et al., 2005; Mao & Witt, 2012), a detection
or upper limit on its brightness can place constraints on SMBH mass, distance
from center, and lens galaxy core size.

Figure 1: Summary of our statistical analysis plan. We consider two
assumptions for inference: (1) $r_{c}\not\propto\sigma$ and (2)
$r_{c}\propto\sigma$. In the first case (1), we constrain the SMBH offset and
galaxy core radius for individual quads depending on a hypothetical central
image non-detection (blue and Figure 6) or detection (yellow and Figure 7).
The cases in the red box, and Figures 10 and 11 make an assumption that all
galaxy lenses have the same SMBH offset and core radius, which is unlikely to
be true. Therefore these estimates should be treated only as "back-of-
envelope" for this general type of galaxies. In the second case (2), we only
make an inference on the $r_{c}-\sigma$ proportionality constant $b$, which is
assumed to be the same for all galaxies. The constraints on $b$ from
individual quads are given for each detection scenario (green and Figure 8).
Combining these results for a global constraint on $b$ is similarly given for
each detection scenario (purple and Figure 9).

Using lens system CLASS B1030+074, where no central image was detected, Quinn
et al. (2016) finds that $\sim$45% of galaxies should yield observable central
images, assuming demagnifications $\leq$10 magnitudes relative to the
brightest image. Hezaveh et al. (2015) propose that 10-hr ALMA observations in
the $mm$ band can detect a central image at high significance for a lens
galaxy core size $\geq$0.2 kpc, allowing for strong constraints on central
density slope, core size, and mass of the central SMBH. Recent ALMA
observations have yielded upper limits on the central image flux and SMBH mass
(Wong et al., 2015; Tamura et al., 2015; Wong et al., 2017b).

If measured or constrained through observations, properties of SMBH, such as
their mass and distance from the host galaxy center will provide invaluable
clues for the understanding of the central regions of galaxies, and SMBH
merger rates. The masses of SMBH in many nearby galaxies have been measured
using a range of other methods (Dullo et al., 2021; Gültekin et al., 2009).
However, their distances from host galaxy centers are less well known (Skipper
& Browne, 2018).

A recent examination of $z=0.3-0.8$ Active Galactic Nuclei (AGN) led to a
determination of the characteristic probability density function of their
offsets’ upper limits from the host galaxy center (Shen et al., 2019). Their
results and sample are given the name of VODKA, which we accordingly adopt.
The employed methodology, known as "varstrometry", is described lucidly in
Hwang et al. (2020). (This technique can be applied to single off-nucleus AGN,
or generalized to dual AGN.)

Observationally, a sub-kpc off-nucleus AGN and the center of its host galaxy
appear as single-source photocenter (photometric center) for $z>0.5$ (Hwang et
al., 2020). The vast majority of AGN exhibit aperiodic photometric variability
on day-year timescales of $\gtrsim 0.03$ mag (Sesar et al., 2007), with
$\sim$30% varying $\geq$0.1 mag (Sesar et al., 2007; Rengstorf et al., 2006).
For a single off-nucleus variable AGN and constant-flux host galaxy, it is
expected that the AGN variability will lead to astrometric variability of the
photocenter of the AGN-galaxy system, which will be strongly correlated with
the total detected flux (Hwang et al., 2020). This allows for a measurement of
the distance separation between the AGN and host galaxy center with linear
regression (see eq. 4 of Hwang et al., 2020). When applied to low redshift
($0.3<z<0.8$) AGN and host galaxy pairs observed with Gaia DR2, it is found
that there are strong constraints on AGN separation. Nearly all AGN are at $<$
1 kpc, 90% are at $<$ 500 pc, and 40% are at $<$ 100 pc (Shen et al., 2019).
While this result is a significant improvement over previous determinations,
it is important to realize that the distribution of separations of AGN may
differ from that of SMBH. All SMBH, not just AGN affect the central regions of
galaxies by reshaping their central mass distributions, and also lead to the
formation of gravitational wave emitting inspirals.

Offset SMBH are of particular importance for studying and potentially
detecting the central lensed image. Ideally, one wants the distribution of
host galaxy-SMBH separations, but such information does not exist. Instead, we
take the AGN-galaxy host separation distribution from VODKA and use it as a
prior in our analysis. The present paper is the first to incorporate the
distribution of AGN offsets.

Gravitational lensing theory predicts that a sufficiently offset SMBH can
produce extra central images. This is displayed in Figure 2. The inclusion of
an offset SMBH creates two additional stationary points (and thus two new
images) in the time-delay surface: a steep maximum very near to the location
of the SMBH and a saddle between the original central maximum and the offset
SMBH. The formation of the two new images and their properties depend on the
offset distance of the central SMBH, the azimuthal location of the SMBH in the
case of non-circular lenses, and the lens galaxy core size. If the offset
distance is too short, then only one image forms at the SMBH maximum. The SMBH
(maximum) image is always strongly demagnified, regardless of whether the two
other central images are formed. When additional central images form in the
case of a large offset distance, the central image near the maximum of the
lens density profile, and the saddle image are not as demagnified, and can
potentially be bright enough to be observed. It is important to emphasize that
if the SMBH is not offset sufficiently, then no extra central images will
form.

Here we present a new modelling framework that will allow constraints to be
made on various galaxy properties with future observations. To make
predictions about these properties, we proceed in three sequential steps. In
the first step we create the macro-model of the galaxy using the positions and
time delays of the four images of its quad. In the second step, we simulate
central images for each of the macro-models. We treat the SMBH as a point mass
and sample its offset according to the aforementioned VODKA distribution.
Lastly, we use Bayesian inference to simulate constraints on galaxy parameters
for specific observation scenarios.

Our analysis is the first one to study several (6) quad lens systems. We
restrict our sample to quads because they provide more constraints for the
galaxy macro-model, which affects the central region through its ellipticity
and the location of the source. We assume two main possible scenarios for each
quad: (i) a non-detection of central images (ii) detection of at least one of
the central images. This allows us to place statistical constraints on the
SMBH offset ($\Delta r$) and core radius ($r_{c}$) of individual lenses.

Additionally, we want to constrain the global distribution of these parameters
for our galaxies. However, the distributions of core radii and SMBH distances
may depend on the galaxy’s mass, or velocity dispersion $\sigma$ (and probably
other parameters), and will not be the same for all galaxies in our sample. If
one wants a ballpark value of a parameter of typical lens galaxy, one can
ignore these differences, and estimate the parameter by combining the results
for individual quads. We also consider a separate scenario of combination of
detection in some, and non-detection in other quads in our sample. This
analysis assumes $r_{c}\not\propto\sigma$. Alternatively, if we assume that
$r_{c}=b\sigma$, then $b$ is the same for all galaxies, and can be estimated
using all quads combined. This analysis plan is outlined in Figure 1.

Importantly, each scenario yields independent statistical constraints on the
SMBH offset distribution and central density profile based on prior
distributions informed from the VODKA analysis and our lens macro-model
sample. Assuming that lenses from ongoing and future surveys will be analogous
to the sample we use, our analysis can be used to forecast how well lens
galaxy properties can be constrained with future data.

In Section 2 we describe our sample of lensed QSOs; in Section 3 we discuss
our lens macro-modelling, simulation of the central images, and statistical
inference analysis; in Section 4 we present the results of our analysis; in
Section 5 we describe a novel technique to detect central images; and in
Section 6 we discuss the implications of our constraints.

Figure 2: Visualization of the Fermat lens potential distortion by the SMBH.
In both panels, the view is restricted to the central 256 mas2 of the lens
galaxy such that the 4 quad images are out of view. The left panel shows the
resulting Fermat potential contours when a $3\times 10^{8}M_{\odot}$ SMBH is
introduced at $(\Delta{\rm RA},\Delta{\rm Dec})=(6.0,-6.0)$ mas. In this case
3 images form (for a total of 7 in the system): a central image at the central
maximum, a saddle image between the central maximum and the SMBH, and a SMBH
image very close to the SMBH maximum. The spacing of contours shows the
steepness of the potential increase due to SMBH, indicating that the SMBH
image is strongly demagnified. The right panel shows the lens potential
contours when the same SMBH is offset a shorter distance at $(\Delta{\rm
RA},\Delta{\rm Dec})=(3.6,-3.6)$ mas. The same 3 images form, with the saddle
images forming closer to the midpoint of the central and SMBH images due to
the less offset SMBH. The images in this case (right panel) are all magnified
relative to the images in the first case (left panel), illustrating the
dependence of central image magnification $\mu$ on SMBH offset $\Delta r$.

##  2 Data

Table 1: Sample of gravitationally lensed QSOs

QSO | zQSO | zlens |  $\sigma$ [km s-1] | MSMBH [M⊙] | mo [AB], Filter | Flux [$10^{-17}f_{\lambda}$] | Time Delays | References
---|---|---|---|---|---|---|---|---
HE0435-1223 | 1.693 | 0.454 | 222 $\pm$ 15 | 3.7 $\pm$ 1.2 $\times 10^{8}$ | 20.33, F275W | 10.36 $\pm$ 0.07 | $\Delta\tau_{12},\Delta\tau_{13},\Delta\tau_{14}$ | 1, 8
PG1115+080 | 1.735 | 0.311 | 287 $\pm$ 18 | 1.3 $\pm$ 0.4 $\times 10^{9}$ | 18.60, F218W | 101.30 $\pm$ 0.50 | $\Delta\tau_{14},\Delta\tau_{\\{23\\}4},\Delta\tau_{1\\{23\\}}$ | 2, 8, 10
RXJ1131-1231 | 0.657 | 0.295 | 323 $\pm$ 20 | 2.3 $\pm$ 0.7 $\times 10^{9}$ | 19.62, F218W | 32.96 $\pm$ 0.29 | $\Delta\tau_{14},\Delta\tau_{24},\Delta\tau_{23}$ | 3, 8
SDSSJ0924+0219 | 1.523 | 0.393 | 215 $\pm$ 21.5 | 3.2 $\pm$ 1.6 $\times 10^{8}$ | 19.41, F275W | 25.40 $\pm$ 0.12 | — | 4
WFI2033-4723 | 1.662 | 0.658 | 250 $\pm$ ${}^{15}_{21}$ | 6.7 $\pm$ 2.3 $\times 10^{8}$ | 18.56, F467M | 18.27 $\pm$ 0.23 | $\Delta\tau_{1\\{23\\}},\Delta\tau_{14},\Delta\tau_{\\{23\\}4}$ | 5, 8, 11
MG0414+0534 | 2.640 | 0.958 | N/A | 8.0 $\times 10^{8\pm 1}$ | 23.58, F621M | 0.11 $\pm$ 0.01 | — | 6, See Section 2
RXJ0911+0551 | 2.763 | 0.769 | N/A | 8.0 $\times 10^{8\pm 1}$ | 19.45, F547M | 6.11 $\pm$ 0.10 | $\Delta\tau_{1\\{234\\}}$ | 7, 9, See Section 2

The columns list QSO identifier, source and lens redshifts, lens velocity
dispersion, SMBH mass inferred from the $M-\sigma$ relation, estimated
apparent magnitude mo for the brightest quad image as observed in the
corresponding HST filter, specific flux of the brightest quad image in the
same HST filter ($f_{\lambda}$ in units of ergs s-1 cm-2 Å-1), the time delays
used, and references from which the data were taken. The subscripts on the
time delays $\Delta\tau$ indicate the arrival order of the two images. Curly
brackets enclose images assumed to have the same arrive time. Two systems did
not have published image arrival times. References: (1) Courbin et al. (2011);
Wong et al. (2017a), (2) Tonry (1998); Treu & Koopmans (2002), (3) Sluse et
al. (2007); Suyu et al. (2013), (4) MacLeod et al. (2015), (5) Sluse et al.
(2019), (6) Hewitt et al. (1992), (7) Bade et al. (1997), (8) Millon et al.
(2020), (9) Hjorth et al. (2002), (10) Bonvin et al. (2018), (11) Bonvin et
al. (2019)

Table 1 presents our sample of gravitationally lensed QSOs. This sample was
chosen to roughly span the redshift range of the VODKA population,
$0.3\\!<\\!z\\!<\\!0.8$ (Shen et al., 2019), allowing direct comparison of our
results with their AGN offset distribution. Most of our systems also have time
delay information, which aids in lens macro-modeling. The SMBH masses we use
are found by applying the $M-\sigma$ relation from Dullo et al. (2021), as
explained in further detail in Section 3.2. Four of the lens galaxies have
measured central velocity dispersions. For SDSSJ0924+0219, MacLeod et al.
(2015) does not directly measure $\sigma$, and instead estimates it with lens
modelling. We assume 10% uncertainty on their result, and record that value in
Table 1. MG0414+0534 and RXJ0911+0551 have no measured $\sigma$ in the
literature. Therefore, we assume a $\sigma$ for these two systems equal to the
average of the measured $\sigma$ in the rest of the sample. From this we
obtain their $M_{\rm SMBH}\approx 8.0\times 10^{8}M_{\odot}$ from the same
$M-\sigma$ relation. For the uncertainty on this $M_{\rm SMBH}$, we assume a
factor of 10, so as to span a wide enough range for possible $M_{\rm SMBH}$.

RXJ0911+0551 presents itself as an exception in our sample. The radial
distribution of the 4 observed quad images is highly asymmetric, probably due
to the presence of a nearby galaxy cluster in the field, providing external
shear (Burud et al., 1998; Tortora et al., 2004). As a result, the core size
distribution of our lens macro-models skews to larger values beyond the range
of the other 6 lenses in our sample. Therefore, we treat this lens separately
from the rest of the sample of 6 and derive independent constraints from it
(see Table 2 and Section A). Henceforth, we refer to "our sample" as those
QSOs mentioned in Table 1 not including RXJ0911+0551.

##  3 Analysis

The goal of the analysis is to obtain constraints on the galaxy core radii and
SMBH offset using observational constraints on the central image.

Our analysis can be summarized in three sequential distinct steps: (1) We
constrain galaxy macro-models based on the 4 non-central images. The result of
this step are $1000$ macro-models per lens, each described by 16 parameters,
including the QSO source $(x,y)$ position. The lens galaxy core radius is
derived from the 14 lens parameters. These macro-models form a prior for the
next step of the analysis. (2) For each macro-model, we simulate the central,
saddle, and SMBH image locations and magnifications using a large sample of
SMBH offsets from the VODKA distribution and SMBH masses from the $M-\sigma$
relation. Image properties depend primarily on the SMBH masses ($M_{\rm
SMBH}$), offsets ($\Delta r$), and galaxy’s core radius ($r_{c}$). (3) The
final step is the statistical inference for the two main observational
scenarios we consider: (i) non-detections of the central image, and (ii)
detection of at least one central image. Additionally, when we combine these
constraints for all galaxies, we also consider a third case of some detections
and some non-detections. With this framework, we can write the posterior
probability as:

$\begin{split}P(\Delta r,r_{c}|D)\propto\int P(D|M_{\rm SMBH},\Delta
r,r_{c})\,P(\Delta r)\,P(r_{c})\\\ \times\,P(M_{\rm SMBH})\,dM_{\rm
SMBH}\end{split}$ (1)

where $D$ are the input image positions and time delays (where available; see
Table 1).

As we explain later in this section, the 16 parameter priors going into the
macro-modeling in step (1) are flat. The resulting distributions of these
parameters emerging from this step are not flat anymore, as the galaxy
properties have been constrained by the quad images. These distributions
become the priors for step (2), from which we obtain the prior for the galaxy
core radius $P(r_{c})$. The distributions of galaxy parameters ($\Delta r$ and
$r_{c}$) get further constrained as a result of step (3), and the constraints
are different for non-detections vs. detections. The priors on SMBH properties
going into step (2) are described in more detail in Section 3.2. Step (3)
further constrains these properties, which we present as the main result of
the paper.

###  3.1 Constraining Galaxy Macro-Model

We generate galaxy-scale macro-models based on the 7 lens systems we consider
in this paper. Since the central images are currently not observed in any of
these, and we envision that our analysis can be extended to future lens
systems, our modelling need not be tailored exactly to these systems. Instead,
we use these 7 as approximate examples of realistic lenses. Because of that,
the many macro-models we generate per lens fit the observed lensing data
(image positions and time delays, where available) only to $\chi^{2}\leq 9$,
and in the systems where satellite galaxies are detected in addition to the
main lens galaxy, their positions are not fixed at the observed position.
Loosening the criteria to accept macro-models with $\chi^{2}\leq 9$ results in
lens plane image rms between $<0.005"$ and $\sim 0.035"$.

We represent lens galaxies by a superposition of two softened power-law
ellipsoid potentials, called alphapot (Keeton, 2011):

$\Psi_{\rm
gal}=b\left(s^{2}+x^{2}+\frac{y^{2}}{q^{2}}+K^{2}xy\right)^{\frac{\alpha}{2}}$
(2)

where $b$ is the normalization, $\alpha$ is the power-law exponent111Not to be
confused with the deflection angle., $s$ is the core radius, $K$ and $q$
determine the ellipticity and the position angle of the ellipsoid (Ghosh et
al., 2020; Barrera et al., 2021). In addition to being a reasonably good
representation of elliptical galaxies, it allows for analytical calculations
of the deflection angles, normalized projected surface density $\kappa$, and
shear $\gamma$ from the first and second derivatives of $\Psi_{\rm gal}$. We
do not include the SMBH at this point because it does not affect the QSO image
positions, and therefore does not affect the lens galaxy macro-model.

The two mass components were allowed to have a non-zero offset between their
centers, $(x_{\rm off},y_{\rm off})$. The reason for the offset is that a
single component, or two co-centered components sometimes cannot reproduce QSO
images to astrometric precision or other observables of the system (Rusu et
al., 2020). Offsets break the elliptical symmetry of the lens and result in
lopsided galaxy mass distributions, which apparently help to model some lenses
and populations of lenses (Bruderer et al., 2016; Gomer & Williams, 2018;
Nightingale et al., 2019; Williams & Zegeye, 2020; Barrera et al., 2021).

The two alphapots have a total of 10 parameters, i.e., two sets of $b$,
$\alpha$, $s$, $q$, and $K$. Combined with $(x_{\rm off},y_{\rm off})$,
external shear amplitude and direction, and the QSO source position the total
number of macro-model parameters is 16. We use downhill simplex to find
solutions. The starting ranges for the macro-model parameters are the same for
all 7 lenses. The density slopes $\alpha=1.0\pm 0.2$, core radii $s=200\pm
199$pc, external shear amplitude, $\gamma=0.1\pm 0.1$, and offsets of the
secondary mass component, $x_{\rm off}=y_{\rm off}=0\pm 40$pc. We made
exceptions for 3 systems that have visible nearby satellite galaxies: WFI
2033, MG 0414, and RXJ 0911, where the second mass components were given
larger initial offsets to represent the satellite. For RXJ 0911+0551 we used a
larger initial external shear, $\gamma=0.1+0.1\pm 0.1$, to account for the
nearby galaxy cluster.

Downhill simplex search is free to modify these initial values, so many of the
final values of all parameters were outside of the starting ranges. This is
also true of the location of the satellite galaxy; we did not fix it at the
observed position but allowed simplex to find different solutions for each
run. Because the number of macro-model parameters exceeds the data
constraints, we generate many models for each observed quad. From these we
reject macro-models that have two density peaks, and those where the single
density peak is not coincident with the light peak, i.e., the center of the
main lens. We also restrict the ellipticity of the final macro-models to
disallow very elliptical or unphysically shaped galaxies: the ellipticity of
the lens potential of each of the two mass components were restricted to have
axis ratios $\geq 0.7$. The surviving macro-models—about 1000 per lens
system—sample the model space allowed by our assumptions and lensing
degeneracies.

For each lens, the fitted macro-models have a range of central density profile
slopes, and hence core radii. In general, steeper central density profiles
imply smaller core radii. For each lens macro-model we calculate the galaxy
density core size as the radius where the log-log density becomes steeper than
$-0.425$. While this value is somewhat arbitrary, it works well to estimate
core sizes of mass distributions that do not have a core radius explicitly
incorporated in their analytical form, as is the case with our two-component
macro-models. Since we generate a range of macro-models for each observed
lens, we also have a range of core sizes. Six of the 7 systems have core radii
distributions that peak at $\lesssim 100$ pc. RXJ 0911 is an exception: its
distribution is broad with a peak at $\sim 400$ pc. These value are consistent
with the range determined based on local ellipticals, 50-500 pc (Ferrarese et
al., 2006; Hezaveh et al., 2015).

The 16 parameters for every one of the 1000 macro-models per system are passed
to the next step of our analysis. These parameters do not have flat
distributions, as they have been constrained by the quad images. Galaxy core
size is calculated based on these parameters; its distribution is shown as
yellow curves in the Figures 10 (for the 6 in our sample) and 15 (for RXJ0911)
of this paper.

###  3.2 Simulating Central Images

With the galaxy macro-models generated, the next focus is to find the
locations and magnifications of the corresponding central images. Since we are
most concerned with the central region, we re-scale the galaxy window to
0.3"$\times$0.3" about the center of the lens. We generate lensing potentials
$\Psi_{\rm gal}$ from the galaxy macro-models parameters (see Section 3.1) and
SMBH parameters.

For the SMBH, we assume it to be a point mass with lens potential:

$\Psi_{\rm SMBH}=\theta_{E}^{2}\ln{\sqrt{x^{2}+y^{2}}},$ (3)

where $\theta_{E}$ is the Einstein radius of the SMBH:

$\displaystyle\theta_{E}$ $\displaystyle=\sqrt{\frac{4GM_{\rm
SMBH}}{c^{2}}\frac{D_{ds}}{D_{s}D_{d}}}.$ (4)

As mentioned earlier, the perturbation of the outer 4 images in each system is
negligible since the SMBH only adds $\sim$10-3 times the mass within the
Einstein radius to the galaxy. From this setup we simply sum the two lens
potential components to get the total lensing potential of a particular
system: $\Psi_{\rm tot}=\Psi_{\rm gal}+\Psi_{\rm SMBH}$. We then scan the re-
scaled window to find the locations and magnifications of the images.

Within this framework, we displace the SMBH from the galaxy center according
to the results from VODKA. To do this we generate a large sample of SMBH
offset positions according to the probability distribution obtained by VODKA.
From the right panel of Figure 1 of Shen et al. (2019), we transform the
presented CDF into a PDF and fit to a Gaussian profile to obtain the best fit
average and standard deviation of their sample. We find $0.131\pm 0.008$ kpc,
and $0.163\pm 0.007$ kpc for the average and standard deviation, respectively.
In practice, this represents a truncated Gaussian since the distribution
extends into the regime of negative SMBH offsets. Using a Box-Muller
transform, we use these results to generate a distribution of SMBH-galaxy
center offsets. The azimuthal locations of SMBH with respect to the galaxy
center were picked randomly. It is important to note that VODKA provides the
distribution of SMBH offset upper limits, and we assume that the offset
positions are the same as these upper limits. This is supported by recent
simulation results which find that 60% of SMBH in brightest cluster galaxies
at $z=0$ are offset by $>0.1\,$kpc, and at $z=2$ about 80% are offset by
$>1\,$kpc (Chu et al., 2022). While the lens galaxies in our sample are not
BCGs, their stellar masses are comparable to those in that study.

Using the $M-\sigma$ relation for Sersic plus Core-Sersic galaxies from Dullo
et al. (2021), and the velocity dispersion, $\sigma$ (Table 1) for each QSO
lens galaxy, we obtain a SMBH mass range for each lens. This range is
determined by propagating the measured uncertainty in $\sigma$ in the
$M-\sigma$ relation. For each macro-model generated, we pick the SMBH mass
randomly from within this range.

Finally, for each of the 1000 macro-models per system we generate images based
on the following parameters described above: SMBH mass, SMBH offset distance,
SMBH azimuthal position, and lens galaxy macro-model (each with a
corresponding galaxy core size). In general, the central image is always
produced and is always demagnified. However, given certain SMBH masses,
offsets, and macro-models, two additional images can be produced: a similarly
demagnified saddle image, and a more strongly demagnified image very near the
SMBH (see Figure 2). The saddle image always forms between the central and
SMBH image. Using the results of magnification and SMBH offset for each macro-
model, we can place statistical constraints on the galaxy lens core size and
SMBH offset.

###  3.3 Statistical Inference

In this section we model the analysis one would do in case of the likely
scenario of a non-detection of the central image. Additionally, we also
consider a less likely but more interesting case of a central image detection.
Deep observations, whether or not they result in a central image detection,
can constrain the galaxy core size and SMBH offset.

To place constraints on the galaxy core size $r_{c}$, and SMBH offset $\Delta
r$, we employ Bayesian inference based on the results from our analysis
described in Sections 3.1 and 3.2. Our first analysis assumes that $r_{c}$ is
independent of the galaxy’s measured velocity dispersion $\sigma$. Later, when
we combine all the quads in our sample for a general constraint, we assume
that $r_{c}$ is proportional to $\sigma$, which is the simplest, yet
physically plausible, relation that these two parameters can have. For each of
these two analyses we consider two cases to inform our likelihood function:
non-detection of all central images and detection of at least one of the
central images, for example, using the technique outlined in Section 5. In the
latter analysis, we consider a case where 2 systems have detections, and 4
have non-detections.

Focusing first on results for individual QSO quads with no $r_{c}-\sigma$
relation, we assume truncated Gaussian prior probabilities on SMBH offset
$\Delta r$ according to the VODKA distribution (Shen et al., 2019), and prior
probability distribution of the galaxy core sizes, $r_{c}$, from the galaxy
lens macro-models, which resemble Gaussians, described in Section 3.1. For the
SMBH mass $M_{\rm SMBH}$ we assume a uniform prior $P(M_{\rm SMBH})$ within
the defined $M_{\rm SMBH}$ range described in Section 3.2. With this, we can
solve for the 2D posterior probability using equation 1.

To get a posterior on $r_{c}$ for each quad, we can marginalize $P(\Delta
r,r_{c}|D)$ over $\Delta r$:

$\begin{split}P(r_{c}|D)\propto\int_{\Delta r}P(\Delta r,r_{c}|D)d\Delta
r\end{split}$ (5)

We can repeat vice versa for a posterior on $\Delta r$. This intrinsically
assumes that there is no relation between $r_{c}$ and $\sigma$. We apply this
analysis individually to each quad to estimate their $\Delta r$ and $r_{c}$.

Alternatively, we can expand the above method to account for a potential
correlation between velocity dispersion and core size. As stated above, we
assume that $r_{c}$ and $\sigma$ are proportional, therefore:

$r_{c}=b\sigma$ (6)

where $b$ is a fit parameter we wish to extract. Under this model, the
parameter $b$ is assumed to be the same for all quads, and we do not attempt
to constrain $\Delta r$. From this, we assume priors for $b$, $P(b)$, with the
same shape as $P(r_{c})$. We continue to use the same priors for $M_{\rm
SMBH}$ and $\Delta r$ as before. Therefore, we can write a separate 2D
posterior:

$\begin{split}P(\Delta r,b|D)\propto\int P(D|M_{\rm SMBH},\Delta
r,b)\,P(\Delta r)\,P(b)\\\ \times\,P(M_{\rm SMBH})\,dM_{\rm SMBH}\end{split}$
(7)

Repeating the same marginalization as in equation 5 will give us individual
posteriors for each lens for $b$. To get a global constraint on $b$, we can
multiply all the resulting marginal posteriors for each quad together:

$\begin{split}P(b|D)\propto\prod_{\rm quads}\left(\int_{\Delta r}P(\Delta
r,b|D)d\Delta r\right)\end{split}$ (8)

This allows us to constrain the $r_{c}-\sigma$ relation.

We emphasize that the constraints for $\Delta r$ and $r_{c}$ are results
applicable to individual lenses, thus providing a robust method to derive
these properties for observed quads. For $b$, however, the constraint is
intended to be a single result for all lenses assuming that equation 6 is
obeyed.

Combining the distributions for $\Delta r$ and $r_{c}$ presents a separate
challenge. As a ballpark estimate, we can use the simplistic assumption that
$\Delta r$ and $r_{c}$ are the same for all lenses. With this, we can derive
these constraints in the following way:

$\begin{split}P(\Delta r|D)\propto\int_{r_{c}}\left(\prod_{\rm quads}P(\Delta
r,r_{c}|D)\right)dr_{c}\end{split}$ (9)

and:

$\begin{split}P(r_{c}|D)\propto\int_{\Delta r}\left(\prod_{\rm quads}P(\Delta
r,r_{c}|D)\right)d\Delta r\end{split}$ (10)

For the scenarios of a central image non-detection and detection, we follow
the outlined inference procedure and simply define different likelihood
functions $P(D|M_{\rm SMBH},\Delta r,r_{c})$ for each case, as explained
below.

####  3.3.1 Central Image Non-Detection

For the likely case of a non-detection, we define a limiting specific flux
$f_{\rm crit}$ below which we assume any central image regardless of its
position in the lens plane will not be detected. While this flux limit can be
varied at one’s wish, we set it to be equivalent to 10 magnitudes fainter than
the brightest quad image in each lens system, corresponding to a
magnification222Magnifications can be converted into magnitudes using:
$\displaystyle m_{i}$ $\displaystyle=-2.5\log_{10}(\mu_{i})+m_{\rm BI}$ where
$m_{i}$, $\mu_{i}$, and $m_{\rm BI}$ are the image magnitude, image
magnification, and brightest quad image magnitude, respectively. $\mu_{\rm
crit}$ of 10-4 relative to the brightest quad image (same condition as used in
Quinn et al., 2016). Since, by definition, magnification $\mu=f/f_{\rm BI}$
where $f_{\rm BI}$ is the flux of the brightest quad image (see Table 1), the
limiting flux $f_{\rm crit}=\mu_{\rm crit}f_{\rm BI}$. Under this assumption,
we define the likelihood of non-detection $P(D|M_{\rm SMBH},\Delta r,r_{c})$
as the integral of a Gaussian probability density centered around each
predicted image flux ($f_{i}$) up to $f_{\rm crit}$:

$P(D|M_{\rm SMBH},\Delta r,r_{c})\propto\int_{-\infty}^{f_{\rm
crit}}\exp(-\frac{1}{2}\left(\frac{f-f_{i}}{\sigma_{f}}\right)^{2})df$ (11)

The dependence of the likelihood on $M_{\rm SMBH}$ enters through $f_{i}$. We
define the image flux uncertainty $\sigma_{f}$ to be $0.4f_{i}$ to allow for a
range of uncertainty in the image flux. Furthermore, we note that the
predicted image positions $\vec{r}_{i}$ are absent in this likelihood as this
information would be unknown in the case of a central image non-detection.

####  3.3.2 Central Image Detection

For the case of a detection of the central image, we define a hypothetical
observed central image with specific flux $f_{o}$ and position in the lens
plane $\vec{r}_{o}$. For each image formed for a given $M_{\rm SMBH}$, $\Delta
r$, and $r_{c}$, the likelihood of detection is assumed to be Gaussian
relative to the hypothetical observation:

$P(D|M_{\rm SMBH},\Delta
r,r_{c})\propto\exp(-\frac{1}{2}\left(\frac{(f_{i}-f_{o})^{2}}{\sigma_{f}^{2}}+\frac{(\vec{r}_{i}-\vec{r}_{o})^{2}}{\sigma_{r}^{2}}\right))$
(12)

Here, $f_{i}$ and $\vec{r}_{i}$ are the predicted image flux and position for
a given $M_{\rm SMBH}$, $\Delta r$, and $r_{c}$, while $\sigma_{f}$ and
$\sigma_{r}$ are the uncertainties on $f_{o}$ and $\vec{r}_{o}$. From our
method, this requires us to choose the location and flux of detection in each
system’s lens plane. This can theoretically be anywhere in the lens plane;
however, for simplicity, we choose the predicted central image location
$\vec{r}_{o}$ to be at the center $(0,0)$ of each galaxy lens macro-model. We
set $\sigma_{r}$ based on the resolution of an image of this detection.
Assuming an HST image, we adopt the astrometric precision of $\sim$0.03". We
convert this $\sigma_{r}$ into kpc units with the measured $D_{d}$ for each
lens. It turns out that this choice for $\sigma_{r}$ is relatively large,
therefore minimizing the importance of the positional dependence of the
likelihood of detection in equation 12. However, since a hypothetical
detection of the central image depends more on brightness rather than
position333This assumption is valid only because our detection technique
(Section 5) assumes the lens galaxy is invisible in the observing filters.,
and our choice of $\sigma_{f}$ is more restrictive, we justify this choice as
realistic. For $f_{o}$ we simply choose the equivalent of 9 magnitudes fainter
than the estimated brightest quad image for each lens system. This choice is
consistent with our choice of 10 magnitudes fainter for non-detections. We set
$\sigma_{f}=0.4f_{o}$, equivalent to $\sim$1 magnitude in uncertainty. This
allows us to encompass a wide range of produced image magnitudes and to
account for magnification dispersion from stellar microlensing (Dobler et al.,
2007).

Using this likelihood function in equation 1 and marginalizing over $M_{\rm
SMBH}$ gives us the posterior probability density function for $\Delta r$ and
$r_{c}$.

##  4 Results

###  4.1 Lens Macro-model Maps

Figure 3: Central image modelling for RXJ1131-1231 quad lens. We used 1576
SMBHs (not shown) randomly offset from the galaxy host center according to the
distribution in Shen et al. (2019). In this figure all SMBH have $M_{\rm
SMBH}\approx 2.11\times 10^{9}M_{\odot}$. (Other maps for this system use
different $M_{\rm SMBH}$ randomly chosen from its range shown in Table 1).
Offset SMBHs that produce central images (green stars), and saddle images
(inverted parity; orange and red squares) are plotted on the mass density
(left panel), and Fermat potential $\Psi_{\rm Fermat}$ (right panel) of one
model of RXJ1131-1231 ($\Psi_{\rm SMBH}$ is not included so as to avoid
cluttering the plot). Cases where only one central image is produced at the
location of the offset SMBH are excluded from this plot, but included in
subsequent calculations. The left and top colorbars indicate the magnitude of
the central and saddle images, respectively. The maps have scales of 0.0035
kpc pix-1. All our maps have 1250 pix arcsec-1. The asymmetry in the image
distribution visible in both the panels comes from modeling of the four quad
images. Figure 4: Histogram of relative magnifications in log space of central
and saddle images for a SMBH of $M_{\rm SMBH}\approx 2.11\times
10^{9}M_{\odot}$ for the case of RXJ1131-1231 shown in Figure 3. $\mu$ and
$\mu_{\rm BI}$ are the magnifications of the modeled central or saddle image,
and the brightest image in RXJ1131-1231, respectively. The two vertical dashed
lines correspond to image magnitudes of $m_{\rm AB}=30$ and $34$.

In this section we present examples of maps for some individual galaxy macro-
models (see Section 3.1), and properties of the corresponding central images
(see Section 3.2), before obtaining statistical results in Section 4.2.

Following the modelling procedure outlined in Section 3, we obtain maps for
each lens galaxy macro-model depicting the locations of SMBH that allow extra
central images to form. For illustrative purposes, the left panel of Figure 3
shows contours of one of these mass maps computed for the case for
RXJ1131-1231 with a $2.11\times 10^{9}M_{\odot}$ SMBH444As shown in Table 1,
this SMBH mass is roughly near the center of the $M_{\rm SMBH}$ range found
from the $M-\sigma$ relation. displaced according to VODKA distribution, and
with random azimuthal positions. The central and saddle images are presented
in these maps, as greenish and reddish points, respectively. The SMBH image is
typically extremely demagnified, far below observability, so these are not
included in the map. Similarly, to avoid cluttering the map, the SMBH
locations are not shown. Our probability calculations include SMBH locations
where no extra images form.

It can be seen from the map that there exists an elliptical ‘barrier’ between
the saddle point images (reddish points) and central images (greenish points),
of radius $\sim$0.02". This barrier is the critical curve formed by the
presence of SMBH. The closest images to this barrier are the brightest images
generated from the model, while the further away the images become dimmer. The
exact radius of this critical curve depends on the galaxy macro-model used.
The right panel of Figure 3 shows the same image locations plotted onto the
same galaxy macro-model’s Fermat potential $\Psi_{\rm
Fermat}=\frac{1}{2}r^{2}-\Psi_{\rm gal}$. The asymmetric Fermat potential of
this galaxy macro-model prefers that for the saddle and SMBH images to form,
SMBH positions should be in the right of the field, as indicated by the
density of saddle images in that region. In general, the shape of the Fermat
potential determines which regions of the field are discriminated against
having saddle and SMBH images, with more elliptical potentials strongly
favoring SMBH locations for central images along the semi-major axis. The
asymmetry of the mass distribution arising from quad-scale macro-modeling is
important for central image properties, yet is often ignored in the
literature.

For this model, the mean SMBH offset distance that produces extra central
images is 0.282 $\pm$ 0.113 kpc. In general, the production of extra central
images requires SMBH offsets greater than the mean $\Delta r$ from VODKA. This
trend is important because the failure to create extra central images means
the singular central image at the maximum of $\Psi_{\rm Fermat}$ will be
strongly demagnified.

Magnification statistics for this model are shown in Figure 4. While in this
case, $M_{\rm SMBH}$ is the same for all the models with the purpose of
showcasing the observed magnification variability for a given $M_{\rm SMBH}$,
our general procedure varies $M_{\rm SMBH}$ randomly for each model according
to the distribution outlined in Section 3.2. As with our map in the left
panel, SMBH locations where no central images form are excluded. In general,
SMBH locations with no central images make up $\sim$50% of all offset
distances, meaning roughly half of all SMBH offsets will not yield observable
images. The vertical dashed lines correspond to image magnitudes of $m_{\rm
AB}=30$ and $34$ in the F218W HST filter. For this particular model we find
that most ($\sim$62%) of the central images are brighter than $m_{\rm
AB}=30.00$. The corresponding saddle images have a much broader distribution
and are generally much fainter than the central images. Combining each of
these maps and magnification histograms for all the models in each QSO allows
us to place statistical constraints on $r_{c}$ and $\Delta r$.

Figure 5: Left: SMBH (all have $M_{\rm SMBH}\approx 3.73\times
10^{8}M_{\odot}$) locations (black dots) corresponding to central images
(green stars) and saddle images (orange squares) generated from a subset of
100 mass models for HE0435-1223 falling within the $r=0.0050"$ region (blue
circle) centered at $(0.0020",0.0108")$ from the mass density peak. Orange
squares depict simultaneously produced saddle images. The right and top
colorbars indicate the magnitude of the central and saddle images,
respectively. Right: Central (green stars) and saddle (orange squares) image
distances from the mass density peak as a function of their corresponding SMBH
offset distance for central images falling within the $r=0.0050"$ region of
the top panel. All our maps have 1250 pix arcsec-1. The vertical streak
pattern indicates the range of image locations for each SMBH distance given by
the 100 model subset. The right and top colorbars show the magnitude of the
central and saddle images, respectively.

Our statistical inference method (see Section 3.3) requires a likelihood
function dependent on the magnitude and position of the predicted central
image. As a visualization of a predicted detection, we define a small circular
region within the lens plane such that any central images falling within the
region are considered to be detected with magnitudes according to their
respective magnification. The left panel of Figure 5 depicts one example
region for a subset of 100 HE0435-1223 galaxy macro-models, centered at
$\vec{r}_{o}=(0.0020",0.0108")$ with a radius of 0.0050". For this lens, this
radius corresponds to $\sim$0.03 kpc. The radius of the region can be thought
of as the uncertainty $\sigma_{r}$ of a detection at $\vec{r}_{o}$. All
central images found within the region and their corresponding saddle images
and SMBH locations are shown. The streak pattern is due to the fact that
saddle images form on the line connecting the center of the lensing potential
and the SMBH, shown here as black points. The SMBH offset distribution average
is measured to be 0.187 kpc with a corresponding standard deviation of 0.017
kpc.

The right panel of Figure 5 presents the information for the same circular
region shown in the left panel, but plots the central and saddle image
distances from the lensing potential peak against the corresponding SMBH
distance. It is clear from the figure that the central and saddle images that
form closer together, and hence closer to the critical curve, are brighter
than those that form further away. Additionally, there is a general positive
correlation between the SMBH position and its corresponding central and saddle
image distances.

###  4.2 Constraints on Lens Galaxy Parameters

The previous section highlights how individual lens galaxy macro-models from
our central image analysis can yield constraints on $\Delta r$ and $r_{c}$.
This is largely the outcome of applying the analysis described in Section 3.2.
The main goal of this work is to combine all models for each QSO in our sample
to obtain general constraints on each lens galaxy’s parameters. For this we
follow the analysis procedure described in Section 3.3. For each QSO, 1000
lens galaxy macro-models were generated. For each of these, and a randomly
picked SMBH $\Delta r$ offset from the VODKA distribution, and randomly picked
SMBH mass from the appropriate distribution, we determined how many central
images were produced, and each image’s position and magnification. This allows
us to create contour probability maps in $(\Delta r,r_{c})$ space based on a
predicted non-detection (see equation 11) or detection (see equation 12). We
extend our statistical inference to 2 outcomes: (i) Non-detections and (ii)
Detections. A summary of our constraints on lens galaxy parameter
distributions for each quad is presented in Table 2. The posterior PDFs
$P(\Delta r,r_{c}|D)$ for the cases of a non-detection and detection in each
quad are shown in Figures 6 and 7. Individual marginal posterior PDFs for each
lens are shown in Figure 8.

Table 2: Constraints on Individual Lens Galaxy Parameters

QSO (Det./Non-det.) |  $\widehat{\Delta r}$ [kpc] | 95% CI ($\Delta r$) [kpc] |  $\widehat{r_{c}}$ [kpc] | 95% CI ($r_{c}$) [kpc] |  $\widehat{b}$ [10-4 kpc km-1 s] | 95% CI ($b$) [10-4 kpc km-1 s] |
---|---|---|---|---|---|---|---
HE0435-1223 (Non-det.) | 0.066 | $0.015<\Delta r<0.367$ | 0.042 | $0.010<r_{c}<0.179$ | 1.90 | $0.459<b<8.069$ |
HE0435-1223 (Det.) | 0.154 | $0.081<\Delta r<0.397$ | 0.095 | $0.042<r_{c}<0.165$ | 4.26 | $1.90<b<7.41$ |
PG1115+080 (Non-det.) | 0.118 | $0.015<\Delta r<0.353$ | 0.044 | $0.010<r_{c}<0.167$ | 1.54 | $0.365<b<5.832$ |
PG1115+080 (Det.) | 0.242 | $0.081<\Delta r<0.426$ | 0.081 | $0.025<r_{c}<0.145$ | 2.84 | $0.886<b<5.051$ |
RXJ1131-1231 (Non-det.) | 0.110 | $0.022<\Delta r<0.367$ | 0.056 | $0.013<r_{c}<0.181$ | 1.73 | $0.407<b<5.616$ |
RXJ1131-1231 (Det.) | 0.162 | $0.059<\Delta r<0.419$ | 0.126 | $0.022<r_{c}<0.206$ | 3.91 | $0.692<b<6.374$ |
SDSSJ0924+0219 (Non-det.) | 0.073 | $0.015<\Delta r<0.367$ | 0.037 | $0.008<r_{c}<0.213$ | 1.72 | $0.349<b<9.928$ |
SDSSJ0924+0219 (Det.) | 0.176 | $0.088<\Delta r<0.404$ | 0.125 | $0.032<r_{c}<0.174$ | 5.82 | $1.49<b<8.10$ |
WFI2033-4723 (Non-det.) | 0.073 | $0.015<\Delta r<0.360$ | 0.071 | $0.013<r_{c}<0.230$ | 2.86 | $0.518<b<9.041$ |
WFI2033-4723 (Det.) | 0.162 | $0.081<\Delta r<0.404$ | 0.117 | $0.034<r_{c}<0.209$ | 4.70 | $1.35<b<8.37$ |
MG0414+0534 (Non-det.) | 0.147 | $0.022<\Delta r<0.375$ | 0.089 | $0.019<r_{c}<0.234$ | 3.45 | $0.728<b<9.020$ |
MG0414+0534 (Det.) | 0.184 | $0.037<\Delta r<0.404$ | 0.132 | $0.033<r_{c}<0.273$ | 5.08 | $1.27<b<10.5$ |
RXJ0911+0551 (Non-det.) | 0.162 | $0.022<\Delta r<0.375$ | 0.485 | $0.084<r_{c}<0.868$ | 18.68 | $3.23<b<33.46$ |
RXJ0911+0551 (Det.) | 0.184 | $0.022<\Delta r<0.389$ | 0.310 | $0.119<r_{c}<0.641$ | 11.96 | $4.57<b<24.72$ |

The columns list the QSO and constraint scenario (based on whether or not a
central image is detected), modes for SMBH offset $\widehat{\Delta r}$, core
size $\widehat{r_{c}}$, and $b$ parameter $\widehat{b}$, and 95% credible
intervals for SMBH offset, core size, and $b$ parameter. All modes and
credible intervals presented are for the posterior PDF distributions for each
parameter. RXJ0911+0551 is treated independently since its core size prior
distribution is skewed to large core sizes due to the presence of large
external shear.

Figure 6: Individual Joint Posterior PDFs $P(\Delta r,r_{c}|D)$ for the
hypothetical case of a non-detection of the central image in all QSOs in our
sample. Marginalizing over each axis gives individual posterior PDFs for core
size $P(r_{c}|D)$ and SMBH offset distance $P(\Delta r|D)$. The vertical red
dashed lines indicate the percentage of SMBH offset by $<$0.1 and $<$0.5 kpc
from VODKA (Shen et al., 2019).

Figure 7: Individual Joint Posterior PDFs $P(\Delta r,r_{c}|D)$ for the
hypothetical case of a central image detection in all QSOs in our sample. The
central image detection is assumed to be located at the center of each lens
galaxy macro-model with a brightness of 9 magnitudes fainter than the
brightest quad image in each lens. Marginalizing over each axis gives
individual posterior PDFs for core size $P(r_{c}|D)$ and SMBH offset distance
$P(\Delta r|D)$. The vertical red dashed lines indicate the percentage of SMBH
offset by $<$0.1 and $<$0.5 kpc from VODKA (Shen et al., 2019).

Figure 8: Individual posterior PDFs for SMBH offset distance $\Delta r$ (top
row), core size $r_{c}$ (middle row)), and $r_{c}-\sigma$ proportionality
constant $b$ (bottom row) for all QSOs in our sample. The case of a central
image non-detection (left column) and detection (right column) yield
independent constraints on lens galaxy parameters. The prior distributions for
each parameter are shown in yellow.  Table 3: Constraints on Lens Galaxy
Parameters

Constraint Scenario |  $\widehat{\Delta r}$ [kpc] | 95% CI ($\Delta r$) [kpc] |  $\widehat{r_{c}}$ [kpc] | 95% CI ($r_{c}$) [kpc] |  $\widehat{b}$ [10-4 kpc km-1 s] | 95% CI ($b$) [10-4 kpc km-1 s] |
---|---|---|---|---|---|---|---
All Detections | 0.242 | $0.154<\Delta r<0.286$ | 0.107 | $0.083<r_{c}<0.136$ | 6.28 | $4.47<b<7.87$ |
All Non-detections | 0.132 | $0.037<\Delta r<0.228$ | 0.049 | $0.019<r_{c}<0.107$ | 3.11 | $0.849<b<5.830$ |
Some Det./Non-det. | 0.191 | $0.110<\Delta r<0.264$ | 0.068 | $0.033<r_{c}<0.102$ | 4.47 | $2.21<b<6.96$ |

The columns list the constraint scenario (based on whether or not a central
image is detected), modes for SMBH offset $\widehat{\Delta r}$, core size
$\widehat{r_{c}}$, and $b$ parameter $\widehat{b}$, and 95% credible intervals
for SMBH offset, core size, and $b$ parameter. All modes and 95% credible
intervals presented are for the marginalized posterior PDF distributions for
each parameter assuming the entire sample is governed by a single constraint
value for $\Delta r$, $r_{c}$, and $b$. The scenario for Some Det./Non-det.
assumes a central image detection in PG1115+080 and WFI2033-4723 and non-
detections in the rest.

Let us first consider the most likely case of central image non-detections.
Our inference method in this case simply checks the probability of an image
being brighter than the $f_{\rm crit}$ threshold (see equation 11). The
posterior PDFs (see equation 1) for each lens in this case is shown in Figure
6. The case of non-detections prefers smaller $\Delta r$ values than the case
of detections. In general, we can say that non-detections of central images
imply that the central SMBH is very likely to be well centered on the galaxy’s
lensing potential. As with $\Delta r$, smaller values of $r_{c}$ and $b$ are
strongly favored for the case of non-detections. Qualitatively, this is a
similar conclusion to those of previous works (e.g. Quinn et al., 2016).

Next, we consider the unlikely, but more interesting, hypothetical outcome
that the central image was detected in the 6 lens systems in our sample. The
resulting posterior PDFs are shown in Figure 7. We can see that in the event
of a detection at the chosen $\vec{r}_{o}$ with $f_{o}$ in each lens system,
$\Delta r$ is most likely to be $>$0.1 kpc. The case of detections favors
larger $\Delta r$. Similarly for $r_{c}$, the case of detections constrains to
larger core sizes and $b$ values.

Figure 9:  Combined posterior PDFs for $r_{c}-\sigma$ proportionality constant
$b$. The constraints from a hypothetical central image detection in each QSO
are shown in blue. The constraints from non-detections of the central image in
each QSO are shown in red. Constraints from a central image detection in
PG1115+080 and WFI2033-4723 but non-detections in all others are shown in
black. The prior distribution is shown in yellow. In general, non-detections
of the central image favor smaller $b$ values.

Lastly, since the $b$ parameter is assumed to be a global parameter for all
lenses, we can multiply the 2D posteriors together and marginalize according
to equation 8 to get a general constraint on $b$. We do this for three
scenarios: (i) All Non-detections in the sample, (ii) All detections in the
sample, and (iii) A combination of Non-detections and detections. For this
third case, we assumed that the lenses with the top two brightest observed
quad images, PG1115+080 and WFI2033-4723, had hypothetical central image
detections. With this assumption, we multiplied together individual non-
detection posteriors for the rest of the lens systems in our sample with the
detection posteriors for the two chosen lens systems that we promoted to
detections. Each scenario yields independent constraints on $b$. The summary
of these constraints is detailed in Table 3 and Figure 9. In general, central
image detections favor steeper $r_{c}-\sigma$ ($b=6.28^{+1.59}_{-1.81}\times
10^{-4}$ kpc km-1 s), while the inverse is true for non-detections
($b=3.11^{+2.72}_{-2.26}\times 10^{-4}$ kpc km-1 s). The combination case
yields $b=4.47^{+2.49}_{-2.26}\pm 1.40\times 10^{-4}$ kpc km-1 s, between the
detection and non-detection case constraint. In fact, in each of the three
scenarios, the resulting distributions for $b$ have credible intervals smaller
than that of the prior, indicating that the search for central images with
this analysis can yield stronger constraints regardless of outcome.

###  4.3 Rough Constraints for $\Delta r$ and $r_{c}$

Figure 10:  Combined posterior PDFs for SMBH offset distance $\Delta r$ (top
panel) and core size $r_{c}$ (bottom panel) assuming these parameters are the
same for all lens galaxies

. The constraints from a hypothetical central image detection in each QSO are
shown in blue. The constraints from non-detections of the central image in
each QSO are shown in red. Constraints from a central image detection in
PG1115+080 and WFI2033-4723 but non-detections in all others are shown in
black. The prior distributions for each parameter are shown in yellow. In
general, non-detections of the central image favor smaller SMBH offset
distance and core sizes.

Figure 11:  Posterior PDFs $P(\Delta r,r_{c}|D)$ for the hypothetical case of
central image non-detections in all QSOs in our sample (top panel), central
image detections in all QSOs in our sample (middle panel), and a central image
detection in PG1115+080 and WFI2033-4723, and non-detections in all other QSOs
in our sample (bottom panel). For detections, the central image detection is
assumed to be located at the center of each lens galaxy macro-model with a
brightness of 9 magnitudes fainter than the brightest quad image in each lens.
Marginalizing over each axis gives individual posterior PDFs for core size
$P(r_{c}|D)$ and SMBH offset distance $P(\Delta r|D)$, assuming $r_{c}$ and
$\Delta r$ are the same for all lens galaxies. The vertical red dashed lines
indicate the percentage of SMBH offset by $<$0.1 and $<$0.5 kpc from VODKA
(Shen et al., 2019).

Our results so far have shown that the non-detection or detection of a central
image can constrain the lens galaxy SMBH offset and core radius. To place a
general constraint on these parameters, we assumed a proportionality with
$r_{c}$ and $\sigma$ and marginalized over $\Delta r$ to obtain the posterior
for this proportionality constant $b$. We are able to do this since $b$ is a
parameter that is assumed to be the same for all lens galaxies. In this
section, we similarly assume that all lens galaxies share the same SMBH offset
and core radius in order to directly constrain $\Delta r$ and $r_{c}$. Since
this is an unrealistic assumption, we treat these results as a "back of
envelope" estimation rather than a strict constraint. In fact, since the
estimated $\Delta r$ and $r_{c}$ for individual lens galaxies are
approximately similar to one another (see Table 2), this exercise serves as a
useful estimate of the rough scale of SMBH offsets and core radii that one can
expect in future analyses.

With this assumption established, we follow equations 9 and 10 to derive
posterior PDFs for $\Delta r$ and $r_{c}$, respectively. We do this for the
same 3 scenarios as for the constraint on $b$ (All Non-detections, All
Detections, and a Combination of Detections and Non-detections. As before, we
assume a detection in PG1115+080 and WFI2033-4723 but non-detections in all
others for the third case.). A summary of all the posterior PDF constraints
for $\Delta r$ and $r_{c}$ is shown in Figure 10 and Table 3. Similarly, the
total multiplied together 2D posteriors for each case are shown in Figure 11.

For all non-detections in our sample, the mode of the posterior $P(\Delta
r|D)$ is 0.132 kpc with a 95% credible interval of $0.037<\Delta r<0.228$ kpc.
This distribution is skewed to smaller $\Delta r$, consistent with the trend
that non-detections of central images imply that the central SMBH is very
likely to be well centered on the galaxy’s lensing potential. With core size,
the $r_{c}$ posterior mode is 0.049 kpc. For all detections, the $\Delta r$
posterior mode is 0.242 kpc with a 95% credible interval of $0.154<\Delta
r<0.286$ kpc. This is a stronger constraint than that of the case of all non-
detections. the $r_{c}$ posterior mode for all detections is 0.107 kpc. In the
combination case, the constraints on $\Delta r$, $r_{c}$, and $b$ lie in
between those found in the cases of non-detections and detections. This
intuitively makes sense as overall, central image detections favor larger SMBH
offsets and core sizes and steeper $r_{c}-\sigma$, while the inverse is true
for non-detections. In addition, the posterior means for each parameter are
similar to the prior means, but with smaller credible intervals.

##  5 A New Detection Technique

The results of our analyses can be applied to future observations of the
elusive central image in hopes of constraining various galaxy properties. Here
we outline a novel observing strategy that can help improve the chances of
central image detection.

In addition to being significantly demagnified, the central image is often
superimposed by the lens galaxy in optical wavelengths. Therefore, to increase
the likelihood of detecting the central image, it is important to increase the
brightness contrast between it and the lens galaxy, that is, to maximize the
central image flux and minimize the lens galaxy flux. To achieve that, the
observed wavelength of the galaxy should be shorter than the rest frame
wavelength of the Balmer break ($3646\,$Å). The lens galaxies will be largely
invisible in these wavelengths, because they are predominantly old, massive
ellipticals, and thus have spectral energy densities that peak in red
wavelengths and drop off toward the blue and ultraviolet. Furthermore, QSO
sources have a blue power-law continuum, with high flux densities in bluer
wavelengths. Additionally, to avoid the wavelength range absorbed by
intervening intergalactic medium (Ly-$\alpha$ forest), the observed wavelength
of the QSO should be greater than the rest frame wavelength of $\sim 1200\,$Å,
which means the source QSOs should have redshifts above $\sim 1.5$ (Francis et
al., 1991; Hewett & Wild, 2010). In practice, this only becomes important for
sources at redshifts greater than $\sim$2\. This implies that the ideal
scenario for which this technique can be utilized is the highest possible
redshift for the lens galaxy for any given source redshift.

In order to accomplish this, we suggest observing in the bluest wavelengths
achievable with the UV filters on the Hubble Space Telescope (HST). In blue
filters, the contrast of the lens galaxy and central image can be increased to
the point where the central image can become detectable.

In addition to brightness considerations, the positioning of the central image
is better constrained for compact lens sources, which is another reason why
QSOs are ideal sources.

As a preliminary illustration of this technique, we consider the HST UV filter
F275W and Optical filter F555W, which have effective wavelengths555Effective
wavelength values taken from the SVO Filter Profile Service:
http://svo2.cab.inta-csic.es/theory/fps/ of $2713.86\,$Å and $5326.96\,$Å,
respectively. Observing the gravitationally lensed QSO quad HE0435-1223
(properties presented in Table 1) in F555W corresponds to rest frame QSO and
lens emission at $\sim\\!1980\,$Å and $\sim\\!3664\,$Å , respectively. In this
filter, the lens emission is longward of the Balmer break limit, so the lens
galaxy would obscure the central image. However, observing in F275W would
correspond to rest frame emission at $\sim\\!1008\,$Å and $\sim\\!1866\,$Å for
the QSO and lens, respectively. This wavelength corresponds to a regime where
the QSO continuum flux is large. Similarly, the lens rest frame is in the
regime that should have minimal flux from the lens galaxy. Therefore, our
central image detection technique can be potentially viable for this system,
as illustrated in Figure 12.

Figure 12:  HE0435-1223 in HST optical filter F555W (left panel) and UV filter
F275W (right panel). Lens and source redshifts are $z_{l}=0.455$ and
$z_{s}=1.693$ for this system. In F555W, the lensing galaxy is visible in the
center of the image since its rest frame emission in this filter lies just
longward of Balmer break. In F275W, the lensing galaxy is not visible, while
the QSO images are bright, showing how UV searches for the central image can
be viable. $1"=5.8\,$kpc at the redshift of the lens.

##  6 Discussion and Conclusions

Using a sample of 7 quad lenses, we modelled the central image and placed
constraints on lens galaxy parameters for various scenarios of central image
detection. We estimated 3 lens galaxy parameters: $b$ (the $r_{c}-\sigma$
proportionality constant), SMBH offset $\Delta r$, and galaxy core radius
$r_{c}$. Constraints on these parameters for individual lenses are presented
in Table 2.

Here we concentrate on the combined constraints for all lenses, assuming
roughly that $\Delta r$ and $r_{c}$ are the same for all galaxies, while $b$
is defined to be this way. Our main results and constraints are presented in
Table 3, and are as follows:

  * •

All the cases we considered are hypothetical. However, the one that comes
closest to the current observational constraints is the case with no
detections of the central image in any of the lenses (red curves in Figure
10). For these, the SMBH offset from the center of the galaxy host is
132${}^{+96}_{-95}$ pc, the galaxy core radius is 49${}^{+58}_{-30}$ pc, and
the constant of proportionality $b$ (eq. 6), relating the core radius and the
galaxy line-of-sight velocity dispersion is $3.11^{+2.72}_{-2.26}\times
10^{-4}$ kpc km-1 s. Each of the constraint distributions yield much tighter
constraints for each parameter than their respective priors, and favor smaller
values.

  * •

The more exciting, but unlikely, case of a central image detection in each
lens also yields tight constraints on lens galaxy parameters (blue curves in
Figure 10). For these, the SMBH offset from the galaxy host center is
242${}^{+44}_{-88}$ pc, the galaxy core radius is 107${}^{+29}_{-24}$ pc, and
the average constant of proportionality $b$, relating the core radius and the
galaxy line-of-sight velocity dispersion is $6.28^{+1.59}_{-1.81}\times
10^{-4}$ kpc km-1s. A central image detection will therefore imply that the
galaxy core sizes and SMBH offsets are larger than the average of their
respective priors.

  * •

The case of a combination of detections and non-detections in our sample
yields independent constraints on lens galaxy parameters. With our assumed
hypothetical detection in PG1115+080 and WFI2033-4723, justified by their
bright quad images, the average SMBH offset, galaxy core radius, and constant
of proportionality relating the core radius and the galaxy line-of-sight
velocity dispersion are 191${}^{+73}_{-81}$ pc, 68${}^{+34}_{-35}$ pc, and
$4.47^{+2.49}_{-2.26}\times 10^{-4}$ kpc km-1s, respectively. In general, even
just a single detection of the central image will strongly influence the
constraints to larger SMBH offset, core radius, and $b$.

  * •

Regardless of whether or not the central image is detected in a lens or not,
tighter constraints can be placed on SMBH offset, galaxy core radius, and the
$r_{c}-\sigma$ proportionality constant. A notable exception to this result is
exemplified in RXJ0911+0551 (see Section A). Due to the presence of a nearby
galaxy cluster which leads to a wide radial range of its image distribution,
and requires a large external shear, RXJ0911+0551 does not constrain any lens
galaxy parameter regardless of detection. From this, it is likely that lenses
that require large external shear due to the presence of nearby external mass
will not be useful in constraining lens galaxy parameters in future studies.

  * •

Our analysis demonstrates that the quad-scale lens macro-model is important
for deriving properties of the central images. This is especially true in the
case of very asymmetric systems, like RXJ0911+0551.

  * •

While our results are not restricted to any particular observational filter,
we recommend a novel observing strategy utilizing UV wavelength filters (see
Section 5). Given that the lens galaxy peaks in red wavelengths and the QSO
source peaks in blue wavelengths, UV filters have a potential use in searches
for central images since they will obscure the lens galaxy and possibly allow
for easier detection of the elusive demagnified blue central image.

In this paper we considered how the fluxes of central QSO images are affected
by the presence of SMBH, which can be offset from the center of its galaxy
host. Because the central galaxy region is dominated by stars, the central
image flux can also be affected by stellar microlensing (Dobler et al., 2007),
which would broaden the predicted flux distribution by $\sim$1 mag. Future
studies should continue to account for this in central image modelling.

While our rough calculations for a single value of $\Delta r$ and $r_{c}$ for
all galaxies (see Section 4.3) are useful in understanding the general scale
of their underlying distribution, a more insightful result would be direct
constraints on the distributions of $\Delta r$ and $r_{c}$. This would require
a more complicated theoretical and statistical model on top of our outlined
framework, and would be the subject of future study.

Future observations with the Vera C. Rubin Observatory, Gaia, and Pan-Starrs
will discover many new lensed systems (Marshall et al., 2010; Lemon et al.,
2019; Cañameras et al., 2020) allowing the extension of our analysis to a
larger sample, giving tighter constraints, and increasing the probability of
detecting a central image. Similarly, combining our constraints with recent
results from gravitational lensing of source AGN at $z>1$ (Millon et al.,
2022; Spingola et al., 2022) can extend our constraints to the inner regions
of higher redshift galaxies. Furthermore, our results can be compared with
simulations (Tremmel et al., 2018; Volonteri et al., 2020; Katz et al., 2020).

In this work we have shown that searches for the central image in quads can
yield tight constraints on lens galaxy parameters regardless of detection or
non-detection. Therefore, we recommend commencing new observation programs for
the central image to formalize new constraints for SMBH offset, galaxy core
radius, and $r_{c}-\sigma$ proportionality constant.

## Acknowledgements

The authors would like to thank Lindsey Gordon, Galin Jones, John Hamilton
Miller Jr., and Sarah Taft for useful suggestions and discussions.

## Data Availability

Data generated from this article will be shared upon reasonable request to the
corresponding author.

## References

  * Bade et al. (1997) Bade N., Siebert J., Lopez S., Voges W., Reimers D., 1997, A&A, 317, L13
  * Barrera et al. (2021) Barrera B., Williams L. L. R., Coles J. P., Denzel P., 2021, The Open Journal of Astrophysics, 4, 12
  * Begelman et al. (1980) Begelman M. C., Blandford R. D., Rees M. J., 1980, Nature, 287, 307
  * Bonvin et al. (2018) Bonvin V., et al., 2018, A&A, 616, A183
  * Bonvin et al. (2019) Bonvin V., et al., 2019, A&A, 629, A97
  * Bruderer et al. (2016) Bruderer C., Read J. I., Coles J. P., Leier D., Falco E. E., Ferreras I., Saha P., 2016, MNRAS, 456, 870
  * Burud et al. (1998) Burud I., et al., 1998, ApJ, 501, L5
  * Cañameras et al. (2020) Cañameras R., et al., 2020, A&A, 644, A163
  * Chen et al. (2022) Chen N., Ni Y., Tremmel M., Di Matteo T., Bird S., DeGraf C., Feng Y., 2022, MNRAS, 510, 531
  * Chu et al. (2022) Chu A., Boldrini P., Silk J., 2022, arXiv e-prints, p. arXiv:2212.13277
  * Courbin et al. (2011) Courbin F., et al., 2011, A&A, 536, A53
  * Di Matteo et al. (2005) Di Matteo T., Springel V., Hernquist L., 2005, Nature, 433, 604
  * Dobler et al. (2007) Dobler G., Keeton C. R., Wambsganss J., 2007, MNRAS, 377, 977
  * Dullo et al. (2021) Dullo B. T., Gil de Paz A., Knapen J. H., 2021, ApJ, 908, 134
  * Ferrarese & Merritt (2000) Ferrarese L., Merritt D., 2000, ApJ, 539, L9
  * Ferrarese et al. (2006) Ferrarese L., et al., 2006, ApJS, 164, 334
  * Francis et al. (1991) Francis P. J., Hewett P. C., Foltz C. B., Chaffee F. H., Weymann R. J., Morris S. L., 1991, ApJ, 373, 465
  * Ghosh et al. (2020) Ghosh A., Williams L. L. R., Liesenborgs J., 2020, MNRAS, 494, 3998
  * Gomer & Williams (2018) Gomer M. R., Williams L. L. R., 2018, MNRAS, 475, 1987
  * Gültekin et al. (2009) Gültekin K., et al., 2009, ApJ, 698, 198
  * Hewett & Wild (2010) Hewett P. C., Wild V., 2010, MNRAS, 405, 2302
  * Hewitt et al. (1992) Hewitt J. N., Turner E. L., Lawrence C. R., Schneider D. P., Brody J. P., 1992, AJ, 104, 968
  * Hezaveh et al. (2015) Hezaveh Y. D., Marshall P. J., Blandford R. D., 2015, ApJ, 799, L22
  * Hjorth et al. (2002) Hjorth J., et al., 2002, ApJ, 572, L11
  * Hwang et al. (2020) Hwang H.-C., Shen Y., Zakamska N., Liu X., 2020, ApJ, 888, 73
  * Katz et al. (2020) Katz M. L., Kelley L. Z., Dosopoulou F., Berry S., Blecha L., Larson S. L., 2020, MNRAS, 491, 2301
  * Keeton (2011) Keeton C. R., 2011, ascl, pp ascl–1102
  * Kormendy & Richstone (1995) Kormendy J., Richstone D., 1995, ARA&A, 33, 581
  * Koss et al. (2018) Koss M. J., et al., 2018, Nature, 563, 214
  * Lemon et al. (2019) Lemon C. A., Auger M. W., McMahon R. G., 2019, MNRAS, 483, 4242
  * Lemon et al. (2022) Lemon C., et al., 2022, arXiv e-prints, p. arXiv:2206.07714
  * MacLeod et al. (2015) MacLeod C. L., et al., 2015, ApJ, 806, 258
  * Mao & Witt (2012) Mao S., Witt H. J., 2012, MNRAS, 420, 792
  * Mao et al. (2001) Mao S., Witt H. J., Koopmans L. V. E., 2001, MNRAS, 323, 301
  * Marshall et al. (2010) Marshall P. J., et al., 2010, in American Astronomical Society Meeting Abstracts #215. p. 401.15
  * Millon et al. (2020) Millon M., et al., 2020, A&A, 640, A105
  * Millon et al. (2022) Millon M., Dalang C., Lemon C., Sluse D., Paic E., Chan J. H. H., Courbin F., 2022, arXiv e-prints, p. arXiv:2207.00598
  * Muller et al. (2020) Muller S., Jaswanth S., Horellou C., Martí-Vidal I., 2020, A&A, 641, L2
  * Nasim et al. (2021) Nasim I. T., Gualandris A., Read J. I., Antonini F., Dehnen W., Delorme M., 2021, MNRAS, 502, 4794
  * Nightingale et al. (2019) Nightingale J. W., Massey R. J., Harvey D. R., Cooper A. P., Etherington A., Tam S.-I., Hayes R. G., 2019, MNRAS, 489, 2049
  * Quinn et al. (2016) Quinn J., et al., 2016, MNRAS, 459, 2394
  * Rantala et al. (2018) Rantala A., Johansson P. H., Naab T., Thomas J., Frigo M., 2018, ApJ, 864, 113
  * Rengstorf et al. (2006) Rengstorf A. W., Brunner R. J., Wilhite B. C., 2006, AJ, 131, 1923
  * Rusin et al. (2005) Rusin D., Keeton C. R., Winn J. N., 2005, ApJ, 627, L93
  * Rusu et al. (2020) Rusu C. E., et al., 2020, MNRAS, 498, 1440
  * Sesar et al. (2007) Sesar B., et al., 2007, AJ, 134, 2236
  * Shen et al. (2019) Shen Y., Hwang H.-C., Zakamska N., Liu X., 2019, ApJ, 885, L4
  * Skipper & Browne (2018) Skipper C. J., Browne I. W. A., 2018, MNRAS, 475, 5179
  * Sluse et al. (2007) Sluse D., Claeskens J. F., Hutsemékers D., Surdej J., 2007, A&A, 468, 885
  * Sluse et al. (2019) Sluse D., et al., 2019, MNRAS, 490, 613
  * Spingola et al. (2022) Spingola C., Schwartz D., Barnacka A., 2022, ApJ, 931, 68
  * Suyu et al. (2013) Suyu S. H., et al., 2013, ApJ, 766, 70
  * Tamura et al. (2015) Tamura Y., Oguri M., Iono D., Hatsukade B., Matsuda Y., Hayashi M., 2015, PASJ, 67, 72
  * Thomas et al. (2016) Thomas J., Ma C.-P., McConnell N. J., Greene J. E., Blakeslee J. P., Janish R., 2016, Nature, 532, 340
  * Tonry (1998) Tonry J. L., 1998, AJ, 115, 1
  * Tortora et al. (2004) Tortora C., Piedipalumbo E., Cardone V. F., 2004, MNRAS, 354, 343
  * Tremmel et al. (2018) Tremmel M., Governato F., Volonteri M., Quinn T. R., Pontzen A., 2018, MNRAS, 475, 4967
  * Treu & Koopmans (2002) Treu T., Koopmans L. V. E., 2002, MNRAS, 337, L6
  * Volonteri et al. (2020) Volonteri M., et al., 2020, MNRAS, 498, 2219
  * Williams & Zegeye (2020) Williams L. L. R., Zegeye D., 2020, The Open Journal of Astrophysics, 3, 10
  * Winn et al. (2004) Winn J. N., Rusin D., Kochanek C. S., 2004, Nature, 427, 613
  * Wong et al. (2015) Wong K. C., Suyu S. H., Matsushita S., 2015, ApJ, 811, 115
  * Wong et al. (2017a) Wong K. C., et al., 2017a, MNRAS, 465, 4895
  * Wong et al. (2017b) Wong K. C., Ishida T., Tamura Y., Suyu S. H., Oguri M., Matsushita S., 2017b, ApJ, 843, L35

##  Appendix A The Special Case of RXJ0911+0551

Figure 13:  Central image modelling for RXJ0911+0551 quad lens. We used 1576
SMBHs (not shown) randomly offset from the galaxy host center according to the
distribution in Shen et al. (2019). In this figure all SMBH have $M_{\rm
SMBH}\approx 9.03\times 10^{8}M_{\odot}$. (Other maps for this system use
different $M_{\rm SMBH}$ randomly chosen from its range shown in Table 1).
Offset SMBHs that produce central images (green stars) and saddle images
(inverted parity; orange squares) are plotted on the mass density (left panel)
and Fermat potential $\Psi_{\rm Fermat}$ (right panel) of one model of
RXJ0911+0551 ($\Psi_{\rm SMBH}$ is not included so as to avoid cluttering the
plot). Cases where only one central image is produced at the location of the
offset SMBH are excluded from this plot, but included in subsequent
calculations. The left and top colorbars indicate the magnitude of the central
and saddle images, respectively. The maps have scales of 0.0059 kpc pix-1. All
our maps have 1250 pix arcsec-1.

Figure 14:  Posterior PDF, $P(\Delta r,r_{c}|D)$ for the hypothetical case of
a central image non-detection (left panel) and detection (right panel) in
RXJ0911+0551. Note that the vertical axis extends to larger core radii than in
Figures 6 and 7. The detected central image is assumed to be located at the
center of each lens galaxy macro-model, with a brightness of 9 magnitudes
fainter than the brightest quad image in each lens. Marginalizing over each
axis gives individual posterior PDFs for core size $P(r_{c}|D)$ and SMBH
offset distance $P(\Delta r|D)$. The vertical red dashed lines indicate the
percentage of SMBH offset by $<$0.1 and $<$0.5 kpc from VODKA (Shen et al.,
2019).

Figure 15:  RXJ0911+0551 posterior PDFs for SMBH offset distance $\Delta r$
(top panel), core size $r_{c}$ (middle panel), and $r_{c}-\sigma$
proportionality constant $b$ (bottom panel). The constraints from a
hypothetical central image detection (non-detection) in this lens are shown in
blue (red). The prior distributions for each parameter are shown in yellow.

RXJ0911+0551 is influenced by a large external shear from a nearby cluster.
Applying the same initial range of quad-scale galaxy parameters as for the
other 6 quads, but a somewhat larger external shear (see Section 3.1) results
in the core size distribution $P(r_{c})$ to be skewed to much larger values
than the rest of the sample. To compensate, we analyze it independently to get
separate constraints on $\Delta r$, $r_{c}$, and $b$. The influence of the
external shear can be seen in the lens mass-model maps in Figure 13. Compared
with the lens maps for RXJ1131-1231 (see Figure 3), which is the system most
similar to RXJ0911+0551 of the 6 we use, we can easily see that the region of
the central images (green stars) is displaced from the center of the galaxy
density map further in RXJ0911+0551 than in RXJ1131-1231. Since the VODKA
distribution is with respect to the galaxy center, the large external shear
and offset position of the source displace the maximum of $\Psi_{\rm Fermat}$
further from the center, leading to small offset SMBH being able to produce
extra images (this did not occur in the other lenses with centrally symmetric
$\Psi_{\rm Fermat}$).

The constraints for this system are listed in Table 2. The posterior PDFs
$P(\Delta r,r_{c}|D)$ for the cases of a non-detection and detection are shown
in Figure 14. For a central image non-detection, the posterior PDF appears
similar with the detection posterior, except with smaller core sizes allowed.
The non-detection posterior PDF for $\Delta r$ has a mode of 0.162 kpc and 95%
credible interval of $0.022<\Delta r<0.375$ kpc.

It is important to note that with RXJ0911+0551 there is significant overlap in
constraints from detections and non-detections for $\Delta r$, something that
is absent in the main sample (see Figure 8). The two regions overlap because
the parameter that differentiates between detections and non-detections is not
one of the plotted parameters: it is the azimuthal angle of the SMBH with
respect to the galaxy center. Since the maximum of the arrival time surface,
$\Psi_{\rm Fermat}$, given the large external shear and, and the positions of
the sources in the lens macro-models, is significantly displaced from the
center of the galaxy, it is the azimuthal angle of the SMBH, rather than its
radial offset distance, that determines whether that SMBH will produce extra
images. Therefore, RXJ0911+0551 does not provide constraints on $\Delta r$
significantly different from VODKA regardless of a central image detection or
non-detection.

The non-detection posterior PDF for $r_{c}$ has a mode and 95% credible
interval of 0.485 kpc and $0.084<r_{c}<0.868$ kpc, respectively. The mode is
very similar to that of the detection posterior distribution, however, the
credible interval is larger in this case as a non-detection does not rule out
small $r_{c}$.

For the unlikely case of a detection in RXJ0911+0551, we assume the central
image is detected with the same properties as the ones in our main sample,
with $f_{o}$ equivalent to 9 magnitudes fainter than the estimated brightest
quad image, $\vec{r}_{o}$ at the center of the lens macro-model, $\sigma_{r}$
defined by the HST astrometric precision, and $\sigma_{f}=0.4f_{o}$. The
detection posterior PDF for $\Delta r$ is very similar to the VODKA prior,
with mode and 95% credible interval of 0.184 kpc and $0.022<\Delta r<0.389$
kpc, respectively. A detection of the central image in RXJ0911+0551 does not
improve SMBH offset constraints from VODKA. For core size $r_{c}$ and $b$, a
detection provides slightly tighter constraints. The prior distribution of
$r_{c}$ is skewed to larger values than the rest of the QSO samples. Likewise,
for a detection, the posterior PDF has a mode of 0.310 kpc and a 95% credible
interval of of $0.119<r_{c}<0.641$. A central image detection in this QSO
would lead to a tighter constraint on $r_{c}$. It is worth noting that there
is no measured velocity dispersion for RXJ0911+0551, and we have adopted a
$\sigma$ measurement equal to the average of the measured $\sigma$ in the main
sample. The constraints on $b$ with this QSO, therefore, could change
significantly with a velocity dispersion measurement.

For both $b$ and $r_{c}$, we see that the detection posterior PDFs yield
slightly tighter constraints than the prior, while the non-detection case does
not. Furthermore, both cases have very similar estimates, something also
observed with $\Delta r$. In general, we see that RXJ0911+0551 does not
provide any strong constraints on $\Delta r$, $r_{c}$, and $b$. Likewise, a
detection or non-detection does not significantly improve the posterior PDFs.
We summarize our constraints for RXJ0911+0551 in Figure 15.

Given that the wide radial distribution of images in this lens was the primary
cause for the large external shear, which lead to the lack of constraints on
$\Delta r$, we can tentatively conclude that quads with a very wide radial
image distribution and/or large external shear are unlikely to provide any
constraints on $\Delta r$, $r_{c}$, and $b$, regardless of central image
detection or non-detection.