11institutetext: Institute of Theoretical Physics, Faculty of Mathematics and
Physics, Charles University,
V Holešovičkách 2, 18000 Praha 8, Czech Republic
11email<EMAIL_ADDRESS><EMAIL_ADDRESS>
# Gravitational lensing by an ellipsoidal Navarro–Frenk–White dark-matter
halo: an analytic solution and its properties

David Heyrovský  Michal Karamazov

(Received XXX / Accepted XXX)

###### Abstract

Context. The analysis of gravitational lensing by galaxies and galaxy clusters
typically relies on ellipsoidal lens models to describe the deflection of
light by the involved dark-matter halos. These models are most often based on
the isothermal density profile—not an optimal description of the halo, but
easy to use because it leads to an analytic deflection-angle formula.

Aims. Dark-matter halos are better described by the Navarro–Frenk–White
(hereafter NFW) density profile. We set out to study lensing by a general
triaxial ellipsoidal NFW halo, with the aim of providing an analytic model
that would be more consistent with the current understanding of dark-matter
halos.

Methods. We computed the conversion between the properties of a triaxial
ellipsoidal lens model and its elliptical surface-density profile. In the case
of the NFW lens model, its angular scale is defined by the projected scale
semi-major axis of the halo, while its lensing regime depends on two
parameters: the projected eccentricity $e$ and the convergence parameter
$\kappa_{\text{s}}$. We employed the Bourassa & Kantowski formalism to compute
the complex scattering function of the model, which yields the deflection-
angle components when separated into its real and imaginary parts.

Results. We present the obtained closed-form expressions for the deflection-
angle components, valid for an arbitrary eccentricity of the surface-density
profile. We use them to compute and describe the lensing properties of the
model, including: the shear, its components, and the phase; the critical
curves, caustics, and the parameter-space mapping of their different
geometries; the deformations and orientations of images.

Conclusions. The analytically solved ellipsoidal NFW lens model is available
for implementation in gravitational lensing software. The techniques
introduced here such as the image-plane analysis can prove to be useful for
understanding the properties of other lens models as well.

###### Key Words.:

Gravitational lensing: strong – Galaxies: clusters: general – Galaxies: halos
– Cosmology: dark matter

##  1 Introduction

Gravitational lensing provides a unique tool for studying the distribution of
matter in galaxies and galaxy clusters. In the case of clusters, detailed maps
of their surface density have been obtained by combined analyses of strong
lensing in their inner and weak lensing in their outer regions (e.g., Bradač
et al., 2005; Limousin et al., 2007; Finney et al., 2018; Jauzac et al.,
2018). In the strong-lensing analysis, the mass distribution is typically
reconstructed using building-block lens models describing the cluster halo and
the subhalos corresponding to cluster-member galaxies (Natarajan et al.,
2024). The analysis of strong lensing by galaxies follows a similar modular
approach, albeit with fewer components. The stellar and dark-matter components
may be modeled separately or by a combined model; these are usually
supplemented by external shear (Treu, 2010; Shajib et al., 2022) to account
for the lensing influence of other nearby galaxies or larger structures.
Subhalos due to satellite galaxies or other substructure may be added to
explain asymmetries or flux-ratio anomalies.

The individual halos, subhalos, or different matter components are generally
modeled by simple density or surface-density profiles. In practice, the range
of possible profiles is narrowed down by the requirement of an analytic
expression for the light-deflection angle. Without it the angle has to be
integrated numerically over the surface-density distribution for each light
ray, which is computationally prohibitive. As a result, density profiles that
cannot be integrated analytically are shunned, even though they may provide a
better description of the mass distribution of the lens. Analytic expressions
for the deflection angle are available for a number of spherically symmetric
lens models (Keeton, 2001). However, in the case of ellipsoidal models, one
geometric step closer to reality, analytic formulae for the deflection-angle
components have been obtained only for a few.

All widely used ellipsoidal models are based on the isothermal density
profile. The basic versions described by Kassiola & Kovner (1993) and Kormann
et al. (1994) are the singular isothermal ellipsoid (SIE, SIEMD), with a
centrally divergent density, and the non-singular isothermal
(pseudoisothermal) ellipsoid (NIE, PIEMD), which is softened by introducing a
central core. Subtracting a second NIE with a larger core leads to an
asymptotically attenuated density profile with a finite total mass. This model
appears in the literature under different names: PIEMD with smooth cutoff
(Kassiola & Kovner, 1993), superposition of PIEMD components (Natarajan &
Kneib, 1997), pseudo-Jaffe ellipsoid (Keeton, 2001), or dual pseudoisothermal
elliptical mass distribution (dPIE, Elíasdóttir et al., 2007). At present this
is the model of choice in cluster-lensing analysis, while galaxy-lensing
analysis often uses a simple SIE with external shear (Shajib et al., 2022). As
shown by Etherington et al. (2023), rather than describing the influence of
nearby matter distributions, the included “external shear” tends to play the
role of a prop compensating for an inadequate lens model.

Another class of analytic models is based on the power-law density profile:
the singular power-law ellipsoid (Grogin & Narayan, 1996; Tessore & Metcalf,
2015; O’Riordan et al., 2020), the broken-power-law-density ellipsoid (Du et
al., 2020), the broken power-law-surface-density elliptical model (O’Riordan
et al., 2021), and the truncated power-law-surface-density elliptical model
(O’Riordan et al., 2021). In all these variants the deflection angle is
expressed in terms of the hypergeometric function of a complex argument, which
requires a suitable numerical technique for efficient evaluation (Tessore &
Metcalf, 2015). The character of some of these models is approaching our
current understanding of the structure of dark-matter-halo density profiles.
In particular, Du et al. (2020) showed that the broken power-law-density
profile can be used to approximate the Navarro–Frenk–White (hereafter NFW,
Navarro et al., 1996) and the Einasto (Einasto, 1965) density profiles.

The NFW profile was originally introduced as an empirical law describing the
density profile of dark-matter halos arising from cosmological structure-
formation simulations. Observational support has built up gradually over the
years, showing that the profile described adequately the continuous matter
(i.e., dark matter + baryonic gas) in galaxy clusters (e.g., Newman et al.,
2013), and agreed with the analyses of cluster lensing (Okabe et al., 2013;
Umetsu & Diemer, 2017) and cluster X-ray emission (Ettori et al., 2013). On
the scale of individual galaxies, it is suitable for describing their halos in
cosmological simulations (Ludlow et al., 2013), and it has been
observationally supported as an approximation for elliptical-galaxy (Shajib et
al., 2021) or massive-spiral-galaxy (Rodrigues et al., 2017) halos.

The spherical NFW lens model was described by Bartelmann (1996) and by Wright
& Brainerd (2000) shortly after the discovery of the profile. The ellipsoidal
NFW model proved to be a harder challenge; its sought analytic solution has
been elusive despite its simple homoeoidal symmetry. In the studies that
implemented the model the deflection angle was either integrated numerically
(e.g., Suyu et al., 2012), or it was constructed from analytic solutions of
various approximate models. One common approximation is based on elliptical
distortions of the spherical NFW lens potential (e.g., Golse & Kneib, 2002;
Meneghetti et al., 2003; Meena & Bagla, 2023), an approach that leads to
unphysical dumbbell-shaped surface-density contours for all but the lowest
ellipticities (Kassiola & Kovner, 1993). The approximation can be improved by
adjusting the ellipticity of the potential to yield the desired ellipticity of
the innermost surface-density contours (Golse & Kneib, 2002). However, the
ellipticity of the surface-density contours changes outward right from the
halo center even in the case of low ellipticities, already breaking the
homoeoidal symmetry before deviations from elliptical symmetry become
prominent (Gomer et al., 2023). Another method of approximation is based on
emulating the ellipsoidal NFW by a sum or series of simpler analytic models.
Recently Oguri (2021) proposed to decompose the NFW surface-density profile
into a series of cored-steep-ellipsoid components and presented sums of 13 and
44 components approximating the lensing characteristics of the ellipsoidal NFW
model to relative accuracy $10^{-2}$ and $10^{-4}$, respectively.

In this work we derive analytic expressions for the deflection angle of the
ellipsoidal NFW model and describe its lensing properties. We start in Sect. 2
from a triaxial homoeoidally symmetric mass distribution with an arbitrary
spatial orientation and compute the properties of its corresponding
elliptically symmetric surface-density distribution. In Sect. 3 we review the
complex formalism of Bourassa & Kantowski (1975) for computing the deflection
angle of a lens with an elliptically symmetric surface-density distribution.
In Sect. 4 we apply the methodology introduced in the previous sections to the
ellipsoidal NFW lens, describing: the convergence profile (Sect. 4.1); the
complex scattering function and the deflection angle (Sect. 4.2); the shear
and phase (Sect. 4.3) with particular attention to the halo center (Sect.
4.3.2) and the zero-shear points (Sect. 4.3.3); the critical curves and
caustics and their variation in the lens-parameter space (Sect. 4.4); the
properties of images (Sect. 4.5). We discuss the observational relevance and
possible further modifications of the model in Sect. 5, and summarize the main
results in Sect. 6. Of the three appendices we would like to draw attention to
Appendix C in which we introduce the technique of image-plane analysis of
critical curves and caustics, a powerful visual tool for tracking the
variations of these curves in the parameter space of the lens.

##  2 Surface density of a triaxial ellipsoidal mass distribution

We consider a model of an astrophysical object with a triaxial ellipsoidal
mass distribution and homoeoidal symmetry: its density is constant on
concentric nestled ellipsoids of a fixed shape. The density

$\rho=\rho(\hat{a})$ (1)

can thus be written as a function of the semi-major axis of these ellipsoids,

$\hat{a}(\hat{x}_{1},\hat{x}_{2},\hat{x}_{3})=\sqrt{\hat{x}_{1}^{2}+\hat{x}_{2}^{2}/(1-e_{1}^{2})+\hat{x}_{3}^{2}/(1-e_{2}^{2})}\,,$
(2)

where the constant eccentricities $e_{1}$, $e_{2}$ satisfy $0\leq e_{1}\leq
e_{2}<1$. In this notation the Cartesian coordinates $\hat{x}_{1}$,
$\hat{x}_{2}$, and $\hat{x}_{3}$ correspond to the principal-axes frame with
orientations along the major, median, and minor axis of the ellipsoids,
respectively.

Figure 1: Reference frames of a triaxially ellipsoidally symmetric mass
distribution and its projected surface density. The transparent shaded gray
ellipsoid represents a surface of constant density, with principal major
$\hat{x}_{1}$, median $\hat{x}_{2}$, and minor $\hat{x}_{3}$ axes marked red,
green, and blue, respectively. The observer’s-frame axes are marked purple,
with $x_{1}^{\prime}$ and $x_{2}^{\prime}$ in the plane of the sky and
$x_{3}^{\prime}$ oriented along the line of sight to the observer. The
perimeter of the dark-gray ellipse in the background represents a
corresponding contour of constant surface density, with plane-of-the-sky axes
marked brown: $x_{1}$ along the major, and $x_{2}$ along the minor axis of the
ellipse.

In order to study gravitational lensing by such an object in the thin-lens
approximation, we need to integrate the density profile along the line of
sight to obtain its surface density. In doing so, we have to account for the
general spatial orientation of the object (see also Oguri et al., 2003; van de
Ven et al., 2009). We present the necessary formulae in Appendix A, where we
use Euler angles $\varphi,\vartheta,\psi$ to transform from the principal-axes
frame to an observer’s frame defined by plane-of-the-sky and line-of-sight
coordinates.

Taking these transformations into account, the surface density can be written
as

$\Sigma=\int_{x_{3-}^{\prime}(\hat{a}_{\text{\tiny{T}}})}^{x_{3+}^{\prime}(\hat{a}_{\text{\tiny{T}}})}\,\rho(\hat{a})\,\text{d}x_{3}^{\prime}\,,$
(3)

where the integration is carried out over the interval defined by the line-of-
sight intersections with the outermost ellipsoid enclosing the mass
distribution. Its semi-major axis $\hat{a}_{\text{\tiny{T}}}$ can be replaced
by $\infty$ for mass-distribution models ranging to infinity. The
intersections are given explicitly by setting
$\hat{a}=\hat{a}_{\text{\tiny{T}}}$ in Eq. (72), expressed as a function of
plane-of-the-sky position, eccentricities $e_{1}$, $e_{2}$, and the Euler
angles. Due to symmetry it is sufficient to evaluate the integral over half of
the interval, from the point at which the line of sight is tangent to the
smallest ellipsoid,

$\Sigma=2\,\int_{[x_{3-}^{\prime}(\hat{a}_{\text{\tiny{T}}})+x_{3+}^{\prime}(\hat{a}_{\text{\tiny{T}}})]/2}^{x_{3+}^{\prime}(\hat{a}_{\text{\tiny{T}}})}\,\rho(\hat{a})\,\text{d}x_{3+}^{\prime}(\hat{a})\,.$
(4)

Next, we use Eq. (76) to change the integration variable to the semi-major
axis $\hat{a}$, yielding

$\Sigma(\tilde{x}_{1},\tilde{x}_{2})=2\,\sqrt{\frac{(1-e_{1}^{2})(1-e_{2}^{2})}{\lambda_{1}\,\lambda_{2}}}\,\int_{\sqrt{\tilde{x}_{1}^{2}/\lambda_{1}+{\tilde{x}_{2}^{2}/\lambda_{2}}}}^{\hat{a}_{\text{\tiny{T}}}}\,\frac{\rho(\hat{a})\,\hat{a}\,\text{d}\hat{a}}{\sqrt{\hat{a}^{2}-\frac{\tilde{x}_{1}^{2}}{\lambda_{1}}-\frac{\tilde{x}_{2}^{2}}{\lambda_{2}}}}\,,$
(5)

where the parameters $\lambda_{1}$ and $\lambda_{2}$ are defined in Eq. (73).
A closer inspection of Eq. (5) reveals the elliptical symmetry of the surface
density which can be written as a function of the semi-major axis $\tilde{a}$
in the plane of the sky,

$\tilde{a}(\tilde{x}_{1},\tilde{x}_{2})=\sqrt{\tilde{x}_{1}^{2}+\frac{\tilde{x}_{2}^{2}}{1-e^{2}}}\,.$
(6)

The major axes of the constant-surface-density ellipses lie along the
$\tilde{x}_{1}$ axis, and the minor axes along the $\tilde{x}_{2}$ axis. These
axes have the same orientation as $x_{1}$ and $x_{2}$, respectively,
illustrated in Fig. 1.

The ellipses have a fixed shape, with eccentricity $e$ given by

$e=\sqrt{1-\frac{\lambda_{2}}{\lambda_{1}}}=\sqrt{2}\,\left\\{1+\frac{2-e_{2}^{2}\,\sin^{2}\vartheta-(1-\cos^{2}\psi\,\sin^{2}\vartheta)\,e_{1}^{2}}{\sqrt{\left[e_{2}^{2}\,\sin^{2}\vartheta+(1-\cos^{2}\psi\,\sin^{2}\vartheta)\,e_{1}^{2}\,\right]^{2}-4\,e_{1}^{2}\,e_{2}^{2}\,\sin^{2}\psi\,\sin^{2}\vartheta}}\right\\}^{-1/2}\,.$
(7)

As a function of general spatial orientation, for any fully triaxial object
with $e_{2}>e_{1}>0$ the eccentricity of the surface-density contours ranges
from a minimum value of $0$ for Euler angles $\vartheta=\arcsin(e_{1}/e_{2})$
and $\psi=\pm\pi/2$ (tilted orientation with the median axis in the plane of
the sky) to a maximum value of $e_{2}$ for $\vartheta=\pi/2$ and $\psi=0$ or
$\pi$ (orientation with the median axis along the line of sight). For an
axially symmetric prolate object with $e_{2}=e_{1}>0$ Eq. (7) simplifies to

$e_{\text{prol}}=e_{1}\,\left\\{\frac{1-\sin^{2}\psi\,\sin^{2}\vartheta}{1-e_{1}^{2}\,\sin^{2}\psi\,\sin^{2}\vartheta}\right\\}^{1/2}\,,$
(8)

while for an axially symmetric oblate object with $e_{2}>e_{1}=0$ we get

$e_{\text{obl}}=e_{2}\,\sin\vartheta\,.$ (9)

In view of its elliptical symmetry the surface density of a general triaxial
ellipsoidal mass distribution can be written as

