# Drewnowski’s index to measure lifespan variation: Revisiting the Gini
coefficient of the life table

José Manuel Aburto  Ugofilippo Basellini  Annette Baudisch  Francisco
Villavicencio  Leverhulme Centre for Demographic Science, Department of
Sociology and Nuffield College, University of Oxford, Oxford OX1 1JD, UK.
Interdisciplinary Centre on Population Dynamics, University of Southern
Denmark, 5230 Odense, Denmark.  Laboratory of Digital and Computational
Demography, Max Planck Institute for Demographic Research, 18057 Rostock,
Germany.  Mortality, Health and Epidemiology Unit, Institut national d’études
démographiques (INED), 93322 Aubervilliers, France  Centre for Demographic
Studies (CED), Universitat Autònoma de Barcelona, 08193 Bellaterra, Spain.
Department of International Health, Johns Hopkins Bloomberg School of Public
Health, Baltimore, MD 21205, USA.

###### Abstract

The Gini coefficient of the life table is a concentration index that provides
information on lifespan variation. Originally proposed by economists to
measure income and wealth inequalities, it has been widely used in population
studies to investigate variation in ages at death. We focus on a complementary
indicator, Drewnowski’s index, which is as a measure of equality. We study its
mathematical properties and analyze how changes over time relate to changes in
life expectancy. Further, we identify the threshold age below which mortality
improvements are translated into decreasing lifespan variation and above which
these improvements translate into increasing lifespan inequality. We
illustrate our theoretical findings simulating scenarios of mortality
improvement in the Gompertz model. Our experiments demonstrate how
Drewnowski’s index can serve as an indicator of the shape of mortality
patterns. These properties, along with our analytical findings, support
studying lifespan variation alongside life expectancy trends in multiple
species.

###### keywords:

concentration index , Gompertz , life expectancy , lifespan inequality ,
mortality , threshold age

††journal: Theoretical Population Biology

##  1 Introduction

The life table is an essential tool in mortality studies. It represents the
current mortality experience of a population and it is usually summarized by
life expectancy at birth $(e_{o})$. Life expectancy at birth is the average
years a synthetic cohort of newborns is expected to live if individuals were
to experience the current mortality conditions throughout their lives. It is,
however, an average indicator that masks variability in ages at death, which
can be substantial. This variability is often referred to as lifespan
variation or lifespan inequality, and has received increasing attention over
the past two decades (see, to mention a few, Wilmoth and Horiuchi 1999;
Edwards and Tuljapurkar 2005; Smits and Monden 2009; Baudisch 2011; Vaupel et
al. 2011; Fernández and Beltrán-Sánchez 2015; Colchero et al. 2016; Ebeling et
al. 2018; van Raalte et al. 2018; Permanyer and Scholl 2019; Aburto et al.
2020; Vaupel et al. 2021). Lifespan variation reveals the uncertainty about
the eventual age at death at the individual level, and measures how evenly
mortality conditions are shared at the population level. There exist several
indicators to measure lifespan variation (for an overview, see Shkolnikov et
al. 2003; van Raalte and Caswell 2013), such as the entropy of the life table
(Leser, 1955; Keyfitz, 1977; Demetrius, 1978), the standard deviation or
variance of the age-at-death distribution (as applied in Tuljapurkar and
Edwards 2011), the coefficient of variation (as applied in Aburto and van
Raalte 2018; Aburto, Wensink et al. 2018), years of life lost $(e^{\dagger})$
(Goldman and Lord, 1986; Vaupel, 1986; Hakkert, 1987; Vaupel and Canudas Romo,
2003), or the Gini coefficient (Hanada, 1983).

Here we study in greater detail the Gini coefficient of the life table $(G)$
and its complementary measure, Drewnowski’s index $(D)$ (Drewnowski, 1982;
Hanada, 1983), from a formal demographic perspective. We additionally aim to
understand how changes in age-specific mortality underpin trends in lifespan
variation. We focus on how changes over time in $D$ relate to changes in
$e_{o}$, and highlight a new measure of absolute variation related to
perturbation theory, named $\vartheta$. We provide the mathematical foundation
of how Drewnowski’s index (and analogously the Gini coefficient) evolves over
time, and give analytical formulae to find the threshold age below which
mortality improvements are translated into decreasing lifespan variation and
above which these improvements translate into increasing lifespan inequality.

##  2 The Gini coefficient and Drewnowski’s index

The Gini coefficient is one of the most popular indices employed in social
sciences to measure concentration in the distribution of a non-negative random
variable (Gini, 1912, 1914). Originally proposed by economists to measure
income or wealth inequality, this coefficient has been applied in demography
and survival analysis to investigate within-group inequality in terms of ages
at death (see, for instance, Hanada 1983; Shkolnikov et al. 2003; Bonetti et
al. 2009; Gigliarano et al. 2017; Barthold Jones et al. 2018; Diaz et al.
2018; Basellini and Camarda 2019; Aburto et al. 2020).

###  2.1 Definition

As thoroughly discussed by Yitzhaki and Schechtman (2013), there are several
equivalent ways to define the Gini coefficient. Let $X$ be a non-negative
random variable with probability density function $f(x)$ and expected value
$\operatorname{\mathbb{E}}[X]$, a common definition is

$G=\frac{1}{2\operatorname{\mathbb{E}}[X]}\int_{0}^{\infty}\int_{0}^{\infty}|x_{1}-x_{2}|\,f(x_{1})\,f(x_{2})\,dx_{1}\,dx_{2}\;.$

Accordingly, if $X$ is a random variable of the ages at death in a population,
the Gini coefficient expresses the average of absolute differences in
individual lifespans relative to the mean length of life
$\operatorname{\mathbb{E}}[X]$.

Michetti and Dall’Aglio (1957), and later Hanada (1983), suggested a re-
formulation of the Gini coefficient in terms of the life table functions,
given by

$G=1-\frac{\int_{0}^{\infty}\ell(x,t)^{2}\,dx}{\int_{0}^{\infty}\ell(x,t)\,dx}=1-\frac{\vartheta}{e_{o}}\;,$
(1)

where $\ell(x,t)$ is the life table survival function at time $t$,
$e_{o}=\int_{0}^{\infty}\ell(x,t)\,dx$ the life expectancy at birth at time
$t$, and $\vartheta=\int_{0}^{\infty}\ell(x,t)^{2}\,dx$ is the resulting life
expectancy at birth of doubling the hazard at all ages. Barthold Jones et al.
(2018) interpret $\vartheta$ as a measure of _shared life expectancy_ , that
is, the average time two newborns at time $t$ are expected to survive
together. For the purposes of this article, the definition of the Gini
coefficient in (1) will be used in the following.

###  2.2 Main properties

The Gini coefficient takes values between 0 and 1, and can be interpreted as a
_measure of inequality_. A value of 0 denotes equality in ages at death, i.e.
when every individual in the population has the exact same length of life. The
index increases approaching 1 as lifespans become more spread and unequal in
the population. This makes the interpretation clear and intuitive: higher
values correspond to greater within-group inequality in ages at death.

An additional attractive feature of the Gini coefficient is that it fulfills
three important properties for inequality indices (Sen, 1973; Anand, 1983):
(i) it does not change if the number of individuals at each age at death is
changed by the same proportion (_population-size independence_); (ii) it does
not change if each individual lifespan is changed by the same proportion
(_scale independence_); and (iii) it decreases if years of life are
transferred from a longer to a shorter lived individual (_Pigou-Dalton
condition_). Note that property (i) enables straightforward comparison between
populations, including comparisons between different species (Wrycza et al.,
2015). Furthermore, the coefficient is not too sensitive to redistributions at
early ages of life, and it well reflects changes at adult ages (Shkolnikov et
al., 2003).

Being bounded between 0 and 1, the Gini coefficient can be readily transformed
from a _measure of inequality_ into a _measure of equality_ of lifespans. From
(1), a complementary measure of lifespan equality immediately derives, the
Drewnowski’s index, defined as

$D=1-G=\frac{\vartheta}{e_{o}}=\frac{\int_{0}^{\infty}\ell(x,t)^{2}\,dx}{\int_{0}^{\infty}\ell(x,t)\,dx}\;,$
(2)

and sharing all the important properties of $G$. According to Hanada (1983),
it was first proposed by Jan Drewnowski on a working group on health
indicators at the World Health Organization in the early 1980s (Drewnowski,
1982).

##  3 Changes over time in Drewnowski’s index

In order to analyze changes over time in Drewnowski’s index—or, equivalently,
the Gini coefficient—we aim to find an analytical expression for the time
derivative of $D$. In the following, a dot over a function will denote the
partial derivative with respect to time, but variable $t$ will be omitted for
simplicity.

###  3.1 Relative derivative of $D$

######  Proposition 1.

Let $D=\vartheta\,/\,e_{o}$ be Drewnowski’s index, where
$\vartheta=\int_{0}^{\infty}\ell(x)^{2}\,dx$,
$e_{o}=\int_{0}^{\infty}\ell(x)\,dx$ is the life expectancy at birth, and
$\ell(x)$ the probability of surviving from birth to age $x$. Then, relative
changes over time in $D$ are given by

$\frac{\dot{D}}{D}=\frac{\dot{\vartheta}}{\vartheta}-\frac{\dot{e}_{o}}{e_{o}}\;.$
(3)

###### Proof.

Note that $D=\vartheta\,/\,e_{o}$ implies that $D\,e_{o}-\vartheta=0$.
Differentiating with respect to time yields

$\displaystyle\dot{D}\,e_{o}+D\,\dot{e}_{o}-\dot{\vartheta}=0\;.$

Solving for $\dot{D}$ and dividing both sides by $D$, we get (3). ∎

Equation (3) decomposes relative changes in $D$ into relative changes of the
shared life expectancy between two individuals $\vartheta$, and relative
changes in the life expectancy at birth $e_{o}$.

###  3.2 Time derivatives of $e_{o}$ and $\vartheta$

Vaupel and Canudas Romo (2003) showed that changes over time in life
expectancy at birth are a weighted average of the total rates of mortality
improvements, expressed as

$\dot{e}_{o}=\int_{0}^{\infty}\rho(x)\,w(x)\,dx\;.$ (4)

Function $\rho(x)=-\dot{\mu}(x)\,/\,\mu(x)$ stands for the age-specific rates
of mortality improvement, where $\mu(x)$ is the force of mortality (hazard
rate) at age $x$. The weights $w(x)=\mu(x)\,\ell(x)\,e(x)$ are a measure of
the importance of death at age $x$, where
$e(x)=\int_{x}^{\infty}\ell(a)\,da\,/\,\ell(x)$ is the remaining life
expectancy at age $x$. Following a similar approach, we aim to express the
time derivative of $\vartheta$ as a weighted average of mortality
improvements, but with different weights.

######  Definition 1.

Let $\ell(x)$ be the probability of surviving from birth to age $x$. A measure
of lifespan equality at age $x$ is Drewnowski’s index conditional upon
survival up to age $x$, defined as

$D(x)=\frac{1}{\ell(x)}\,\frac{\int_{x}^{\infty}\ell(a)^{2}\,da}{\int_{x}^{\infty}\ell(a)\,da}\;.$
(5)

