In [1]:
%matplotlib inline
from astropy import units,constants
import numpy as np
import evillens as evil
For this demo, we will use two SIE lenses to show the extent at which finding a numerical solution to one affects the bias in the estimated mass of the other. To do this, we will find the numeric solution to the deflection angles for the combination of the two lenses, and subtract the analytic solution to the deflection angles for the larger lens.
In [2]:
Main_lens = evil.AnalyticSIELens(0.25,2.85)
Subhalo = evil.AnalyticSIELens(0.25,2.85)
Main_lens.setup_grid(NX=400,NY=400,pixscale=0.05, n=4)
Subhalo.setup_grid(NX=400,NY=400,pixscale=0.05, n=4)
We will use a central velocity dispersion of 200 km / s for the main lens, and a central velocity dispersion of 40 km / s for the subhalo. This corresponds to masses of approximately _ and _ respectively. Also we use a moderately elliptical mass distribution for the main lens, and one that is nearly an SIE for the subhalo.
In [3]:
Main_lens.build_kappa_map(sigma=200.0 , q3 = 0.75 , inclination = 0.0 , rotation = 0.0 , centroid = [0.025,0.025], r_c=0.1)
Subhalo.build_kappa_map(sigma = 40.0, q3 = 0.99 , inclination = 0.0 , rotation = 0.0 , centroid = [0.525,0.625])
Main_lens.plot('kappa')
Subhalo.plot('kappa')
We then create a new lens called 'composite' which has kappa equal to the sum of the kappas from the main lens and subhalo, and numerically calculate its deflection angles.
In [4]:
Composite = Main_lens+Subhalo
Composite.plot('kappa')
Composite.deflect()
Now we can calculate the deflection angles from the Main_lens and subhalo (analytically and individually). Below are the deflection angles calculated for all three lenses.
In [5]:
Main_lens.deflect()
Subhalo.deflect()
Main_lens.plot('alpha')
Subhalo.plot('alpha')
Composite.plot('alpha')
Now create a 4th lens called Subhalo_numerical which is formed by subtracting Main_lens from Composite. This result will have the same kappa map as the subhalo, and should have the same alpha maps as the subhalo also. It will not however, since the numeric solution can only ever be approximate.
But since the deflection angles SHOULD have an analytic solution of the form
$\alpha_x = \frac{b}{\sqrt{1-q^2}} \tan{(....)}$
and
$\alpha_y = \frac{b}{\sqrt{1-q^2}} \tanh{( ..... )}$
and have the same deflection angle along any "radial" axis (we will use the x axis for simplicity), we can use the average of the x component of the deflection angles along the axis that intersects the center of the subhalo parallel to the x axis (for this particular example, since the pixel size is 0.02, the lower limit of y is -2.0 and the center of the subhalo is at y=0.61 this will be on the 130th x axis) to estimate the einstein radius of Subhalo_numerical. This can be easily translated into the Einstein radius mass using the relation:
$M_{Einstein} = \pi\Sigma_{crit}r_{Einstein}^2$ (Hezaveh et al., 2013, ApJ, 767, 9)
In [6]:
Subhalo_numerical = Composite - Main_lens
Subhalo_numerical.plot('kappa')
Subhalo_numerical.plot('alpha')
In [7]:
axis = np.argmin(abs(Subhalo_numerical.image_y[:,0]-Subhalo.centroid[1]))
axis2 = np.argmin(abs(Subhalo_numerical.image_x[0,:]-Subhalo.centroid[0]))
b = np.average(abs(Subhalo_numerical.alpha_x[axis,:]))/np.arctan(np.sqrt((1-Subhalo.q3**2)/Subhalo.q3**2))*np.sqrt(1-Subhalo.q3**2)
bsis = b/(np.sqrt(1-Subhalo.q3**2)/np.arcsin(np.sqrt(1-Subhalo.q3**2)))
M_actual = (((Subhalo.b*units.arcsec).to('', equivalencies=units.dimensionless_angles())*Subhalo.Dd)**2 * Subhalo.SigmaCrit *np.pi)
M_estimated = (((b*units.arcsec).to('', equivalencies=units.dimensionless_angles())*Subhalo_numerical.Dd)**2 * Subhalo_numerical.SigmaCrit *np.pi)
In [8]:
print('Estimated Mass:', M_estimated)
print('Actual Mass:', M_actual)
print('Percent error:', 100.0*(M_estimated-M_actual)/M_actual)
This suggests that the estimated mass is ~40 percent larger than would be expected from the deflection angles calculated. However if we instead plot the value of alpha_x along that axis, we see that most of the bias is due the impact of an increasing magnitude of the deflction angle far away from the subhalo.
In [9]:
import matplotlib.pyplot as plt
plt.plot((Subhalo_numerical.alpha_x[axis,:]),'b-')
Out[9]:
If we instead examine only the innermost 25 pixels on either side, we see that the measured bias is much smaller.
In [10]:
b = np.average(abs(Subhalo_numerical.alpha_x[axis,(axis2-25):(axis2+26)]))/np.arctan(np.sqrt((1-Subhalo.q3**2)/Subhalo.q3**2))*np.sqrt(1-Subhalo.q3**2)
bsis = b/(np.sqrt(1-Subhalo.q3**2)/np.arcsin(np.sqrt(1-Subhalo.q3**2)))
M_actual = (((Subhalo.b*units.arcsec).to('', equivalencies=units.dimensionless_angles())*Subhalo.Dd)**2 * Subhalo.SigmaCrit *np.pi)
M_estimated = (((b*units.arcsec).to('', equivalencies=units.dimensionless_angles())*Subhalo_numerical.Dd)**2 * Subhalo_numerical.SigmaCrit *np.pi)
print('Estimated Mass:', M_estimated)
print('Actual Mass:', M_actual)
print('percent error:', 100.0*(M_estimated-M_actual)/M_estimated)
It still overestimates the mass of the lens, but less significantly (~10%). So lets again examine closer to the halo. Now the innermost 15 pixels on either side.
In [11]:
b = np.average(abs(Subhalo_numerical.alpha_x[axis,(axis2-15):(axis2+16)]))/np.arctan(np.sqrt((1-Subhalo.q3**2)/Subhalo.q3**2))*np.sqrt(1-Subhalo.q3**2)
bsis = b/(np.sqrt(1-Subhalo.q3**2)/np.arcsin(np.sqrt(1-Subhalo.q3**2)))
M_actual = (((Subhalo.b*units.arcsec).to('', equivalencies=units.dimensionless_angles())*Subhalo.Dd)**2 * Subhalo.SigmaCrit *np.pi)
M_estimated = (((b*units.arcsec).to('', equivalencies=units.dimensionless_angles())*Subhalo_numerical.Dd)**2 * Subhalo_numerical.SigmaCrit *np.pi)
print('Estimated Mass:', M_estimated)
print('Actual Mass:', M_actual)
print('Percent Error:', 100.0*(M_estimated-M_actual)/M_estimated)
Now it underestimates the mass by only three percent. This is because the value of kappa is diverging close to the center. However, the error in the measured mass is not significantly reduced by ignoring the innermost region.
In [12]:
b = np.average([np.average(abs(Subhalo_numerical.alpha_x[axis,(axis2-15):axis2])),np.average(abs(Subhalo_numerical.alpha_x[axis,(axis2+1):(axis2+16)]))])/np.arctan(np.sqrt((1-Subhalo.q3**2)/Subhalo.q3**2))*np.sqrt(1-Subhalo.q3**2)
M_actual = (((Subhalo.b*units.arcsec).to('', equivalencies=units.dimensionless_angles())*Subhalo.Dd)**2 * Subhalo.SigmaCrit *np.pi)
M_estimated = (((b*units.arcsec).to('', equivalencies=units.dimensionless_angles())*Subhalo_numerical.Dd)**2 * Subhalo_numerical.SigmaCrit *np.pi)
print('Estimated Mass:', M_estimated)
print('Actual Mass:', M_actual)
print('Percent Error:', 100.0*(M_estimated-M_actual)/M_estimated)
Thus, it seems that ~3% error is about the best that can be done on a 400x400 kappa grid and a 100x100 image grid for a non-singular isothermal ellipsoid, which is pretty good given the mass ratio between the subhalos and the main objects (approximately 1:1000).
For the singular isothermal ellipsoid, possibly due to poorer sampling of the gradiant of the density distribution, we find that errors are approximately 20-25% for this method. Interestingly, computing the average over the whole array yields good accuracy (3%), but it is unclear whether this is because it is intrinsically more accurate, or because the errors near the ends coincidentally average against the errors in the middle. Future testing will enlighten us.
Fortunately or unfortunately (depending on perspective), The same accuracy will be recovered regardless of the mass of the perturber (since the biases are caused by errors in alpha from the main halo). Thus, larger biases will be found for larger main lenses, and better accuracy is achieved for smaller central objects.
Also, the magnitude of the errors on the returned mass of a substructure seems to be independant of position in the image. Since the uncertainties of the x direction deflection angles increase (slowly) with increasing x position, the error is expected to become significant at some point. However, this is not expected to become relevant until significantly outside of the Einstein radius, again in areas where we do not expect to find lensed images anyways.