$\Sigma(\tilde{a})=2\,\sqrt{\frac{(1-e_{1}^{2})(1-e_{2}^{2})}{(1-e_{1}^{2})\,\cos^{2}\vartheta+(1-e_{2}^{2})(1-e_{1}^{2}\,\sin^{2}\psi)\,\sin^{2}\vartheta}}\,\int_{\tilde{a}/\sqrt{\lambda_{1}}}^{\hat{a}_{\text{\tiny{T}}}}\,\frac{\rho(\hat{a})\,\hat{a}\,\text{d}\hat{a}}{\sqrt{\hat{a}^{2}-\tilde{a}^{2}/\lambda_{1}}}\,.$
(10)

In the special case of an axially symmetric oblate object (with $e_{1}=0$ and
$\lambda_{1}=1$) Eq. (10) reduces to the simpler form found in equation (14)
of Bourassa & Kantowski (1975).

##  3 Deflection angle in the Bourassa & Kantowski formalism

For a general mass distribution with surface density $\Sigma(\bm{w})$, the
gravitational deflection angle at lens-plane position $\bm{\tilde{}}{x}$ is
given by

$\bm{\alpha}(\bm{\tilde{}}{x})=\frac{4\,G}{c^{2}}\,\int\,\Sigma(\bm{w})\,\frac{\bm{\tilde{}}{x}-\bm{w}}{|\bm{\tilde{}}{x}-\bm{w}|^{2}}\,\text{d}^{2}\bm{w}\,,$
(11)

where $G$ is the gravitational constant and $c$ is the speed of light
(Schneider et al., 1992). For surface densities with elliptical homoeoidal
symmetry the two-dimensional integrals can be analytically reduced to one
dimension, as shown for example by Schramm (1990).

We follow another method based on complex analysis, described by Bourassa &
Kantowski (1975) and fine-tuned by Bray (1984). The components of the
deflection angle are computed as the real and imaginary parts from the
equation

$\alpha_{1}(\tilde{x}_{1},\tilde{x}_{2})+\text{i}\,\alpha_{2}(\tilde{x}_{1},\tilde{x}_{2})=\frac{4\,G}{c^{2}}\,I^{*}(\tilde{x}_{1},\tilde{x}_{2})\,,$
(12)

where i is the imaginary unit and $I^{*}$ is the complex conjugate of the
scattering function

$I(\tilde{x}_{1},\tilde{x}_{2})=\frac{2\,\pi\,\sqrt{1-e^{2}}}{\tilde{x}_{1}+\text{i}\,\tilde{x}_{2}}\,\int_{0}^{\tilde{a}(\tilde{x}_{1},\tilde{x}_{2})}\frac{\Sigma(a^{\prime})\,a^{\prime}\,\text{d}a^{\prime}}{\sqrt{1-\frac{e^{2}}{(\tilde{x}_{1}+\text{i}\,\tilde{x}_{2})^{2}}\,a^{\prime
2}}}\,.$ (13)

Here $\tilde{x}_{1}$ and $\tilde{x}_{2}$ are coordinates along the major and
minor axis of the ellipses, $e$ is their eccentricity, and the semi-major axis
$\tilde{a}(\tilde{x}_{1},\tilde{x}_{2})$ is defined in Eq. (6). The upper
limit of integration is technically valid for a surface-density distribution
extending to infinity. In the case of distributions with a finite extent, it
should be replaced by the limiting semi-major axis for positions $\bm{x}$
outside the mass distribution. Note that while it may appear tempting to
transfer the first denominator into the square-root expression inside the
integral, doing so would cancel natural symmetries of the result such as
$I(-\bm{\tilde{}}{x})=-I(\bm{\tilde{}}{x})$ and others (Bray, 1984).

Distributions for which the integral in Eq. (13) can be solved analytically
yield the scattering function $I$ in the form of a complex function. In order
to obtain real formulae for the deflection-angle components, this function has
to be separated into its real and imaginary parts.

##  4 Lensing by an ellipsoidal NFW halo

###  4.1 Surface density and convergence

The NFW density profile (Navarro et al., 1996) describes the mass distribution
in a dark-matter halo of galaxy clusters as well as galaxies. While the model
in its original form describes a spherical halo, we study its generalization
to triaxial ellipsoidal symmetry with a density profile

$\rho(\hat{a})=\rho_{\text{s}}\,\left(\frac{\hat{a}}{\hat{a}_{\text{s}}}\right)^{-1}\,\left(1+\frac{\hat{a}}{\hat{a}_{\text{s}}}\right)^{-2}\,,$
(14)

where $\hat{a}$ is the semi-major axis of a constant-density ellipsoid passing
through a given point as defined in Eq. (2), $\hat{a}_{\text{s}}$ is the scale
semi-major axis at which $\text{d}\ln{\rho}/\text{d}\ln{\hat{a}}=-2$, and
$\rho_{\text{s}}$ is a characteristic density such that
$\rho(\hat{a}_{\text{s}})=\rho_{\text{s}}/4$.

Substituting the density profile $\rho(\hat{a})$ in Eq. (10) and setting the
upper limit $\hat{a}_{\text{\tiny{T}}}\rightarrow\infty$ leads to an integral
of the same form as for the spherical NFW profile, yielding the surface
density

$\Sigma(\tilde{a})=2\,\rho_{\text{s}}\,\hat{a}_{\text{s}}\,\sqrt{\frac{(1-e_{1}^{2})(1-e_{2}^{2})}{(1-e_{1}^{2})\,\cos^{2}\vartheta+(1-e_{2}^{2})(1-e_{1}^{2}\,\sin^{2}\psi)\,\sin^{2}\vartheta}}\;\left[1-\mathcal{F}\left(\frac{\tilde{a}}{\hat{a}_{\text{s}}\sqrt{\lambda_{1}}}\right)\right]\,\left(\frac{\tilde{a}^{2}}{\hat{a}_{\text{s}}^{2}\lambda_{1}}-1\right)^{-1}\,,$
(15)

where the monotonically decreasing function $\mathcal{F}(\xi)$ is defined
piecewise as follows:

$\mathcal{F}(\xi)=\begin{cases}\cfrac{\displaystyle\operatorname{arctanh}{\sqrt{1-\xi^{2}}}}{\displaystyle\sqrt{1-\xi^{2}}}&\text{for
$\xi<1$}\,,\\\\[9.0pt] \hfil 1&\text{for $\xi=1$}\,,\\\\[6.0pt]
\cfrac{\displaystyle\arctan{\sqrt{\xi^{2}-1}}}{\displaystyle\sqrt{\xi^{2}-1}}&\text{for
$\xi>1$}\,.\end{cases}$ (16)

Note that for a given triaxial mass distribution in a fixed orientation the
entire expression before the square brackets is a constant; only the square
brackets and the final parentheses are functions of the plane-of-the-sky
position via the semi-major axis $\tilde{a}$. In addition, note that it is
advantageous to write the plane-of-the-sky separations and positions in units
of the projected scale semi-major axis of the halo,

$\tilde{a}_{\text{s}}=\hat{a}_{\text{s}}\sqrt{\lambda_{1}}\,,$ (17)

as shown in Appendix A. In these units

$\\{a,x_{1},x_{2}\\}=\\{\tilde{a},\tilde{x}_{1},\tilde{x}_{2}\\}/\tilde{a}_{\text{s}}$
(18)

denote the semi-major axis and the coordinates along the major and minor axes
of the ellipse in the plane of the sky, satisfying

$a(x_{1},x_{2})=\sqrt{x_{1}^{2}+\frac{x_{2}^{2}}{1-e^{2}}}\,.$ (19)

We account for these properties in the expression for the convergence, which
is obtained by dividing the surface density by the critical surface density

$\Sigma_{\text{cr}}=\frac{c^{2}}{4\pi
G}\frac{D_{\text{s}}}{D_{\text{l}}\,D_{\text{ls}}}\,,$ (20)

where $D_{\text{s}}$, $D_{\text{l}}$, and $D_{\text{ls}}$ are the angular
diameter distances between the observer and the source, between the observer
and the lens, and between the lens and the source, respectively. The
convergence of an NFW halo with triaxial ellipsoidal symmetry can be written
as

$\kappa(a)=\begin{cases}2\,\kappa_{\text{s}}\;\cfrac{\displaystyle
1-\mathcal{F}(a)}{\displaystyle a^{2}-1}&\text{for $a\neq 1$}\,,\\\\[9.0pt]
\hfil\cfrac{2}{3}\,\kappa_{\text{s}}&\text{for $a=1$}\,.\end{cases}$ (21)

While the functional form is exactly the same as for a spherical NFW halo
(e.g., Bartelmann, 1996; Keeton, 2001), the variable is generalized from the
radial position in units of halo scale radius to the semi-major axis in units
of projected halo scale semi-major axis. In addition, the halo convergence
parameter

$\kappa_{\text{s}}=\frac{\rho_{\text{s}}\,\hat{a}_{\text{s}}}{\Sigma_{\text{cr}}}\,\sqrt{\frac{(1-e_{1}^{2})(1-e_{2}^{2})}{(1-e_{1}^{2})\,\cos^{2}\vartheta+(1-e_{2}^{2})(1-e_{1}^{2}\,\sin^{2}\psi)\,\sin^{2}\vartheta}}$
(22)

depends not only on the characteristic density and scale semi-major axis of
the halo, but also on its shape and spatial orientation through the square-
root term which reduces to unity in the spherical case ($e_{1}=e_{2}=0$).

As a function of orientation, for a given triaxially ellipsoidal halo the
square root reaches its maximum value of 1 when viewed down the major axis
($\vartheta=\pi/2$ and $\psi=\pi/2$ or $3\pi/2$), its minimum value of
$\sqrt{1-e_{2}^{2}}$ down the minor axis ($\vartheta=0$ or $\pi$), and the
intermediate value of $\sqrt{1-e_{1}^{2}}$ down the median axis
($\vartheta=\pi/2$ and $\psi=0$ or $\pi$). For a more intuitive understanding
of the term, note that the volume of an ellipsoid is proportional to the
square root of the numerator, while the area of its elliptical projection is
proportional to the square root of the denominator, as indicated by Eq. (74).
We may thus re-write Eq. (22) in a form that is easier to interpret,

$\kappa_{\text{s}}=\frac{3}{4}\,\frac{\rho_{\text{s}}}{\Sigma_{\text{cr}}}\,\frac{V(\hat{a}_{\text{s}})}{S_{\text{proj}}(\hat{a}_{\text{s}})}\,,$
(23)

where
$V(\hat{a}_{\text{s}})=4\pi\hat{a}_{\text{s}}^{3}\sqrt{(1-e_{1}^{2})(1-e_{2}^{2})}/3$
is the volume of an ellipsoid with semi-major axis $\hat{a}_{\text{s}}$, and
$S_{\text{proj}}(\hat{a}_{\text{s}})=\pi\hat{a}_{\text{s}}^{2}\sqrt{\lambda_{1}\lambda_{2}}$
is the area of its plane-of-the-sky projection (see Appendix A and Fig. 1).

Figure 2: Lens-plane plots of the deflection angle
$\bm{\alpha}_{0}(\bm{x})/\kappa_{\text{s}}$ of the ellipsoidal NFW lens with
components given by Eq. (30) and Eq. (31). The three panels correspond to
eccentricities $e\in\\{0,0.6,0.9\\}$, as marked above the plots. The arrows
mark the deflection-angle vectors at the position of their central dots. The
amplitude of each vector is indicated by its color (with a separate color bar
for each panel), and by constant-amplitude contours at values marked
vertically across the color bars. Positions with maximum deflection amplitude
are marked by the large bold black circle in the left panel, and by the two
solid black dots along the minor axis in the second and third panels. The two
small black circles in these panels indicate saddle points with maximum
deflection along the major axis. Note that in the case of non-zero
eccentricity the deflection at positions off the major and minor axes
typically has a non-radial orientation.

###  4.2 Deflection angle

In order to obtain the deflection angle we first re-write Eq. (13) in terms of
a dimensionless integration variable $\bar{a}=a^{\prime}/\tilde{a}_{\text{s}}$
and the dimensionless quantities from Eq. (18),

$I(x_{1},x_{2})=\frac{2\,\pi\,\tilde{a}_{\text{s}}\,\sqrt{1-e^{2}}}{x_{1}+\text{i}\,x_{2}}\,\int_{0}^{a(x_{1},x_{2})}\frac{\Sigma(\bar{a})\,\bar{a}\,\text{d}\bar{a}}{\sqrt{1-\frac{e^{2}}{(x_{1}+\text{i}\,x_{2})^{2}}\,\bar{a}^{2}}}\,.$
(24)

In our case, $\Sigma(\bar{a})=\Sigma_{\text{cr}}\,\kappa(\bar{a})$ with the
convergence given by Eq. (21). The integral can be computed analytically,
yielding the complex scattering function

$I(x_{1},x_{2})=\frac{4\,\pi\,\Sigma_{\text{cr}}\,\kappa_{\text{s}}\,\tilde{a}_{\text{s}}\,\sqrt{1-e^{2}}}{(x_{1}+\text{i}\,x_{2})^{2}-e^{2}}\,\left[\sqrt{(x_{1}+\text{i}\,x_{2})^{2}-e^{2}\,\left(x_{1}^{2}+\frac{x_{2}^{2}}{1-e^{2}}\right)}\quad\mathcal{F}\left(\sqrt{x_{1}^{2}+\frac{x_{2}^{2}}{1-e^{2}}}\right)\right.\\\
\left.+(x_{1}+\text{i}\,x_{2})\,\ln{\frac{\sqrt{x_{1}^{2}+\frac{x_{2}^{2}}{1-e^{2}}}}{1+\sqrt{1-\frac{e^{2}}{(x_{1}+\text{i}\,x_{2})^{2}}\left(x_{1}^{2}+\frac{x_{2}^{2}}{1-e^{2}}\right)}}}\,\right]\,.$
(25)

This result can be substituted on the r.h.s. of Eq. (12) to obtain the
components of the deflection angle. Checking the result in the limiting $e=0$
case, we get

$\alpha_{1}(x_{1},x_{2})+\text{i}\,\alpha_{2}(x_{1},x_{2})=\frac{4\,D_{\text{s}}\,\tilde{a}_{\text{s}}\,\kappa_{\text{s}}}{D_{\text{l}}\,D_{\text{ls}}}\,\frac{x_{1}+\text{i}\,x_{2}}{x_{1}^{2}+x_{2}^{2}}\,\left[\mathcal{F}\left(\sqrt{x_{1}^{2}+x_{2}^{2}}\right)+\ln{\frac{\sqrt{x_{1}^{2}+x_{2}^{2}}}{2}}\,\right]\,,$
(26)

the expression for the deflection angle of a spherical NFW halo (e.g.,
Karamazov et al., 2021). For simplicity, we present the remaining results in
terms of the reduced deflection angle expressed in units of the angular
projected scale semi-major axis of the halo,

$\bm{\alpha}_{0}(\bm{x})=\frac{D_{\text{ls}}\,D_{\text{l}}}{D_{\text{s}}\,\tilde{a}_{\text{s}}}\,\bm{\alpha}(\bm{x})\,.$
(27)

For an arbitrary projected eccentricity $e$, Eq. (25) also directly yields
simple formulae for the deflection angle along the major axis

$\bm{\alpha}_{0}(x_{1},0)=\left(4\,\kappa_{\text{s}}\,\frac{x_{1}}{x_{1}^{2}-e^{2}}\,\left\\{(1-e^{2})\,\mathcal{F}(x_{1})+\sqrt{1-e^{2}}\,\ln{\frac{|x_{1}|}{1+\sqrt{1-e^{2}}}}\,\right\\},\,0\right)\,,$
(28)

and along the minor axis

$\bm{\alpha}_{0}(0,x_{2})=\left(0,\,4\,\kappa_{\text{s}}\,\frac{x_{2}}{x_{2}^{2}+e^{2}}\,\left\\{\mathcal{F}\left(\frac{x_{2}}{\sqrt{1-e^{2}}}\right)+\sqrt{1-e^{2}}\,\ln{\frac{|x_{2}|}{1+\sqrt{1-e^{2}}}}\,\right\\}\right)\,.$
(29)

In order to obtain real expressions for the deflection angle at an arbitrary
lens-plane position $\bm{x}$, Eq. (25) has to be decomposed into its real and
imaginary parts. This procedure leads to the following general results for the
deflection-angle component parallel to the major axis:

$\alpha_{01}(\bm{x})=\frac{4\,\kappa_{\text{s}}\,\sqrt{1-e^{2}}}{\left[(x_{1}-e)^{2}+x_{2}^{2}\right]\left[(x_{1}+e)^{2}+x_{2}^{2}\right]}\,\left\\{x_{1}\left[(x_{1}^{2}-e^{2})(1-e^{2})+x_{2}^{2}(1+e^{2})\right]\,f_{1}(x_{1},x_{2})\qquad\qquad\qquad\right.\\\
\left.+x_{1}(x_{1}^{2}+x_{2}^{2}-e^{2})\,f_{2}(x_{1},x_{2})-x_{2}(x_{1}^{2}+x_{2}^{2}+e^{2})\,f_{3}(x_{1},x_{2})\right\\}\,,$
(30)