######  Proposition 2.

Let $\vartheta=\int_{0}^{\infty}\ell(x)^{2}\,dx$, where $\ell(x)$ is the
probability of surviving from birth to age $x$. Then, its partial derivative
with respect to time can be expressed as

$\dot{\vartheta}=\int_{0}^{\infty}\rho(x)\,w(x)\,2\,\ell(x)\,D(x)\,dx\;,$ (6)

where $\rho(x)$ are the age-specific rates of mortality improvement, $w(x)$
the same weights defined in (4), and $D(x)$ as defined in (5).

###### Proof.

Applying the chain rule, the derivative of $\vartheta$ with respect to time is
simply

$\dot{\vartheta}=\int_{0}^{\infty}2\,\ell(x)\,\dot{\ell}(x)\,dx\;.$

Using that $\dot{\ell}(x)=-\ell(x)\int_{0}^{x}\dot{\mu}(a)\,da$, and reversing
the order of integration, we get

$\begin{split}\dot{\vartheta}&=-2\int_{0}^{\infty}\ell(x)^{2}\,\int_{0}^{x}\dot{\mu}(a)\,da\,dx=-2\int_{0}^{\infty}\dot{\mu}(a)\int_{a}^{\infty}\ell(x)^{2}\,dx\,da\\\
&=2\int_{0}^{\infty}\rho(x)\,\mu(x)\,\ell(x)\,e(x)\,\frac{\int_{x}^{\infty}\ell(a)^{2}\,da}{\int_{x}^{\infty}\ell(a)\,da}\,dx\\\
&=\int_{0}^{\infty}\rho(x)\,w(x)\,2\,\ell(x)\,D(x)\,dx\;,\end{split}$

where $w(x)=\mu(x)\,\ell(x)\,e(x)$, which proves (6). ∎

###  3.3 Changes over time in $D$ in terms of mortality improvements

Equations (4) and (6) enable us to express changes over time in $D$ in terms
of mortality improvements. Using (2) and subsequently replacing (4) and (6) in
(3) yields

$\displaystyle\dot{D}$
$\displaystyle=D\,\left(\frac{\dot{\vartheta}}{\vartheta}-\frac{\dot{e}_{o}}{e_{o}}\right)=\frac{\dot{\vartheta}-\dot{e}_{o}D}{e_{o}}$
$\displaystyle=\int_{0}^{\infty}\rho(x)\,w(x)\,\frac{2\,\ell(x)\,D(x)-D}{e_{o}}\,dx$
$\displaystyle=\int_{0}^{\infty}\rho(x)\,w(x)\,W(x)\,dx\;.$ (7)

This result shows that changes over time in $D$ (and analogously in $G$) are a
total average of mortality improvements weighted by $w(x)\,W(x)$, where
$w(x)=\mu(x)\,\ell(x)\,e(x)$ are the same weights as in (4) and

$W(x)=\frac{2\,\ell(x)\,D(x)-D}{e_{o}}\;.$ (8)

##  4 The threshold age

###  4.1 Positive and negative contributions to lifespan equality

Because Drewnowski’s index is a measure of equality, $\dot{D}>0$ indicates
that lifespan variation decreases over time, whereas $\dot{D}<0$ implies that
lifespan variation increases over time. Equation (3.3) can then be used to
analyze the existence of a threshold age that separates _positive_ from
_negative_ contributions to lifespan equality as a result of mortality
improvements.

Note that in the assumption that mortality improvements occur at all ages,
$\rho(x)=-\dot{\mu}(x)\,/\,\mu(x)>0$ is a strictly positive function.
Therefore, from (3.3),

  1. 1.

Those ages $x$ for which $w(x)\,W(x)>0$ will contribute _positively_ to
Drewnowski’s index $D$ and increase lifespan equality;

  2. 2.

Those ages $x$ for which $w(x)\,W(x)<0$ will contribute _negatively_ to
Drewnowski’s index $D$ and favor lifespan inequality;

  3. 3.

Those ages $x$ for which $w(x)\,W(x)=0$ will have no effect on the variation
over time of $D$.

Any existing threshold age that separates positive from negative contributions
to lifespan equality will occur whenever $w(x)\,W(x)=0$. Since $\mu(x)$,
$\ell(x)$, and $e(x)$ are all positive functions, so are $w(x)$ and $e_{o}$.
Hence,

$w(x)\,W(x)=0\quad\Longleftrightarrow\quad 2\,\ell(x)\,D(x)-D=0\;.$ (9)

###  4.2 Existence and uniqueness of the threshold age

By means of the following two propositions and one theorem, we aim to prove
that in a scenario in which mortality improvements occur at all ages and
$\rho(x)>0$ for all $x\geq 0$, there exists a unique threshold age $a^{D}$
that separates positive from negative contributions to lifespan equality
(measured by $D$) as a result of those improvements.

######  Remark.

Following (2), Drewnowski’s index $D$ is bounded between 0 and 1, reaching a
value of 1 at complete equality in the ages at death within a population. A
score of 0 would express maximum inequality in the ages at death. By
definition, however, this value can never be reached:

$D=0\,\Longleftrightarrow\,\frac{\int_{0}^{\infty}\ell(x)^{2}\,dx}{\int_{0}^{\infty}\ell(x)\,dx}=0\,\Longleftrightarrow\,\int_{0}^{\infty}\ell(x)^{2}\,dx=0\,\Longleftrightarrow\,\ell(x)=0$
(10)

for all ages $x\geq 0$. This implies that the denominator in (10) equals 0
because $\ell(x)\geq 0$ is always positive and, consequently, $D$ would be
undefined. Hence, $0<D\leq 1$.

######  Proposition 3.

Let $\ell(x)$ be the probability of surviving from birth to age $x$, $D$
Drewnowski’s index as defined in (2), and $D(x)$ as defined in (5). Define the
function $g(x):=2\,\ell(x)\,D(x)-D$. Then, there exists at least one age
$a^{D}$ such that $g(a^{D})=0$.

###### Proof.

At age $x=0$,

$g(0)=2\,\ell(0)\,D(0)-D=2\,D-D=D>0$ (11)

by definition, since $0<D\leq 1$.

When ages become arbitrarily large,

$\lim_{x\to\infty}g(x)=\lim_{x\to\infty}(2\,\ell(x)\,D(x)-D)=2\lim_{x\to\infty}\ell(x)\,D(x)-D\;,$

which only depends on the behavior of $\ell(x)\,D(x)$. Because
$\ell(x)\in[0,1]$ for all ages $x\geq 0$, we have that
$0\leq\ell(x)^{2}\leq\ell(x)$ and

$0\leq\lim_{x\to\infty}\int_{x}^{\infty}\ell(a)^{2}\,da\leq\lim_{x\to\infty}\int_{x}^{\infty}\ell(a)\,da=\lim_{x\to\infty}e(x)\,\ell(x)=0\;,$

where $e(x)$ is the remaining life expectancy at age $x$. This proves that
both integrals tend to 0 as $x$ approaches $\infty$. Consequently,

$\lim_{x\to\infty}\ell(x)\,D(x)=\lim_{x\to\infty}\frac{\int_{x}^{\infty}\ell(a)^{2}\,da}{\int_{x}^{\infty}\ell(a)\,da}$

is indeterminate, but applying L’Hôpital’s rule, we get

$\lim_{x\to\infty}\ell(x)\,D(x)=\lim_{x\to\infty}\frac{\frac{\partial}{\partial
x}\,\int_{x}^{\infty}\ell(a)^{2}\,da}{\frac{\partial}{\partial
x}\,\int_{x}^{\infty}\ell(a)\,da}=\lim_{x\to\infty}\frac{-\ell(x)^{2}}{-\ell(x)}=\lim_{x\to\infty}\ell(x)=0\;.$

As a result,

$\lim_{x\to\infty}g(x)=2\lim_{x\to\infty}\ell(x)\,D(x)-D=0-D<0\;.$ (12)

Finally, using (11) and (12), on a continuous framework the intermediate value
theorem guarantees the existence of at least one positive age $a^{D}$ at which
$g(a^{D})=0$. ∎

######  Proposition 4.

Let $\ell(x)$ be the probability of surviving from birth to age $x$, $D$
Drewnowski’s index as defined in (2), and $D(x)$ as defined in (5). Then,
$g(x):=2\,\ell(x)\,D(x)-D$ is a strictly decreasing function.

###### Proof.

In order to demonstrate that $g(x)$ is a strictly decreasing function it
suffices to show that its first derivative is negative for all ages $x\geq 0$.
Note that since $D$ does not depend on age,

$\frac{\partial}{\partial
x}\,g(x)<0\quad\Longleftrightarrow\quad\frac{\partial}{\partial
x}\,\big{(}\ell(x)\,D(x)\big{)}<0\;.$

Applying the quotient rule together with the fundamental theorem of calculus,
we get

$\begin{split}\frac{\partial}{\partial
x}\,\big{(}\ell(x)\,D(x)\big{)}&=\frac{\partial}{\partial
x}\,\left(\frac{\int_{x}^{\infty}\ell(a)^{2}\,da}{\int_{x}^{\infty}\ell(a)\,da}\right)\\\
&=\frac{\int_{x}^{\infty}\ell(a)\,da\,\frac{\partial}{\partial
x}\,\left(\int_{x}^{\infty}\ell(a)^{2}\,da\right)-\int_{x}^{\infty}\ell(a)^{2}\,da\,\frac{\partial}{\partial
x}\,\left(\int_{x}^{\infty}\ell(a)\,da\right)}{\left(\int_{x}^{\infty}\ell(a)\,da\right)^{2}}\\\
&=\frac{\int_{x}^{\infty}\ell(a)\,da\,\left(-\ell(x)^{2}\right)-\int_{x}^{\infty}\ell(a)^{2}\,da\,\left(-\ell(x)\right)}{\left(\int_{x}^{\infty}\ell(a)\,da\right)^{2}}\,.\end{split}$

Hence,

$\begin{split}\frac{\partial}{\partial
x}\,g(x)<0&\,\Longleftrightarrow\,\ell(x)\int_{x}^{\infty}\ell(a)^{2}\,da-\ell(x)^{2}\int_{x}^{\infty}\ell(a)\,da<0\\\
&\,\Longleftrightarrow\,\frac{1}{\ell(x)^{2}}\int_{x}^{\infty}\ell(a)^{2}\,da<\frac{1}{\ell(x)}\int_{x}^{\infty}\ell(a)\,da\;.\end{split}$

Note that $\ell(x)=\exp\left[-\int_{0}^{x}\mu(a)\,da\right]$ for a given age-
specific hazard function $\mu(x)$. Therefore,
$\ell(x)^{2}=\exp\left[-\int_{0}^{x}2\,\mu(a)\,da\right]$ can be interpreted
as the survival schedule with doubling hazard $2\,\mu(x)$ at all ages $x\geq
0$. We can then define
$\tilde{e}(x)=\int_{x}^{\infty}\ell(a)^{2}\,da\,/\,\ell(x)^{2}$ as the
remaining life expectancy at age $x$ of a population with survival schedule
$\ell(x)^{2}$ and age-specific force of mortality $2\,\mu(x)$. Accordingly,

