15.1. The lensing equation¶
15.1.1. The deflection angle¶
As discussed above, the study of the deflection of light in gravitational fields is a phenomenon that requires the general theory of relativity. In the limit of weak gravitational fields—\(|\Phi/c^2| \ll 1\) for \(\Phi \rightarrow 0\) as \(r\rightarrow \infty\), always the case on galactic scales—luckily for us, the GR calculation of the deflection experienced by a light ray traversing a gravitational potential \(\Phi\) is simply given by \begin{align}\label{eq-gravlens-light-bend-alpha} \hat{\boldsymbol{\alpha}} & = 2\int \mathrm{d} s\, \nabla_\perp \left({\Phi \over c^2}\right)\,, \end{align} where \(\hat{\boldsymbol{\alpha}}\) is the two-dimensional deflection angle perpendicular to the propagation direction, \(\nabla_\perp f= \nabla f - k^{-2}\,(\vec{k}\cdot \nabla f)\,\vec{k}\) is the gradient perpendicular to the unperturbed photon direction \(\vec{k}\), and the integral is over the unperturbed trajectory. This equation is derived in detail in Appendix C.3.
The simplest application of Equation (15.1) is to the case of light deflection by a point mass potential \(\Phi = -GM/r\). This describes, for example, the deflection of light grazing the surface of the Sun as measured in the 1919 eclipse. In the formalism below, we define the coordinate system such that the unperturbed light ray travels along the \(z\) direction and the vector \(\vec{x} = (x,y)\) describes the plane perpendicular to this direction. Assuming that the light ray originates far from the point mass at \(z=-\infty\) and is also observed far from the point mass at \(z=+\infty\), we then have that \begin{align} \hat{\boldsymbol{\alpha}} & = {2\over c^2}\int_{-\infty}^{+\infty} \mathrm{d} z\, \begin{pmatrix} \partial \Phi / \partial x \\ \partial \Phi / \partial y\end{pmatrix} = {2GM\over c^2}\,\vec{x}\,\int_{-\infty}^{+\infty} \mathrm{d} z\, {1\over (b^2 + z^2)^{3/2}}= {4GM\over b c^2}\,\vec{\hat{x}}\,,\label{eq-gravlens-deflect-point} \end{align}
where \(b = \sqrt{x^2+y^2}\) and \(\vec{\hat{x}}\) is the unit vector pointing towards \((x,y)\). Thus, the deflection points towards \(\vec{x}\): we observe a background object further from the point-mass at \((x,y) = (0,0)\) than it would be without the gravitational light deflection. The magnitude of \(\hat{\boldsymbol{\alpha}}\) for \(M=M_\odot\) and \(b = R_\odot\) is the value predicted by Einstein (1916) and confirmed by Dyson et al. (1920):
[6]:
from astropy.constants import G, c
print((4.*G*u.Msun/u.Rsun/c**2.)\
.to(u.arcsec,equivalencies=u.dimensionless_angles()))
1.751190325559984 arcsec
Thus, we find that \(|\hat{\boldsymbol{\alpha}}| = 1.75''\).
From the calculation above, it is immediately clear that if \(\Phi(x,y,z)\) does not depend on the angle \(\phi\) in the \((x,y)\) plane [that is, if \(\Phi(x,y,z) \equiv \Phi(R,z)\)], then \(\hat{\boldsymbol{\alpha}} \parallel \vec{x}\) and, thus, for axially-symmetric lenses, the deflection is in the direction pointing from the lens’ center to the object and is given by \begin{align}\label{eq-gravlens-light-bend-alpha-axi} \hat{\boldsymbol{\alpha}} & = {2\over c^2}\,\vec{x}\int_{-\infty}^{+\infty} \mathrm{d} z\, {1\over R}{\partial \Phi \over \partial R}\,, \end{align} where \(\partial \Phi / \partial R\) is the derivative of the gravitational potential with respect to the cylindrical radius in the \((x,y)\) plane. If the potential is spherically symmetric, we can also compute this as \(\hat{\boldsymbol{\alpha}} = {2\over c^2}\,\vec{x}\int_{-\infty}^{+\infty} \mathrm{d} z\, {1\over r}{\partial \Phi \over \partial r}\).
Galaxies are extended objects and for the purpose of lensing are not well-modeled as point masses. As discussed in Chapter 2.4.5, because galaxies are observed to have flat rotation curves (Chapter 8), a simple one-parameter model for a galaxy is that of a logarithmic potential \(\Phi = v_c^2 \ln r\) or, equivalently, the potential of the singular isothermal sphere \(\Phi = 2\sigma^2 \ln r\) (see Chapter 5.6.2) ; in the context of lensing, the logarithmic potential is typically referred to as the singular isothermal sphere. For this potential, Equation (15.3) gives a deflection of \begin{align} \hat{\boldsymbol{\alpha}} & = {2 v_c^2\over c^2}\,\vec{x}\,\,\int_{-\infty}^{+\infty} \mathrm{d} z\, {1\over r^2} = {2 v_c^2\over c^2}\,\vec{x}\,\,\int_{-\infty}^{+\infty} \mathrm{d} z\, {1\over (b^2 + z^2)} = 2\pi \left({v_c \over c}\right)^2\,\vec{\hat{x}}= 4\pi \left({\sigma \over c}\right)^2\,\vec{\hat{x}}\,,\label{eq-gravlens-deflect-sis} \end{align}
where we have also written the result in terms of the velocity dispersion of the singular isothermal sphere. Massive galaxies like the Milky Way or elliptical galaxies have \(\sigma \approx 200\, \mathrm{km\,s}^{-1}\) and therefore give rise to deflections with typical sizes of \(\approx 1''\):
[7]:
print((4.*numpy.pi*(200.*u.km/u.s/c)**2.)\
.to(u.arcsec,equivalencies=u.dimensionless_angles()))
1.1535955781163916 arcsec
Such angular deflections are at the edge of the range that large ground-based telescopes can resolve; space telescopes such as the Hubble Space Telescope can resolve these with relative ease. Note that the deflection angle for the singular isothermal sphere does not depend on the position of the source!
Calculating deflection angles using Equation (15.1) is simple for any gravitational potential and we implement the general expression in the following function:
[8]:
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())
For example, for an NFW halo with virial mass \(M_\mathrm{vir} = 10^{13}\,M_\odot\)—a typical mass for an elliptical lens—and a concentration of 6.63 following the Dutton & Macciò (2014) relation as in Chapter 2.4.6, we find the deflection angle as a function of separation shown in Figure 15.2.
[9]:
from galpy.potential import NFWPotential
np= NFWPotential(mvir=10,conc=6.63,
H=67.1,overdens=200.,wrtcrit=True,
ro=8.,vo=220.)
figure(figsize=(6,4))
Rs= numpy.linspace(0.01*u.kpc,np.rvir(H=67.1,overdens=200.,wrtcrit=True),101)
plot(Rs,u.Quantity([deflection(np,R)[0].to(u.arcsec) for R in Rs]))
xlabel(r'$R\,(\mathrm{kpc})$')
ylabel(r'$\widehat{\boldsymbol{\alpha}}\,(\mathrm{arcsec})$');
Figure 15.2: Deflection angle for a \(M_\mathrm{vir} = 10^{13}\,M_\odot\) NFW halo, a typical massive elliptical.
The value of the deflection angle is similar to that for the singular-isothermal-sphere example above, because the typical velocity dispersion of a \(M_\mathrm{vir} = 10^{13}\,M_\odot\) is \(\approx 200\,\mathrm{km\,s}^{-1}\).
15.1.2. The lensing equation¶
Knowledge of the deflection angle \(\hat{\boldsymbol{\alpha}}\) allows us to relate the unperturbed position on the sky of a source to the observed, lensed position. Assuming that light only travels through a single gravitational field (e.g., that of one star or one galaxy) strong enough to observably deflect it, the geometry is demonstrated in Figure 15.3.
[10]:
def angle_plot(line1,line2,offset=1,color='k',origin=(0, 0),
len_x_axis=1,len_y_axis=1):
# Draw angle arc between two lines
# Edited from https://stackoverflow.com/a/25228427 and
# https://gist.github.com/battlecook/0c0bdb7097ec7c8fa160e342b1bf51ef
from matplotlib.patches import Arc
# Angle between line1 and x-axis
l1xy= line1.get_xydata()
slope1= (l1xy[1][1]-l1xy[0][1])/float(l1xy[1][0]-l1xy[0][0])
angle1= numpy.degrees(numpy.arctan(slope1))
# Angle between line2 and x-axis
l2xy= line2.get_xydata()
slope2= (l2xy[1][1]-l2xy[0][1])/float(l2xy[1][0]-l2xy[0][0])
angle2= numpy.degrees(numpy.arctan(slope2))
# Angle between them
theta1= numpy.amin([angle1,angle2])
theta2= numpy.amax([angle1,angle2])
angle= theta2-theta1
return Arc(origin,len_x_axis*offset,len_y_axis*offset,
angle=0,theta1=theta1,theta2=theta2,
color=color,label=str(angle)+u"\u00b0")
def line_between_points(o1,o2,*args,**kwargs):
return plot([o1[0],o2[0]],[o1[1],o2[1]],*args,**kwargs)
def extend_line(o1,o2,xnew,*args,**kwargs):
"""Extend the line between o1 and o2 to xnew"""
slope= (o2[1]-o1[1])/(o2[0]-o1[0])
ynew= slope*(xnew-o1[0])+o1[1]
return line_between_points(o1,[xnew,ynew],*args,**kwargs)
figure(figsize=(11,5.5))
# Source, lens, observer position
source= [0.,.5]
lens= [1.1,.75]
observer= [2.5,0.]
line_os= line_between_points(source,observer,'k--')[0]
line_sl= line_between_points(source,lens,'k-')[0]
line_ol= line_between_points(observer,lens,'k-')[0]
gca().add_patch(angle_plot(line_ol,line_os,-1.5,origin=observer))
line_beta= line_between_points([0,0],observer,'k:',lw=2.)[0]
gca().add_patch(angle_plot(line_beta,line_os,-2.,origin=observer))
gca().add_patch(angle_plot(line_beta,line_ol,-.75,origin=observer))
line_ex= extend_line(lens,observer,0.65,'k:',lw=2.)[0]
gca().add_patch(angle_plot(line_ex,line_sl,-0.5,origin=lens))
# Angle labels
galpy_plot.text(0.8,0.8,r'$\hat{\alpha}$',ha='center',va='center',fontsize=18.)
galpy_plot.text(1.74,0.3,r'$\alpha$',ha='center',va='center',fontsize=18.)
galpy_plot.text(1.45,0.1,r'$\beta$',ha='center',va='center',fontsize=18.)
galpy_plot.text(2.1,0.13,r'$\theta$',ha='center',va='center',fontsize=18.)
# Some dots
plot([source[0],lens[0],observer[0],lens[0]],[source[1],lens[1],observer[1],0.],
'ko',ms=10.)
# Impact parameter
annotate(text='',xy=(lens[0],0.02),xytext=(lens[0],lens[1]-0.02),
arrowprops=dict(arrowstyle='<|-|>',color='k'))
galpy_plot.text(lens[0],lens[1]/2.,r'$\xi$',
ha='center',va='center',fontsize=18.,backgroundcolor='w')
# Distance arrows
annotate(text='',xy=(source[0],-0.1),xytext=(lens[0],-0.1),
arrowprops=dict(arrowstyle='<|-|>',color='k'))
galpy_plot.text(0.5*(source[0]+lens[0]),-0.1,r'$D_{\mathrm{LS}}$',
ha='center',va='center',fontsize=18.,backgroundcolor='w')
annotate(text='',xy=(observer[0],-0.1),xytext=(lens[0],-0.1),
arrowprops=dict(arrowstyle='<|-|>',color='k'))
galpy_plot.text(0.5*(observer[0]+lens[0]),-0.1,r'$D_{\mathrm{L}}$',
ha='center',va='center',fontsize=18.,backgroundcolor='w')
annotate(text='',xy=(source[0],-0.175),xytext=(observer[0],-0.175),
arrowprops=dict(arrowstyle='<|-|>',color='k'))
galpy_plot.text(0.5*(source[0]+observer[0]),-0.175,r'$D_{\mathrm{S}}$',
ha='center',va='center',fontsize=18.,backgroundcolor='w')
# Observer, source, lens annotations
galpy_plot.text(source[0],source[1]+0.05,r'$\mathrm{Source}$',
ha='center',va='bottom',fontsize=18.)
galpy_plot.text(lens[0]-0.1,0.05,r'$\mathrm{Lens}$',
ha='center',va='bottom',fontsize=18.)
galpy_plot.text(observer[0],observer[1]+0.1,r'$\mathrm{Observer}$',
ha='center',va='bottom',fontsize=18.)
xlim(-0.28,2.6)
ylim(-.3,1.1)
gca().axis('off');
Figure 15.3: The basic geometry of gravitational lensing.
In this figure, the source is on the left, the lens in the middle, and the observer is on the right. In all galactic applications of gravitational lensing, the extent of the lens is much smaller than the distances between the source, lens, and observer. We can therefore approximate the lens as being infinitely thin in the direction perpendicular to the light ray. This is the thin screen approximation. We can then approximate the perturbed light path in the manner of Figure 15.3, with light propagating in a straight line until it reaches the lens at a closest approach \(\boldsymbol{\xi}\) where its direction is changed by the deflection angle; after this, the ray again propagates on a straight line until it is observed. In the absence of the lens, we would observe the source at an angular position \(\boldsymbol{\beta}\) on the source plane, but because of the deflection, we actually observe it at \(\boldsymbol{\theta}\). Figure 15.3 demonstrates how these angles are related in a two dimensional slice of the full three-dimensional picture. Because the angles involved are small, the relation between the angles is \begin{equation}\label{eq-gravlens-lensing-eq} \boldsymbol{\beta} = \boldsymbol{\theta}-{D_{\mathrm{LS}}\over D_{\mathrm{S}}}\,\hat{\boldsymbol{\alpha}}(\boldsymbol{\xi})\,. \end{equation} (Figure 15.3 shows that this equation holds for both the \(x\) and \(y\) components of the angles and therefore the vector form follows). This is the lensing equation (Bourassa & Kantowski 1975). Because of the way the relative distances appear in this equation, it is common to define the reduced deflection angle \(\boldsymbol{\alpha}\) as \begin{equation}\label{eq-gravlens-reduced-deflect} \boldsymbol{\alpha}(\boldsymbol{\theta}) = {D_{\mathrm{LS}}\over D_{\mathrm{S}}}\,\hat{\boldsymbol{\alpha}}(\boldsymbol{\theta})\,, \end{equation} where \(\hat{\boldsymbol{\alpha}}\) is a function of \(\boldsymbol{\theta}\) through \(\boldsymbol{\xi} = D_\mathrm{L}\boldsymbol{\theta}\). Then, the lensing equation can be written as \begin{equation}\label{eq-gravlens-lensing-eq-reduced} \boldsymbol{\beta} = \boldsymbol{\theta}-\boldsymbol{\alpha}(\boldsymbol{\theta})\,. \end{equation} The distances that appear in the lensing equation are the regular Euclidean distance in the case of lensing within our own Galaxy (or within the local Universe), but for cosmological lenses (such as a high-redshift background quasar lensed by an intermediate-redshift galaxy), the distances that appear are angular diameter distances. The angular diameter distance depends on the energy and matter content of the Universe and on its spatial curvature. The angular diameter distance is discussed in more detail in Appendix C.4.3. Note that on cosmological scales, it is not the case that \(D_{\mathrm{S}} = D_{\mathrm{LS}} + D_{\mathrm{L}}\).
For strong-lensing systems, it is generally a good assumption that there is only a single deflector along the line of sight (the lens), but in the weak-lensing regime (see below) or for certain more detailed applications in strong lenses, it is necessary to consider the effect of additional deflectors in front of or behind the main lens. In that case, we have to deal with multi-plane lensing. The lensing equation in this case can be straightforwardly derived by looking at the single-plane lensing equation (15.5) as mapping an image onto its source location; in multi-plane lensing the source obtained by a single application of Equation (15.5) applied to the closest lens to the observer, becomes the image from the second-closest lens and so on. Thus, the lensing equation is replaced by a set of equations (e.g., Blandford & Narayan 1986; Kovner 1987) \begin{equation} \boldsymbol{\beta}_{i+1} = \boldsymbol{\theta}-{1\over D_{\mathrm{S}}}\sum_{k=1}^i{D_{\mathrm{kS}}\,\hat{\boldsymbol{\alpha}}(D_k\,\boldsymbol{\beta}_k)}\,, \end{equation} where \(D_k\) is the (angular diameter) distance to lens \(k\) and \(D_{\mathrm{kS}}\) that between lens \(k\) and the source; the index \(i\) runs from \(i=1\) to \(N\) for \(N\) lenses and we set \(\boldsymbol{\beta}_1 = \boldsymbol{\theta}\) (thus, the index runs backwards and the source is at \(\boldsymbol{\beta}_{N+1}\)). While this equation looks simple, solving it is much more difficult than solving the single-plane lensing equation. In this chapter, we will only consider the single-plane lensing equation, but a large literature on solving the multi-plane lensing equation in various regimes exists (e.g., Keeton et al. 1997; Schneider 1997; Birrer et al. 2017; Schneider 2019).
15.1.3. Examples: point mass and singular isothermal sphere¶
In the case of a point-mass lens, we have from Equation (15.2) that \begin{equation} \boldsymbol{\alpha}(\boldsymbol{\theta}) = {D_{\mathrm{LS}}\over D_{\mathrm{S}}D_{\mathrm{L}}}\,{4GM\over c^2}\,{\boldsymbol{\theta}\over |\boldsymbol{\theta}|^2}\,, \end{equation} where we have substituted \(b = D_\mathrm{L}|\boldsymbol{\theta}|\) and \(\vec{\hat{x}} = \boldsymbol{\theta}/|\boldsymbol{\theta}|\). From the discussion above, we know that for an axially-symmetric lens, the deflection is along the line connecting \(\boldsymbol{\beta}\) and the center of the lens, and the lensing equation for axially-symmetric lenses therefore reduces to a one-dimensional equation for the position \(\theta\) on this line (at least when \(\boldsymbol{\beta} \neq 0\)). For the point-mass lens this becomes \begin{equation}\label{eq-gravlens-lensing-eq-reduced-pm} \beta = \theta-{D_{\mathrm{LS}}\over D_{\mathrm{S}}D_{\mathrm{L}}}\,{4GM\over c^2}\,{1\over \theta}\,. \end{equation} To solve this equation, we first define the Einstein angle \(\theta_\mathrm{E}\) as \begin{equation}\label{eq-gravlens-einstein-angle-point-mass} \theta_\mathrm{E} = \sqrt{{D_{\mathrm{LS}}\over D_{\mathrm{S}}D_{\mathrm{L}}}\,{4GM\over c^2}}\,, \end{equation} because then we have that \(0= \theta^2-\beta\theta -\theta_\mathrm{E}^2\). This has the solutions \begin{equation}\label{eq-gravlens-lensing-eq-sol-pm} \theta = {1\over 2}\,\left(\beta\pm\sqrt{\beta^2+4\theta_\mathrm{E}^2}\right)\,. \end{equation} When \(\boldsymbol{\beta} = 0\), this gives \(\theta = \theta_\mathrm{E}\). However, in this case, our reduction of the problem to a one-dimensional equation does not hold, and in fact the entire ring \(|\boldsymbol{\theta}| = \theta_\mathrm{E}\) solves the lensing equation. This is an Einstein ring with \(\theta_\mathrm{E}\) in this context known as the Einstein radius: a source directly behind a point-like object gets lensed into a ring.
We can similarly solve the lensing equation for lensing by a singular isothermal sphere. From Equation (15.4), we have that \begin{equation}\label{eq-gravlens-reduceddeflect-sis} \boldsymbol{\alpha}(\boldsymbol{\theta}) = 4\pi {D_{\mathrm{LS}}\over D_{\mathrm{S}}}\,\left({\sigma \over c}\right)^2\,{\boldsymbol{\theta}\over |\boldsymbol{\theta}|}\,. \end{equation} Again, as long as \(\boldsymbol{\beta} \neq 0\), the lensing equation reduces to one-dimensional form and can now be written as \(\beta = \theta-\theta_\mathrm{E}\,\theta/|\theta|\), where the Einstein angle is given by \begin{equation}\label{eq-gravlens-sis-einstein-angle} \theta_\mathrm{E} = 4\pi {D_{\mathrm{LS}}\over D_{\mathrm{S}}}\,\left({\sigma \over c}\right)^2\,. \end{equation} The lensing equation has a single solution \(\theta = \beta + \mathrm{sign}(\beta)\,\theta_\mathrm{E}\) if \(|\beta| \geq \theta_\mathrm{E}\), but when \(|\beta| < \theta_\mathrm{E}\) it has two solutions \(\theta = \beta \pm\theta_\mathrm{E}\). When \(\boldsymbol{\beta} = 0\), we again find that the entire ring \(|\boldsymbol{\theta}| = \theta_\mathrm{E}\) satisfies the lensing equation and we again get an Einstein ring.
The point-mass and singular-isothermal-sphere solutions illustrate several important features of gravitational lensing. The first is that gravitational lensing can lead to the appearance of multiple images: the point-mass lens always has two solutions corresponding to two apparent positions for every source, while the singular isothermal sphere multiply-images sources that fall within its Einstein ring. The second is that gravitational lensing can lens a compact source into a one-dimensional Einstein ring if the source and the lens are aligned. More generally, the fact that the lensing equation for axially-symmetric lenses reduces to a one-dimensional equation (except when an Einstein ring forms) means that quadrupolar, four-image configurations (“quads”) such as the “Einstein cross” (Huchra et al. 1985) cannot result from axially-symmetric lenses. We also see that every source in the sky is technically lensed by a single point-like lens or by a single singular isothermal sphere: No matter how distant, the gravitational field of a singular isothermal sphere moves the image position by \(\theta_\mathrm{E}\). However, at large distances from the centers of galaxies or clusters of galaxies, the density distribution drops off faster and a point-mass potential is a better model. At large separations \(|\beta| \gg |\theta_\mathrm{E}|\), Equation (15.12) shows that the two images approach \(\theta_1 \rightarrow \beta\) and \(\theta_2 \rightarrow 0\) and one of the images therefore overlaps with the lens. Moreover, as we will see in an exercise, this image is also strongly demagnified and therefore almost impossible to observe.
.
.
.