and for the deflection-angle component parallel to the minor axis:

$\alpha_{02}(\bm{x})=\frac{4\,\kappa_{\text{s}}\,\sqrt{1-e^{2}}}{\left[(x_{1}-e)^{2}+x_{2}^{2}\right]\left[(x_{1}+e)^{2}+x_{2}^{2}\right]}\,\left\\{x_{2}\left[x_{1}^{2}(1-2\,e^{2})+x_{2}^{2}+e^{2}\right]\,f_{1}(x_{1},x_{2})\qquad\qquad\qquad\qquad\qquad\right.\\\
\left.+x_{2}(x_{1}^{2}+x_{2}^{2}+e^{2})\,f_{2}(x_{1},x_{2})+x_{1}(x_{1}^{2}+x_{2}^{2}-e^{2})\,f_{3}(x_{1},x_{2})\right\\}\,.$
(31)

In both expressions we introduced the auxiliary functions

$\displaystyle f_{1}(x_{1},x_{2})$ $\displaystyle=$
$\displaystyle(1-e^{2})^{-1/2}\,\mathcal{F}\left(\sqrt{x_{1}^{2}+x_{2}^{2}\,/(1-e^{2})}\right)$
$\displaystyle f_{2}(x_{1},x_{2})$ $\displaystyle=$
$\displaystyle\ln{\frac{\sqrt{x_{1}^{2}+x_{2}^{2}}}{1+\sqrt{1-e^{2}}}}$ (32)
$\displaystyle f_{3}(x_{1},x_{2})$ $\displaystyle=$
$\displaystyle\arctan{\frac{x_{1}x_{2}(1-\sqrt{1-e^{2}})}{x_{1}^{2}\sqrt{1-e^{2}}+x_{2}^{2}}}\,.$

Figure 3: Left panel: Position of deflection-angle maximum $x_{2\text{max}}$
(dashed line) and position of saddle point $x_{1\text{sad}}$ (solid line) as a
function of ellipticity $1-\sqrt{1-e^{2}}$ (linear bottom axis) and
eccentricity $e$ (non-linear top axis). Right panel: Deflection-angle
amplitudes $|\bm{\alpha}_{0}|$ at these positions in units of convergence
parameter $\kappa_{\text{s}}$ as a function of ellipticity and eccentricity.

The formulae for both components have the same prefactor before the curly
braces. Its denominator is equal to the product of the squared distances from
the two foci of the $a=1$ ellipse. Nevertheless, the deflection angle is not
divergent at these points; it is continuous and smooth there due to the
vanishing expressions in the curly braces. This can be seen by expanding Eq.
(28) near $x_{1}=\pm e$, as shown also in Appendix B.

Both expressions in the curly braces are formed by combinations of the same
three transcendental functions defined in Eq. (4.2), each weighted by a
polynomial in terms of $x_{i}$. The formulae for the two components thus
differ only in these polynomials. Note also that $x_{1}$ cannot be factored
from $\alpha_{01}$, just as $x_{2}$ cannot be factored from $\alpha_{02}$ — in
both cases due to the $\arctan$ term. This implies that the general
orientation of $\bm{\alpha}$ is not radial. An analysis of Eqs. (30) and (31)
shows that the deflection is oriented radially only along the major and minor
axes.

The pattern of the deflection angle $\bm{\alpha}_{0}(\bm{x})$ as a function of
lens-plane position is shown in Fig. 2 for eccentricities $e=0,0.6,0.9$. The
arrows indicate the orientation of the deflection angle at their centers
(marked by dots); their color and the contours indicate its amplitude. The
color bars with values in units of the convergence parameter
$\kappa_{\text{s}}$ are set for each panel separately. The values
corresponding to the contours are marked by black lines in the bars.

For all values of the eccentricity $e$, the deflection amplitude is zero at
the center, from where it increases outward and drops back to zero
asymptotically. In the first panel of Fig. 2 illustrating the $e=0$ case, the
ring corresponding to maximum deflection is marked by the bold contour. For
non-zero eccentricity the symmetry is broken and the global maximum of
deflection occurs along the minor (here vertical) axis, at points marked by
black dots in the two right panels. The maximum value along the major axis
occurs at saddle points which are marked in the two right panels by small
black circles.

The change of the deflection plots with eccentricity is further illustrated by
Fig. 3, with the curves in both panels plotted as a function of the
ellipticity (flattening) $1-\sqrt{1-e^{2}}$ with the eccentricity $e$ marked
along the non-linear top axis. The solid curve in the left panel shows the
position $x_{2\text{max}}$ of the global maximum along the minor axis: it
starts at $x_{2\text{max}}\approx 1.32$ for $e=0$ and decreases with
eccentricity, reaching $x_{2\text{max}}=0$ at $e=1$. The dashed curve in the
same panel shows the position $x_{1\text{sad}}$ of the saddle point along the
major axis: it starts at the same distance of $x_{1\text{sad}}\approx 1.32$
for $e=0$ and decreases with eccentricity. However, unlike in the case of the
maximum, at $e=1$ it reaches $x_{1\text{sad}}=1$, the point along the major
axis at a distance of the projected scale semi-major axis from the halo
center.

The right panel of Fig. 3 shows the deflection-angle amplitude at the global
maximum (solid curve) and at the saddle point (dashed curve). Both curves
start at $\alpha_{0}\approx 1.24\,\kappa_{\text{s}}$ for $e=0$ and decrease to
$\alpha_{0}=0$ for $e=1$, with the decline particularly pronounced at higher
projected eccentricities.

From the perspective of the three-dimensional halo, it can be shown from Eq.
(7) that the extreme value $e=1$ can be achieved only for a halo with
$e_{2}=1$ and either $\vartheta=\pi/2$ (flat halo viewed edge-on) or $e_{1}=1$
(extremely prolate or filament-like halo). Except in the case of an exactly
edge-on orientation, for a halo with $e_{2}\rightarrow 1$ the convergence
parameter given by Eq. (22) drops to zero. The decline of the deflection-angle
amplitude at high projected eccentricities is thus even more pronounced in
absolute terms.

###  4.3 Shear and phase

The formulae for the deflection angle can be used to compute the two
components $\gamma_{1}$, $\gamma_{2}$ of the shear $\gamma$. These are given
in terms of the lens-plane-position derivatives of the deflection angle as
follows:

$\gamma_{1}\equiv\gamma\cos{2\phi}=\frac{1}{2}\left(\frac{\partial\,\alpha_{01}}{\partial\,x_{1}}-\frac{\partial\,\alpha_{02}}{\partial\,x_{2}}\right)$
(33)

and

$\gamma_{2}\equiv\gamma\sin{2\phi}=\frac{\partial\,\alpha_{01}}{\partial\,x_{2}}=\frac{\partial\,\alpha_{02}}{\partial\,x_{1}}\,,$
(34)

where $\phi$ is the lens phase. Using the deflection angle components from
Eqs. (30) and (31) yields the expressions

$\gamma_{1}(\bm{x})=\frac{4\,\kappa_{\text{s}}\,\sqrt{1-e^{2}}}{\left[(x_{1}-e)^{2}+x_{2}^{2}\right]^{2}\left[(x_{1}+e)^{2}+x_{2}^{2}\right]^{2}}\,\left\\{\left[(x_{1}-e)^{2}+x_{2}^{2}\right]\left[(x_{1}+e)^{2}+x_{2}^{2}\right](x_{1}^{2}-x_{2}^{2}-e^{2})\,f_{0}(x_{1},x_{2})\right.\\\\[9.0pt]
\left.+\left[2\,e^{2}(x_{1}^{2}-1)\left[(x_{1}^{2}-x_{2}^{2}-e^{2})^{2}-4\,x_{1}^{2}x_{2}^{2}\right]-(x_{1}^{2}-x_{2}^{2}-e^{2})^{3}(3+e^{2})/2+6\,x_{1}^{2}x_{2}^{2}(e^{2}-1)(x_{1}^{2}-x_{2}^{2}-e^{2})\right]\,f_{1}(x_{1},x_{2})\,\right.\\\\[9.0pt]
\left.-\left[(x_{1}^{2}-x_{2}^{2}-e^{2})\left([x_{1}^{2}+x_{2}^{2}]^{2}-e^{4}\right)-8\,e^{2}x_{1}^{2}x_{2}^{2}\right]\,f_{2}(x_{1},x_{2})+2\,x_{1}x_{2}\left[(x_{1}^{2}+x_{2}^{2}+e^{2})^{2}-4\,e^{2}(x_{2}^{2}+e^{2})\right]\,f_{3}(x_{1},x_{2})\right\\}$
(35)

and

$\gamma_{2}(\bm{x})=\frac{4\,\kappa_{\text{s}}\,\sqrt{1-e^{2}}}{\left[(x_{1}-e)^{2}+x_{2}^{2}\right]^{2}\left[(x_{1}+e)^{2}+x_{2}^{2}\right]^{2}}\,\left\\{2\,x_{1}x_{2}\,\left[(x_{1}-e)^{2}+x_{2}^{2}\right]\left[(x_{1}+e)^{2}+x_{2}^{2}\right]\,f_{0}(x_{1},x_{2})\right.\\\\[9.0pt]
\left.+x_{1}x_{2}\,\left([x_{1}^{2}+x_{2}^{2}+e^{2}]\left[(5\,e^{2}-3)\,x_{1}^{2}-3\,(1+e^{2})\,x_{2}^{2}+(5-3\,e^{2})\,e^{2}\right]-4\,e^{2}x_{1}^{2}[1+e^{2}]\right)\,f_{1}(x_{1},x_{2})\right.\\\\[9.0pt]
\left.-2\,x_{1}x_{2}\,\left[(x_{1}^{2}+x_{2}^{2}+e^{2})^{2}-4\,e^{2}(x_{2}^{2}+e^{2})\right]\,f_{2}(x_{1},x_{2})-\left[(x_{1}^{2}-x_{2}^{2}-e^{2})\left([x_{1}^{2}+x_{2}^{2}]^{2}-e^{4}\right)-8\,e^{2}x_{1}^{2}x_{2}^{2}\right]\,f_{3}(x_{1},x_{2})\right\\}\,,$
(36)

where the auxiliary functions $f_{1},f_{2},f_{3}$ are defined in Eq. (4.2) and
we newly introduce

$f_{0}(x_{1},x_{2})=1+\cfrac{1}{2\sqrt{1-e^{2}}}\,\cfrac{x_{1}^{2}+x_{2}^{2}+e^{2}-2+(1-e^{2}x_{1}^{2})\,\mathcal{F}\left(\sqrt{x_{1}^{2}+x_{2}^{2}\,/(1-e^{2})}\right)}{1-x_{1}^{2}-x_{2}^{2}/(1-e^{2})}\,.$
(37)

The expressions in Eqs. (35) and (36) can be combined to yield the shear

$\gamma(\bm{x})=\sqrt{\gamma_{1}^{2}(\bm{x})+\gamma_{2}^{2}(\bm{x})}$ (38)

and the phase $\phi(\bm{x})$ can then be computed from Eqs. (33) and (34).

Figure 4: Large-scale lens-plane maps of the shear component $\gamma_{1}$ (top
row), shear component $\gamma_{2}$ (second row), shear $\gamma$ (third row),
and phase $\phi$ (bottom row) of the ellipsoidal NFW lens for eccentricities
$e\in\\{0,0.6,0.9\\}$, as marked above the three columns. Colors and contour
values for each row are marked in the color bars on the right side; shear and
its components are plotted in units of $\kappa_{\text{s}}$.

####  4.3.1 Large-scale properties

Lens-plane maps of the shear components, the shear, and the phase are
presented in Fig. 4 for projected eccentricities $e=0$, 0.6, and 0.9.

In the zero-eccentricity case shown in the left column, the shear component
$\gamma_{1}$ (top row) is zero along both diagonals, and it reaches its
azimuthal positive maxima and negative minima along the vertical and
horizontal axes, respectively. The component $\gamma_{2}$ (second row) has the
same pattern rotated counter-clockwise by $\pi/4$: it is zero along both axes,
it reaches its azimuthal extremes along the diagonals, and it alternates in
the plot quadrants between negative (first and third) and positive values
(second and fourth). The shear (third row) has a purely radial dependence,
reaching its maximum value $\gamma=\kappa_{\text{s}}$ at the halo center and
decreasing outward asymptotically to zero (e.g., Karamazov & Heyrovský, 2022).
The phase $\phi$ (bottom row) has a purely azimuthal dependence: it is equal
to the azimuthal angle measured from the horizontal axis shifted by $-\pi/2$
and limited here to the interval $(-\pi/2,\pi/2]$. The phase is thus zero
along the vertical axis, $\pi/2$ along the horizontal axis, and it alternates
in the plot quadrants between negative and positive values in the same manner
as $\gamma_{2}$.

For non-zero eccentricities shown in the two right columns, the character of
the plots changes significantly. Least so in the case of $\gamma_{2}$ (second
row): the azimuthal extremes occur along lines closer to the major axis, and
the absolute values in units of $\kappa_{\text{s}}$ decrease, as indicated by
the shrinking contours. In the case of $\gamma_{1}$ (top row), the absolute
values decrease in the positive areas along the minor axis, but they increase
in the negative areas along the major axis. However, while the inner negative
contours of $\gamma_{1}$ expand, the outer ones shrink. Note also that while
the zero contours of $\gamma_{2}$ lie always along the coordinate axes
irrespective of the eccentricity, the zero-contours of $\gamma_{1}$ do not
reach the halo center when $e>0$. The pattern of the shear $\gamma$ (third
row) changes dramatically: the outer oval contours give way to the central
dipole-like pattern dominated by the character of $\gamma_{1}$. Along the
minor axis just above and below the center there are localized low-shear
regions. The phase (bottom row) preserves the large-scale pattern in the
quadrants, but the contours are tilted toward the horizontal axis. More
specifically, the $\gamma_{2}=0$ contours correspond to $\phi=0$ and
$\phi=\pi/2$ contours in all the phase plots in Fig. 4 (as well as in Fig. 5),
while the $\gamma_{1}=0$ contours correspond to $|\phi|=\pi/4$ contours. Note
in the bottom right panel of Fig. 4 that non-zero eccentricities lead to
deviations from the large-scale pattern in the vicinity of the halo center,
which are explored in the following Sect. 4.3.2.

Figure 5: Halo-center details of the lens-plane maps from Fig. 4, with the
same color bars and additional shear contours in the first and third rows. In
the maps of the shear on the third row: the central dot in the first panel
corresponds to $\gamma=\kappa_{\text{s}}$; due to the very low shear variation
in the first panel the dotted $\gamma=0.9\,\kappa_{\text{s}}$ contour is added
in all three panels; the cyan dots in the second and third panels mark the
positions of the zero-shear points. Remaining notation as in Fig. 4.

####  4.3.2 Halo-center features

In order to understand the central patterns, we present in Fig. 5 the same set
of plots zoomed in ten times, and interpret them using expansions of the
plotted quantities in the vicinity of the halo center. In the zero-
eccentricity case, the expansions of the shear and its components to leading
order for $x\ll 1$ yield

$\gamma(\bm{x})=\kappa_{\text{s}}+\mathcal{O}(x^{2}\,\ln{x})\,,$ (39)
$\gamma_{1}(\bm{x})=\kappa_{\text{s}}\,\frac{x_{2}^{2}-x_{1}^{2}}{x_{1}^{2}+x_{2}^{2}}+\mathcal{O}(x^{2}\,\ln{x})\,,$
(40)

and

$\gamma_{2}(\bm{x})=-\,\kappa_{\text{s}}\,\frac{2\,x_{1}\,x_{2}}{x_{1}^{2}+x_{2}^{2}}+\mathcal{O}(x^{2}\,\ln{x})\,,$
(41)

while the phase at any distance $x$ from the center fulfills

$\left(\,\cos{2\phi(\bm{x})},\,\sin{2\phi(\bm{x})}\,\right)=\left(\frac{x_{2}^{2}-x_{1}^{2}}{x_{1}^{2}+x_{2}^{2}},\,-\,\frac{2\,x_{1}\,x_{2}}{x_{1}^{2}+x_{2}^{2}}\right)\,.$
(42)

