15.4. The occurrence of images and the critical and caustic curves¶
15.4.1. The topology of the Fermat potential¶
We have already solved the lensing equation (15.7) for the special cases of a point mass or singular-isothermal sphere and we saw that under certain conditions, multiple images of a single source can form. We have also shown that for axially-symmetric lenses, the images lie on a single line, the line that connects the source to the center of the lens. Now, we turn to a more general discussion of where images form, how many images form, and the properties of these images.
To broadly understand when and where images form, we can use the fact that images form at extrema of the Fermat potential, as we demonstrated in Section 15.2.2 above. The Fermat potential, given in Equation (15.32), is a simple sum of (i) a quadratic function \(\left|\boldsymbol{\theta}-\boldsymbol{\beta}\right|^2/2\) and (ii) minus the lensing potential \(-\psi(\boldsymbol{\theta})\). The quadratic function has a single extremum (a minimum) at \(\boldsymbol{\theta} = \boldsymbol{\beta}\), while \(-\psi(\boldsymbol{\theta})\) generally has a single extremum (a maximum) at the center of the lens for realistic galactic mass distributions. Where and how many images form is then a function of to what extent the lens \(-\psi(\boldsymbol{\theta})\) can push up the \(\left|\boldsymbol{\theta}-\boldsymbol{\beta}\right|^2/2\) surface to create additional extrema in addition to the \(\boldsymbol{\theta} \approx \boldsymbol{\beta}\) one.
For axially-symmetric lenses, we displayed an example Fermat potential in Figure 15.5 that illustrates the different possibilities. If the source position \(\boldsymbol{\beta}\) is far from the lens’ center, \(|\psi(\boldsymbol{\theta})| \ll \left|\boldsymbol{\theta}-\boldsymbol{\beta}\right|^2/2\) and there is a single image with \(\boldsymbol{\theta} \approx \boldsymbol{\beta}\). As the source gets close to the lens’ center, \(-\psi(\boldsymbol{\theta})\) increases and, if the surface density of the lens is high enough, at some point becomes large enough to create a maximum in the Fermat potential, where an additional image forms. Because for galaxy lenses, we can assume that \(\psi(\boldsymbol{\theta}) \rightarrow 0\) as \(\boldsymbol{\theta} \rightarrow \infty\), at large \(|\boldsymbol{\theta}|\) the Fermat potential approaches \(\left|\boldsymbol{\theta}-\boldsymbol{\beta}\right|^2/2\), which increases as \(\boldsymbol{\theta} \rightarrow \infty\). Therefore, aside from the minimum at \(\boldsymbol{\theta} \approx \boldsymbol{\beta}\), another minimum has to be present on the other side of the maximum and a third image has to be present. An axially-symmetric galaxy lens therefore produces either one or three images. However, as we discuss below, we may not be able to observe all three images and for a lens with a central singularity, one of the images may be absent altogether.
For lenses that do not possess axial symmetry, more complex image configurations are possible. One possible scenario is shown in Figure 15.9, where we move a source along the line \(\beta_1 = 2\beta_2\) from an initial wide separation between source and lens center to a close alignment. We show contours of the Fermat potential for a massive elliptical galaxy centered at the origin (again modeled as a perfect ellipsoid), with the unlensed position of the source indicated by a plus sign.
[17]:
# First we define the lensing potential and fermat potentials again,
# using code from earlier sections in this chapter
from scipy import integrate
from astropy.cosmology import Planck18 as cosmo
from astropy.constants import c
from galpy.potential import evaluatePotentials
def lensing_potential(Pot,theta,zlens,zsource,zmax=numpy.inf):
DL= cosmo.angular_diameter_distance(zlens)
cosmo_fac= cosmo.angular_diameter_distance_z1z2(zlens,zsource)\
/cosmo.angular_diameter_distance(zsource)/DL
# Some gymnastics to deal with units
R= (numpy.sqrt(theta[0]**2.+theta[1]**2.)*DL)\
.to(u.kpc,equivalencies=u.dimensionless_angles())
R= (R/(Pot[0]._ro if isinstance(Pot,list) else Pot._ro)/u.kpc)\
.to_value(u.dimensionless_unscaled)
phi= numpy.arctan2(theta[1],theta[0])
# Result includes more gymnastics to deal with units
return (cosmo_fac*2./c**2.
*integrate.quad(lambda z: evaluatePotentials(Pot,R,z,phi=phi,
use_physical=False),
-zmax,zmax)[0]
*(Pot[0]._ro*Pot[0]._vo**2. if isinstance(Pot,list)
else Pot._ro*Pot._vo**2.)*u.kpc*u.km**2/u.s**2)\
.to(u.arcsec**2,equivalencies=u.dimensionless_angles())
def fermat_potential(Pot,beta,theta,zlens,zsource,zmax=numpy.inf):
return 0.5*numpy.sum((beta-theta)**2.)\
-lensing_potential(Pot,theta,zlens,zsource,zmax=zmax)
from galpy.potential import PerfectEllipsoidPotential
# Twice as massive as axially-symmetric example, could
# also reduce zlens to 0.3-ish
pep= PerfectEllipsoidPotential(amp=10.**12*u.Msun,
a=4*u.kpc,b=0.5,c=1.)
zlens, zsource= 0.5, 1.
betas= [[2,1],[0.4,0.2],[0.2,0.1]]
levelss= [[21.8,22.,22.2,22.6,23.,23.4,23.8,24.2,24.6,25.,
25.5,26.5,27.5,28.5,30.,32.,34.,36.,38.,40.,43.,46.,49.],
[23.4,23.6,23.8,23.9,24.,24.2,24.4,24.6,24.669,
24.8,25.,25.4,25.8,26.2,26.6,27.,27.5,28.,28.5,29.,
29.5,30.,30.5,31.,31.5],
[23.6,23.65,23.7,23.75,23.8,23.8225,23.83,23.85,
23.9,24.,24.1,24.2,24.3,24.354,24.4,24.5,24.75,25.,
25.25,25.5,26.,26.5,27.,27.5,28.,28.5,29.,29.5,30.,30.5]]
saddle_contourss= [[],[24.669],[23.83,24.354]]
figure(figsize=(12,4))
for ii,(beta,levels,saddle_contours) \
in enumerate(zip(betas,levelss,saddle_contourss)):
subplot(1,3,ii+1)
if numpy.linalg.norm(beta) < 0.4:
thetas= numpy.linspace(-2.51,2.51,50)*u.arcsec
elif numpy.linalg.norm(beta) < 0.5:
thetas= numpy.linspace(-3.01,3.01,50)*u.arcsec
else:
thetas= numpy.linspace(-4.01,4.01,50)*u.arcsec
thetaxs= thetas
thetays= thetas
plotz= u.Quantity([u.Quantity(\
[fermat_potential(pep,u.Quantity(beta)*u.arcsec,
u.Quantity([thetax,thetay]),
zlens,zsource).to(u.arcsec**2)
for thetay in thetays]) for thetax in thetaxs])
levels= numpy.array(levels)
cntrlw= numpy.array([1. for l in levels])
for saddle in saddle_contours:
cntrlw[levels==saddle]= 2.
galpy_plot.dens2d(plotz.to_value(u.arcsec**2).T,
origin='lower',
xrange=[thetaxs[0].to_value(u.arcsec),
thetaxs[-1].to_value(u.arcsec)],
yrange=[thetays[0].to_value(u.arcsec),
thetays[-1].to_value(u.arcsec)],
xlabel=r'$\theta_1\,(\mathrm{arcsec})$',
ylabel=r'$\theta_2\,(\mathrm{arcsec})$' if ii == 0 else None,
justcontours=True,
levels=levels,
cntrlw=cntrlw,
gcf=True)
plot([beta[0]],[beta[1]],'+',color='#1f77b4',ms=15,mew=3)
tight_layout();
Figure 15.9: The Fermat potential for a massive elliptical galaxy without axial symmetry for three different source positions (the cross).
Again, when the source is far from the center of the lens in the left panel, there is a single image with \(\boldsymbol{\theta} \approx \boldsymbol{\beta}\). For a large-enough lens surface density, as the source gets closer to the lens center in the middle panel, an additional extremum (a maximum) appears near the center of the lens. Additionally, there is again a third extremum introduced by this, a saddle point that occurs where the thick contour intersects itself. In this case, there are therefore three images. As the source gets even closer to the lens center in the right panel, a further two extrema are introduced when the maximum is distorted into two maxima separated by another saddle point (at the intersection of the second thick contour). Thus, five images are produced.
While the exact locations and properties of the lensed images of a source depend on the details of the lensing potential, the general properties of images can be studied using methods from the mathematical field of topology. This field is concerned with properties of spaces that are invariant under, e.g., continuous deformations of the space. For example, the genus of a surface is the number of holes it has and is a topological quantity. A sphere has genus zero, as has any continuous deformation of the sphere, while a donut has genus one. The number of images that form in gravitational lensing is similarly a topological quantity. Using topology, one can prove a theorem (a theorem!) about the number of images produced by a gravitational lens:
The odd-number theorem (Burke 1981; McKenzie 1985): The number of images produced by a transparent, non-singular lens is odd.
This is a very general theorem that does not, for example, require that the lensing potential increases monotonically as the distance to the lens center increases (as we assumed above for realistic galactic mass distributions). All that the theorem requires is that the lens is transparent, that is, that it does not obscure any of the images through, e.g., dust attenuation. It further requires that the lens is non-singular, meaning that the deflection angle is finite everywhere.
The odd-number theorem is the best kind of theorem in astrophysics: a mathematically-proven statement that does not seem to hold in nature! Recall the examples given in Figure 15.1, the “twin quasar” 0957+561 and the RXJ1131-1231 system, or other famous gravitational lenses like the Einstein cross. These all have an even number of images, two or four typically. While some lenses with three images are known (e.g., Lawrence et al. 1984; Hewitt et al. 1987), they are rare and some lenses that are originally thought to have three images turn out to have four, with one image pair hard to resolve (e.g., the second lens system to be discovered: PG1115+08; Weymann et al. 1980; Hege et al. 1981). Thus, the assumptions behind the odd-number theorem have to be broken by most real gravitational lenses. We will come back to this later.
Other general statements about the lensed images that we can make relate to their parity—whether or not the image is flipped along two perpendicular directions—and the image magnification. These are statements about extended images and therefore require the magnification matrix from Section 15.3 above. From Equation (15.59), we can also write the inverse magnification matrix \(\mathcal{A}\) as the Hessian of the Fermat potential \begin{equation} \mathcal{A}_{ij} = {\partial^2 \tau(\boldsymbol{\theta};\boldsymbol{\beta})\over \partial \theta_i\partial \theta_j}\,, \end{equation} and the parity and magnification of images is therefore directly determined by the curvature of the Fermat potential or, equivalently, the time-delay surface \(\Delta t_\mathrm{lensing} = D_{\Delta t}\tau(\boldsymbol{\theta};\boldsymbol{\beta})/c\) (Equation 15.36). The parity and magnification are both determined by the eigenvalues \(\lambda_r\) and \(\lambda_t\) of the magnification matrix. We say that an image has positive partial parity in one of the eigenvectors if the corresponding eigenvalue is positive and negative partial parity otherwise. If the product of the eigenvalues is positive, we say that the images has positive total parity and similar for negative total parity. Minima and maxima have positive total parity, while saddle points have negative total parity, but only a minimum has two positive partial parities. These different parity combinations are also known as Type I (for two positive partial parities, a minimum), Type II (for opposite sign partial parities, a saddle point), and Type III (for two negative partial parities, a maximum). The following then hold (Blandford & Narayan 1986):
For a non-singular, transparent lens, \(n+1\) of the \(2n+1\) images have positive total parity and \(n\) have negative total parity. For example, in the example in Figure 15.9, we have either (i) a single minimum, (ii) a minimum, a maximum, and a saddle point, or (iii) two minima, a maximum, and two saddle points. The addition of an image pair when the surface density increases always leads to the introduction of a saddle point and either a minimum or a maximum. The sum of the total parities of all of the images is always one.
The global minimum has positive total and partial parities. For a variable source, the global minimum can be determined as that which is observed first (that is, whose variable time series is not repeated in the other images), because it is the minimum of the time delay surface.
Images that are minima have magnifications \(\mu > 1\) (Schneider 1984). This follows from the fact that they have two positive partial parities and thus \(\lambda_{_{t}^{r}}= 1-\kappa\pm\gamma > 1\) and \(|\gamma| < 1-\kappa < 1\) (because we assume that \(\Sigma > 0\)). The magnification from Equation (15.66) is then \(\mu = 1/|(1-\kappa)^2-\gamma^2| > 1\), because \((1-\kappa)^2-\gamma^2 < 1-\gamma^2 < 1\). Because the global minimum has to exist, there is always an image with more flux than the unlensed source. This again demonstrates that gravitational lensing always increases the total flux received from background sources.
Images that are closer together tend to be more highly magnified. Qualitatively, this is because images that are closer together cannot have too much curvature in the Fermat potential to go between two nearby extrema. As images get so close together that they merge, the Fermat potential has to become flat between them and the magnification goes to infinity. We will come back to this when we discuss the importance of the critical curve below.
Because we cannot see the unlensed source for lensing by galaxies, we can of course not directly observe the parity of each image, but we can determine the relative parities of all observed images. With time-delay observations, we can determine the global minimum and because we know that it has two positive partial parities, we can use it to determine the absolute parity of all images.
From these considerations, we can make the following general statements about possible image configurations:
If there is a single image, then it has to be a minimum, with positive total and partial parities and increased size/flux. This is the case for lenses that have too low surface density to allow multiple images to form, or far away from the center of a massive lens that can produce multiple images. Thus, lensing keeps the orientation of the single image intact, but it generically magnifies it.
If there are three images, one has to be a minimum and one has to be a saddle point. The remaining image can either be a minimum or a maximum. The contour through the saddle point can then have different topologies, with the first (two minima, one saddle) having the topology of the lemniscate and the second one that of the limaçon. The limaçon is the topology of the three-image saddle-contour in the example above. The lemniscate is a horizontal figure eight shape. From the relative parities, we can determine which topology a given observed lens has.
If there are five images, there are six combinations possible (Blandford & Narayan 1986). These are built from the lemniscate/limaçon building blocks, with the additional images forming from the splitting of one of the minima or maxima into three images with lemniscate or limaçon topology. There are six rather than eight topological possibilities, because the splitting of the minima in the two-minima (lemniscate) three-image configuration are equivalent to each other. In the example above, the topology is that of the minimum in the three-image limaçon splitting into a lemniscate topology with two minima and one saddle point.
More than five images is possible as well, but requires increasingly non-typical mass distributions.
15.4.2. Critical curves and caustics¶
For a massive galaxy that is capable of producing multiple images, an interesting question is in which regions of the source and image planes multiple images are possible. For sources at large distances from the lens center, a single image is formed at the minimum of the Fermat potential. For a non-singular lens, the lensing potential smoothly distorts the \(\left|\boldsymbol{\theta}-\boldsymbol{\beta}\right|^2/2\) contribution to the Fermat potential and as a source gets closer to the lens center, the Fermat potential (or the time-delay surface) distorts smoothly. At some point, the source will enter the region where three images form, with either lemniscate or limaçon topology. Typically, this involves a non-extremal part of the Fermat potential surface transitioning into a pair of extrema; looked at in the other direction of going from three to one image, this process is known as the merging of two images, after which they disappear (as we discuss below, it is also possible for all three images to merge simultaneously, after which a single image remains). Because the transformation of the surface into two additional extrema is smooth, it is clear that at the transition, the Fermat potential surface has to become flat in at least one direction, the direction along which the split between the two new images happens. Because the curvature matrix (the Hessian) of the Fermat potential surface is the inverse magnification matrix, these loci of zero curvature correspond to infinite magnification. The same happens when going from three to five images and so on: the part of the Fermat potential surface where two additional images form, always has to go through a configuration where the Fermat potential surface is flat along the splitting direction and the magnification diverges there (note that infinite magnification does not actually occur in reality, because of the finite size of the source and the breakdown of the geometrical-optics approximation near caustics; Ohanian 1983). Thus, curves of (formally) infinite magnification separate regions with different numbers of images. In the image plane, these curves are called critical curves, while their locations in the source plane are known as caustics (caustics is also more generally used for both critical curves and caustics). Image pairs that appear when a source crosses a caustic are always on opposite sides of the corresponding critical curve, because they always have to have opposite total parity from each other.
As we saw in Section 15.3.3 above, the magnification is \(\mu = 1/|\lambda_t\lambda_r|\), where \(\lambda_{_{t}^{r}}= 1-\kappa\pm\gamma\) are the radial and tangential eigenvalues. The critical curves are therefore the loci in the image plane where either the radial or tangential eigenvalue of the inverse magnification matrix is zero. When \(\lambda_r = 0\), we speak of the radial critical curve, while when \(\lambda_t = 0\), we call this the tangential critical curve. When images form near the radial critical curve, they are strongly deformed in the radial direction, while images that form near the tangential critical curve are strongly deformed in the tangential direction (see Equation 15.65). Because gravitational-lens systems are easier to spot when the magnification is high, observed lens systems typically have images near the critical curve and if the lensed source is an extended object, they are therefore highly distorted (see, for example, the example of RXJ1131-1231 at the beginning of this chapter).
A general method for computing the critical curves and caustics of a given lens is to compute the inverse of the magnification in the source plane and solve for the locus where it is zero. Separately determining the radial and tangential critical curves is possible by solving for where the radial or tangential eigenvalue is zero, respectively. The corresponding caustic can then be obtained using the lensing equation (15.7), which directly allows the position \(\boldsymbol{\beta}\) of the caustic to be computed from the \(\boldsymbol{\theta}\) location of the critical curve. For some lens mass distributions, we can compute the critical curve and the caustic analytically and this is useful to gain intuition about what to expect for the critical curves and caustics of realistic galaxy lenses. For example, as we demonstrated in Section 15.3.3, we can compute the eigenvalues of the inverse magnification matrix analytically for the singular-isothermal sphere and we have that \(\lambda_r = 1\) and \(\lambda_t = 1-2\kappa\). Only the tangential critical curve therefore exists, because \(\lambda_r\) never equals zero. In general, the radial critical curve is inside of the tangential critical curve. This is easy to see from the eigenvalues: \(\lambda_t = 0\) requires \(1-\kappa = \gamma > 0\) and thus, \(\kappa < 1\), while \(\lambda_r = 0\) needs \(1-\kappa = -\gamma < 0\), and therefore \(\kappa > 1\). Assuming that the surface density, and thus the convergence, increases towards the center of the lens, \(\lambda_r = 0\) therefore occurs closer to the lens center than \(\lambda_t=0\). However, for caustics the order is generally opposite, with the radial caustic occurring outside of the tangential caustic.
To further illustrate the critical curves and caustics, we can compute them for the lens mass distribution for which we displayed the Fermat potential in Figure 15.9. For this mass distribution, we need to calculate the critical curves and caustics numerically, and to do this we need to compute \(\kappa\) and \(\gamma\) from the second derivatives of the lensing potential. The second derivatives can be calculated by numerically integrating derivatives of the gravitational potential itself over \(z\) similar to how we computed the deflection angle in Section 15.2.1 above. With the numerically calculated values of \(\kappa(\boldsymbol{\theta})\) and \(\gamma(\boldsymbol{\theta})\), we can plot contours of where \(\lambda_{_{t}^{r}} = 0\) to show the critical curves and transform these curves to the source plane using the lensing equation to display the caustics. The result is shown in Figure 15.10.
[18]:
from scipy import integrate
from astropy.cosmology import Planck18 as cosmo
from astropy.constants import c
from galpy.potential import evaluateRforces, evaluatephitorques, \
evaluateR2derivs, evaluatephi2derivs, evaluateRphiderivs
def convergence_and_shear(Pot,theta,zlens,zsource,zmax=numpy.inf):
DL= cosmo.angular_diameter_distance(zlens)
cosmo_fac= DL*cosmo.angular_diameter_distance_z1z2(zlens,zsource)\
/cosmo.angular_diameter_distance(zsource)
R= (numpy.sqrt(theta[0]**2.+theta[1]**2.)*DL)\
.to(u.kpc,equivalencies=u.dimensionless_angles())
R= (R/(Pot[0]._ro if isinstance(Pot,list) else Pot._ro)/u.kpc)\
.to_value(u.dimensionless_unscaled)
phi= numpy.arctan2(theta[1],theta[0])
cp, sp= numpy.cos(phi), numpy.sin(phi)
cp2, sp2, cpsp= cp**2., sp**2., cp*sp
# Need to integrate Rforce, phiforce, R2deriv, phi2deriv, and Rphideriv
int_Rforce= integrate.quad(lambda z: evaluateRforces(Pot,R,z,phi=phi,
use_physical=False),
-zmax,zmax)[0]
int_phiforce= integrate.quad(lambda z: evaluatephitorques(Pot,R,z,phi=phi,
use_physical=False),
-zmax,zmax)[0]
int_R2deriv= integrate.quad(lambda z: evaluateR2derivs(Pot,R,z,phi=phi,
use_physical=False),
-zmax,zmax)[0]
int_phi2deriv= integrate.quad(lambda z: evaluatephi2derivs(Pot,R,z,phi=phi,
use_physical=False),
-zmax,zmax)[0]
int_Rphideriv= integrate.quad(lambda z: evaluateRphiderivs(Pot,R,z,phi=phi,
use_physical=False),
-zmax,zmax)[0]
# Now compute x2deriv, y2deriv, and xyderiv
int_x2deriv= int_R2deriv*cp2-2.*int_Rphideriv*cpsp/R+int_phi2deriv*sp2/R**2.\
-int_Rforce*sp2/R-2.*int_phiforce*cpsp/R**2.
int_y2deriv= int_R2deriv*sp2+2.*int_Rphideriv*cpsp/R+int_phi2deriv*cp2/R**2.\
-int_Rforce*cp2/R+2.*int_phiforce*cpsp/R**2.
int_xyderiv= int_R2deriv*cpsp+int_Rphideriv*(cp2-sp2)/R-int_phi2deriv*cpsp/R**2.\
+int_Rforce*cpsp/R+int_phiforce*(cp2-sp2)/R**2.
# and finally the convergence and shear
kappa= 0.5*(int_x2deriv+int_y2deriv)
gamma1= 0.5*(int_x2deriv-int_y2deriv)
gamma2= int_xyderiv
gamma= numpy.sqrt(gamma1**2.+gamma2**2.)
return ((cosmo_fac*2./c**2.*kappa
*(Pot[0]._vo**2./Pot[0]._ro if isinstance(Pot,list)
else Pot._vo**2./Pot._ro)/u.kpc*u.km**2/u.s**2)\
.to_value(u.dimensionless_unscaled),
(cosmo_fac*2./c**2.*gamma
*(Pot[0]._vo**2./Pot[0]._ro if isinstance(Pot,list)
else Pot._vo**2./Pot._ro)/u.kpc*u.km**2/u.s**2)\
.to_value(u.dimensionless_unscaled))
# Deflection angle code from Section 15.1.1
from scipy import integrate
from astropy.constants import c
from galpy.potential import evaluateRforces, evaluatephitorques
def deflection(Pot,R,phi=0.):
# Some gymnastics to deal with units
if isinstance(R,u.Quantity):
R= (R/(Pot[0]._ro if isinstance(Pot,list) else Pot._ro)/u.kpc)\
.to_value(u.dimensionless_unscaled)
cp,sp= numpy.cos(phi), numpy.sin(phi)
int_Rforce= integrate.quad(lambda z: evaluateRforces(Pot,R,z,phi=phi,
use_physical=False),
-numpy.inf,numpy.inf)[0]
int_phiforce= integrate.quad(lambda z: evaluatephitorques(Pot,R,z,phi=phi,
use_physical=False),
-numpy.inf,numpy.inf)[0]
# Result includes some gymnastics to deal with units
return (-2./c**2.*numpy.array([int_Rforce*cp-int_phiforce*sp/R,
int_Rforce*sp+int_phiforce*cp/R])\
*(Pot[0]._vo**2. if isinstance(Pot,list) else Pot._vo**2.)*u.km**2/u.s**2)\
.to(u.arcsec,equivalencies=u.dimensionless_angles())
def source_from_image(Pot,theta,zlens,zsource):
"""Compute source position from image position
using the lensing equation"""
DL= cosmo.angular_diameter_distance(zlens)
cosmo_fac= cosmo.angular_diameter_distance_z1z2(zlens,zsource)\
/cosmo.angular_diameter_distance(zsource)
R= (numpy.sqrt(theta[0]**2.+theta[1]**2.)*DL)\
.to(u.kpc,equivalencies=u.dimensionless_angles())
R= (R/(Pot[0]._ro if isinstance(Pot,list) else Pot._ro)/u.kpc)\
.to_value(u.dimensionless_unscaled)
phi= numpy.arctan2(theta[1],theta[0])
alpha= cosmo_fac*deflection(Pot,R,phi=phi)
return theta-alpha
from galpy.potential import PerfectEllipsoidPotential
# Twice as massive as axially-symmetric example, could
# also reduce zlens to 0.3-ish
pep= PerfectEllipsoidPotential(amp=10.**12*u.Msun,
a=4*u.kpc,b=0.5,c=1.)
zlens, zsource= 0.5, 1.
# Compute convergence and shear on a grid in theta1 and theta2
thetas= numpy.linspace(-2.21,2.21,44)*u.arcsec
thetaxs= thetas
thetays= thetas
kappagamma= numpy.array([[convergence_and_shear(pep,
u.Quantity([thetax,thetay]),
zlens,zsource)
for thetay in thetays] for thetax in thetaxs])
kappa= kappagamma[:,:,0]
gamma= kappagamma[:,:,1]
lambdar= 1.-kappa+gamma
lambdat= 1.-kappa-gamma
figure(figsize=(7,4))
# Critical curves
subplot(1,2,1)
ltcont= galpy_plot.dens2d(lambdat.T,origin='lower',
xrange=[thetaxs[0].to_value(u.arcsec),
thetaxs[-1].to_value(u.arcsec)],
yrange=[thetays[0].to_value(u.arcsec),
thetays[-1].to_value(u.arcsec)],
xlabel=r'$\theta_1\,(\mathrm{arcsec})$',
ylabel=r'$\theta_2\,(\mathrm{arcsec})$',
justcontours=True,cntrlw=2.,
levels=[0.],gcf=True,
retCont=True)
lrcont= galpy_plot.dens2d(lambdar.T,origin='lower',
xrange=[thetaxs[0].to_value(u.arcsec),
thetaxs[-1].to_value(u.arcsec)],
yrange=[thetays[0].to_value(u.arcsec),
thetays[-1].to_value(u.arcsec)],
xlabel=r'$\theta_1\,(\mathrm{arcsec})$',
ylabel=r'$\theta_2\,(\mathrm{arcsec})$',
justcontours=True,cntrlw=2.,cntrls=':',
levels=[0.],overplot=True,
retCont=True)
galpy_plot.text(-2.,1.6,r'$\mathrm{Tangential\ critical\ curve}$',
fontsize=15.)
plot([-1.6,-1.],[1.4,0.85],'k-',lw=0.8)
galpy_plot.text(-2.,-1.9,r'$\mathrm{Radial\ critical\ curve}$',
fontsize=15.)
plot([-1.2,-.35],[-1.5,-.35],'k-',lw=0.8)
annotate(text=r'$\mathrm{Image\ plane}$',xy=(0.5,1.01),xycoords='axes fraction',
ha='center',va='bottom',fontsize=18.)
# Corresponding caustics
subplot(1,2,2)
from packaging.version import parse as parse_version
import matplotlib
for ii,cont in enumerate([ltcont,lrcont]):
# Extract critical curve contour
if parse_version(matplotlib.__version__) < parse_version('3.8'):
critx= cont.collections[0].get_paths()[0].vertices[:,0]
crity= cont.collections[0].get_paths()[0].vertices[:,1]
# Contour may consist of multiple paths
for further in cont.collections[0].get_paths()[1:]:
critx= numpy.concatenate((critx,further.vertices[:,0]))
crity= numpy.concatenate((crity,further.vertices[:,1]))
else:
critx= cont.get_paths()[0].vertices[:,0]
crity= cont.get_paths()[0].vertices[:,1]
# Contour may consist of multiple paths
for further in cont.get_paths()[1:]:
critx= numpy.concatenate((critx,further.vertices[:,0]))
crity= numpy.concatenate((crity,further.vertices[:,1]))
# Map to source plane using the lensing equation
beta= numpy.array([source_from_image(pep,u.Quantity([tcritx,tcrity])*u.arcsec,zlens,zsource)
for tcritx,tcrity in zip(critx,crity)])
galpy_plot.plot(beta[:,0],beta[:,1],'k-' if ii == 0 else 'k:',lw=2.,
xrange=[thetaxs[0].to_value(u.arcsec),thetaxs[-1].to_value(u.arcsec)],
yrange=[thetays[0].to_value(u.arcsec),thetays[-1].to_value(u.arcsec)],
xlabel=r'$\theta_1\,(\mathrm{arcsec})$',
ylabel=r'$\theta_2\,(\mathrm{arcsec})$',
gcf=ii==0,overplot=ii>0)
galpy_plot.text(-1.8,1.8,r'$\mathrm{Tangential\ caustic}$',
fontsize=15.)
plot([-1.2,-.25],[1.65,0.25],'k-',lw=0.8)
galpy_plot.text(-1.8,-1.9,r'$\mathrm{Radial\ caustic}$',
fontsize=15.)
plot([-1.2,-0.8],[-1.6,-1.2],'k-',lw=0.8)
annotate(text=r'$\mathrm{Source\ plane}$',xy=(0.5,1.01),xycoords='axes fraction',
ha='center',va='bottom',fontsize=18.)
tight_layout();
Figure 15.10: Critical curves and caustics.
On the left we see the expected behavior in the image plane: the tangential critical curve occurs outside of the radial critical curve. However, the caustics in the source plane on the right demonstrate that the curves are reversed in the source plane: the tangential caustic is inside the radial caustic. Thus, generically what happens when we move a source closer to the lens center is the following: a source outside the radial caustic produces a single image; as the source crosses the radial caustic, two additional images are produced and these additional images appear close to the center of the lens near the radial critical curve; as the source moves closer and crosses the tangential caustic, two more images are created and these images appear further from the lens center, near the tangential critical curve. This is indeed the behavior seen in the Fermat potential surfaces in Figure 15.9 when we move the source closer to the center. All critical curves and caustics are roughly elliptically shaped, except for the tangential caustic, which instead looks like an astroid; this is typically the case and the tangential caustic is therefore also known as the astroid caustic.
The critical curves and caustics depend strongly on the shape of the lens mass distribution. Taking the same mass distribution as above, but compressing it to an axis ratio of \(0.2\) or rounding it so it is almost axially-symmetric with an axis ratio of \(0.99\), we get the critical curves and caustics displayed in Figure 15.11.
[19]:
pep2= PerfectEllipsoidPotential(amp=10.**12*u.Msun,
a=4*u.kpc,b=0.2,c=1.)
pep99= PerfectEllipsoidPotential(amp=10.**12*u.Msun,
a=4*u.kpc,b=0.99,c=1.)
zlens, zsource= 0.5, 1.
# Compute convergence and shear on a grid in theta1 and theta2
thetas= numpy.linspace(-3.01,3.01,44)*u.arcsec
thetaxs= thetas
thetays= thetas
gcf= True
figure(figsize=(7,4))
for pep,color in zip([pep2,pep99],['#1f77b4','#ff7f0e']):
kappagamma= numpy.array([[convergence_and_shear(pep,
u.Quantity([thetax,thetay]),
zlens,zsource)
for thetay in thetays] for thetax in thetaxs])
kappa= kappagamma[:,:,0]
gamma= kappagamma[:,:,1]
lambdar= 1.-kappa+gamma
lambdat= 1.-kappa-gamma
# Critical curves
subplot(1,2,1)
ltcont= galpy_plot.dens2d(lambdat.T,origin='lower',
xrange=[thetaxs[0].to_value(u.arcsec),
thetaxs[-1].to_value(u.arcsec)],
yrange=[thetays[0].to_value(u.arcsec),
thetays[-1].to_value(u.arcsec)],
xlabel=r'$\theta_1\,(\mathrm{arcsec})$',
ylabel=r'$\theta_2\,(\mathrm{arcsec})$',
justcontours=True,cntrlw=2.,cntrcolors=color,
levels=[0.],gcf=gcf,overplot=not gcf,
retCont=True)
lrcont= galpy_plot.dens2d(lambdar.T,origin='lower',
xrange=[thetaxs[0].to_value(u.arcsec),
thetaxs[-1].to_value(u.arcsec)],
yrange=[thetays[0].to_value(u.arcsec),
thetays[-1].to_value(u.arcsec)],
xlabel=r'$\theta_1\,(\mathrm{arcsec})$',
ylabel=r'$\theta_2\,(\mathrm{arcsec})$',
justcontours=True,cntrlw=2.,cntrls=':',
levels=[0.],overplot=True,cntrcolors=color,
retCont=True)
annotate(text=r'$\mathrm{Image\ plane}$',xy=(0.5,1.01),xycoords='axes fraction',
ha='center',va='bottom',fontsize=18.)
# Corresponding caustics
subplot(1,2,2)
for ii,cont in enumerate([ltcont,lrcont]):
# Extract critical curve contour
if parse_version(matplotlib.__version__) < parse_version('3.8'):
critx= cont.collections[0].get_paths()[0].vertices[:,0]
crity= cont.collections[0].get_paths()[0].vertices[:,1]
# Contour may consist of multiple paths
for further in cont.collections[0].get_paths()[1:]:
critx= numpy.concatenate((critx,further.vertices[:,0]))
crity= numpy.concatenate((crity,further.vertices[:,1]))
else:
critx= cont.get_paths()[0].vertices[:,0]
crity= cont.get_paths()[0].vertices[:,1]
# Contour may consist of multiple paths
for further in cont.get_paths()[1:]:
critx= numpy.concatenate((critx,further.vertices[:,0]))
crity= numpy.concatenate((crity,further.vertices[:,1]))
# Map to source plane using the lensing equation
beta= numpy.array([source_from_image(pep,u.Quantity([tcritx,tcrity])*u.arcsec,zlens,zsource)
for tcritx,tcrity in zip(critx,crity)])
galpy_plot.plot(beta[:,0],beta[:,1],ls='-' if ii == 0 else ':',lw=2.,color=color,
xrange=[thetaxs[0].to_value(u.arcsec),thetaxs[-1].to_value(u.arcsec)],
yrange=[thetays[0].to_value(u.arcsec),thetays[-1].to_value(u.arcsec)],
xlabel=r'$\theta_1\,(\mathrm{arcsec})$',
ylabel=r'$\theta_2\,(\mathrm{arcsec})$',
gcf=gcf,overplot=not gcf)
gcf= False
subplot(1,2,2)
annotate(text=r'$\mathrm{Source\ plane}$',xy=(0.5,1.01),xycoords='axes fraction',
ha='center',va='bottom',fontsize=18.)
tight_layout();
Figure 15.11: The effect of a lens’ shape on the critical curves and caustics. The blue curves correspond to a highly flattened lens.
For the near-axially-symmetric case, we see that all curves become near-circles; this makes intuitive sense, because there is no preferred direction in an axially-symmetric lens. It remains the case that the radial critical curve is inside of the tangential critical curve no matter the shape of the lens, but in the axially-symmetric case, the tangential caustic becomes a point. This is a general property of axially-symmetric lenses. Thus, sources never cross the tangential caustic for axially-symmetric lenses and only three images can be produced (as we already argued above). However, the relevance of the tangential critical curve is that it is the Einstein ring for a source that is right behind the lens. Such a source is lensed into an entire ring consisting of the tangential critical curve. For the strongly-flattened lens, the critical curves are also highly flattened, in the same direction as the mass distribution (we flattened it along \(\theta_2\)). The tangential caustic also stretches in the same direction as the mass distribution, but the radial caustic elongates in the perpendicular direction. The tangential caustic also starts to encroach on the radial one and they even meet at two points. This means that along directions that cross either of these points, there is a transition from a single image to five images! For highly-flattened lenses, the tangential caustic can even extend beyond the radial caustic, a configuration that is known as a naked cusp (Keeton & Kochanek 1998). The significance of naked-cusp lenses is that these can produce three images of a source on the same side of the lens; when a source crosses the radial caustic first, the two additional images generally lie on the opposite side of the lens as the original image (as in Figure 15.9).
From the discussion so far, it is clear that what happens to sources near the caustics is highly important in studies of gravitational lensing, because observational biases cause high-magnification lenses to be preferentially found in observational data. Again, it is possible to make very general statements about the behavior near caustics, because the universal aspects of this behavior are described by the mathematical theory of catastrophes. We will only briefly describe important aspects of this for the study of gravitational lensing. Catastrophe theory is the general theory of dynamical systems that display qualitative changes in their behavior under small changes to the input parameters, called “control parameters”. For a gravitational lens of a fixed mass distribution and with the source and lens at fixed redshifts, there are two control parameters, the two-dimensional location \(\boldsymbol{\beta}\) of the source in the source plane. Fundamentally, there are only two basic catastrophes for a system with two control parameters: the fold and the cusp. The fold is the form of both the radial and astroid caustic, except in the corners of the astroid, where the catastrophe has the form of a cusp. The difference between a fold and a cusp is that when we move across a fold caustic going from the region of \(n\) to \(n-2\) images (e.g., \(n=3\) or \(n=5\) in the examples here), two images approach the corresponding critical curve and merge and disappear as the source crosses the fold caustic. The remaining images are far removed from this location in the image plane. On the other hand, when a source crosses the cusp caustic going from the region of \(n\) to \(n-2\) images, three images approach each other along the critical curve. When they merge as the source crosses the cusp caustic, they combine into a single image (but note that there may be other pairs of images in the system, if \(n>3\)). Because this single image lies close to the critical curve, it is highly magnified. A source near a cusp caustic is therefore a way to create highly magnified isolated images (highly magnified images near a fold come in pairs).
The example in Figure 15.9 displayed a source crossing fold caustics. It is clear from that figure that pairs of images are created where there was no image before as the source moves across the radial and astroid caustics closer to the center. As an example of a source crossing a cusp, we can take the same massive elliptical galaxy centered at the origin, but now moving the source closer and closer to the lens center along the \(\theta_2 = 0\) direction, where the astroid caustic has a cusp. This gives the Fermat potentials shown in Figure 15.12.
[20]:
from galpy.potential import PerfectEllipsoidPotential
# Twice as massive as axially-symmetric example, could
# also reduce zlens to 0.3-ish
pep= PerfectEllipsoidPotential(amp=10.**12*u.Msun,
a=4*u.kpc,b=0.5,c=1.)
zlens, zsource= 0.5, 1.
betas= [[0.4,0.],[0.3,0.],[0.075,0.]]
levelss= [[23.58,23.7,23.9,24.05,24.2,24.4,24.63,24.75,25.1,25.45,25.8,
26.15,26.5,26.85,27.2,27.55,27.9,28.25,28.6,28.95,
29.3,29.65,30.,30.35,30.7],
[23.5,23.65,23.6817,23.8,23.95,24.1,24.25,24.4,24.486,24.55,24.7,
24.85,25.,25.15,25.3,25.45,25.6,25.75,25.9,26.05,
26.2,26.35,26.5,26.65,26.8,26.95,27.1,27.25,27.4,
27.55,27.7,27.85],
[23.775,23.875,23.96912,24.05,24.125,24.1705,24.25,24.375,24.5001,
24.75,25.,25.25,25.5,25.75,26.,26.25,26.5,27.,27.4]]
saddle_contourss= [[24.63],[23.6817,24.486],[23.96912,24.1705]]
figure(figsize=(12,4))
for ii,(beta,levels,saddle_contours) \
in enumerate(zip(betas,levelss,saddle_contourss)):
subplot(1,3,ii+1)
if numpy.linalg.norm(beta) < 0.4:
thetas= numpy.linspace(-2.251,2.251,50)*u.arcsec
elif numpy.linalg.norm(beta) < 0.5:
thetas= numpy.linspace(-2.751,2.751,50)*u.arcsec
else:
thetas= numpy.linspace(-3.51,3.51,50)*u.arcsec
thetaxs= thetas
thetays= thetas
plotz= u.Quantity([u.Quantity(\
[fermat_potential(pep,u.Quantity(beta)*u.arcsec,
u.Quantity([thetax,thetay]),
zlens,zsource).to(u.arcsec**2)
for thetay in thetays]) for thetax in thetaxs])
levels= numpy.array(levels)
cntrlw= numpy.array([1. for l in levels])
for saddle in saddle_contours:
cntrlw[levels==saddle]= 2.
galpy_plot.dens2d(plotz.to_value(u.arcsec**2).T,origin='lower',
xrange=[thetaxs[0].to_value(u.arcsec),
thetaxs[-1].to_value(u.arcsec)],
yrange=[thetays[0].to_value(u.arcsec),
thetays[-1].to_value(u.arcsec)],
xlabel=r'$\theta_1\,(\mathrm{arcsec})$',
ylabel=r'$\theta_2\,(\mathrm{arcsec})$' if ii == 0 else None,
justcontours=True,
levels=levels,
cntrlw=cntrlw,
gcf=True)
plot([beta[0]],[beta[1]],'+',color='#1f77b4',ms=15,mew=3)
tight_layout();
Figure 15.12: Fermat potentials like in Figure 15.9, but moving the source along \(\theta_2 = 0\) such that it crosses the caustic at a cusp.
Again, when the source is far from the lens center, a single image is produced (not shown) and as the source crosses the fold radial caustic in the left panel, a pair of images appears on the opposite side of the lens. However, when the source crosses the cusp astroid caustic in the middle panel, the single image in the left panel splits into three images (with limaçon topology, two minima and a saddle point; bottom left panel). Thus, we have an image configuration with three images on one side of the lens, an image on the other side of the lens, and a final image near the center of the lens. This is similar to the image morphology of the RXJ1131-1231 lens system shown in Figure 15.1 (at least if the central image is ignored; see below). As the source gets closer to the center, the two minima move towards larger \(|\theta_2|\) and smaller \(|\theta_1|\) and as the source gets very close to center in the right panel, the images appear on four sides of the lens, a configuration similar to that of the famous “Einstein cross” (again, ignoring the central image; Huchra et al. 1985). Lensing close to the center of the astroid caustic is sometimes referred to as the cross configuration, with the four basic strong lensing configurations then being fold, cusp, cross, and Einstein ring.
Because lensing near a fold or a cusp caustic produces the highly-magnified images that are preferentially detected in flux-limited surveys of galaxies (Turner et al. 1984), it is useful to understand lensing near these catastrophes. One general set of results that can be derived for lensing near folds or cusps has to do with the relative magnifications of the three images involved in the crossing of the caustic near a fold or cusp. For smooth lens mass distributions (that is, such that the lens does not significantly change around the caustic), the following relations should hold (Blandford & Narayan 1986; Schneider & Weiss 1992):
Lensing near a fold: the magnifications of the two images \(A\) and \(B\) that appear as a source crosses a fold caustic are equal: \(\mu_A = \mu_B\);
Lensing near a cusp: the magnification of the central image \(A\) is equal to the sum of the magnifications of the flanking images \(B\) and \(C\): \(\mu_A = \mu_B + \mu_C\).
These two relations can be succinctly expressed if we give the magnification the sign of the total parity (such that the magnification can be positive or negative). In that case, both relations become \(\sum_i \mu_i = 0\): the sum of the signed magnifications equals zero near a fold or a cusp. Even though it is difficult to obtain absolute magnifications, we can determine relative magnifications by comparing the fluxes of different images. As with the odd-number theorem, these relations are also not typically observed to hold for observed lenses. This could be due to intrinsic variability of the source, because then the relative time delays mean that we are not observing the same intrinsic flux for all images. Or the assumption that the lens is smooth could be broken, for example, because of microlensing by stars in the lens (if the source is so small that it is susceptible to microlensing, see Section 15.2.3) or if the mass distribution is not smooth due to substructure in the dark matter distribution. If all other possibilities can be excluded, so-called flux-ratio anomalies become a sensitive and robust method for determining the clustering of dark matter on small scales and constraining the microphysics of dark matter (Mao & Schneider 1998; Metcalf & Zhao 2002; Dalal & Kochanek 2002; Gilman et al. 2020).
15.4.3. Lensing by cuspy mass distributions¶
We now understand many of the general characteristics of where lensed images appear and what their basic properties are. But we haven’t yet addressed why the most basic general result in lensing does not seem to hold: the odd-number theorem appears to be in conflict with almost all observed lens systems! The basic reason for this (and what it helps demonstrate) is that the cores of galaxies are very dense and indeed, the density distribution in the cores of galaxies may be singular. For example, as we discussed in Chapter 12.1, the cores of elliptical galaxies range from near-isothermal, with stellar densities \(\rho(r) \propto r^{-\alpha}\) with \(0 \lesssim \alpha \lesssim 1\), to strong cusps with \(1 \lesssim \alpha \lesssim 2\). The dark-matter halo is similarly believed to have \(\rho(r) \propto 1/r\). As we discussed above, for singular lens mass distributions, the odd-number theorem does not hold, because for singular lenses, the image that would appear near the lens center is instead a defect in the Fermat potential. The relevant singularity is that of the deflection angle rather than the density: the odd-number theorem holds for singular densities that produce finite deflection angles everywhere, but not for singular deflection angles. Thus, for \(\rho(r) \propto r^{-\alpha}\), shallow cusps have \(0 \lesssim \alpha \lesssim 1\) and the odd-number theorem holds, while for strong cusps with \(1 \lesssim \alpha \lesssim 2\) the odd-number theorem fails and we expect an even number of images (this is the same distinction between shallow and non-shallow cusps that was relevant for the orbital structure and self-consistent shape of triaxial models in Chapter 13.4.2).
If lenses aren’t actually singular, but just very dense (as must be the case for real systems), the odd-number theorem should hold, but it will appear to be broken. This is because for a very dense lens, the image that forms near the center of the lens is highly demagnified. For a dense lens, the lensing potential is deep and highly curved and one of the images that forms beyond the first is a maximum near the center of the lens with high curvature of the Fermat potential surface. Because the magnification is the inverse of the curvature, this image is highly demagnified and therefore difficult to observe. Moreover, because it near-coincides with the lens, the lensing galaxy itself can obscure the image further (thus breaking the other assumption in the odd-number theorem: the transparency of the lens). Thus, whether because of a true singularity or simply a dense galaxy core, we expect to not see the central image and thus we expect to observe an even number of images.
Beyond the central image, lensing in singular or high-density lenses is similar to that in non-singular, lower-density lenses such as the examples that we considered above. For example, we show the Fermat potentials when moving a source towards the center of a lens modeled as a flattened, cuspy Hernquist density along a fold astroid caustic in Figure 15.13 and this produces a similar sequence of images as in the cored perfect-ellipsoid example in Figure 15.12. The Hernquist’s central behavior is such that the deflection angle is finite everywhere and there are therefore an odd number of images.
[21]:
from galpy.potential import TriaxialHernquistPotential
pep= TriaxialHernquistPotential(amp=3*10.**12*u.Msun,
a=4*u.kpc,b=0.5,c=1.)
zlens, zsource= 0.5, 1.
betas= [[0.6,0.],[0.2,0.],[0.05,0.]]
levelss= [sorted(numpy.concatenate(([18.6566,18.75],
numpy.arange(17.5,18.,.3),
numpy.arange(18.,33.,1.)))),
sorted(numpy.concatenate(([17.735,17.7657,18.18675],
numpy.arange(17.7,18.,.05),
numpy.arange(18.,22.4,.3)))),
sorted(numpy.concatenate(([17.9172,18.02335],
numpy.arange(17.8,17.9,.05),
numpy.arange(18.1,22.,.3))))]
saddle_contourss= [[18.6566],[17.7657,18.18675],[17.9172,18.02335]]
figure(figsize=(12,4))
for ii,(beta,levels,saddle_contours) \
in enumerate(zip(betas,levelss,saddle_contourss)):
subplot(1,3,ii+1)
if numpy.linalg.norm(beta) < 0.4:
thetas= numpy.linspace(-2.51,2.51,50)*u.arcsec
elif numpy.linalg.norm(beta) < 0.5:
thetas= numpy.linspace(-3.01,3.01,50)*u.arcsec
else:
thetas= numpy.linspace(-4.01,4.01,50)*u.arcsec
thetaxs= thetas
thetays= thetas
plotz= u.Quantity([u.Quantity(\
[fermat_potential(pep,u.Quantity(beta)*u.arcsec,
u.Quantity([thetax,thetay]),
zlens,zsource).to(u.arcsec**2)
for thetay in thetays]) for thetax in thetaxs])
levels= numpy.array(levels)
cntrlw= numpy.array([1. for l in levels])
for saddle in saddle_contours:
cntrlw[levels==saddle]= 2.
galpy_plot.dens2d(plotz.to_value(u.arcsec**2).T,origin='lower',
xrange=[thetaxs[0].to_value(u.arcsec),
thetaxs[-1].to_value(u.arcsec)],
yrange=[thetays[0].to_value(u.arcsec),
thetays[-1].to_value(u.arcsec)],
xlabel=r'$\theta_1\,(\mathrm{arcsec})$',
ylabel=r'$\theta_2\,(\mathrm{arcsec})$' if ii == 0 else None,
justcontours=True,
levels=levels,
cntrlw=cntrlw,
gcf=True)
plot([beta[0]],[beta[1]],'+',color='#1f77b4',ms=15,mew=3)
tight_layout();
Figure 15.13: Like Figure 15.12, but for the singular Hernquist profile.
Considering a flattened version of the singular isothermal sphere instead, which has \(\rho(r) \propto 1/r^2\), we still get a similar sequence, but the image that appears in the Hernquist and perfect-ellipsoid lenses near the center is absent. To illustrate this, the Fermat potentials in Figure 15.14 use a set of source positions that cross the fold caustic as in the example in Figure 15.9.
[22]:
from galpy.potential import PowerTriaxialPotential
pep= PowerTriaxialPotential(amp=5*10.**10*u.Msun,alpha=2.01,b=0.5,c=1.)
zlens, zsource= 0.5, 1.
betas= [[1,0.5],[0.3,0.15],[0.08,0.04]]
levelss= [numpy.arange(4982,5000,0.5),
sorted(numpy.concatenate(([4984.26,4984.265,4984.983],
numpy.arange(4983.9,4984.5,0.1),
numpy.arange(4984.6,4990.0,0.2)))),
sorted(numpy.concatenate(([4984.4908,4984.69],
numpy.arange(4984.1,4984.5,0.075),
numpy.arange(4984.6,4989.,0.2))))]
saddle_contourss= [[],[4984.265,4984.983],[4984.4908,4984.69]]
figure(figsize=(12,4))
for ii,(beta,levels,saddle_contours) \
in enumerate(zip(betas,levelss,saddle_contourss)):
subplot(1,3,ii+1)
if numpy.linalg.norm(beta) < 0.4:
thetas= numpy.linspace(-3.01,3.01,50)*u.arcsec
else:
thetas= numpy.linspace(-4.01,4.01,50)*u.arcsec
thetaxs= thetas
thetays= thetas
plotz= u.Quantity([u.Quantity(\
[fermat_potential(pep,u.Quantity(beta)*u.arcsec,
u.Quantity([thetax,thetay]),
zlens,zsource,zmax=50.).to(u.arcsec**2)
for thetay in thetays]) for thetax in thetaxs])
levels= numpy.array(levels)
cntrlw= numpy.array([1. for l in levels])
for saddle in saddle_contours:
cntrlw[levels==saddle]= 2.
galpy_plot.dens2d(plotz.to_value(u.arcsec**2).T,origin='lower',
xrange=[thetaxs[0].to_value(u.arcsec),
thetaxs[-1].to_value(u.arcsec)],
yrange=[thetays[0].to_value(u.arcsec),
thetays[-1].to_value(u.arcsec)],
xlabel=r'$\theta_1\,(\mathrm{arcsec})$',
ylabel=r'$\theta_2\,(\mathrm{arcsec})$' if ii == 0 else None,
justcontours=True,
levels=levels,
cntrlw=cntrlw,
gcf=True)
plot([beta[0]],[beta[1]],'+',color='#1f77b4',ms=15,mew=3)
tight_layout();
Figure 15.14: Like Figure 15.9, but for the singular isothermal sphere.
Note that what appears to be an extremum near the center in the middle and left panels is in fact a defect without an image. The kink in the contour near the center in the left panel shows the effect of the defect.
While for the flattened Hernquist lens, the critical curves and caustics are qualitatively similar to that of the perfect ellipsoid, for the flattened singular-isothermal sphere the situation is more complex. As we saw in Equation (15.71), for the axially-symmetric singular isothermal sphere, \(\lambda_r = 1\) and it therefore does not have a radial critical curve or caustic. This continues to be the case for the flattened singular-isothermal sphere (Kormann et al. 1994). The tangential critical curve and caustic, however, are similar as before, as can be seen in Figure 15.15.
[23]:
from galpy.potential import PowerTriaxialPotential
# Similar surface density as previous models
pep= PowerTriaxialPotential(amp=5*10.**10*u.Msun,alpha=2.01,b=0.5,c=1.)
zlens, zsource= 0.5, 1.
# Compute convergence and shear on a grid in theta1 and theta2
thetas= numpy.linspace(-3.01,3.01,44)*u.arcsec
thetaxs= thetas
thetays= thetas
kappagamma= numpy.array([[convergence_and_shear(pep,
u.Quantity([thetax,thetay]),
zlens,zsource)
for thetay in thetays] for thetax in thetaxs])
kappa= kappagamma[:,:,0]
gamma= kappagamma[:,:,1]
lambdar= 1.-kappa+gamma
lambdat= 1.-kappa-gamma
figure(figsize=(7,4))
# Critical curves
subplot(1,2,1)
ltcont= galpy_plot.dens2d(lambdat.T,origin='lower',
xrange=[thetaxs[0].to_value(u.arcsec),
thetaxs[-1].to_value(u.arcsec)],
yrange=[thetays[0].to_value(u.arcsec),
thetays[-1].to_value(u.arcsec)],
xlabel=r'$\theta_1\,(\mathrm{arcsec})$',
ylabel=r'$\theta_2\,(\mathrm{arcsec})$',
justcontours=True,cntrlw=2.,
levels=[0.],gcf=True,
retCont=True)
galpy_plot.text(-2.3,2.,r'$\mathrm{Tangential\ critical\ curve}$',
fontsize=15.)
plot([-1.8,-1.05],[1.8,1.05],'k-',lw=0.8)
annotate(text=r'$\mathrm{Image\ plane}$',xy=(0.5,1.01),xycoords='axes fraction',
ha='center',va='bottom',fontsize=18.)
# Corresponding caustics
from packaging.version import parse as parse_version
import matplotlib
subplot(1,2,2)
for ii,cont in enumerate([ltcont]):
# Extract critical curve contour
if parse_version(matplotlib.__version__) < parse_version('3.8'):
critx= cont.collections[0].get_paths()[0].vertices[:,0]
crity= cont.collections[0].get_paths()[0].vertices[:,1]
# Contour may consist of multiple paths
for further in cont.collections[0].get_paths()[1:]:
critx= numpy.concatenate((critx,further.vertices[:,0]))
crity= numpy.concatenate((crity,further.vertices[:,1]))
else:
critx= cont.get_paths()[0].vertices[:,0]
crity= cont.get_paths()[0].vertices[:,1]
# Contour may consist of multiple paths
for further in cont.get_paths()[1:]:
critx= numpy.concatenate((critx,further.vertices[:,0]))
crity= numpy.concatenate((crity,further.vertices[:,1]))
# Map to source plane using the lensing equation
beta= numpy.array([source_from_image(pep,u.Quantity([tcritx,tcrity])*u.arcsec,zlens,zsource)
for tcritx,tcrity in zip(critx,crity)])
galpy_plot.plot(beta[:,0],beta[:,1],'k-' if ii == 0 else 'k:',lw=2.,
xrange=[thetaxs[0].to_value(u.arcsec),thetaxs[-1].to_value(u.arcsec)],
yrange=[thetays[0].to_value(u.arcsec),thetays[-1].to_value(u.arcsec)],
xlabel=r'$\theta_1\,(\mathrm{arcsec})$',
ylabel=r'$\theta_2\,(\mathrm{arcsec})$',
gcf=ii==0,overplot=ii>0)
galpy_plot.text(-2.3,2.,r'$\mathrm{Tangential\ caustic}$',
fontsize=15.)
plot([-1.8,-.35],[1.85,0.35],'k-',lw=0.8)
annotate(text=r'$\mathrm{Source\ plane}$',xy=(0.5,1.01),xycoords='axes fraction',
ha='center',va='bottom',fontsize=18.)
tight_layout();
Figure 15.15: The critical curve and the caustic for the singular isothermal sphere.
Because there is only a single caustic, you might think that there is only a transition between one and two images for the flattened singular-isothermal sphere (two, because the odd-number theorem does not hold). However, in singular lenses, there is in fact a region similar to that encompassed by the radial caustic in the previous examples where there are two images, and passing the tangential caustic then creates four images (this is what happens in the Fermat-potential example above). The curve that separates the single-image region from the two-image region is known as the cut (Kovner 1987). Softening the singularity, the cut becomes the radial caustic. The cut can therefore be computed by softening the singularity (e.g., by replacing \(x^2 + y^2/b^2\) by \(x^2 + y^2/b^2 + R_c^2\)) and obtaining the cut as the radial caustic in the limit of zero softening (\(R_c \rightarrow 0\)). It is clear from this discussion that aside from the absence of the central image, lensing by singular mass distributions is essentially the same as that in non-singular lenses.
.
.
.