$\frac{\partial}{\partial
x}\,g(x)<0\,\Longleftrightarrow\,\frac{1}{\ell(x)^{2}}\int_{x}^{\infty}\ell(a)^{2}\,da<\frac{1}{\ell(x)}\int_{x}^{\infty}\ell(a)\,da\,\Longleftrightarrow\,\tilde{e}(x)<e(x)$

for all $x\geq 0$, which holds true since doubling the hazard corresponds to a
lower remaining life expectancy, in the reasonable assumption that $\mu(x)>0$
for all ages. ∎

######  Theorem.

Let $D=\vartheta\,/\,e_{o}$ be Drewnowski’s index, where
$\vartheta=\int_{0}^{\infty}\ell(x)^{2}\,dx$,
$e_{o}=\int_{0}^{\infty}\ell(x)\,dx$ is the life expectancy at birth, and
$\ell(x)$ the probability of surviving from birth to age $x$. Assume mortality
improvements over time occur at all ages. Then, there exists a unique
threshold age $a^{D}$ that separates positive from negative contributions to
lifespan equality, measured by $D$, as a result of those improvements.

###### Proof.

Following (3.3), changes over time in $D$ can be expressed as a weighted
average of mortality improvements, given by

$\dot{D}=\int_{0}^{\infty}\rho(x)\,w(x)\,W(x)\,dx\;,$

where $\rho(x)$ are the age-specific rates of mortality improvement over time,
and $w(x)\,W(x)$ the weights. By assumption, $\rho(x)>0$ for all ages $x\geq
0$. Therefore, any threshold age that separates positive from negative
contributions to lifespan equality as a result of mortality improvements will
occur whenever $w(x)\,W(x)=0$. From (9),

$w(x)\,W(x)=0\quad\Longleftrightarrow\quad 2\,\ell(x)\,D(x)-D=0\;,$

where $D(x)$ is as defined in (5). Proposition 3 guarantees the existence of
at least one positive age $a^{D}$ such that
$2\,\ell\left(a^{D}\right)\,D\left(a^{D}\right)-D=0$. In addition, from
Proposition 4 the function $g(x):=2\,\ell(x)\,D(x)-D$ is strictly decreasing.
Hence, assuming continuity, $g(x)$ is a one-to-one function and therefore the
threshold age $a^{D}$ is unique. ∎

######  Corollary.

Let $G$ be the Gini coefficient as defined in (1). Provided that $G=1-D$,
following (3.3)

$\dot{G}=-\dot{D}=-\int_{0}^{\infty}\rho(x)\,w(x)\,W(x)\,dx\;.$

Hence, $G$ and $D$ have the same threshold age, which is unique in the
assumption that mortality improvements occur at all ages. Improvements below
the threshold age will reduce lifespan variation (lowering $G$, but enlarging
$D$), and improvements above will increase lifespan inequality (enlarging $G$,
but lowering $D$).

##  5 Application

We illustrate our theoretical findings simulating scenarios of mortality
improvement depending on level versus rate changes in the Gompertz mortality
model (Gompertz, 1825). This model captures an exponential change of the force
of mortality $\mu(x)=\alpha\cdot\mathrm{e}^{\beta\,x}$ over age $x$, where
$\alpha$ is the baseline level of mortality and $\beta$ the rate of aging. The
R code (R Core Team, 2021) to reproduce all the figures and results discussed
in the following will be publicly available upon acceptance.

###  5.1 Gompertz force of mortality with positive aging: level vs rate change

Our first scenario demonstrates how improvements in mortality affect our
outcome variables when progress results from reducing the level of mortality,
but not from slowing the rate of aging. Figure 1 depicts the Gompertz
mortality function on a log-scale (Panel A), the corresponding age-at-death
distribution (Panel B), and associated weights $W(x)$ (Eq. 8) (Panel C), for a
fixed rate of aging $\beta=0.1$ and changing levels of baseline mortality
$\alpha$. Values of Drewnowski’s index $D$ and corresponding threshold ages
$a^{D}$ are reported in Panels B) and C).

Figure 1: Gompertz mortality model with positive aging $(\beta>0)$ for
different levels of baseline mortality $\alpha$ and a fixed rate of aging
$\beta=0.1$

Lower levels of baseline mortality $\alpha$ translate into vertical downward
shifts in the log-hazard (Panel A). Correspondingly, lifespan variation falls
as deaths concentrate at older ages, as captured by larger values of $D$
(Panel B). Lower levels of mortality raise the threshold age $a^{D}$ at which
the weights $W(x)$ switch sign from positive to negative (Panel C). Positive
weights at some age imply that saving lives at that age increases equality,
i.e. enlarging $D$; negative weights imply that saving lives decreases
equality, i.e. lessening $D$. The threshold age $a^{D}$ marks the boundary
between these positive and negative effects of life saving on lifespan
equality and is indicated by the dashed vertical lines. Panel C) shows that
under the assumption of mortality improvements over all ages, this threshold
age is unique and that the magnitude of the weights $W(x)$ over age changes
depending on the different shapes of the distribution of deaths, supporting
our theoretical results.

Complementing these findings, Fig. 2 analyzes how slowing the rate of aging
$\beta$ affects lifespan variation. For a fixed baseline level of mortality
$\alpha=\mathrm{e}^{-5}\approx 0.0067$, Panel A) illustrates how lower values
of $\beta$ translate, as expected, into a slower rise in the log-hazard. In
contrast to the first scenario, Panel B) reveals that slowing the rate of
aging $\beta$ reduces $D$ (lighter colors) and widens the distribution of ages
at death. That is, extending lifespan by slowing the rate of aging increases
lifespan inequality, which is consistent with previous research by Colchero et
al. (2021). The threshold age $a^{D}$ increases as lifespan equality
decreases. While life saving has opposite effects on lifespan variability and
the level of the threshold age in the two scenarios, the general dynamics of
the weights $W(x)$ remain similar. Saving lives before the threshold age
reduces inequality; saving lives later increases inequality. The threshold age
$a^{D}$ remains unique, independent of specific values of the mortality
parameters, further supporting our theoretical predictions.

Figure 2: Gompertz mortality model with positive aging for different rates of
aging $\beta>0$ and a fixed level of mortality $\alpha=\mathrm{e}^{-5}$

###  5.2 Zero rate of aging

Figure 3 captures a scenario in which the rate of aging equals zero. Under
varying but constant levels of mortality (Panel A), lifespan variation remains
unchanged at $D=0.5$ and the largest number of deaths occur at the earliest
ages (Panel B). The lower the level of mortality, the longer the length of
life, and the later the threshold age (Panel C).

This can also be verified analytically. When $\beta=0$, the survival function
for the Gompertz mortality model simplifies to
$\ell(x)=\mathrm{e}^{-\alpha\,x}$. Then, life expectancy as the integral of
the survivorship over all ages becomes $1/\alpha$, which is the inverse of the
force of mortality. Following Definition 1, in this scenario the conditional
Drewnowski index $D(x)$ is

$D(x)=\frac{1}{\mathrm{e}^{-\alpha\,x}}\,\frac{\mathrm{e}^{-2\,\alpha\,x}\,/\,2\,\alpha}{\mathrm{e}^{-\alpha\,x}\,/\,\alpha}=0.5\;.$

In other words, when mortality is constant over age, $D(x)=0.5$ for all ages
and it is independent of the baseline mortality level $\alpha$. Accordingly,
from Proposition 3

$g(x)=0\Longleftrightarrow
2\,\mathrm{e}^{-\alpha\,x}\,0.5-0.5=0\Longleftrightarrow
x=\ln(2)\,/\,\alpha\;,$

which implies that the threshold age is $a^{D}=\ln(2)\,\,/\alpha$ and
$\ell(a^{D})=\mathrm{e}^{-\alpha\,a^{D}}=0.5$. This means that, when mortality
is constant, the threshold age and the median age at death are the same.

Figure 3: Gompertz mortality model with zero rate of aging $(\beta=0)$ and
different levels of baseline mortality $\alpha$

###  5.3 Negative rate of aging

To complete the analyses, Figs. 4 and 5 capture scenarios of decreasing
mortality, which are consistently marked by values of $D$ below one half
(Panels B). Different to the scenario with positive aging rates, here effects
of both level and rate changes point in the same direction: higher mortality
results in less lifespan variation (higher values of $D$ in Panels B). Even
though the value of the Drewnowski’s index varies substantially (from $0.03$
to $0.48$), the threshold age barely changes and remains close to the earliest
age (Panels C). This is because in these scenarios most individuals die young,
and life expectancy is extremely short for all cases.

The two negative aging scenarios substantially differ, however, in the maximum
age for the survivors (not shown). This reflects the exceeding disparity among
the majority of individuals who die upon—or shortly after—birth, and those few
surviving to manifold higher ages as they benefit from the improvements in
mortality over age. These differences are also reflected in the density
functions (Panels B). While varying the rate of aging leaves no
distinguishable effect on the density function (Fig. 5, Panel B), varying the
baseline level of mortality $\alpha$ distinctly affects the shape of the
density function with lower mortality levels yielding an increasing
rectangularization.

Figure 4: Gompertz mortality model with negative rate of aging $(\beta<0)$ for
different levels of baseline mortality $\alpha$ and a fixed rate of aging
$\beta=-0.1$ Figure 5: Gompertz mortality model with negative rate of aging
($\beta<0$) for different rates of aging $\beta$ and a fixed level of
mortality $\alpha=\mathrm{e}^{-0.5}$

##  6 Conclusion

In this article we uncovered how age-specific patterns of mortality underpin
trends in lifespan variation as measured by Drewnowski’s index $D$, a
complementary indicator to the Gini coefficient of the life table, by means of
formal demography. We contribute to the literature by disentangling how
changes in age patterns of mortality affect lifespan variation and provide an
analytical proof for the existence of a threshold age $a^{D}$ below which
mortality improvements are translated into increasing lifespan equality and
above which these improvements translate into increasing lifespan inequality.
Previous research determined such age for the life table entropy (Aburto et
al., 2019), the variance (Gillespie et al., 2014) and $e^{\dagger}$ (Zhang and
Vaupel, 2009).

We test and illustrate our results under multiple scenarios of mortality
changes with age, including positive, zero and negative rates of aging. This
is important because shapes of mortality patterns across the tree of life vary
substantially (Jones, Scheuerlein et al. 2014). For example, while some
species such as humans exhibit increasing mortality over adult ages, other
species experience an unchanged force of mortality over age, such as _Hydra_
(Schaible, Scheuerlein et al. 2015), or even declining adult mortality, called
‘negative senescence’ (Vaupel et al., 2004; Baudisch, 2008; Cayuela et al.,
2019). Our experiments demonstrate how Drewnowski’s index $D$ can serve as an
indicator of the ‘shape’ of mortality patterns (Baudisch, 2011; Wrycza et al.,
2015), distinguishing between increasing mortality $(D>0.5)$, constant
mortality $(D=0.5)$, and decreasing mortality $(D<0.5)$ with age. These
properties, along with our analytical findings, support studying lifespan
variation alongside life expectancy trends in multiple species.