The expression for the shear explains the featureless bottom left plot showing
the slow radial decrease of the shear. The expressions for the shear
components explain their alternating sign around the center. In addition, they
demonstrate that both components are undefined directly at the halo center:
approaching the center from different directions leads to all values in the
interval $[-\kappa_{\text{s}},\kappa_{\text{s}}]$. Nevertheless, their
combination given by Eq. (38) that yields the shear is well defined and
analytic at the halo center. The phase given by Eq. (42) is always oriented
perpendicular to the position vector $\bm{x}$ (see Karamazov & Heyrovský,
2022); at the center it is undefined.

In the case of non-zero eccentricity, the central expansions require the
analogous condition $a\ll 1$ plus $|x_{1}|,|x_{2}|\ll e$. Hence, the following
results are valid well within the unit ellipse and at the same time well
within the circle passing through its foci. In this region the shear and its
components can be approximated by

$\gamma(\bm{x})=\frac{2\,\kappa_{\text{s}}}{e^{2}}\left\\{\left[\left(1-\sqrt{1-e^{2}}\right)^{2}\,\left(\ln{\frac{\sqrt{x_{1}^{2}+x_{2}^{2}/(1-e^{2})}}{2}}+1\right)+2\,\sqrt{1-e^{2}}\,\ln{\left(\frac{1+\sqrt{1-e^{2}}}{2}\,\sqrt{1+\frac{e^{2}}{1-e^{2}}\,\frac{x_{2}^{2}}{x_{1}^{2}+x_{2}^{2}}}\right)}\;\right]^{2}\right.\\\
\left.+4\,(1-e^{2})\,\arctan^{2}{\frac{x_{1}x_{2}(1-\sqrt{1-e^{2}})}{x_{1}^{2}\sqrt{1-e^{2}}+x_{2}^{2}}}\right\\}^{1/2}+\mathcal{O}(x^{2}\,\ln{x})\,,$
(43)
$\gamma_{1}(\bm{x})=\frac{2\,\kappa_{\text{s}}}{e^{2}}\left[\left(1-\sqrt{1-e^{2}}\right)^{2}\,\left(\ln{\frac{\sqrt{x_{1}^{2}+x_{2}^{2}/(1-e^{2})}}{2}}+1\right)+2\,\sqrt{1-e^{2}}\,\ln{\left(\frac{1+\sqrt{1-e^{2}}}{2}\,\sqrt{1+\frac{e^{2}}{1-e^{2}}\,\frac{x_{2}^{2}}{x_{1}^{2}+x_{2}^{2}}}\right)}\;\right]+\mathcal{O}(x^{2}\,\ln{x})\,,$
(44)

and

$\gamma_{2}(\bm{x})=-4\,\kappa_{\text{s}}\,\frac{\sqrt{1-e^{2}}}{e^{2}}\,\arctan{\frac{x_{1}x_{2}(1-\sqrt{1-e^{2}})}{x_{1}^{2}\sqrt{1-e^{2}}+x_{2}^{2}}}+\mathcal{O}(x^{2}\,\ln{x})\,,$
(45)

where the expression “$x^{2}\,\ln{x}$” in the next higher order stands for
both $x_{i}^{2}\,\ln{a}$ and $x_{i}^{2}\,\ln{x}$ terms. The phase has to be
computed from these approximations using Eq. (33) and Eq. (34).

The central pattern of $\gamma_{2}$ (second row in Fig. 5) reveals the same
character as in the case of zero eccentricity, illustrating the two systematic
changes with increasing $e$ noted above. First, the axes of maximum negative
and positive values tilt from the diagonals toward the major axis. The
position angle of the negative axis in the first quadrant is given by
$\varphi_{\text{neg}}=\arctan{(1-e^{2})^{1/4}}$, which can be shown from Eq.
(45). Second, the amplitude of $\gamma_{2}$ at the center when approached from
different directions decreases from $\kappa_{\text{s}}$: at first slowly, then
dropping to zero for $e\to 1$. The component $\gamma_{2}$ is thus undefined at
the halo center, just as in the case of zero eccentricity.

The central pattern of $\gamma_{1}$ (top row) reveals that for non-zero
eccentricity the zero contours do not follow the diagonals to the center, as
indicated above. Instead, they cross the minor axis at points that shift away
from the center as the eccentricity increases. Directly at the center,
$\gamma_{1}\to-\infty$ from all directions. Equation (44) shows that
$\gamma_{1}\propto\ln{a}$ as $a\to 0$, so that the extreme negative region is
more extended along the major axis. The divergent term vanishes for $e\to 0$,
so that the next term in Eq. (44) determines the overall pattern for low
eccentricities. Nevertheless, despite its divergence, $\gamma_{1}$ is defined
at the halo center in the case of non-zero eccentricity.

The central pattern of $\gamma$ (third row) reveals the $\gamma\propto-\ln{a}$
divergence driven by $\gamma_{1}$, extending primarily along the major axis.
However, when progressing outward along the minor axis the shear drops to
zero, then increases and, as seen in Fig. 4, decreases asymptotically to zero.
The immediate vicinity of the halo center thus exhibits a dramatic variation
of the shear, which ranges there from $0$ to $\infty$ for any non-zero
eccentricity. This behavior is in stark contrast to the near-constant shear of
a halo with zero projected eccentricity (left column).

In the bottom row, note that the phase $\phi$ switches along the minor axis
from $0$ to $\pi/2$ between the zero-shear points. In the immediate vicinity
of the center the phase is thus vertical along both axes and near-vertical off
the axes. All contours meet at the zero-shear points, with $\phi=0$ pointing
outward along the minor axis, $\phi=\pi/2$ inward along the minor axis and all
other values along contours pointing sideways under different angles. Contours
pointing in opposite directions have perpendicular phases ($\pm\pi/2$
difference). We conclude that in the case of non-zero eccentricity the phase
is defined at the halo center, but it is undefined at the zero-shear points.

The phase pattern in combination with the position relative to the unit-
convergence contour determines the orientation of image distortions. The
interplay in the case of the studied ellipsoidal NFW model is demonstrated
further in Sect. 4.5.

Figure 6: Position of the zero-shear point $x_{0\,\displaystyle{\gamma}}$
along the minor axis as a function of ellipticity $1-\sqrt{1-e^{2}}$ (linear
bottom axis) and eccentricity $e$ (non-linear top axis). Solid black curve:
exact value obtained by numerically solving Eq. (46); dashed gray curve: low-
eccentricity approximation given by Eq. (47).

####  4.3.3 Zero-shear points

As shown in Karamazov & Heyrovský (2022), zero-shear points mark the lens-
plane positions of undistorted positive-parity images, while undistorted
negative-parity images occur along the unit-convergence contour. In the
special case when the zero-shear points lie on the unit-convergence contour,
they form umbilic points belonging to the critical curve (see Sect. 4.4). In
the case of our lens model the zero-shear points occur for any non-zero
eccentricity at the intersections of the $\gamma_{1}=0$ and $\gamma_{2}=0$
contours. The pattern of these contours in Fig.s 4 and 5 shows that the zero-
shear points are located at
$\bm{x}_{0\,\displaystyle{\gamma}}=(0,\pm\,x_{0\,\displaystyle{\gamma}})$,
where $x_{0\,\displaystyle{\gamma}}$ is the positive root of
$\gamma_{1}(0,x_{0\,\displaystyle{\gamma}})=0$. Using Eq. (35), this condition
can be simplified to

$-2\,(x_{0\,\displaystyle{\gamma}}^{2}+e^{2})\,f_{0}(0,x_{0\,\displaystyle{\gamma}})+\left[x_{0\,\displaystyle{\gamma}}^{2}(3+e^{2})-e^{2}(1-e^{2})\right]\,f_{1}(0,x_{0\,\displaystyle{\gamma}})+2\,(x_{0\,\displaystyle{\gamma}}^{2}-e^{2})\,f_{2}(0,x_{0\,\displaystyle{\gamma}})=0\,,$
(46)

which can be solved numerically, taking into account the definitions of
$f_{1}$ and $f_{2}$ in Eq. (4.2), and $f_{0}$ in Eq. (37). Zero-shear points
identified in this way are marked by the cyan dots in the $e=0.6$ and $e=0.9$
columns of the third row of Fig. 5.

The position $x_{0\,\displaystyle{\gamma}}$ of the zero-shear points on the
minor axis is shown in Fig. 6. For low eccentricities (marked along the top
axis) the zero-shear points lie exponentially close to the halo center. They
start to deviate from the center noticeably for $e\approx 0.5$, reaching a
peak separation of $x_{0\,\displaystyle{\gamma}}\approx 0.132$ at $e\approx
0.973$ and dropping back to zero for $e\to 1$. Resolving any effects
associated with zero-shear points thus requires halos with projected
eccentricities $e>0.5$.

For a more analytical insight into the low-eccentricity behavior, the solution
of the implicit Eq. (46) can be approximated by the expression

$x_{0\,\displaystyle{\gamma}}\simeq
2\,\exp{\left[\,-2\,e^{-2}-\frac{1}{4}-\frac{7}{24}\,e^{2}\,\right]}\,,$ (47)

which is accurate to $1\%$ for $e\leq 0.5$ and to $5\%$ for $e\leq 0.96$. This
expression was obtained by expanding Eq. (46) for $x\ll 1$, setting the
combination of the two highest-order terms equal to zero, and expanding this
equation in terms of $e\ll 1$. The approximation is marked by the dashed gray
curve in Fig. 6. It traces the exact solution excellently, deviating just
before the high-eccentricity peak.

###  4.4 Critical curves and caustics

In the present notation it is convenient to write the angular position of the
source $\bm{\beta}$ in units of the angular projected scale semi-major axis of
the halo, i.e.,

$\bm{y}=\frac{D_{\text{l}}}{\tilde{a}_{\text{s}}}\,\bm{\beta}\,.$ (48)

Utilizing the expression for the reduced deflection angle in the same units
from Eq. (27), we can write the lens equation in the simple form

$\bm{y}=\bm{x}-\bm{\alpha}_{0}(\bm{x})\,.$ (49)

In order to compute the lens-equation Jacobian we combine the expressions for
the convergence from Eq. (21) and for the shear components from Eq.s (35) and
(36) in the general formula

$\mathrm{det}\,J(\bm{x})=\frac{\partial\,\bm{y}}{\partial\,\bm{x}}=\left[\,1-\kappa\left(a[\bm{x}]\right)\,\right]^{2}-\gamma_{1}^{2}(\bm{x})-\gamma_{2}^{2}(\bm{x})\,.$
(50)

Lens-plane points with zero Jacobian form the critical curve, which may have
tangential and radial components defined by the factorization of the
expression on the right-hand side of Eq. (50). The tangential critical curve
is restricted to regions with sub-critical surface density and its points
$\bm{x}_{\text{T}}$ satisfy

$\gamma(\bm{x}_{\text{T}})=1-\kappa\left(a[\bm{x}_{\text{T}}]\right)\,.$ (51)

The radial critical curve is restricted to regions with super-critical surface
density and its points $\bm{x}_{\text{R}}$ satisfy

$\gamma(\bm{x}_{\text{R}})=\kappa\left(a[\bm{x}_{\text{R}}]\right)-1\,.$ (52)

The corresponding tangential and radial components of the caustic are obtained
by mapping the critical curve to the source plane using Eq. (49), i.e.,
$\bm{y}_{\text{T}}=\bm{y}(\bm{x}_{\text{T}})$ and
$\bm{y}_{\text{R}}=\bm{y}(\bm{x}_{\text{R}})$, respectively.

Figure 7: Tangential (green) and radial (purple) critical curves (top row) and
caustics (bottom row) of ellipsoidal NFW lenses with convergence parameter
$\kappa_{\text{s}}=0.25$ and eccentricities $e\in\\{0,0.3,0.6,0.811,0.9\\}$,
as marked above the top panels. Background: lens-plane color maps of the
Jacobian $\mathrm{det}\,J(\bm{x})$ in the top row; source-plane gray-scale
maps of total point-source magnification $\mu_{\text{tot}}(\bm{y})$ in the
bottom row.

####  4.4.1 Variation of geometries with eccentricity

Figure 7 shows a typical sequence of critical-curve (top row) and caustic
geometries (bottom row) as a function of eccentricity increasing from the left
to the right column, in this example for convergence parameter
$\kappa_{\text{s}}=0.25$. The tangential components of the curves are marked
in green; the radial components in purple.

The tangential critical curve for $e=0$ is the outer circular component. At
$e=0.3$ it is oval-shaped with its major axis horizontal, along the major axis
of the projected mass distribution. At $e=0.6$ the oval is more compressed,
locally concave in the vicinity of the vertical axis. Eccentricity $e=0.811$
corresponds to the hyperbolic-umbilic configuration for
$\kappa_{\text{s}}=0.25$. The tangential and radial critical curves come into
contact at the two umbilic points along the vertical axis. These are defined
by zero shear and unit convergence, as indicated by Eq. (51) and Eq. (52).
Note that the tangential critical curve has sharp points with an undefined
tangent at the umbilic points, but no inflection points. At $e=0.9$ the
tangential critical curve is peanut-shaped, very elongated horizontally and
locally concave close to the vertical axis.

The radial critical curve for $e=0$ is the inner circular component. At
$e=0.3$ and $e=0.6$ it is oval-shaped with its major axis vertical, along the
minor axis of the projected mass distribution. At the $e=0.811$ hyperbolic-
umbilic configuration, the radial critical curve has sharp points at the two
umbilic points along the vertical axis. At these points the radial critical
curve has an undefined tangent; at the same time it is compressed and locally
concave in the vicinity of the horizontal axis. At $e=0.9$ the radial critical
curve is small and peanut-shaped, elongated vertically and locally concave
close to the horizontal axis.

The color-map background of the critical-curve plots shows the Jacobian, with
shades of red and blue indicating positive and negative values, respectively.
Since image magnification is given by the absolute inverse of the Jacobian,
paler regions correspond to higher magnification. Note that while these are
fairly extended at low eccentricity, at high eccentricity (and the same
$\kappa_{\text{s}}$) high magnification is limited to a narrow region along
the critical curve.

In the bottom row of Fig. 7, the tangential caustic is limited to the central
point in the source plane at $e=0$. At $e=0.3$ it is a small four-cusped
caustic with nearly identical extent along the horizontal and vertical axes.
At $e=0.6$ the size is larger and the asymmetry is apparent, with the
horizontal extent larger than the vertical and the cusps along the horizontal
axis piercing the radial caustic. At the $e=0.811$ hyperbolic-umbilic
configuration the tangential caustic fully encloses the radial caustic. Here
the cusps along the vertical axis are replaced by hyperbolic umbilic points,
shared by both caustics. At $e=0.9$ the tangential caustic is large and
horizontally spindle-shaped, locally convex near the vertical axis and with
only two cusps, which are located along the horizontal axis. The cusps along
the vertical axis were transferred to the radial caustic in the hyperbolic-
umbilic metamorphosis.

The radial caustic at $e=0$ is circular. At $e=0.3$ and $e=0.6$ it is oval-
shaped with its major axis vertical, fully enclosing the tangential caustic at
$e=0.3$, but pierced by its cusps along the horizontal axis at $e=0.6$. At the
$e=0.811$ hyperbolic-umbilic configuration the radial caustic extends along
the vertical axis to the umbilic points, at which it has an undefined tangent.
At $e=0.9$ the radial caustic is small and vertically spindle-shaped, convex
near the horizontal axis with two cusps along the vertical axis.

The gray-scale maps forming the background of the caustic plots show the total
point-source magnification,

$\mu_{\text{tot}}(\bm{y})=\sum_{i}|\mathrm{det}\,J(\bm{x}_{i})|^{-1}\,,$ (53)

where $\bm{x}_{i}$ are the positions of all images formed by a point source at
$\bm{y}$ obtained by solving the lens Eq. (49). In the
$\kappa_{\text{s}}=0.25$ sequence presented in Fig. 7, sources located outside
both, in between, and inside both caustics have one, three, and five images,
respectively. Features to notice in the magnification maps are: the
particularly high magnification inside the tangential caustic in the $e=0.3$
and $e=0.6$ maps (including within the protruding cusps); the high-
magnification regions outside the cusps truncated by the radial caustic along
the vertical axis in the $e=0.6$ map; or the low magnification increase from
the center inside the radial caustic in the $e=0.811$ and $e=0.9$ maps.