## Acknowledgements

This work was partially funded by the European Research Council (Grant
ERC-2019-AdG-884328) (J.M.A. and F.V.). J.M.A. also acknowledges funding from
the British Academy’s Newton International Fellowship (Grant NIFBA19/190679)
and the Leverhulme Trust Large Centre, and F.V. from the Juan de la Cierva
Programme, Spanish Ministry of Science and Innovation (Grant
IJC2019-039144-I).

## References

  * Aburto et al. (2019) Aburto, J.M., Alvarez, J.A., Villavicencio, F., Vaupel, J.W., 2019\.  The threshold age of the lifetable entropy.  Demographic Research 41, 83–102.
  * Aburto and van Raalte (2018) Aburto, J.M., van Raalte, A.A., 2018\.  Lifespan dispersion in times of life expectancy fluctuation: The case of Central and Eastern Europe.  Demography 55, 2071–2096.
  * Aburto et al. (2020) Aburto, J.M., Villavicencio, F., Basellini, U., Kjærgaard, S., Vaupel, J.W., 2020.  Dynamics of life expectancy and life span equality.  Proceedings of the National Academy of Sciences of USA 117, 5250–5259.
  * Aburto et al. (2018) Aburto, J.M., Wensink, M., van Raalte, A., Lindahl-Jacobsen, R., 2018\.  Potential gains in life expectancy by reducing inequality of lifespans in Denmark: An international comparison and cause-of-death analysis.  BMC Public Health 18, 831\.
  * Anand (1983) Anand, S., 1983.  Inequality and poverty in Malaysia: Measurement and decomposition.  Oxford University Press for the World Bank, New York, NY.
  * Barthold Jones et al. (2018) Barthold Jones, J.A., Lenart, A., Baudisch, A., 2018.  Complexity of the relationship between life expectancy and overlap of lifespans.  PLoS ONE 13, e0197985.
  * Basellini and Camarda (2019) Basellini, U., Camarda, C.G., 2019\.  Modelling and forecasting adult age-at-death distributions.  Population Studies 73, 119–138.
  * Baudisch (2008) Baudisch, A., 2008.  Inevitable aging? Contributions to evolutionary-demographic theory.  Springer, Berlin/Heidelberg.
  * Baudisch (2011) Baudisch, A., 2011.  The pace and shape of ageing.  Methods in Ecology and Evolution 2, 375–382.
  * Bonetti et al. (2009) Bonetti, M., Gigliarano, C., Muliere, P., 2009.  The Gini concentration test for survival data.  Lifetime Data Analysis 15, 493–518.
  * Cayuela et al. (2019) Cayuela, H., et al., 2019.  Life history and age-dependent mortality processes in tropical reptiles.  Biological Journal of the Linnean Society 128, 251–262.
  * Colchero et al. (2016) Colchero, F., et al., 2016.  The emergence of longevous populations.  Proceedings of the National Academy of Sciences of USA 113, E7681–E7690.
  * Colchero et al. (2021) Colchero, F., et al., 2021.  The long lives of primates and the ‘invariant rate of ageing’ hypothesis.  Nature Communications 12, 3666\.
  * Demetrius (1978) Demetrius, L., 1978.  Adaptive value, entropy and survivorship curves.  Nature 275, 213–214.
  * Diaz et al. (2018) Diaz, G., Debón, A., Giner-Bosch, V., 2018.  Mortality forecasting in Colombia from abridged life tables by sex.  Genus 74, 15.
  * Drewnowski (1982) Drewnowski, J., 1982.  Social indicators and their application for measuring the progress of development, in: Working Group on Indicators for Monitoring and Evaluation of Strategies for Health For All by the Year 2000, WHO Regional Office for the Western Pacific, Manila, Philippines, 2–9 March 1982. pp. 103–130.  URL: https://apps.who.int/iris/handle/10665/207170.
  * Ebeling et al. (2018) Ebeling, M., Rau, R., Baudisch, A., 2018.  Rectangularization of the survival curve reconsidered: The maximum inner rectangle approach.  Population Studies 72, 369–379.
  * Edwards and Tuljapurkar (2005) Edwards, R.D., Tuljapurkar, S., 2005\.  Inequality in life spans and a new perspective on mortality convergence across industrialized countries.  Population and Development Review 31, 645–674.
  * Fernández and Beltrán-Sánchez (2015) Fernández, O.E., Beltrán-Sánchez, H., 2015\.  The entropy of the life table: A reappraisal.  Theoretical Population Biology 104, 26–45.
  * Gigliarano et al. (2017) Gigliarano, C., Basellini, U., Bonetti, M., 2017.  Longevity and concentration in survival times: The log-scale-location family of failure time models.  Lifetime Data Analysis 23, 254–274.
  * Gillespie et al. (2014) Gillespie, D.O., Trotter, M.V., Tuljapurkar, S.D., 2014.  Divergence in age patterns of mortality change drives international divergence in lifespan inequality.  Demography 51, 1003–1017.
  * Gini (1912) Gini, C., 1912.  Variabilità e mutabilità: Contributi allo studio delle distribuzioni delle relazioni statistiche.  Studi Economico-Giuridici della Università di Cagliari 3.
  * Gini (1914) Gini, C., 1914.  Sulla misura della concentrazione e della variabilità dei caratteri.  Atti del Reale Istituto Veneto di Scienze, Lettere ed Arti 73, 1203–1248.
  * Goldman and Lord (1986) Goldman, N., Lord, G., 1986\.  A new look at entropy and the life table.  Demography 23, 275–282.
  * Gompertz (1825) Gompertz, B., 1825.  On the nature of the function expressive of the law of human mortality, and on a new mode of determining the value of life contingencies. In a letter to Francis Baily, Esq. F.R.S. &c.  Philosophical Transactions of the Royal Society of London 115, 513–583.
  * Hakkert (1987) Hakkert, R., 1987.  Life table transformations and inequality measures: Some noteworthy formal relationships.  Demography 24, 615–622.
  * Hanada (1983) Hanada, K., 1983.  A formula of Gini’s concentration ratio and its application to life tables.  Journal of Japanese Statistical Society 13, 95–98.
  * Jones et al. (2014) Jones, O.R., Scheuerlein, A., et al., 2014\.  Diversity of ageing across the tree of life.  Nature 505, 169–173.
  * Keyfitz (1977) Keyfitz, N., 1977.  What difference would it make if cancer were eradicated? An examination of the Tæuber paradox.  Demography 14, 411–418.
  * Leser (1955) Leser, C.E.V., 1955.  Variations in mortality and life expectation.  Population Studies 9, 67–71.
  * Michetti and Dall’Aglio (1957) Michetti, B., Dall’Aglio, G., 1957\.  La differenza semplice media.  Statistica 7, 159–255.
  * Permanyer and Scholl (2019) Permanyer, I., Scholl, N., 2019\.  Global trends in lifespan inequality: 1950–2015.  PLoS ONE 14, e0215742.
  * R Core Team (2021) R Core Team, 2021.  R: A Language and Environment for Statistical Computing.  R Foundation for Statistical Computing. Vienna, Austria.  URL: https://www.R-project.org/.
  * van Raalte and Caswell (2013) van Raalte, A.A., Caswell, H., 2013\.  Perturbation analysis of indices of lifespan variability.  Demography 50, 1615–1640.
  * van Raalte et al. (2018) van Raalte, A.A., Sasson, I., Martikainen, P., 2018.  The case for monitoring life-span inequality.  Science 362, 1002–1004.
  * Schaible et al. (2015) Schaible, R., Scheuerlein, A., Dańko, M.J., Gampe, J., Martínez, D.E., Vaupel, J.W., 2015\.  Constant mortality and fertility over age in _Hydra_.  Proceedings of the National Academy of Sciences of USA 112, 15701–15706.
  * Sen (1973) Sen, A.K., 1973.  On economic inequality.  Clarendon Press, Oxford.
  * Shkolnikov et al. (2003) Shkolnikov, V.M., Andreev, E.E., Begun, A.Z., 2003.  Gini coefficient as a life table function: Computation from discrete data, decomposition of differences and empirical examples.  Demographic Research 8, 305–358.
  * Smits and Monden (2009) Smits, J., Monden, C., 2009\.  Length of life inequality around the globe.  Social Science & Medicine 68, 1114–1123.
  * Tuljapurkar and Edwards (2011) Tuljapurkar, S., Edwards, R.D., 2011\.  Variance in death and its implications for modeling and forecasting mortality.  Demographic Research 24, 497–526.
  * Vaupel (1986) Vaupel, J.W., 1986.  How change in age-specific mortality affects life expectancy.  Population Studies 40, 147–157.
  * Vaupel et al. (2004) Vaupel, J.W., Baudisch, A., Dölling, M., Roach, D.A., Gampe, J., 2004.  The case for negative senescence.  Theoretical Population Biology 65, 339–351.
  * Vaupel and Canudas Romo (2003) Vaupel, J.W., Canudas Romo, V., 2003\.  Decomposing change in life expectancy: A bouquet of formulas in honor of Nathan Keyfitz’s 90th birthday.  Demography 40, 201–216.
  * Vaupel et al. (2021) Vaupel, J.W., Villavicencio, F., Bergeron-Boucher, M.P., 2021.  Demographic perspectives on the rise of longevity.  Proceedings of the National Academy of Sciences of USA 118, e2019536118.
  * Vaupel et al. (2011) Vaupel, J.W., Zhang, Z., van Raalte, A.A., 2011.  Life expectancy and disparity: An international comparison of life table data.  BMJ Open , bmjopen–2011–000128.
  * Wilmoth and Horiuchi (1999) Wilmoth, J.R., Horiuchi, S., 1999\.  Rectangularization revisited: Variability of age at death within human populations.  Demography 36, 475–495.
  * Wrycza et al. (2015) Wrycza, T.F., Missov, T.I., Baudisch, A., 2015.  Quantifying the shape of aging.  PLoS ONE 10, e0119163.
  * Yitzhaki and Schechtman (2013) Yitzhaki, S., Schechtman, E., 2013\.  The Gini Methodology: A primer on a statistical methodology.  Springer, New York, NY.
  * Zhang and Vaupel (2009) Zhang, Z., Vaupel, J.W., 2009\.  The age separating early deaths from late deaths.  Demographic Research 20, 721–730.


# Quasar Microlensing Statistics and Flux-Ratio Anomalies in Lens Models

E. Mediavilla Instituto de Astrofísica de Canarias, Vía Láctea S/N, La Laguna
38200, Tenerife, Spain Departamento de Astrofísica, Universidad de la Laguna,
La Laguna 38200, Tenerife, Spain J. Jiménez-Vicente Departamento de Física
Teórica y del Cosmos, Universidad de Granada, Campus de Fuentenueva, 18071
Granada, Spain Instituto Carlos I de Física Teórica y Computacional,
Universidad de Granada, 18071 Granada, Spain V. Motta European Southern
Observatory, Alonso de Córdova 3107, Vitacura, Santiago, Chile. Instituto de
Física y Astronomía, Facultad de Ciencias, Universidad de Valparaíso, Avda.
Gran Bretaña 1111, 2360102 Valparaíso, Chile