Figure 8: Low- and intermediate-$\kappa_{\text{s}}$ lensing regimes of the
ellipsoidal NFW lens mapped in its parameter space, as described in Sect.
4.4.2 and Sect. 4.4.3. Eccentricity $e$ is marked on the right vertical axes;
ellipticity $1-\sqrt{1-e^{2}}$ on the left vertical axes; convergence
parameter on the horizontal axis up to $\kappa_{\text{s}}=2.1$ in the left
panel and up to $\kappa_{\text{s}}=8$ in the right panel. Plotted boundaries
mark the following transitions: cusp piercing (bold solid black); tangential
peanut (thin dashed green); radial peanut (thin solid purple); hyperbolic
umbilic (dashed black); tangential bulge (thin dotted green); butterfly
metamorphosis (bold dashed green). Parameter combinations of the critical-
curve and caustic examples presented in Fig. 7 and Fig. 9 are marked by red
and blue dots, respectively. Cyan dotted circles mark combinations for which
image-orientation plots are presented in Fig. 17; the green dot marks the
branching point of the tangential-bulge boundary from the tangential-peanut
boundary. Vertical dashed lines: $\kappa_{\text{s}}$ values for which
critical-curve and caustic dimensions are plotted in Fig. 12 and Fig. 14,
respectively. Horizontal dashed lines: $e$ values for which the image-plane
analysis is demonstrated in Fig. 19 and Fig. 20.

####  4.4.2 Lensing regimes in parameter space: low convergence parameter

In the previous section we demonstrated the main lensing regimes of the
ellipsoidal NFW lens, as defined by the corresponding critical curves and
caustics, for a single value of the convergence parameter
$\kappa_{\text{s}}=0.25$ and a sequence of eccentricities. Here we describe
the parameter-space map indicating the occurrence of the different regimes for
the full range of eccentricities $e\in[0,1)$ and two convergence-parameter
intervals, as shown in Fig. 8. In the left panel we present a detail of the
map in the low-convergence regime, shown here for
$\kappa_{\text{s}}\in[0,2.1]$. For convenience, we mark the ellipticity
($1-\sqrt{1-e^{2}}\,$) along the left axis and show the corresponding
eccentricity values along the non-linear right axis.

The plotted curves mark boundaries corresponding to changes in the geometry of
the critical curves and of the caustics. Green lines indicate changes on the
tangential components, purple lines on the radial components, and black lines
involve both. Thin lines indicate changes only on the critical curve; bold
lines indicate changes on the caustic, with the top-most dashed bold line
involving the critical curve as well. Solid lines indicate changes occurring
along the major axis of the mass distribution, dashed lines indicate changes
along the minor axis, and the dotted line indicates changes occurring off
either axis. The five red dots along the $\kappa_{\text{s}}=0.25$ vertical
line mark the configurations illustrated in Fig. 7.

The four primary boundaries seen in the left panel all extend further beyond
its range, intersecting all vertical and horizontal lines, i.e., the
corresponding transitions occur for any fixed value of $\kappa_{\text{s}}$ or
$e\in(0,1)$. Note that the boundaries do not intersect within the parameter
limits of the left panel, so that the sequence of the corresponding
transitions remains the same. In order of increasing eccentricity, the lowest
bold solid black boundary marks the piercing of the radial caustic by the
tangential caustic along the major axis. The equation identifying this
condition can be derived from the lens Eq. (49). The piercing occurs when the
two caustic points along the major axis each generate one image on the
tangential critical curve (on the same side from the center) and another image
on the radial critical curve (on the opposite side):

$x_{\text{T}}-\alpha_{0}(x_{\text{T}},0)=-x_{\text{R}}-\alpha_{0}(-x_{\text{R}},0)\,,$
(54)

where $\alpha_{0}$ is the first component of the deflection angle from Eq.
(28). In our choice $x_{\text{T}}>0$ and $x_{\text{R}}>0$ are the major-axis
intercepts of the tangential and radial critical curves, respectively. Using
Eqs. (51) and (52) and taking into account the properties of the shear
components, the equations for the intercepts are

$1-\kappa(x_{\text{T}})+\,\gamma_{1}(x_{\text{T}},0)=0$ (55)

and

$\kappa(x_{\text{R}})-1+\,\gamma_{1}(x_{\text{R}},0)=0\,.$ (56)

This set of three equations can be solved iteratively for a given value of
$\kappa_{\text{s}}$ to obtain the value of $e$ for which Eq. (54) is
satisfied.

Crossing the thin dashed green boundary to higher eccentricities causes a
“peanut” transition on the tangential critical curve. Two symmetric pairs of
inflection points appear along the minor axis, leading to a horizontal peanut-
like shape. The conditions for this transition are thus given by the common
solution of the tangential-critical-curve Eq. (51) along the minor axis

$\kappa(x_{2}/\sqrt{1-e^{2}})-1+\,\gamma_{1}(0,x_{2})=0$ (57)

and the equation for its inflection point along the minor axis

$\frac{\sqrt{1-e^{2}}}{|x_{2}|}\;\frac{{\rm d}\,\kappa}{{\rm
d}\,a}(x_{2}/\sqrt{1-e^{2}})+\frac{\partial^{2}\gamma_{1}}{\partial\,x_{1}^{2}}\,(0,x_{2})+\frac{1}{\gamma_{1}(0,x_{2})}\left(\frac{\partial\,\gamma_{2}}{\partial\,x_{1}}\right)^{\\!\\!2}(0,x_{2})=0\,,$
(58)

obtained by differentiating Eq. (51). Note that while Eq. (57) depends on the
inflection-point coordinate $x_{2}$, $e$, and $\kappa_{\text{s}}$, in Eq. (58)
$\kappa_{\text{s}}$ cancels out. It is thus advantageous to choose a value of
$e$, solve this equation numerically to obtain $x_{2}$, and use both to
directly compute the corresponding value of $\kappa_{\text{s}}$ from Eq. (57).

The thin solid purple boundary marks a similar “peanut” transition on the
radial critical curve when crossed to higher eccentricities. Here two
symmetric pairs of inflection points appear along the major axis, leading to a
vertical peanut-like shape. In this case the conditions for the transition are
given by the common solution of the radial-critical-curve Eq. (52) along the
major axis

$\kappa(x_{1})-1+\,\gamma_{1}(x_{1},0)=0$ (59)

and the equation for its inflection point along the major axis

$\frac{1}{(1-e^{2})\,|x_{1}|}\;\frac{{\rm d}\,\kappa}{{\rm
d}\,a}(x_{1})+\frac{\partial^{2}\gamma_{1}}{\partial\,x_{2}^{2}}\,(x_{1},0)+\frac{1}{\gamma_{1}(x_{1},0)}\left(\frac{\partial\,\gamma_{2}}{\partial\,x_{2}}\right)^{\\!\\!2}(x_{1},0)=0\,,$
(60)

obtained by differentiating Eq. (52). Equation (59) depends on the inflection-
point coordinate $x_{1}$, $e$, and $\kappa_{\text{s}}$, but in Eq. (60)
$\kappa_{\text{s}}$ cancels out. Here it is advantageous to choose a value of
$e$, solve this equation numerically to obtain $x_{1}$, and use both to
directly compute the corresponding value of $\kappa_{\text{s}}$ from Eq. (59).

Figure 9: Critical curves (top row) and caustics (bottom row) of ellipsoidal
NFW lenses with $(\kappa_{\text{s}},e)=(6,0.9881)$ (left column) and
$(\kappa_{\text{s}},e)=(8,0.9986)$ (right column), with tangential and radial
components marked green and purple, respectively. Both tangential critical
curves are centrally bulged; both radial critical curves are flattened along
the major axis. The inset in the bottom right panel shows a detail of the top
part of the tangential caustic, revealing the three-cusped butterfly
structure.

The topmost bold dashed black boundary marks the hyperbolic umbilic, in which
the tangential and radial components come into contact on the critical curve
as well as on the caustic. The conditions for this transition are defined by
unit convergence and zero shear. For a given value of $\kappa_{\text{s}}$ we
can obtain the semi-major axis of the unit-convergence ellipse by solving
numerically the equation

$\cfrac{\displaystyle 1-\mathcal{F}(a)}{\displaystyle
a^{2}-1}=\frac{1}{2\,\kappa_{\text{s}}}\,.$ (61)

Using the obtained value of $a$, we substitute
$x_{0\,\displaystyle{\gamma}}=a\,\sqrt{1-e^{2}}$ and solve the zero-shear-
point Eq. (46) numerically to obtain the corresponding eccentricity $e$.

One important result seen in the left panel of Fig. 8 at low values of the
convergence parameter (e.g., for $\kappa_{\text{s}}\lesssim 0.1$) is the
dominance of the lensing regime illustrated by the last column of Fig. 7, with
an outer two-cusped tangential caustic and a small inner two-cusped radial
caustic. For large values of the convergence parameter (e.g., for
$\kappa_{\text{s}}\gtrsim 1$) the lensing regime in the second column of Fig.
7 dominates, with an outer oval radial caustic and a small inner four-cusped
tangential caustic.

####  4.4.3 Lensing regimes in parameter space: intermediate and high
convergence parameter

The right panel of Fig. 8 illustrates the more complicated situation for a
broader convergence-parameter interval, $\kappa_{\text{s}}\in[0,8]$. The four
boundaries seen in the left panel continue to rise monotonically in the right
panel. However, the order of the two “peanut” transitions on the critical
curve changes at $\kappa_{\text{s}}\approx 3$. For higher values of
$\kappa_{\text{s}}$ the radial “peanut” appears at a lower eccentricity than
the tangential “peanut”.

More interestingly, two new boundaries can be seen in the right panel. The
thin dotted green boundary branching off from the thin dashed green boundary
at the green point at $(\kappa_{\text{s}},e)\approx(2.0487,0.9828)$ marks a
“bulge” transition on the tangential critical curve, in which four symmetric
pairs of inflection points appear off the principal axes. In the parameter-
space region between the dotted and the dashed thin green boundaries the
tangential critical curve thus has eight inflection points forming a central
bulge. An example corresponding to the
$(\kappa_{\text{s}},e)\approx(6,0.9881)$ parameter combination marked by the
left blue dot in Fig. 8 is shown in the left column of Fig. 9. The size and
extent of the bulge varies substantially within its region of parameter space,
from a mere bump near the minor axis to lemon-shaped critical curves. Similar
bulged critical curves have been found in the singular power-law-ellipsoid
lens model (Tessore & Metcalf, 2015) at high ellipticity.

Note that the described “bulge” boundary peaks at
$(\kappa_{\text{s}},e)\approx(2.296,0.9881)$, which also identifies the
highest-eccentricity parameter combination within the range of Fig. 8 that has
an entirely convex tangential critical curve. To the right from the green
point in Fig. 8 the dashed thin green boundary describes a different version
of the “peanut” transition. When crossed from lower to higher eccentricities
it corresponds to the disappearance of two symmetric pairs of inflection
points along the minor axis, turning the bulged critical curve into a peanut-
shaped one.

Equations describing the condition for the “bulge” transition can be obtained
by Taylor-expanding the tangential-critical-curve Eq. (51) around its general
point $\bm{x}_{\text{T}0}$, and setting the quadratic and cubic terms of the
expansion along its tangent equal to zero. We used an alternative approach
based on the image-plane analysis of critical curves, described in Appendix C.
For a given eccentricity this method provides a tool for visualizing the
existence of the transition and it indicates the corresponding value of
$\kappa_{\text{s}}$.

Figure 10: High-$\kappa_{\text{s}}$ lensing-regime boundaries of the
ellipsoidal NFW lens mapped in its parameter space, as described in Sect.
4.4.3, plotted here for ellipticities $1-\sqrt{1-e^{2}}>0.65$ and for
convergence parameters up to $\kappa_{\text{s}}=140$. Notation same as in Fig.
8, with the additional thin dashed purple boundary marking the minor-axis
radial-peanut transition.

The second new boundary, plotted in bold dashed green at the top right corner
of the right panel of Fig. 8 indicates two butterfly metamorphoses occurring
at the minor-axis cusps of the tangential caustic. Inside the region to the
right of the boundary, these cusps are replaced by small three-cusped self-
intersecting butterfly structures. An example corresponding to the
$(\kappa_{\text{s}},e)\approx(8,0.9986)$ parameter combination marked by the
right blue dot is shown in the right column of Fig. 9, with an inset showing a
blow-up of the vicinity of the top cusp revealing the butterfly. These caustic
tangles lead to very localized changes in the number of images: five images in
the butterfly “wings”, and seven images in the butterfly “proboscis” (given
that the structure lies inside the oval radial caustic). The described
butterfly metamorphosis occurs for convergence-parameter values
$\kappa_{\text{s}}\gtrapprox 7.185$. We identified and located this transition
by studying the changes in the number of cusps using the image-plane analysis
described in Appendix C. Note that while the vertical extent of the butterfly
features in the example in Fig. 9 is merely $\approx 0.06\%$ of the vertical
extent of the tangential caustic, their prominence grows with increasing
convergence parameter so that at $\kappa_{\text{s}}=100$ and higher
eccentricities it exceeds $10\%$ of the vertical caustic extent.

The structure of the parameter-space map at higher convergence-parameter
values, further to the right of the range of Fig. 8 is shown in Fig. 10 for
$\kappa_{\text{s}}\in[0,140]$. We show here only the region of higher
ellipticity, $1-\sqrt{1-e^{2}}\geq 0.65$, corresponding to eccentricities
$e\gtrsim 0.937$. No boundaries occur at lower eccentricity for
$\kappa_{\text{s}}>5$ within the range of Fig. 10. The plot illustrates the
changes of the boundaries from Fig. 8 as well as the appearance of another new
boundary. The thin dashed purple boundary extending to the right from
$\kappa_{\text{s}}\approx 41.6$ corresponds to another “peanut” transition on
the radial critical curve. In this one two symmetric pairs of inflection
points appear at the minor axis of the mass distribution, so that to the right
of the boundary the radial critical curve is pinched in vertically.

Regarding the sequence of the boundaries, a closer inspection reveals there
are only two boundaries that gradually intersect others within the range of
Fig. 10; the remaining boundaries do not cross each other. The solid black
cusp-piercing boundary intersects the green dotted “bulge” boundary at
$\kappa_{\text{s}}\approx 7.37$, the bold green dashed butterfly boundary at
$\kappa_{\text{s}}\approx 26.1$, and the purple dashed radial “peanut”
boundary at $\kappa_{\text{s}}\approx 121$. The solid purple radial “peanut”
boundary crosses the thin green dashed “peanut” boundary at
$\kappa_{\text{s}}\approx 2.98$, the bold green dashed butterfly boundary at
$\kappa_{\text{s}}\approx 7.62$, and the dashed purple radial “peanut”
boundary $\kappa_{\text{s}}\approx 78.3$. Over most of the range of Fig. 10
the dotted green “bulge” boundary represents the first transition when
increasing the eccentricity from zero. The boundary drops to a minimum at
$(\kappa_{\text{s}},e)\approx(62.6,0.950)$, then rises slowly with increasing
$\kappa_{\text{s}}$. In general, note the entirely different sequence of nine
transitions with increasing eccentricity at the right edge of Fig. 10 in
comparison with the generic low-convergence sequence of four transitions as
seen in the left panel of Fig. 8.

We haven’t explored the parameter space for higher values of
$\kappa_{\text{s}}$ beyond the range of Fig. 10 comprehensively. Overall, the
boundaries tend to level out, their rate of change declines with increasing
$\kappa_{\text{s}}$. The lower part of the new purple dashed peanut boundary
from Fig. 10 declines monotonically asymptotically to $e\approx 0.9945$. The
declining lower part of the bold green dashed butterfly boundary eventually
crosses the slowly rising dotted green bulge boundary at
$(\kappa_{\text{s}},e)\approx(8000,0.956)$. For higher values of
$\kappa_{\text{s}}$ the butterfly boundary marks the first transition to occur
when increasing the eccentricity from zero. However, at such high values of
the convergence parameter these results purely illustrate the changing
properties of the simple lens model.

Figure 11: Radii of characteristic curves of the spherical NFW lens, plotted
as a function of the convergence parameter up to $\kappa_{\text{s}}=2.1$ (left
panel) and up to $\kappa_{\text{s}}=8$ (right panel). Included are the radii
of the tangential ($x_{\text{T}}$) and radial ($x_{\text{R}}$) critical
curves, of the tangential ($y_{\text{T}}$) and radial ($y_{\text{R}}$)
caustics, and of the unit-convergence circle ($x_{0}$). The radius of the
radial limiting curve ($x_{\text{R}\infty}\approx 1.32$) is marked along the
right edge of the right panel. The five vertical dashed lines mark the values
of $\kappa_{\text{s}}$ for which critical-curve and caustic dimensions are
plotted for the ellipsoidal NFW lens in Fig. 12 and Fig. 14, respectively.

####  4.4.4 Critical-curve and caustic radii in the circularly symmetric case

In Sect. 4.4.2 and Sect. 4.4.3 we described the range of possible shapes of
the critical curves and caustics. In this and the following two sections we
study the linear dimensions and areas of their tangential and radial
components.