###### Abstract

Precise lens modeling is a critical step in time delay studies of multiply
imaged quasars, which are key for measuring some important cosmological
parameters (specially $H_{0}$). However, lens models (in particular those
semi-automatically generated) often show discrepancies with the observed flux-
ratios between the different quasar images. These flux-ratio anomalies are
usually explained through differential effects between images (mainly
microlensing) that alter the intrinsic magnification ratios predicted by the
models. To check this hypothesis, we collect direct measurements of
microlensing to obtain the histogram of microlensing magnifications. We
compare this histogram with recently published model flux-ratio anomalies and
conclude that they cannot be statistically explained by microlensing. The
average value of the model anomalies ($0.74\,$magnitudes) significantly
exceeds the mean impact of microlensing ($0.33\,$magnitudes). Moreover, the
histogram of model anomalies presents a significant tail with high anomalies
($|\Delta m|\geq 0.7$ magnitudes) which is completely unexpected from the
statistics of microlensing observations. Microlensing simulations neither
predict the high mean nor the fat tail of the histogram of model anomalies. We
perform several statistical tests which exclude that microlensing can explain
the observed flux-ratio anomalies (although Kolmogorov-Smirnov, which is less
sensitive to the tail of the distributions, is not always conclusive). Thus,
microlensing cannot statistically explain the bulk of flux-ratio anomalies,
and models may explore different alternatives to try to reduce them. In
particular, we propose to complement photometric observations with accurate
flux ratios of the broad emission lines obtained from integral field
spectroscopy to check and, ideally, constrain lens models.

Gravitational lensing: strong, models — gravitational lensing: micro

##  1 Introduction

Systems of multiple images of distant quasars formed by the gravitational
field of intervening galaxies are one of the most useful ”laboratories” in
astrophysics and cosmology, allowing to study the structure of the quasar
sources, the properties of matter in the lens galaxies, and the cosmological
parameters, among other applications. A necessary step in these studies is the
modeling of the lens, which needs to be very precise and robust, particularly
for cosmographic applications like the prediction of gravitational time delays
between the images, which may be used to solve the current tension in the
determination of the Hubble constant, $H_{0}$, from different methods (see,
e.g., Di Valentino et al. 2023 and references therein). On the other hand, the
number of observed systems which will need to be modeled is expected to
increase considerably in the near future, which will prevent a detailed
individual modeling, calling for (semi-)automated procedures (see, e.g. Shajib
et al. 2019, Schmidt et al. 2023 and references therein).

Common observable photometric/astrometric quantities of the lensed systems are
the positions of the lens and images, and the fluxes of the images. Models
must, therefore, take into account the structure of the lens (often including
secondary lenses), the structure of the source, and a careful modeling of the
point spread function (Koopmans et al. 2003, Suyu et al. 2010, Birrer et al.
2022). To avoid an unmanageable large number of unknowns, the lens mass
distribution is usually parametrized (e.g. as a power law or Navarro-Frenk-
White profile). Spectroscopic information is less frequently used in spite
that it may be crucial to break important degeneracies present in lens
modeling.

Photometry based lens models are usually much more constrained by the
astrometric observables (because they are accurate, with typical uncertainties
of a few milliarcseconds) than by fluxes. Moreover, in the case of quasar
images, broad- or narrow-band fluxes can be, in principle, affected by several
sources of uncertainty like intrinsic variability of the source combined with
time-delays, micro and millilensing, extinction, etc. (see Pooley et al. 2007,
Yonehara et al. 2008, Motta et al. 2012 and references therein) which make
them practically of no use in constraining the models. In fact, calculated
models very often present strong differences between their predicted flux
ratios and the observed ones (Witt et al. 1995, Mao & Schneider 1998, Chiba
2002, Metcalf & Madeau 2001, Dalal & Kochanek 2002, Schechter & Wambsganss
2002, Keeton 2002, Bradac et al. 2002, Metcalf & Zhao 2002, Moustakas &
Metcalf 2003, Metcalf & Amara 2012, Xu et al. 2009, 2015, Gilman et al. 2017).
As a statistically representative case, Shajib et al. (2019) find strong flux-
ratio anomalies, which they attribute mainly to microlensing, in a sample of
13 quadruple imaged quasars studied to devise a general framework to model
multiply imaged quasars (with the aim of processing the large number of
systems to be discovered in deep wide-field surveys like the Wide-Field
Infrared Survey Telescope, LSST or Euclid). If microlensing is the cause of
the anomalies, continuum flux-ratios in the visible are basically useless to
model lens systems, but if the impact of microlensing cannot explain the flux-
ratio departures from the predictions, it would make sense to explore possible
ways to improve the models (like the use of spectroscopic data). Microlensing
magnification of the quasar images can be directly (independently from lens
modeling) measured by using spectroscopic information of the quasar images.
Microlensing is size sensitive (the larger the size, the smaller the impact of
microlensing) and in the spectrum of each image we have information from
different regions in the quasar: the continuum comes from the tiny accretion
disk, which can be strongly affected by microlensing, while the broad emission
lines come from the relatively large broad line region, which is rather
insensitive to this effect (see, e.g, Wisotzki et al. 1993, Mediavilla et al.
2009, 2011 and references therein). Consequently, we can use the flux ratios
corresponding to the emission lines as zero microlensing baseline to measure
the impact of microlensing in the continuum flux ratios at a given epoch
(single epoch microlensing measurements).

Alternatively, microlensing can also be studied from photometric monitoring of
the images of a lensed quasar. Subtracting the (time delay corrected) light
curves of two images, we can obtain microlensing light curves. Owing to the
presence of extinction and to the relatively reduced extension of the
monitoring period we can not fully quantify the total amplitude of
microlensing at a given epoch, but the amplitudes of the peaks of microlensing
events should provide a conservative upper bound to single epoch microlensing.

The first objective of this paper is, therefore, to obtain the histogram of
observed flux-ratio anomalies induced by microlensing. In order to do so, we
calculate the experimental differential microlensing magnifications for a
sample of 44 measurements in 34 image pairs of 23 lens systems with available
spectroscopic information (collected by Esteban-Gutiérrez et al. 2022 from
Rojas et al. 2020, 2014, Motta et al. 2017, 2012, Jiménez-Vicente et al.,
2015, and Mediavilla et al. 2009). We then compare this experimental histogram
with the flux-ratio anomalies inferred from model predictions in two samples
of quadruple lens systems (Shajib et al. 2019, Schmidt et al. 2023), to
illustrate how the hypothesis that the flux-ratio anomalies are caused by
microlensing can be tested. The present work does not intend to make any
general statement about lens modeling. Instead, we just aim at providing some
tools to detect potential problems in some models and to, eventually, improve
them in some cases.

The paper is organized as follows. In §2 we present the histogram of observed
microlensing magnifications (based on spectroscopic data), and compare it with
the statistics of microlensing peak amplitudes (derived from microlensing
light curves) and with the predictions of microlensing simulations. In §3 we
collect histograms of microlensing model anomalies from the samples in Shajib
et al. (2019) and Schmidt et al. (2023) and compare them with the histograms
of observed microlensing magnifications and of microlensing peak amplitudes.
In §4 we discuss possible observational strategies to derive useful
constraints based on the flux ratios either to cross-check the models or to
improve them. In §5 we summarize the main conclusions. Finally, we devote an
Appendix to explore the relationship of flux-ratio anomalies with the
degeneracy of lens models with respect to the radial mass distribution.

##  2 Observed microlensing in multiple imaged quasars

###  2.1 Estimates from emission-lines

To directly measure the impact of microlensing in the images of lensed
quasars, we can take advantage of the sensitivity of microlensing to the size
of the source. Microlensing by a distribution of stars induces strong spatial
changes (”microlensing roughness”) in the otherwise uniform (smooth)
magnification at the source plane. If the size of the source is large enough,
the inhomogeneities of the magnification are spatially averaged and washed
out. The spatial scale of the magnification roughness is related to the
Einstein radius of the microlenses, which for a typical mass of 0.3$M_{\odot}$
and typical values for the redshifts of the lens ($z_{l}=0.5$) and the source
($z_{l}=2$) amounts to approximately $10\,$light-days. Consequently, the
impact of microlensing can be potentially high for the quasar continuum source
(a few light-days in size) but negligible for the Broad Line Region (BLR)
(with sizes above hundred light-days) (e.g., Jiménez-Vicente et al. 2022,
Guerras et al. 2013, Fian et al. 2018). Thus, we can use the broad lines
present in quasar spectra (in particular, the core of the lines) to determine
the zero microlensing baseline. For a pair of images of a lens system, we can
define the relative microlensing magnification between them as the continuum
ratio relative to the zero point defined by the emission line ratio. Expressed
in magnitudes we can write (see, e.g., Mediavilla et al. 2009),

$\Delta m_{ij}=(m_{i}-m_{j})^{cont}-(m_{i}-m_{j})^{line.},$ (1)

Using the sample of 44 microlensing measurements collected by Esteban-
Gutiérrez et al.111Collected from Rojas et al. 2020, 2014, Motta et al. 2017,
2012, Jiménez-Vicente et al., 2015, and Mediavilla et al. 2009 (2022) obtained
according to Eq. 1, we derive the histogram of (unsigned) microlensing
magnifications (shown in Figure 1). The average rest wavelength for which
these measured differential microlensing magnifications are estimated is
$\lambda\sim 1700$ Å (cf. Jiménez-Vicente et al., 2012). The histogram has a
mean of $\langle|\Delta m_{ij}|\rangle=0.33\pm 0.22$. This is, in fact, an
overestimate of the expected impact of microlensing, as all four single
measurements in Eq.1 are also affected by experimental uncertainties, which
will broaden the intrinsic histogram of microlensing magnifications. We have
estimated directly the mean error, $\langle\sigma_{\Delta m}\rangle$ from the
different data sources used: $0.13\pm 0.09$ (Mediavilla et al.222There is no
rms estimates for all the measurements. 2009), $0.11\pm 0.04$ (Motta et al.
2011), $0.23\pm 0.02$ (Rojas et al. 2014), $0.15\pm 0.13$ (Motta et al. 2017),
and $0.11\pm 0.07$ (Rojas et al. 2020). The weighted average of the errors is
$0.13$.

Notice that the experimental $\Delta m_{ij}$ also includes the differences in
flux arising from intrinsic quasar variability combined with the time delay
between images, which are supposed to be small, specially for quads. On the
contrary, as far as the continuum and emission lines are observed at close
wavelengths, we can assume that the $\Delta m_{ij}$, calculated using Eq.1,
are virtually free from extinction.