For initial orientation we plot in Fig. 11 the radii of the critical curves
and caustics in the $e=0$ case for convergence-parameter intervals
$\kappa_{\text{s}}\in[0,2.1]$ (left panel) and $\kappa_{\text{s}}\in[0,8]$
(right panel), matching the horizontal extent of the panels in Fig. 8. Note
that this case includes purely spherical NFW lenses as well as ellipsoidal NFW
lenses that are circular in projection, as discussed in the text following Eq.
(7). The plots are extensions of Fig. 3 from Karamazov et al. (2021) to higher
values of $\kappa_{\text{s}}$.

As described and illustrated in Karamazov et al. (2021), for convergence-
parameter values $\kappa_{s}\lessapprox 0.2$ all the radii shrink
exponentially fast to zero in the limit $\kappa_{\text{s}}\to 0$, and all
except the tangential-caustic radius increase monotonically with increasing
$\kappa_{\text{s}}$. The tangential and radial critical-curve radii maintain
the ordering $x_{\text{T}}>x_{0}>x_{\text{R}}$, where $x_{0}$ is the radius of
the unit-convergence circle defined implicitly by

$\kappa(x_{0})=1\,.$ (62)

The tangential critical-curve radius is larger than the scale radius of the
halo for $\kappa_{\text{s}}\gtrapprox 0.815$. The radial critical-curve radius
exceeds the scale radius for $\kappa_{\text{s}}\gtrapprox 9.44$ (just beyond
the range of the right panel) and for $\kappa_{\text{s}}\to\infty$ it reaches
the asymptotic value $x_{\text{R}\infty}\approx 1.32$ marked by the orange
level mark at the right edge of the right panel of Fig. 11. The existence of
this limiting radius is a special case of the more general result that for any
value of the eccentricity $e$ the radial critical curve approaches a finite
curve in the limit $\kappa_{\text{s}}\to\infty$, as demonstrated in Appendix
C.

The tangential caustic is limited to the single point behind the lens,
$y_{\text{T}}=0$, for any value of $\kappa_{\text{s}}$. The radius of the
radial caustic is initially smaller than both critical-curve radii but it
increases faster with $\kappa_{\text{s}}$, so that it exceeds the radial
critical-curve radius for $\kappa_{\text{s}}\gtrapprox 0.544$, the tangential
critical-curve radius for $\kappa_{\text{s}}\gtrapprox 4.28$, and it grows
linearly with $\kappa_{\text{s}}$ asymptotically.

####  4.4.5 Critical-curve dimensions in the general ellipsoidal case

In Fig. 12 we present the sizes of the tangential and radial critical curves
for five sample values of the convergence parameter,
$\kappa_{\text{s}}\in\\{0.125,0.25,0.5,1,2\\}$. Each of these values
corresponds to a row in Fig. 12 with both panels showing plots as a function
of ellipticity $1-\sqrt{1-e^{2}}\in[0,1]$, with the corresponding eccentricity
values marked along the nonlinear top axis. Plotted in the left column are:
the horizontal (solid) and vertical (dashed) “radii” (defined as the
coordinates of the critical-curve intercepts with the major and minor axes of
the mass distribution); the vertical maximum extent (dotted; defined as half
of the maximum vertical coordinate range of the critical curve); the semi-
minor axis of the unit-convergence ellipse (dashed black) and the semi-major
axis of the 1/2-convergence ellipse (solid black) for orientation. Plotted in
the right column are the areas enclosed by the critical curves, in units of
$\tilde{a}_{\text{s}}^{2}$, squared projected scale semi-major axis of the
halo.

Figure 12: Linear dimensions and areas of critical curves of the ellipsoidal
NFW lens for five values of the convergence parameter $\kappa_{\text{s}}$
indicated at the left side. Plotted as a function of ellipticity (bottom axis)
and eccentricity (top axis) in the left column are: tangential-critical-curve
horizontal (solid green) and vertical (long-dashed green) radii plus maximum
vertical extent (short-dashed green); radial-critical-curve horizontal (solid
purple) and vertical (long-dashed purple) radii; unit-convergence ellipse
semi-minor axis (dashed black); convergence$-1/2$ ellipse semi-major axis
(solid black). The coincidence point of the purple, green, and black dashed
curves marks the hyperbolic umbilic; the branching point of the short-dashed
from the long-dashed green curve marks the tangential-peanut transition. The
right column shows the areas enclosed by the tangential (green) and radial
(purple) critical curves plotted as a function of ellipticity and
eccentricity.

The same five convergence-parameter values are also marked for convenience in
Fig. 11 and Fig. 8 by the thin vertical dashed lines. The former may serve as
a reference for the $e=0$ starting points along the left axes of the plots in
the left column of Fig. 12. The latter is useful for keeping track of the
transitions between different lensing regimes encountered while increasing the
ellipticity from zero to one.

Looking first at the green tangential critical curve at
$\kappa_{\text{s}}=0.125$ in the left column of the first row, at the left
axis its horizontal radius starts from the spherical-case radius
$x_{\text{T}}$ (Fig. 11) and increases monotonically with ellipticity,
extending more than four times before reaching the length of the semi-major
axis of the $\kappa=1/2$ contour at the $e=1$ right axis. This $e\to 1$ limit
is valid for arbitrary $\kappa_{\text{s}}$, as can be shown analytically from
Eq. (51). The vertical radius also starts from the spherical-case radius at
the left axis, decreasing monotonically until it reaches the umbilic point at
ellipticity $1-\sqrt{1-e^{2}}\approx 0.23$ (see the umbilic boundary crossing
in Fig. 8). Here it is equal to the semi-minor axis of the unit-convergence
ellipse and the vertical radius of the radial critical curve. For higher
ellipticities, the vertical radius of the tangential critical curve rises
slowly before declining to zero at $e=1$.

The dotted green curve of the maximum vertical extent branches off from the
green dashed curve at ellipticity $1-\sqrt{1-e^{2}}\approx 0.086$,
corresponding to the tangential peanut boundary in Fig. 8. It remains nearly
constant before dropping to zero at $e=1$. The plot of the enclosed area in
the right column shows an initial mild decrease with ellipticity, followed by
a prominent increase peaking at $1-\sqrt{1-e^{2}}\approx 0.73$ before dropping
to zero at $e=1$. Overall, the tangential critical curve changes with
increasing ellipticity from circular to increasingly horizontally elongated,
shrinking to a line segment spanning the major axis of the $\kappa=1/2$
ellipse in the unit ellipticity / unit eccentricity limit.

At $\kappa_{\text{s}}=0.25$ the overall change in the tangential critical
curve is similar to the first row, with several differences. The increase in
horizontal radius is less dramatic, the vertical radius as well as the maximum
vertical extent both decrease monotonically in the full range of
ellipticities, and the umbilic and peanut transitions occur at higher
ellipticity as seen in Fig. 8. More notably, the area of the tangential
critical curve decreases monotonically, with a mere bump around intermediately
high ellipticities. Similar changes can be seen in the next row, at
$\kappa_{\text{s}}=0.5$.

At $\kappa_{\text{s}}=1$ the horizontal radius stays nearly constant,
increasing slightly before decreasing to the semi-major axis of the
$\kappa=1/2$ contour at $e=1$. The enclosed area decreases nearly linearly
with ellipticity. In the $\kappa_{\text{s}}=2$ case in the bottom row the
horizontal radius reverses the trend seen in the top rows, decreasing
monotonically with ellipticity. The tangential critical curve thus shrinks
monotonically along both axes from the spherical case with an increasingly
horizontally elongated shape.

The purple horizontal and vertical radii of the radial critical curve start
from the spherical-case radius $x_{\text{R}}$ (Fig. 11) at the left axis in
the left column of Fig. 12. The horizontal radius decreases monotonically with
ellipticity to zero at $e=1$ for all values of the convergence parameter
$\kappa_{\text{s}}$. The vertical radius initially increases with ellipticity
in the first two rows until it reaches the umbilic point described above,
after which it decreases to zero at $e=1$. At higher values of
$\kappa_{\text{s}}$, from the third row downward, the vertical radius
decreases monotonically even before reaching the umbilic point.

The maximum horizontal extent of the radial critical curve branches off from
the horizontal radius at ellipticities larger than the radial peanut boundary
in Fig. 8. However, the relative increase is so small that the difference
between the two could not be seen in Fig. 12. The area enclosed by the radial
critical curve plotted in the right column decreases monotonically with
ellipticity to zero at $e=1$ in all rows. The area practically vanishes at the
scales of the plots already at lower ellipticities: for
$1-\sqrt{1-e^{2}}\gtrapprox\\{0.4,0.6,0.7,0.8,0.8\\}$ in the first to fifth
row, respectively.

At lower values of $\kappa_{\text{s}}$ (first three rows), the radial critical
curve changes with increasing ellipticity from circular to increasingly
vertically elongated. The elongation peaks at the umbilic, after which it
decreases as the two radii converge while the radial critical curve shrinks to
the central point in the unit ellipticity limit. At $\kappa_{\text{s}}=1$, the
elongation stays initially negligible and the critical curve maintains nearly
circular shape until $1-\sqrt{1-e^{2}}\approx 0.4$ corresponding to
eccentricity $e\approx 0.8$. At higher convergence parameters, as seen in the
$\kappa_{\text{s}}=2$ example in the last row, the initial orientation is
reversed, with the radial curve becoming increasingly horizontally elongated
with increasing ellipticity. The orientation switches to vertical elongation
just before the umbilic, followed by convergence of the two radii as the curve
shrinks to a point in the unit ellipticity limit.

Figure 13: Linear dimensions of the radial limiting curve of the ellipsoidal
NFW lens: horizontal (solid) and vertical (dashed) radius plotted as a
function of ellipticity and eccentricity. The value of the radius
$x_{\text{R}\infty}\approx 1.32$ in the $e=0$ limit is marked along the left
axis. With increasing $\kappa_{\text{s}}$ the dimensions of the radial
critical curve approach but never exceed these values.

With increasing $\kappa_{\text{s}}$ the growth of the radial critical curve
slows down as it converges to a radial limiting curve, as mentioned above and
described in Appendix C. The horizontal and vertical “radii” of this curve are
plotted in Fig. 13 as a function of ellipticity. The radii start at zero
ellipticity from the value $x_{\text{R}\infty}\approx 1.32$ marked in the
right panel of Fig. 11. Both decrease monotonically with ellipticity: the
horizontal reaches $x_{1}=1$; the vertical $x_{2}=0$ at $e=1$. The radial
limiting curve shrinks from a circle at $e=0$ to an increasingly horizontally
elongated shape, converging at $e=1$ to a horizontal line segment extending
one scale radius to either side of the halo center. Note that although the
horizontal and vertical radii of radial critical curves converge to the lines
plotted in Fig. 13 as $\kappa_{\text{s}}\to\infty$, their horizontal radius
always shrinks to zero in the limit $e\to 1$, as mentioned above.

Notice also the remarkable similarity of Fig. 13 and the left panel of Fig. 3.
The plotted curves are in fact identical, bringing us back to the properties
of the deflection angle described in Sect. 4.2: the points of maximum
deflection $(0,\pm x_{2\text{max}})$ and deflection-saddle points $(\pm
x_{1\text{sad}},0)$ lie on the radial limiting curve and define its vertical
and horizontal radii, respectively. Radial critical curves thus approach these
points along the axes as $\kappa_{\text{s}}\to\infty$, but never extend beyond
them.

Figure 14: Linear dimensions and areas of caustics of the ellipsoidal NFW lens
for five values of the convergence parameter $\kappa_{\text{s}}$ indicated at
the left side. Plotted as a function of ellipticity (bottom axis) and
eccentricity (top axis) in the left column are: tangential-caustic horizontal
(solid green) and vertical (dashed green) radii; radial-caustic horizontal
(solid purple) and vertical (dashed purple) radii; convergence$-1/2$ ellipse
semi-major axis (solid black). The intersection of the purple and green dashed
curves marks the hyperbolic umbilic; the intersection of the purple and green
solid curves marks the cusp-piercing transition. The middle column shows the
areas enclosed by the tangential (green) and radial (purple) caustics and the
total area inside the caustics (gray), which is equal to the radial-caustic
area at eccentricities lower than the cusp piercing and equal to the
tangential-caustic area at eccentricities higher than the hyperbolic umbilic.
The right column includes similar plots showing the caustic areas yielding
three (blue), five (magenta), and at least three (gray) point-source images.

####  4.4.6 Caustic dimensions in the general ellipsoidal case

Figure 14 is the source-plane counterpart of Fig. 12, including similarly
designed plots describing the sizes of the tangential and radial caustics for
the same set of values of the convergence parameter. Plotted in the left
column are: the horizontal (solid) and vertical (dashed) “radii” (defined as
the coordinates of the caustic intercepts with the major and minor axes of the
mass distribution); the semi-major axis of the $\kappa=1/2$ ellipse (solid
black) for orientation. Plotted in the central column are the areas enclosed
by the caustics (including the area enclosed by the combined caustic plotted
in gray), in units of $\tilde{a}_{\text{s}}^{2}$, squared projected scale
semi-major axis of the halo. Plotted in the right column are the areas
separated here into those enclosed by a single caustic, corresponding to
point-source positions yielding three images (blue); those enclosed by both
caustics, yielding five images (magenta); and the same total area as in the
central column, yielding multiple images (gray).

Looking first at the green tangential caustic in any row of the left column,
both radii start at zero at the left axis corresponding to the point-like
nature of the caustic in the spherical case (Fig. 11). Its horizontal radius
increases monotonically with ellipticity, crosses the solid purple horizontal
radius of the radial caustic at the piercing point, and reaches the length of
the semi-major axis of the $\kappa=1/2$ contour at the $e=1$ right axis, the
same limit reached by the horizontal radius of the tangential critical curve,
as seen in Fig. 12. The vertical radius first increases with ellipticity at
the same rate as the horizontal radius, gradually slows down, reaches a peak,
and decreases to zero at $e=1$. En route it crosses the dashed purple vertical
radius of the radial caustic at the umbilic point. In the rows from top to
bottom the piercing and umbilic transitions gradually occur at higher
ellipticities, as seen in Fig. 8.

The plot of the enclosed area in the central column shows a strong increase
from zero, a peak at high ellipticity and a final drop to zero at $e=1$.
Overall, the tangential caustic changes with increasing ellipticity from
point-like to initially symmetrically astroidal, becoming increasingly
horizontally elongated, piercing the radial caustic along the major axis,
losing its two cusps along the minor axis at the umbilic, then shrinking
vertically to a line segment spanning the major axis of the $\kappa=1/2$
ellipse in the unit ellipticity / unit eccentricity limit.

The purple horizontal and vertical radii of the radial caustic start from the
spherical-case radius $y_{\text{R}}$ (Fig. 11) at the left axis in the left
column of Fig. 14. The horizontal radius decreases monotonically with
ellipticity, crossing the solid green horizontal radius of the tangential
caustic at the piercing point, and decreasing to zero at $e=1$ for all values
of the convergence parameter $\kappa_{\text{s}}$. The vertical radius
initially increases weakly with ellipticity in the first three rows, then
decreases to zero at $e=1$, crossing the dashed green vertical radius of the
tangential caustic at the umbilic point en route. At higher values of
$\kappa_{\text{s}}$, in the two lowest rows, the vertical radius decreases
monotonically in the full range of ellipticities.

The area enclosed by the radial caustic plotted in the central column
decreases monotonically with ellipticity to zero at $e=1$ in all rows. The
area practically vanishes at the scales of the plots already at lower
ellipticities: for $1-\sqrt{1-e^{2}}\gtrapprox\\{0.4,0.65,0.8,0.85,0.9\\}$ in
the first to fifth row, respectively. Overall, the radial caustic changes with
increasing ellipticity from circular to increasingly vertically elongated,
gaining two cusps along the minor axis at the umbilic. Both radii shrink to
zero, with the horizontal radius shrinking faster, so that the caustic first
resembles a short vertical line segment that shrinks to zero length in the
unit ellipticity limit.

The gray total area in the central column is equal to the radial caustic area
at low ellipticities. It branches off at the piercing point and merges with
the tangential caustic area at the umbilic point. The central column shows the
dramatic change in dominance of the tangential and radial components by area.
At any value of $\kappa_{\text{s}}$ the radial caustic dominates at low
ellipticities, while the tangential caustic dominates at high ellipticities.
The ellipticity at which the tangential caustic overtakes the radial caustic
shifts dramatically with the convergence parameter from
$1-\sqrt{1-e^{2}}\approx 0.14$ in the top row to $1-\sqrt{1-e^{2}}\approx
0.75$ in the bottom row. At low convergence parameters, the total area of the
caustic is absolutely dominated by the tangential caustic at all but the
lowest ellipticities, while at high convergence parameters, the total area is
absolutely dominated by the radial caustic at all but the highest
ellipticities. This result adds an explicit quantification to the related
comment at the end of the parameter-space analysis at the end of Sect. 4.4.2.