In the histogram of microlensing magnifications (Figure 1) doubles and quads
are mixed. If we separate both groups, we obtain $\langle|\Delta
m_{ij}|\rangle_{quads}=0.27\pm 0.22$ and $\langle|\Delta
m_{ij}|\rangle_{doubles}=0.40\pm 0.19$. The results indicate that doubles
exhibit slightly larger microlensing than quads (likely because doubles have
larger time delays, which combined with intrinsic variability may increase
flux-ratio anomalies; on the other hand in doubles one of the images is also
often close to the lens galaxy and, consequently, more prone to microlensing).
Then, if we restrict to quads, micro-lensing anomalies would be even slightly
smaller. Some of the microlensing measurements in our sample correspond to
different epochs of a same image pair. To avoid the impact of possible
covariance between repeated measurements we have substituted for each image
pair where more than one measurement were available, all the measurements by
its mean, finding neglectable differences in $\langle|\Delta m_{ij}|\rangle$
either for quads or doubles.

###  2.2 Comparison with peak amplitudes of microlensing light-curves

It is interesting to compare the values of microlensing magnifications
obtained using the emission lines as baseline, with the peak amplitudes of
microlensing light curves. In Figure 1, we also present the histogram of
microlensing peak amplitudes taken from Mediavilla et al. (2016). These
amplitudes use the ”flat” part of the light-curve before or after the
microlensing event as baseline. This is an idealization, because due to the
relatively high optical depth, quasar light curves can not generally be
described as isolated microlensing events/peaks over a flat baseline. In fact,
it is common that microlensing light-curves do not present a well defined flat
region. Moreover, the peak or part of the baseline can fall in one of the
(seasonal or incidental) gaps of the light-curves. For this reason, the peak
amplitude defined with respect to the left or right sides from the peak, can
be different in some cases. We have selected always the largest one.

There is a great similarity between the histograms of peak amplitudes and of
microlensing magnifications from emission lines, although an offset towards
larger values of the histogram of microlensing peak amplitudes would be
expected. The coincidence of the means of both histograms can result, in part,
from the above mentioned overestimate in the mean value of microlensing
inferred from the emission lines due to measurement uncertainties. Notice also
that intrinsic variability is contributing to the microlensing magnifications
from emission lines while it is not affecting to microlensing light curves
which are obtained subtracting time delay corrected light curves of two
images. Moreover, an underestimate in peak amplitudes can be produced because
the true microlensing zero-point may fall below the ”flat” regions of the
microlensing light-curves taken as baseline333In not long enough light curves,
caustic clustering can mimic a fake ”flat” region., or because the maximum is
in a seasonal gap.

Finally, note that for source sizes comparable to the Einstein radius of the
microlenses, lensed quasar images can be frequently engaged in microlensing
events with gentle slopes and broad peaks of relatively small amplitude. In
fact, taking into account the total monitoring time (310 years) of the light
curves of the ensemble from Mediavilla et al. (2016), the total number of
microlensing events detected in the ensemble (20) and the mean Einstein radius
crossing time scale444Calculated using Eq. 8 of Mediavilla et al. 2016 for a
typical lens system: $z_{lens}=0.5$, $z_{source}=2$,
$\sigma_{*}=200\rm\,km\,s^{-1}$, and $\sigma_{pec}(z=0.5)=638\rm\,km\,s^{-1}$.
of $9.4$ years, a 61% of the images will be engaged in a microlensing event at
any time (slightly above the 50% estimate by Mosquera & Kochanek, 2011). Thus,
many of the measurements from emission lines likely correspond to images
undergoing a microlensing event with amplitudes not very different from that
of the peak. Although, even with all these considerations in mind, the
coincidence between the means obtained either from the line-emission or from
the light-curves may remain questionable, the absence of a significant high
magnification tail in both histograms is a very robust common result.

###  2.3 Theoretical microlensing estimates using reverberation mapping sizes
for the quasar source

It is also possible to make a theoretical estimate of the expected impact of
microlensing from simulations based on microlensing magnification maps. The
key parameter in the simulations is the size of the continuum quasar source,
which can be estimated around $r_{s}=5\rm\,light-days$ from reverberation
mapping studies (see, e.g., Edelson et al. 2015, Fausnaugh et al. 2016, Jiang
et al. 2017, Esteban Gutiérrez et al. 2022 and references therein). Taking
this value for $r_{s}$, Esteban-Gutiérrez et al. (2022) calculate the
probability distributions of microlensing magnifications corresponding to a
population of stars, for all the objects in the sample used in the present
work. As it can be observed (see the blue lines in their Figure A1), the
impact of microlensing is concentrated around zero, with typically,
$\sigma(\Delta m)\leq 0.4$, and with a negligible tail above $|\Delta m|>1$,
in agreement with the emission-line based measurements.

##  3 Comparison with model flux ratio anomalies

In order to compare the observed microlensing flux ratios with model
predictions, we consider here the work of Shajib et al. 2019 (and its
extension by Schmidt et al. 2023, see below). These authors explicitly
introduce the question of flux-ratio anomalies related to microlensing and
provide a homogeneous data set which has been modeled in a very systematic
way. On the other hand, this work has the interesting perspective of exploring
semi-automated modeling to face the future massive data availability. In
Figure 1 we include the histogram of frequencies of (unsigned) flux-ratio
anomalies obtained from Shajib et al. (2019) corresponding to a sample of 13
quads. We use this sample just as a test bench to illustrate how to check the
impact of microlensing knowing that no general conclusion about lens modeling
can be inferred from a particular set of models that the authors themselves
consider susceptible of refinement in several ways. We hope that the present
work can be one of them.

In each quad we take image A as reference to compute the magnitude
differences. We have used the data in the F475X filter from Shajib et al.
(2019), which have the closest wavelength correspondence with the average rest
wavelength of the microlensing measurements described in Section 2555 The
average rest wavelength (see Jiménez-Vicente et al. 2012) of the used emission
lines is $1736\pm 373$Å, which for a mean redshift $z=1.9$ results in an
observed wavelength of $5034.4$Å, close to the central wavelength, $4940.7$Å,
of the 475X HST filter.. The mean of the histogram, $\langle|\Delta
m_{models\\_shajib}|\rangle=0.74$, greatly exceeds the mean of the
microlensing measurements from emission lines ($\langle|\Delta
m_{lines}|\rangle=0.33$). On the other hand, comparison of the tails of the
histograms, also reveals strong differences. The high magnification tail is
very populated in the case of the model flux anomalies (36% of pairs with
$|\Delta m_{models\\_shajib}|\geq 0.74$ magnitudes) while there is only one
case (2.2%) in the histogram of microlensing from emission lines, and two
(10%) in the sample of microlensing peaks. Finally, we have computed several
statistical tests that reject the hypothesis that this sample of flux ratio
anomalies and the sample of direct microlensing measurements estimated from
the emission lines have the same parent population (see Table 1). The same
negative conclusion is reached using the sample of microlensing estimates from
the light curves peaks (Table 1).

Very recently, the sample by Shajib et al. (2019) has been enlarged with 16
additional systems by Schmidt et al. (2023), who compute new models for all
the quads. In Figure 1 we include the histogram of anomalies corresponding to
Schmidt et al. (2023), which presents an offset mean $\langle|\Delta
m_{models\\_schmidt}|\rangle=0.64$ and a populated high magnification tail
(30% of pairs with $|\Delta m_{models\\_schmidt}|\geq 0.74$ magnitudes).
Although the results confirm those obtained from Shajib et al. (2019) sample,
the differences are slightly smaller. The performed statistical tests (see
Table 1) reject the hypothesis that model flux ratio anomalies and observed
microlensing anomalies come from the same underlying distribution. Only the
Kolmogorov-Smirnov test (known to be less sensitive than the others,
particularly to differences in the tails of the distributions) is inconclusive
for the case of the comparison of the observed sample based on emission lines
with the results by Schmidt et al. (2023).

It is interesting to notice that the average redshifts of the lens galaxies
(present sample: $0.53\pm 0.18$; Shajib et al sample: $0.44\pm 0.14$; Schmidt
et al. sample: $0.47\pm 0.23$) and of the quasar sources (present sample:
$1.90\pm 0.54$; Shajib et al sample: $2.10\pm 0.75$; Schmidt et al. sample:
$2.20\pm 0.82$), are statistically consistent.

Finally, notice that early calculations of flux-ratio anomalies by Pooley et
al. (2007) with a smaller (10 systems) sample result in $\langle|\Delta
m_{models\\_pooley}|\rangle=0.59$ with a 24% of pairs with $|\Delta
m_{models\\_pooley}|\geq 0.74$, consistent with the results corresponding to
Shajib et al. (2019) or even larger taking into account that Pooley et al.
photometric data correspond to $\sim 8000\rm\,\AA$ (observed).

Therefore, from the statistical comparison, we can conclude that the flux
ratio anomalies inferred from the lens models that we have taken as example,
can not be attributed to microlensing, and that, consequently, there is room
to improve the flux-ratio predictions of the models and reduce the bulk of the
anomalies.

##  4 Discussion

###  4.1 Statistical cross-check

According to the previous analysis, the flux-ratio predictions of the
considered lens models do not pass the statistical cross-check based on quasar
images observations. To estimate the true amplitude of the deviations of the
model predictions, we can remove666Considering that the involved parent PDFs
are one sided-normal distributions. the expected average effect of
microlensing, after which we are left with a mean model-predicted flux-ratio
anomaly of 0.66 magnitudes for Shajib et al. (2019) models (0.55 for Schmidt
et al. 2023 models) which can neither be attributed to microlensing nor to
intrinsic source variability777Because, as commented before, the random
effects of quasar variability combined with time delay are also included in
the microlensing estimates based on Eq.1..

A mean deviation of 0.66 magnitudes is so large, that other causes frequently
invoked to explain the flux-ratio anomalies can be confidently ruled out
according to previous estimates from the literature (see Motta et al. 2012,
Pooley et al. 2007). To confirm this directly from the data, we repeat the
histograms removing 12 image-pairs where extinction can play a significant
role finding that the high magnification tail remains. On the other hand, to
check the possible impact of time-delay in the flux-ratio anomalies we have
compared the histograms considering all the image-pairs in the samples with
histograms excluding pairs with time delays greater than 10 or 30 days,
respectively, finding no significant differences among them. In fact,
following the steps described in Yonehara (2008) and assuming a time delay of
30 days and an absolute demagnified magnitude of the sources in the -21 to -23
magnitude range, we find an uncertainty of just 0.26 magnitudes in the worst
case.

Millilensing (by dark matter subhaloes, for instance) may also, in principle,
contribute to flux-ratio anomalies. However, Pooley et al. (2007, 2012) find
that X-ray anomalies are much larger than the optical ones. This result,
confirmed by Jimenez-Vicente et al. (2015), indicates that the effect causing
the anomalies is sensitive to the differences in size between the X-ray and
optical sources, while they should perform as point like under the lensing
action of large mass millilenses (Pooley et al. 2007, 2012). Based on a
similar reasoning, Pooley et al. (2007) also exclude that changes in the
smooth lens model component can explain at once the anomalies in X-ray and in
the optical.