In the right column of Fig. 14, the magenta curve representing the source-
plane area yielding five point-source images is equal to the tangential
caustic area at ellipticities lower than the cusp-piercing transition, and it
is equal to the radial caustic area at ellipticities higher than the umbilic
transition. Thus, it drops to zero at low as well as at high ellipticities.
Between these two regimes the five-image area peaks, with the position of the
peak shifting from low to high ellipticity with increasing convergence
parameter $\kappa_{\text{s}}$. In the first three rows there even is a narrow
interval of ellipticities around the peak at which the five-image area
narrowly exceeds the three-image area marked by the blue curve. At higher
$\kappa_{\text{s}}$ in the bottom two rows the five-image area is smaller than
the three-image area at all ellipticities.

At a more basic level, the gray curve in the central and right columns vividly
illustrates the change in the overall probability of multiple imaging. At low
values of $\kappa_{\text{s}}$ the probability is negligible at low ellipticity
and high at high ellipticity. At high $\kappa_{\text{s}}$ the opposite is
true: the probability is high at low ellipticity and negligible at high
ellipticity. At the same time we point out the varying scale on the vertical
axes marked in the central column, illustrating the dramatic increase of the
probability with the convergence parameter.

###  4.5 Image properties

The analysis of the critical curves and caustics in Sect. 4.4 shows that the
ellipsoidal NFW lens may produce 1, 3, 5, or 7 images of a point source lying
off the caustic. Here we describe briefly the geometrical properties of the
individual images. In Sect. 4.5.1 we utilize convergence–shear diagrams to
study the range of image-magnification and flattening values achievable by a
lens with a given combination of $\kappa_{\text{s}}$ and $e$. In Sect. 4.5.2
we illustrate the patterns of image orientation, as determined by the
combination of phase and convergence. In Sect. 4.5.3 we present a few examples
of interesting image configurations.

####  4.5.1 Magnification and flattening

The convergence–shear diagram (hereafter CS diagram; Karamazov & Heyrovský,
2022) is based on a simple plot with convergence $\kappa$ along the horizontal
axis and shear $\gamma$ along the vertical axis, as seen in the examples in
Fig. 15. In such a plot, lines of constant magnification
$\mu=1/\mathrm{det}\,J$ are hyperbolae

$(1-\kappa)^{2}-\gamma^{2}=1/\mu$ (63)

centered around the point with the umbilic combination
$(\kappa,\gamma)=(1,0)$. Their common diagonal asymptotes (marked solid bold
black in Fig. 15) correspond to $|\mu|=\infty$, and the absolute magnification
declines outward from the diagonals. The bold black hyperbolae and the point
at the origin correspond to $|\mu|=1$; the parts of the diagram above and to
the right of the region between these hyperbolae correspond to demagnified
images. The $(\kappa,\gamma)$ combinations in the regions under the bold
diagonals yield positive magnification (i.e., positive-parity images); in the
complementary region above the bold diagonals the magnification is negative
(i.e., negative-parity images).

Figure 15: Convergence–shear diagrams for the ellipsoidal NFW lens with
$\kappa_{\text{s}}=0.25$ and eccentricities $e\in\\{0,0.6,0.811,0.9\\}$ as
marked in the panels. In the top row with $e=0$ the green curve marks the
$(\kappa,\gamma)$ combinations occurring in the lens plane as a function of
radius $a$ marked by the green tick marks. In each of the three lower rows the
magenta band marks the $(\kappa,\gamma)$ combinations occurring in the lens
plane, with its vertical cross-sections indicating the shear range along the
ellipse with semi-major axis $a$. The red and blue curves mark the shear along
the major and minor axis, respectively. For reading the image properties from
the diagrams, see Sect. 4.5.1.

The flattening, which is defined as the ellipticity of the image of a small
circular source,

$f=1-\min\left({\left|\frac{1-\kappa-\gamma}{1-\kappa+\gamma}\right|,\left|\frac{1-\kappa+\gamma}{1-\kappa-\gamma}\right|}\right)\,,$
(64)

is constant along the four straight lines

$\begin{matrix}f\,\kappa+(2-f)\,\gamma-f\\\ (2-f)\,\kappa+f\,\gamma+f-2\\\
(2-f)\,\kappa-f\,\gamma+f-2\\\
f\,\kappa-(2-f)\,\gamma-f\end{matrix}\;\begin{matrix}=0\\\ =0\\\ =0\\\
=0\end{matrix}$ (65)

extending radially from the point $(\kappa,\gamma)=(1,0)$. The sets of these
lines are shown for $f\in\\{0.25,0.5,0.75\\}$ in Fig. 15, where they are
marked gray. The first two lines appear in the diagram in the $\kappa<1$
region, where the image is extended parallel to the phase. The third and
fourth lines appear in the $\kappa>1$ region, where the image is extended
perpendicular to the phase. In the case of undistorted images ($f=0$) the four
lines reduce to the horizontal axis, $\gamma=0$, and the vertical line at
$\kappa=1$. In the case of extremely distorted images ($f=1$) the four lines
reduce to the two solid black diagonal lines, $\gamma=1-\kappa$ and
$\gamma=\kappa-1$, with the former corresponding to the tangential critical
curve and the latter to the radial critical curve.

Lines of constant image-scaling factor parallel to the phase
$\lambda_{\parallel}$ correspond to the set of lines

$1-\kappa-\gamma=1/\lambda_{\parallel}$ (66)

parallel to the tangential-critical-curve diagonal, along which
$|\lambda_{\parallel}|=\infty$. The absolute value declines on both sides,
with positive values under and negative values above the diagonal. The point
at the origin corresponds to $\lambda_{\parallel}=1$; the dashed black
$\gamma=2-\kappa$ line from the top left of the diagram corresponds to
$\lambda_{\parallel}=-1$. Lines of constant image-scaling factor perpendicular
to the phase $\lambda_{\perp}$ correspond to the set of lines

$1-\kappa+\gamma=1/\lambda_{\perp}\,$ (67)

parallel to the radial-critical-curve diagonal, along which
$|\lambda_{\perp}|=\infty$. The absolute value declines on both sides, with
positive values above and negative values under the diagonal. The dashed black
$\gamma=\kappa$ line from the origin corresponds to $\lambda_{\perp}=1$; the
dashed black $\gamma=\kappa-2$ line at the right side of the diagram
corresponds to $\lambda_{\perp}=-1$.

For any given point $(\kappa,\gamma)$ in the CS diagram one may thus directly
read off the corresponding magnification $\mu$ (including parity as the sign
of $\mu$), flattening $f$, orientation with respect to the phase, and the two
scaling factors $\lambda_{\parallel}$ and $\lambda_{\perp}$ of a small image
appearing in the lens plane at such a convergence–shear combination. For a CS
diagram of a particular gravitational lens we utilize the plot described
above, and mark in it the range of $(\kappa,\gamma)$ combinations occurring in
the full image plane of the lens. In this way we may easily identify and
visualize the above-mentioned properties for the full range of images possibly
formed by the lens.

In Fig. 15 we present four CS diagrams for an ellipsoidal NFW lens with
convergence parameter $\kappa_{\text{s}}=0.25$ and eccentricities
$e\in\\{0,0.6,0.811,0.9\\}$. The same parameter combinations plus the
additional $(\kappa_{\text{s}},e)=(0.25,0.3)$ were used for the critical-curve
and caustic examples in Fig. 7; all five are marked by the red dots in the
parameter-space maps in Fig. 8.

In the zero-eccentricity case (top panel of Fig. 15) the semi-major axis is
equal to the radial distance, $a=x$. The convergence and shear are purely
radial functions, $\kappa(a)$ and $\gamma(a)$, so that all their combinations
trace the green curve parameterized by the semi-major axis / radial distance,
marked along the curve by logarithmically spaced tick marks. Starting from the
asymptotic region $a\to\infty$ at the $(\kappa,\gamma)=(0,0)$ origin and
progressing toward the halo center, the magnification $\mu$ initially
increases from 1; the flattening $f$ increases from 0; $\lambda_{\parallel}$
increases from 1; $\lambda_{\perp}$ first decreases from 1 before reversing
the trend and increasing, crossing the dashed line $\lambda_{\perp}=1$ at
$a=x_{\text{R}\infty}\approx 1.32$. Progressing further inward (to the right
along the green curve), we reach the tangential critical curve at $a\approx
0.17$. As we cross it and enter the negative-parity region of the diagram, the
magnification diverges $\mu\to\infty$ and changes sign to $\mu=-\infty$; the
flattening reaches its maximum $f=1$ and starts to decrease;
$\lambda_{\parallel}$ exhibits the same behavior as $\mu$, diverging
$\lambda_{\parallel}\to\infty$ and changing sign to
$\lambda_{\parallel}=-\infty$; $\lambda_{\perp}\approx 2.14$ at the tangential
critical curve and keeps increasing. At $a\approx 0.103$ the convergence
exceeds $\kappa=1$ and the image distortion switches orientation from parallel
to the phase to perpendicular to the phase; the flattening correspondingly
reaches its minimum $f=0$ and starts to increase.

At $a\approx 0.102$ the absolute magnification reaches $|\mu|\approx 17$, its
minimum value for negative-parity images. As we cross the radial critical
curve at $a\approx 0.062$ and enter the inner positive-parity region, the
magnification diverges $\mu\to-\infty$ and changes sign to $\mu=\infty$; the
flattening reaches its maximum $f=1$ and starts to decrease;
$\lambda_{\parallel}\approx-2.03$ at the radial critical curve and keeps
increasing; $\lambda_{\perp}$ exhibits similar behavior as $\mu$ with a
flipped sign, diverging $\lambda_{\perp}\to\infty$ and changing sign to
$\lambda_{\perp}=-\infty$. At $a\approx 0.022$ the dashed line
$\lambda_{\parallel}=-1$ is crossed, indicating that images closer to the halo
center are shrunk in the direction parallel to the phase. At $a\approx 0.013$
we cross the magnification $\mu=1$ hyperbola, indicating that images closer to
the halo center are demagnified. At $a\approx 0.0082$ the dashed line
$\lambda_{\perp}=-1$ is crossed, indicating that images closer to the halo
center are shrunk even in the direction perpendicular to the phase. Beyond the
right edge of the plot, closer to the halo center, the convergence keeps
increasing $\kappa\to\infty$ and the shear reaches its central value
$\gamma=\kappa_{\text{s}}=0.25$. All the followed quantities tend to zero:
$\mu\searrow 0$, $f\searrow 0$, $\lambda_{\parallel}\nearrow 0$,
$\lambda_{\perp}\nearrow 0$.

In the case of non-zero eccentricity the shear does not have the elliptical
symmetry of the convergence, as described in Sect. 4.3. As a result, the
region in the CS diagram corresponding to such a lens is two-dimensional, as
indicated by the magenta-colored areas in the lower three panels of Fig. 15.
Nevertheless, the connection between the semi-major axis and the convergence
$\kappa(a)$ remains unchanged and independent of eccentricity. Hence, vertical
lines in the diagram correspond to ellipses with semi-major axes $a$ and the
eccentricity $e$ of the model. The vertical span of the colored area is thus
given by the variation of the shear along these ellipses. The interval of the
shear values at semi-major axis $a$ is determined by finding the minimum and
maximum values of $\gamma(a\,\cos{\theta},a\,\sqrt{1-e^{2}}\,\sin{\theta})$.
Due to symmetry it is sufficient to vary the parameter $\theta$ from 0 along
the major axis to $\pi/2$ along the minor axis.

In the lower three panels of Fig. 15, the shear along the major and minor axes
is marked by the red and blue lines, respectively. The green line with the
semi-major-axis ticks corresponds to $\theta=\pi/4$, a line dividing the
quarter-ellipse into equal-area parts. We see that in the inner region of the
halo covering nearly the full extent of the plotted CS diagrams, the shear is
lowest along the minor axis and highest along the major axis. In the outer
region extending asymptotically to infinity, the order is reversed: minimum
shear occurs along the major axis and maximum along the minor axis. In the
plotted examples, the transition between the inner and outer regions occurs
between $a=1$ and $a=2$. An inspection of the shear plots in Fig. 5 and Fig. 4
confirms these results.

The next feature to notice in the non-zero-eccentricity panels of Fig. 15 is
the presence of the zero-shear point along the minor axis, as described in
Sect. 4.3.3. For $e=0.6$ in the second panel the point lies on the $a\approx
0.065$ ellipse; for $e=0.811$ in the third panel on the $a\approx 0.103$ unit-
convergence ellipse; for $e=0.9$ in the fourth panel on the $a\approx 0.24$
ellipse. These semi-major-axis values are equal to
$x_{0\,\displaystyle{\gamma}}/\sqrt{1-e^{2}}$, and can be thus computed using
the zero-shear-point positions along the minor axis plotted in Fig. 6. The
example in the third panel with the zero-shear point lying on the unit-
convergence ellipse illustrates an umbilic combination of lens parameters, at
which the tangential and radial critical curves come into contact along the
minor axis of the $a\approx 0.103$ ellipse.

Unlike in the zero-eccentricity case, the shear in the lower three panels of
Fig. 15 diverges at the halo center. The rate of its divergence is slower than
that of the convergence $\kappa$, but it increases with increasing
eccentricity as indicated by the rising slope of the magenta band. The shear
can be approximated for $a\to 0$ and $e>0$ by the expression in Eq. (43). The
convergence near the halo center can be approximated by

$\kappa(a)=-2\,\kappa_{\text{s}}\,(\,\ln{\frac{a}{2}}+1+\frac{3}{2}\,a^{2}\,\ln{\frac{a}{2}}+\frac{5}{4}\,a^{2}\,)+\mathcal{O}(a^{4}\,\ln{a})\,,$
(68)

as shown by Karamazov et al. (2021)111The expression presented here includes a
correction in the third-order term.. The two approximations can be combined to
yield the dependence of the shear on the convergence and position along the
ellipse $\theta$ for $a\to 0$:

$\gamma\simeq\frac{\left(1-\sqrt{1-e^{2}}\,\right)^{2}}{e^{2}}\,\kappa(a)-4\,\kappa_{\text{s}}\,\frac{\sqrt{1-e^{2}}}{e^{2}}\,\ln{\frac{1+\sqrt{1-e^{2}}}{2\,\sqrt{1-e^{2}\,\sin^{2}{\theta}}}}\,.$
(69)

The factor preceding the convergence explains the slope of the band in Fig.
15: 0 for $e\to 0$, and 1 for $e\to 1$. The second term confirms that maximum
shear $\gamma_{\text{\tiny MAX}}$ occurs along the major axis ($\theta=0$) and
minimum shear $\gamma_{\text{\tiny MIN}}$ occurs along the minor axis
($\theta=\pi/2$). The vertical width of the band (at a constant $a$) is well
approximated for $a\to 0$ and $e>0$ by

$\gamma_{\text{\tiny MAX}}-\gamma_{\text{\tiny
MIN}}\simeq-2\,\kappa_{\text{s}}\,\frac{\sqrt{1-e^{2}}}{e^{2}}\,\ln{(1-e^{2})}\,,$
(70)

a slowly varying function yielding a width of $2\,\kappa_{\text{s}}$ for low
eccentricities and decreasing significantly only for eccentricities $e\gtrsim
0.9$.

Figure 16: Convergence–shear diagrams for $e=0$ NFW lenses with five values of
$\kappa_{\text{s}}$ as marked along the right axis. The conversion between
different values of $\kappa_{\text{s}}$ is obtained simply by scaling the plot
from the origin by the target/original $\kappa_{\text{s}}$ value ratio. The
three dashed orange lines connect points on the different curves corresponding
to radius $a=1$, to radius $a=0.1$, and to radius $a=0.01$.

The CS diagrams in Fig. 15 are all plotted for convergence parameter
$\kappa_{\text{s}}=0.25$. Karamazov & Heyrovský (2022) originally plotted the
CS diagram for a spherical NFW lens with the very similar value
$\kappa_{\text{s}}=0.239$ (compare their Fig. 4 with the top panel of Fig.
15). The structure and properties of the CS diagrams for any particular value
of $\kappa_{\text{s}}$ can be directly inferred from the diagrams for any
other value, such as those in Fig. 15. The expressions for the convergence and
shear are directly proportional to $\kappa_{\text{s}}$, so that
$\kappa(\bm{x})/\kappa_{\text{s}}$ and $\gamma(\bm{x})/\kappa_{\text{s}}$ are
independent of $\kappa_{\text{s}}$. In terms of the CS diagram this means that
increasing or decreasing $\kappa_{\text{s}}$ amounts to expanding the plot in
the diagram radially away from the origin, or shrinking the plot radially
toward the origin, respectively.