Then, the large anomalies are, indeed, an indication that there is room to
improve flux predictions of lens models. In particular, an examination of the
specific procedure followed to fit the model for each system should be
performed to analyze the origin of the discrepancies in the fluxes and their
relationship with possible uncertainties in the time delay estimates for
cosmographic studies. It is possible that with a small effort of
sophistication in the models, the predictions of the flux-ratios improve
drastically (Ertl et al. 2023), though, the impact of these changes on the
time-delays should be, anyway, examined. Although lens modeling of specific
systems is outside of the scope of this work, in Appendix A we explore flux-
ratio anomalies under some simplifying assumptions which, while they may not
reproduce the complexity of the real problem, may still provide some
interesting insight.

###  4.2 Individual model cross-check and modeling using (integral field)
spectroscopic flux-ratios

Perhaps the most interesting reflection to improve lens modeling is that
accurate experimental determinations of intrinsic flux ratios (free from the
effects of microlensing and extinction) obtained from the broad emission lines
of quasars can be used to cross-check individual models888Incidentally,
extinction can be disentangled from microlensing by using emission lines at
different wavelengths.. Moreover, for quads, the effects of variability
combined with time delays between images can be reasonably controlled (in
particular when a determination of the time delays is available). In
principle, narrow emission lines, mid-infrared or radio emission may also be
used to determine the intrinsic flux ratios free from microlensing, but as far
as the emitting regions involved are much larger than the continuum source,
their images might present different flux ratios depending on the shape,
centroid, extension and location of the source respect to the macro caustic.
Notice, also, that only the use of emission lines of wavelengths close to the
continuum as baseline automatically cancels the effects of extinction (see
Eq.1). Accurate flux-ratios free from these systematic effects can be used to
select reliable systems for cosmographic studies, rejecting those systems with
large anomalies.

The emission lines of lensed quasars are relatively bright, and Integral Field
Spectroscopy (IFS) combined with adaptive optics in large telescopes can be a
very reasonable experimental possibility to simultaneously observe the
emission in the continuum and in the lines of a large enough ($\sim 50$)
sample of lensed systems as to estimate $H_{0}$ with a few percent
precision999Assuming that systematic effects are not biasing the estimates..
For instance, a 20 magnitude object (I filter) can be observed with
HARMONI@ELT (Thatte et al. 2016, 2020), with S/N $\sim 100$ with a total time
exposure $\lesssim 10\rm\,min$, with an spatial resolution of $\sim 10\,\rm
mas$ and a spectral resolution of 0.208$\,\rm nm$ (adjacent wavelength slices
could be co-added to further increase the S/N ratio). Thanks to the broad
spectral range covered with HARMONI, several emission lines and continuum
bands can be observed at once, and it is possible to perform a simultaneous
fitting of the lens system images in all of them, significantly increasing the
reliability and robustness of the analysis.

Going a step further, can the flux ratios inferred from the emission lines be
used not only to cross-check individual models but also to effectively
constrain them? Commonly adopted uncertainties for continuum based flux-ratios
between images (of even 20%) make their use irrelevant as compared with that
of the astrometry and, likely for these reason, the use of flux-ratios has
been, in general, considered accessory. However, with IFS based flux-ratios
with a few percents of relative uncertainties, the role of flux-ratios may be
interesting to constrain the lens models101010Of course, other type of
ancillary observations (e.g. data from stellar kinematics, Kochanek, 2020) can
provide the information needed to solve the degeneracies..

In the simple study performed by us to explore the impact of the radial
dependence of the gravitational potential (see Appendix A), we see that, in
most cases, flux ratios are much less sensitive than time delays to changes in
the potential, although the impact in some image-pairs of some specific
systems may be large enough as to help breaking the degeneracy. In any event,
we have used a very simple model which, among other issues, does not take into
account any complexity of the angular part. A more thorough exploration of
lens modeling is needed to ascertain the real usefulness of precise flux
ratios to constrain the models and improve the robustness of theoretical time
delay estimates.

##  5 Summary and conclusions

We use a sample of 44 measurements from 34 image pairs of 23 lensed systems
with spectroscopic observations to obtain the histogram of microlensing
magnifications, using the emission lines to define the non-microlensed
baseline. This histogram can be used to perform a statistical cross-check of
lens models comparing with the deviations of the predicted flux-ratios with
respect to the observed ones. To illustrate this possibility we obtain the
histogram of model flux-ratio anomalies (predicted minus observed flux-ratios)
from Shajib et al. (2019) and Schmidt et al. (2023). The main conclusions are
the following:

1 - The mean value of the model anomalies ($\langle|\Delta
m_{models}|\rangle=0.74$) exceeds significantly the mean impact of
microlensing ($\langle|\Delta m_{lines}|\rangle=0.33$). The histogram of model
anomalies shows a significant tail ($|\Delta m_{models}|\geq$ 0.7 magnitudes)
not present in the histogram of directly measured microlensing magnifications.
The histogram of peak amplitudes of microlensing events obtained from
microlensing light curves, neither presents this extended tail. These results
strongly disfavors the hypothesis that the model flux-ratio anomalies arise
mainly from microlensing.

2 - Consequently, the remaining flux-ratio anomalies (after removing
microlensing and intrinsic variability combined with time delay effects) of
$\langle|\Delta m|\rangle=0.6$ magnitudes may be reduced by further
refinements of the models, which are well outside of the scope of the present
work. Using, nevertheless, an exploratory simple model, we find that the
degeneracy of astrometric model fitting with the radial distribution of mass
in the lens can account only for a relatively small part of the observed flux-
ratio anomalies (departures from ellipticity of lens galaxies can play a more
significant role).

3 - In principle, models can be cross-checked, not only statistically but also
individually, using flux ratios from emission (mid infrared, radio, broad and
narrow emission lines, for instance) coming from regions large enough as to be
insensitive to microlensing. However, if the region is too large (including
the dusty torus, the NLR or the radio jet, for instance), additional modeling
of the source and of its (extended) images is needed. In this sense, the use
of the relatively compact BLR can be less complex. We propose to use
spectroscopic data, specifically based on integral field spectroscopy, to
measure with current (SINFONI, MUSE) or future (HARMONI) instrumentation,
accurate broad emission line fluxes, to obtain flux-ratios free from
microlensing and extinction to check the models just to the experimental
uncertainties of flux photometry.

4 - In most cases, the uncertainties associated to the possible impact of
microlensing in the observed continuum flux ratios, have made them irrelevant
in lens model fitting. The use of very accurate broad emission line flux-
ratios to establish effective constrains in the models should be explored.

Finally, we recommend the consideration of the flux-ratio anomalies as a
quality check for the fitted models, and we advice to discard those systems
with unexpectedly large anomalies. Although a large flux-ratio anomaly in an
individual system is certainly not warranty of a wrong model, it may be a good
warning signal which is worth taking into account. Given that ongoing and
future surveys will produce large numbers of lensed systems suitable of being
used for cosmographic studies, discarding a fraction of suspicious systems
shall not damage the statistical quality of those studies.

This research was supported by grants PID2020-118687GBC33 and
PID2020-118687GB-C31, financed by MCIN/AEI/10.13039/501100011033. J.J.V. is
also financed by projects FQM-108, P20_00334, and A-FQM-510-UGR20/FEDER,
financed by Junta de Andalucía. V.M. acknowledges support from ANID Fondecyt
Regular #1231418 and Centro de Astrofísica de Valparaíso.

##  Appendix A What can be learnt from model flux anomalies?

The degeneracy of models based in astrometry with respect to the law
describing the radial mass distribution in the lens galaxy is an often invoked
difficulty of lens models to provide accurate estimates of the true time
delays (see Kochanek, 2020, and references therein). We can explore here
whether the flux-ratio anomalies can also be related to this
degeneracy111111Because the non geometrical part of the time delays is
proportional to the potential strength, they are strongly affected by the
degeneration in the potential, but the question is not immediate in the case
of flux-ratios, as they depend on the second derivatives of the potential.,
estimating and comparing its impact in both quantities. With this exploratory
aim, we can consider the singular isothermal ellipsoid potential (SIE)
generalized to take into account a power law dependence with size,
$\phi_{SIE}=b\left[(x^{1}\cos{\theta_{\epsilon}}-x^{2}\sin{\theta_{\epsilon}})^{2}/q^{2}+(x^{1}\sin{\theta_{\epsilon}}+x^{2}\cos{\theta_{\epsilon}})^{2}q^{2}\right]^{\beta/2}$,
where $q^{2}$ is the axial ratio. Expanding to first order this potential and
adding an external shear, we can write,

$\phi=b|\vec{r}|^{\beta}-b{|\vec{r}|^{2}\over|\vec{r}|^{2-\beta}}\gamma_{\epsilon}\cos{2(\theta-\theta_{\epsilon})}-{1\over
2}b|\vec{r}|^{2}\gamma_{ext}\ cos{2(\theta-\theta_{ext}}),$ (A1)

where the second term is analogous to the SIE quadrupole (see, e.g., Kochanek
2002). Using this potential, we fit the images and lens positions of the nine
quadruple lens systems that Shajib et al. (2019) modeled with a single lens
mass profile. We consider logarithmic slopes of the power-law in the range
$\beta=0.5$ to 1.5. Eight of the nine systems are very well fitted by our
simple model for all the values of $\beta$, ($\chi^{2}(\beta)\leq 1$),
confirming the degeneracy of models based on astrometric data with respect to
plausible radial dependences121212In the case of PS J0147+4360, a good fit of
the images positions can also be achieved, but the position of the lens galaxy
needs to be offset by more than 40 miliarcsec.. Then, we compute the variation
of flux ratios and time delays between images in the considered range of
$\beta$. In Figure 2 we illustrate, for one of the lens systems (SDSS
J0248+1913), the fractional deviation of both magnitudes (time delays vs. flux
ratios) with respect to a fiducial model that we arbitrarily select for
$\beta=1$, corresponding to SIS+$\gamma_{\epsilon}$+$\gamma_{ext}$. As it is
shown in this Figure, the maximum deviation of the flux ratios ranges from 0%
to $\pm$25% (depending on the pair of images), while a larger maximum
fractional deviation of $\pm$50% is obtained for the time delays for all the
three image pairs. Similar results (maximum fractional variations of the flux
ratios between 0 and $\pm 30\rm\%$ with a mean value of about $10\%$ while the
time-delays exhibit a much larger typical variation of about $50\rm\%$) are
derived considering all the image pairs of the lens systems in the sample when
$\beta$ is changed.

In principle, these results show that the degeneracy in the radial mass
distribution of the lens may be a common source of uncertainties for flux
ratios and time-delays. However, the range of variability of flux ratios with
$\beta$ would only account for a small part of the measured model anomalies
($\langle\Delta m\rangle=0.6$ magnitudes is equivalent to a fractional
deviation of 60%). This indicates that other ingredients of the model (aside
from the radial dependence of the lens potential) might be affecting the flux
ratios131313We are assuming that, on average, observed flux ratios are
reasonably reliable, and not contaminated by extended emission or other
potential sources of uncertainty.. In fact, on top of the degeneracies on the
radial dependence of the mass distribution, the angular dependence may also
induce biasses in the models (e.g. Kochanek, 2021, Van de Vyvere et al. 2022,
Gomer et al. 2022). Several recent works have indeed shown that assuming
elliptical models can bias the estimate of $H_{0}$ up to 10% (see e.g. Gomer &
Williams 2021; Cao et al. 2022) and flux-ratio anomalies have proven to be
indicative of non-elliptical components in the mass distribution in some
systems (e.g. Hsueh et al. 2016, 2017). The anomalies virtually disappeared
when the non-ellipticity is included in the models. Then, more complex models
(and ancillary data) are needed to account for the observed anomalies in the
flux-ratios and to explore possible correlations with uncertainties on the
time-delays.