We illustrate this property in Fig. 16, for simplicity on the example of the
circularly symmetric case from the top panel of Fig. 15. The five green curves
correspond to convergence-parameter values
$\kappa_{\text{s}}\in\\{0.125,0.25,0.5,1,2\\}$ as marked along the right side
of the plot (the same values were used in Fig. 12 and Fig. 14). For
$\kappa_{\text{s}}=0.125$ all points (including their $a$ values) are moved
two times closer to the origin with respect to the $\kappa_{\text{s}}=0.25$
case, for $\kappa_{\text{s}}=0.5$ all points are moved two times further from
the origin, etc. For orientation, we connected points with $a=1$ on all the
curves by the top dashed orange line, and added two more dashed orange lines
similarly for all $a=0.1$ and all $a=0.01$ points.

The changes of imaging properties for different $\kappa_{\text{s}}$ in the
circularly symmetric case can be interpreted from Fig. 16 analogously to the
walk-through performed for $\kappa_{\text{s}}=0.25$ above. Here we point out
just two examples of features readily seen from the figure. First, as
$\kappa_{\text{s}}$ increases, the minimum absolute magnification of negative-
parity images decreases, so that for $\kappa_{\text{s}}=2$ we already see a
zone with demagnified images ($|\mu|<1$) between the tangential and radial
critical curves. Second, the crossing of the dashed line $\gamma=\kappa$
always occurs at $a=x_{\text{R}\infty}\approx 1.32$, independent of
$\kappa_{\text{s}}$. As discussed above, at larger radii $\lambda_{\perp}<1$
so that images are compressed perpendicular to the phase, narrower than the
source. Since the critical-curve radii increase with increasing
$\kappa_{\text{s}}$ (see Fig. 11), for $\kappa_{s}=2$ this zone reaches well
under the tangential critical curve, as seen from Fig. 16.

Changing $\kappa_{\text{s}}$ in the elliptical case (lower three panels in
Fig. 15) analogously amounts to expanding or shrinking the entire magenta area
from or to the origin of the CS diagram, respectively. For example, reducing
the convergence parameter in the second panel $\sim 2.35$ times would lead to
the umbilic configuration for $e=0.6$. Equation (69) confirms that the slope
of the inner-halo band in the CS diagram does not change with $\kappa_{s}$,
while the width of the band scales with $\kappa_{s}$ linearly, as seen from
Eq. (70).

Figure 17: Orientation of image distortion by the ellipsoidal NFW lens,
presented in color maps with contours. The parameter combinations are
indicated above the panels: the first panel shows the circularly symmetric
case with convergence parameter $\kappa_{\text{s}}=0.25$; the three right
panels show the $e=0.9$ elliptical case and a sub-umbilic, umbilic, and super-
umbilic value of $\kappa_{\text{s}}$. The respective tangential (green) and
radial (purple) critical curves are plotted in each panel. The color bar at
the right side indicates the orientation angles, the contour values spaced by
$\pi/8$, plus corresponding tilted-ellipse icons for visual guidance.

####  4.5.2 Orientation

Here we use the term “image orientation” as a shorthand for the orientation of
the image distortion, or, even more specifically, for the orientation of the
major axis of an elliptical image of a small circular source222Based on its
equality to the phase in weak lensing, Karamazov & Heyrovský (2022) use the
term “weak phase” for the same quantity.. As such the orientation is primarily
defined by the lens phase $\phi$, described in Sect. 4.3 and illustrated in
Fig. 4 and Fig. 5. The secondary factor is the convergence: outside the unit-
convergence ellipse where $\kappa<1$, the image orientation is equal to the
phase; inside the unit-convergence ellipse where $\kappa>1$, the image
orientation is perpendicular to the phase. We quantify the orientation by the
polar angle in the lens plane, measured counter-clockwise from the major axis
of the mass distribution modulo $\pi$, using the corresponding value from the
interval $[-\pi/2,\pi/2]$.

We present four examples of image-orientation maps in Fig. 17: one with
circular symmetry and convergence parameter $\kappa_{\text{s}}=0.25$, and
three with eccentricity $e=0.9$ and convergence parameters
$\kappa_{\text{s}}\in\\{0.25,0.398,0.5\\}$. These combinations are marked in
the parameter-space maps in Fig. 8 by the cyan dotted circles. The circularly
symmetric case in the first panel of Fig. 17 illustrates the $\pi/2$ switch
along the unit-convergence circle (marked in gray). Outside the circle
horizontally oriented images appear along the $x_{2}$ axis, vertically
oriented images along the $x_{1}$ axis, and along the green tangential
critical curve images are always oriented tangentially. Inside the circle
horizontally oriented images appear along the $x_{1}$ axis, vertically
oriented images along the $x_{2}$ axis, and along the purple radial critical
curve images are always oriented radially. The radial orientation is
maintained into the halo center, where all constant-orientation contours meet.

The second panel of Fig. 17 illustrates the case when the convergence
parameter is lower than the umbilic value of $\kappa_{\text{s}}$ at the given
eccentricity. In addition to the gray unit-convergence ellipse, the two zero-
shear points outside the ellipse along the minor axis clearly play an
important role here. Outside the unit-convergence ellipse horizontally
oriented images appear only outside the zero-shear points along the minor
axis, while vertically oriented images appear along the major axis plus along
the minor axis inside the zero-shear points. Along the tangential critical
curve image orientation deviates by barely more than $\pm\pi/8$ from the
vertical, so that while the images are oriented tangentially near the major
axis, near the minor axis they are oriented radially. In either half-plane,
all constant-orientation contours meet at the zero-shear points, but the
orientation changes by $\pi/2$ when crossing the point instead of staying the
same as in the case of the halo center in the first panel. Inside the unit-
convergence ellipse horizontally oriented images appear along both the major
and the minor axis, vertically oriented images are entirely absent here and
the image orientation deviates by less than $\pm\pi/8$ from the horizontal.
Along the radial critical curve the images are oriented radially near the
major axis, while near the minor axis they are oriented tangentially.

The third panel of Fig. 17 illustrates the case when the convergence parameter
is equal to the umbilic value of $\kappa_{\text{s}}$ at the given
eccentricity. Here the zero-shear points lie on the unit-convergence ellipse
and thus they form the umbilic points. Outside the unit-convergence ellipse
horizontally oriented images appear along the minor axis, while vertically
oriented images appear along the major axis. Along the tangential critical
curve image orientation deviates by barely more than $\pm\pi/4$ from the
vertical, so that while they are oriented tangentially near the major axis,
when approaching the minor axis they are oriented nearly perpendicularly to
the critical curve. Inside the unit-convergence ellipse horizontally oriented
images appear along both the major and the minor axis, vertically oriented
images are entirely absent here and the image orientation deviates by less
than $\pm\pi/4$ from the horizontal. Along the radial critical curve the
images are oriented radially near the major axis, and when approaching the
minor axis they are oriented nearly perpendicularly to the critical curve. In
either half-plane, constant-orientation contours with orientation in the
interval $[-\pi/4,\pi/4]$ meet at the umbilic points, and there is no change
in orientation when crossing the points. Note that there are no vertically
oriented images in the vicinity of the umbilic points, all images deviate by
no more than $\pm\pi/4$ from the horizontal.

The fourth panel of Fig. 17 illustrates the case when the convergence
parameter is higher than the umbilic value of $\kappa_{\text{s}}$ at the given
eccentricity. In such a case the two zero-shear points appear inside the unit-
convergence ellipse. Outside the ellipse horizontally oriented images appear
along the minor axis, and vertically oriented images appear along the major
axis. Along the tangential critical curve images are oriented tangentially
near the major axis as well as near the minor axis, although between these two
positions image orientation deviates by more than $\pi/4$ from the critical-
curve tangent. Inside the unit-convergence ellipse horizontally oriented
images appear along the major axis as well as along the minor axis inside the
zero-shear points, vertically oriented images appear only along the minor-axis
segments outside the zero-shear points. Only in their vicinity does the
orientation deviate by more than $\pm\pi/4$ from the horizontal. Along the
radial critical curve the images are oriented radially near the major axis as
well as near the minor axis. In either half-plane, all constant-orientation
contours meet at the zero-shear points, with the orientation changing by
$\pi/2$ when crossing either point.

At low eccentricities, the zero-shear points lie very close to the halo
center, as shown in Fig. 6. Even for very small values of $\kappa_{\text{s}}$
such a lens is in the regime of the fourth panel, but the adjacent pair of
zero-shear points generates an image-orientation map resembling a distorted
version of the circular case shown in the first panel.

####  4.5.3 Configuration examples

In the previous parts of Sect. 4.5 we illustrated the diversity of the imaging
properties of the ellipsoidal NFW lens, their variation in the lens plane and
within the parameter space of the lens. Rather than attempting to present an
exhaustive overview of possible image configurations, we show in Fig. 18 four
specific examples of images formed by lenses with four different
$(\kappa_{\text{s}},e)$ parameter combinations. In each pair of panels, the
left one shows the source plane, indicating the position of the source (red
and yellow circles for two source sizes) with respect to the caustics. The
right panel shows the lens plane, indicating the positions of the images
(colored by the corresponding source size) with respect to the critical
curves.

In the top left example, the lens parameters
$(\kappa_{\text{s}},e)=(0.25,0.9)$ correspond to the fourth CS diagram in Fig.
15 and the second image-orientation map in Fig. 17. The circular sources with
radii 0.004 (red) and 0.002 (yellow) are centered at the point
$(y_{1},y_{2})=(0,0.0345)$, so that the yellow source lies fully between the
radial and tangential caustics while the red source crosses both caustics. Due
to its position, the yellow source has three full separate images. The top-
most part of the red source lying outside the tangential caustic has one
image, the part between the caustics has three, and the lowest part lying
within the radial caustic has five images. One full image of the red source
lies along the minor axis above the critical curves, while all other images
are merged into one large bell-shaped macro-image under the halo center. Note
the horizontal orientation of the middle yellow image, which is tangent to the
radial critical curve, and the vertical orientation of the bottom yellow
image, which is perpendicular to the tangential critical curve. Both are in
agreement with the second image-orientation map in Fig. 17, as discussed in
Sect. 4.5.2. Similar distortions in the corresponding parts of the red macro-
image lead to its peculiar shape.

Figure 18: Examples of image configurations formed by four ellipsoidal NFW
lenses. Each of the four pairs of panels shows the position of the source with
respect to the caustics (left), and the positions of the images with respect
to the critical curves (right). Images are plotted for two source sizes (red,
yellow); in the bottom right example a third size is added (blue). Lens
parameters $(\kappa_{\text{s}},e)$ for top left example: $(0.25,0.9)$; top
right: $(0.398,0.9)$; bottom left: $(0.5,0.9)$; bottom right: $(40,0.9986)$.

In the top right example in Fig. 18, the lens parameters
$(\kappa_{\text{s}},e)=(0.398,0.9)$ represent the umbilic combination seen
also in the third image-orientation map of Fig. 17. The circular sources with
radii 0.01 (red) and 0.005 (yellow) are centered at the point
$(y_{1},y_{2})=(0,0.09)$, so that the yellow source lies fully inside both
caustics while the red source crosses both caustics. The yellow source has
five full separate images. The narrow parts of the top edge of the red source
lying outside the tangential caustic have one image, the narrow parts between
the caustics have three, and the rest of the red source lying inside the
radial caustic has five images. One full image of the red source lies along
the minor axis above the critical curves, while all other images are merged
into one large bumpy ring-like macro-image around the umbilic point (rather
than around the halo center). Note the horizontal orientation of the yellow
images above and below the umbilic point and the corresponding dominant
horizontal distortion of the top and bottom parts of the red ring-like image.
The left and right “bumps” on the ring correspond to the abrupt change in
orientation along the unit-convergence ellipse, as seen in the third image-
orientation map in Fig. 17.

In the bottom left example in Fig. 18, the lens parameters
$(\kappa_{\text{s}},e)=(0.5,0.9)$ correspond to the fourth image-orientation
map in Fig. 17. The circular sources with radii 0.012 (red) and 0.006 (yellow)
are centered at the point $(y_{1},y_{2})=(0,0.14)$, so that the yellow source
lies fully inside both caustics while the red source lies inside the radial
caustic but crosses the tangential caustic. The yellow source has five full
separate images. The narrow parts of the top edge of the red source lying
outside the tangential caustic have three images, and the rest of the red
source lying inside both caustics has five images. One full image of the red
source lies along the minor axis above the critical curves, another full image
lies inside the radial critical curve, and all other images are merged into
one large boomerang-like macro-image. Note the horizontal orientation of all
yellow images along the minor axis, with the lowest one oriented tangentially
to the tangential critical curve, and the near-perpendicular orientation of
the two side images. The top of the red boomerang-like macro-image is
distorted and flattened due to its location inside the unit-convergence
ellipse, as seen in the fourth image-orientation map in Fig. 17. In comparison
with the top left example, the image configuration is dramatically different
even though the critical-curve geometry is superficially similar.

We include the bottom right example in Fig. 18 to demonstrate the special
seven-image configuration arising when the tangential caustic has a pair of
butterflies along the minor axis (see the bottom right panel of Fig. 9 and
Sect. 4.4.3). For purposes of illustration, we use the same eccentricity and a
higher convergence parameter, $(\kappa_{\text{s}},e)=(40,0.9986)$. Instead of
plotting the full extent of the caustics, we include only the part of the
tangential caustic near the top butterfly noting that the entire plotted area
lies inside the radial caustic. In the right panel we omit the outer parts of
the tangential critical curve, which extend along the major axis out to
$x_{1}\approx\pm 12.4$. In addition to the two circular sources with radii
0.08 (red) and 0.04 (yellow) we include here a third small blue source with
radius 0.008. All are centered at the point $(y_{1},y_{2})=(0,2.382)$. The
blue source lies fully inside the lower part of the butterfly feature, the red
source encloses the entire butterfly, and the intermediate yellow source
encloses its main part without the three cusps. The blue source has seven full
separate images: one along the minor axis above the tangential critical curve,
a tiny image near the center of the radial critical curve, and five images
strung out along the lower part of the tangential critical curve alternating
on its sides (three inside, two outside). Moving the blue source across any of
the four different segments of the tangential caustic that surround it would
lead to the disappearance of one of the four pairs of neighboring images lying
on opposite sides of the tangential critical curve.

The top, left, and right parts of the yellow source lying outside the
tangential caustic have three images, the parts in the left and right wings
and the lower left and lower right parts inside the tangential caustic but
outside the butterfly have five images, and the central part adjacent to the
blue source has seven images. The red source can be dissected in a similar
manner. For the yellow as well as the red sources, one full image lies along
the minor axis above the critical curves, another tiny full image lies inside
the radial critical curve, and all other images are merged into one long arc-
like macro-image with its center of curvature lying on the minor axis above
the halo center.

##  5 Discussion

###  5.1 Relevance of high-eccentricity results

We first return to the seven-image regions described and illustrated above in
Sect. 4.5.3. While their existence in the source plane of the ellipsoidal NFW
lens is unexpected and very interesting in the context of simple lens models,
their observational significance is rather limited. First, they occur at a
combination of a very high eccentricity $e$ and a very high convergence
parameter $\kappa_{\text{s}}$. Second, even though the butterfly structures
increase in size for higher values of $\kappa_{\text{s}}$, our tests show that
at $\kappa_{\text{s}}=100$ their combined vertical extent reaches at most
$12\%$ of the vertical width of the tangential caustic (at most $7.2\%$ for
the seven-image regions of the butterflies). Due to the small extent of these
regions, a lensed source would rarely fit in entirely. Two such possibilities
are worth mentioning: a perfectly aligned quasar lensed by a galaxy-scale NFW
halo; and an exploding supernova in a galaxy lensed by a cluster- or galaxy-
scale NFW halo. For larger sources extending across the seven-image region,
high-resolution imaging resolving fine structure of the source would be
necessary to distinguish the individual partial images stretched along the
arc.

A final point worth mentioning is the rarity of butterflies. Their occurrence
requires near-perfect symmetry of the mass distribution about the butterfly
axis (in our case the minor axis of the projected ellipse). Even with a highly
symmetric lens, in most astrophysical settings one may expect perturbation by
additional shear – due to lens substructure, due to companion lenses, or due
to larger-scale mass distributions. This influence may easily break the
butterfly structure on the caustic, typically into a cusp and an adjacent
swallow-tail structure (e.g., Schneider et al., 1992). Nevertheless, the tail