## References

  * Bradač et al. (2002) Bradač, M., Schneider, P., Steinmetz, M., et al. 2002, A&A, 388, 373. doi:10.1051/0004-6361:20020559
  * Birrer et al. (2022) Birrer, S., Bhamre, V., Nierenberg, A., et al. 2022, Astrophysics Source Code Library. ascl:2210.005
  * Cao et al. (2022) Cao, X., Li, R., Nightingale, J. W., et al. 2022, Research in Astronomy and Astrophysics, 22, 025014. doi:10.1088/1674-4527/ac3f2b
  * Chiba (2002) Chiba, M. 2002, ApJ, 565, 17. doi:10.1086/324493
  * Dalal & Kochanek (2002) Dalal, N. & Kochanek, C. S. 2002, ApJ, 572, 25. doi:10.1086/340303
  * Di Valentino et al. (2023) Di Valentino, E., Giarè, W., Melchiorri, A., et al. 2023, MNRAS, 520, 210. doi:10.1093/mnras/stad152
  * Edelson et al. (2015) Edelson, R., Gelbord, J. M., Horne, K., et al. 2015, ApJ, 806, 129. doi:10.1088/0004-637X/806/1/129
  * Ertl et al. (2023) Ertl, S., Schuldt, S., Suyu, S. H., et al. 2023, A&A, 672, A2. doi:10.1051/0004-6361/202244909
  * Esteban-Gutiérrez et al. (2022) Esteban-Gutiérrez, A., Agües-Paszkowsky, N., Mediavilla, E., et al. 2022, ApJ, 929, 123. doi:10.3847/1538-4357/ac57c5
  * Fausnaugh et al. (2016) Fausnaugh, M. M., Denney, K. D., Barth, A. J., et al. 2016, ApJ, 821, 56. doi:10.3847/0004-637X/821/1/56
  * Fian et al. (2018) Fian, C., Guerras, E., Mediavilla, E., et al. 2018, ApJ, 859, 50. doi:10.3847/1538-4357/aabc0d
  * Gilman et al. (2017) Gilman, D., Agnello, A., Treu, T., et al. 2017, MNRAS, 467, 3970. doi:10.1093/mnras/stx158
  * Gomer & Williams (2021) Gomer, M. R. & Williams, L. L. R. 2021, MNRAS, 504, 1340. doi:10.1093/mnras/stab930
  * Gomer et al. (2022) Gomer, M. R., Sluse, D., Van de Vyvere, L., et al. 2022, A&A, 667, A86. doi:10.1051/0004-6361/202244324
  * Guerras et al. (2013) Guerras, E., Mediavilla, E., Jimenez-Vicente, J., et al. 2013, ApJ, 764, 160. doi:10.1088/0004-637X/764/2/160
  * Jiang et al. (2016) Jiang, L., Shen, Y., McGreer, I. D., et al. 2016, ApJ, 818, 137. doi:10.3847/0004-637X/818/2/137
  * Jiménez-Vicente et al. (2012) Jiménez-Vicente, J., Mediavilla, E., Muñoz, J. A., et al. 2012, ApJ, 751, 106. doi:10.1088/0004-637X/751/2/106
  * Jiménez-Vicente et al. (2015) Jiménez-Vicente, J., Mediavilla, E., Kochanek, C. S., et al. 2015, ApJ, 806, 251. doi:10.1088/0004-637X/806/2/251
  * Jiménez-Vicente & Mediavilla (2022) Jiménez-Vicente, J. & Mediavilla, E. 2022, ApJ, 941, 80. doi:10.3847/1538-4357/ac9e59
  * Hsueh et al. (2016) Hsueh, J.-W., Fassnacht, C. D., Vegetti, S., et al. 2016, MNRAS, 463, L51. doi:10.1093/mnrasl/slw146
  * Hsueh et al. (2017) Hsueh, J.-W., Oldham, L., Spingola, C., et al. 2017, MNRAS, 469, 3713. doi:10.1093/mnras/stx1082
  * Keeton (2002) Keeton, C. R. 2002, ApJ, 575, L1. doi:10.1086/342558
  * Kochanek (2002) Kochanek, C. S. 2002, ApJ, 578, 25. doi:10.1086/342476
  * Kochanek (2020) Kochanek, C. S. 2020, MNRAS, 493, 1725. doi:10.1093/mnras/staa344
  * Kochanek (2021) Kochanek, C. S. 2021, MNRAS, 501, 5021. doi:10.1093/mnras/staa4033
  * Koopmans et al. (2003) Koopmans, L. V. E., Treu, T., Fassnacht, C. D., et al. 2003, ApJ, 599, 70. doi:10.1086/379226
  * Mao & Schneider (1998) Mao, S. & Schneider, P. 1998, MNRAS, 295, 587. doi:10.1046/j.1365-8711.1998.01319.x
  * Mediavilla et al. (2009) Mediavilla, E., Muñoz, J. A., Falco, E., et al. 2009, ApJ, 706, 1451. doi:10.1088/0004-637X/706/2/1451
  * Mediavilla et al. (2011) Mediavilla, E., Muñoz, J. A., Kochanek, C. S., et al. 2011, ApJ, 730, 16. doi:10.1088/0004-637X/730/1/16
  * Mediavilla et al. (2016) Mediavilla, E., Jiménez-Vicente, J., Muñoz, J. A., et al. 2016, ApJ, 832, 46. doi:10.3847/0004-637X/832/1/46
  * Metcalf & Amara (2012) Metcalf, R. B. & Amara, A. 2012, MNRAS, 419, 3414. doi:10.1111/j.1365-2966.2011.19982.x
  * Metcalf & Madau (2001) Metcalf, R. B. & Madau, P. 2001, ApJ, 563, 9. doi:10.1086/323695
  * Metcalf & Zhao (2002) Metcalf, R. B. & Zhao, H. 2002, ApJ, 567, L5. doi:10.1086/339798
  * Mosquera & Kochanek (2011) Mosquera, A. M. & Kochanek, C. S. 2011, ApJ, 738, 96. doi:10.1088/0004-637X/738/1/96
  * Motta et al. (2012) Motta, V., Mediavilla, E., Falco, E., et al. 2012, ApJ, 755, 82. doi:10.1088/0004-637X/755/1/82
  * Motta et al. (2017) Motta, V., Mediavilla, E., Rojas, K., et al. 2017, ApJ, 835, 132. doi:10.3847/1538-4357/835/2/132
  * Moustakas & Metcalf (2003) Moustakas, L. A. & Metcalf, R. B. 2003, MNRAS, 339, 607. doi:10.1046/j.1365-8711.2003.06055.x
  * Pooley et al. (2007) Pooley, D., Blackburne, J. A., Rappaport, S., et al. 2007, ApJ, 661, 19. doi:10.1086/512115
  * Rojas et al. (2014) Rojas, K., Motta, V., Mediavilla, E., et al. 2014, ApJ, 797, 61. doi:10.1088/0004-637X/797/1/61
  * Rojas et al. (2020) Rojas, K., Motta, V., Mediavilla, E., et al. 2020, ApJ, 890, 3. doi:10.3847/1538-4357/ab63cb
  * Schechter & Wambsganss (2002) Schechter, P. L. & Wambsganss, J. 2002, ApJ, 580, 685. doi:10.1086/343856
  * Schmidt et al. (2023) Schmidt, T., Treu, T., Birrer, S., et al. 2023, MNRAS, 518, 1260. doi:10.1093/mnras/stac2235
  * Shajib et al. (2019) Shajib, A. J., Birrer, S., Treu, T., et al. 2019, MNRAS, 483, 5649. doi:10.1093/mnras/sty3397
  * Suyu et al. (2010) Suyu, S. H., Marshall, P. J., Auger, M. W., et al. 2010, ApJ, 711, 201. doi:10.1088/0004-637X/711/1/201
  * Thatte et al. (2016) Thatte, N. A., Clarke, F., Bryson, I., et al. 2016, Proc. SPIE, 9908, 99081X. doi:10.1117/12.2230629
  * Thatte et al. (2020) Thatte, N. A., Bryson, I., Clarke, F., et al. 2020, Proc. SPIE, 11447, 114471W. doi:10.1117/12.2562144
  * Van de Vyvere et al. (2022) Van de Vyvere, L., Sluse, D., Gomer, M. R., et al. 2022, A&A, 663, A179. doi:10.1051/0004-6361/202243382
  * Wisotzki et al. (1993) Wisotzki, L., Koehler, T., Kayser, R., et al. 1993, A&A, 278, L15
  * Witt et al. (1995) Witt, H. J., Mao, S., & Schechter, P. L. 1995, ApJ, 443, 18. doi:10.1086/175499
  * Xu et al. (2009) Xu, D. D., Mao, S., Wang, J., et al. 2009, MNRAS, 398, 1235. doi:10.1111/j.1365-2966.2009.15230.x
  * Xu et al. (2015) Xu, D., Sluse, D., Gao, L., et al. 2015, MNRAS, 447, 3189. doi:10.1093/mnras/stu2673
  * Yonehara et al. (2008) Yonehara, A., Hirashita, H., & Richter, P. 2008, A&A, 478, 95. doi:10.1051/0004-6361:20067014

Table 1: p-values for several statistical tests

|  SHAJ191 | SHAJ19 |  SCHM232 | SCHM23
---|---|---|---|---
|  vs. PEAKS3 |  vs. MED094 | vs. PEAKS | vs. MED09
Kolmogorov-Smirnov | 0.01 | 0.03 | 0.00 | 0.18
Epps-Singleton | 0.00 | 0.00 | 0.00 | 0.00
Anderson | 0.02 | 0.02 | 0.02 | 0.02
Cramer-von Mises | 0.03 | 0.04 | 0.01 | 0.06

Note. — 1Sample of flux-ratio anomalies from Shajib et al. (2019). 2Sample of
flux-ratio anomalies from Schmidt et al. (2023). 3Sample of microlensing
estimates based on light-curves peaks. 4Sample of single-epoch microlensing
measurements based on emission lines. See Figure 1.

Figure 1: Flux-ratio histograms corresponding to: single epoch microlensing
obtained using the broad emission lines as reference (thick, black line),
peaks in microlensing light curves (thin, blue line), models from Shajib et
al. (2019) (red shaded histogram), and models from Schmidt et al. (2023)
(green shaded histogram). Vertical lines show the mean values. The grey shaded
region corresponds to one standard deviation around the mean for single epoch
microlensing and microlensing peaks.  Figure 2: Fractional variations of time
delays vs. flux ratios from models obtained for SDSS J2048+1913 changing the
exponent of the power law mass distribution of the lens in the $\beta=0.5$ to
$1.5$ range. Each color corresponds to an image pair.