Single lenses
Gravity.jl uses a hierarchical model for lenses. A multiple-lens system is implemented as a LensSystem, which is essentially a list of lens planes at different redshifts. A lens plane is implemented as a LensPlane, and this is essentially a list of single lenses. The various lenses in a LensPlane are concrete types all belonging to the abstract type Gravity.SingleLens.
Single Lens
All lens models are implemented as subtypes of Gravity.SingleLens. This is an important concept in Gravity.jl: a single lens is essentially a type that allows one to compute the time delay, the deflection angle, and the Jacobian of the ray-tracing.
Gravity.AbstractLens — Type
Abstract type AbstractLens, the parent of all lenses.
Gravity.SingleLens — Type
Abstract type SingleLens, the parent of all (single) lens models.
This type is important because it is used as a building block for a lens plane, i.e. a collection of lenses at the same redshift (for example, the galaxies and the DM halo of a galaxy cluster).
One can add more lens models by creating a concrete subtype of SingleLens. For concrete types, one must define at least the potential method, which will return the dimensionless lensing potential. Optionally, one can also define the deflection method (returning the deflection angle for a source at redshift ∞) and the jacobian method (returning the jacobian of the deflection angle). If either of these methods are not defined, they are obtained through automatic differentiation (see LensTrait).
Spherical lenses
Many lens models are based on spherical (i.e., axial) symmetry, so that the lens mass distribution is a function of the distance from the lens center only. In this case, the lens equations can be greatly simplified using Gauss' theorem: moreover, the lens deflection is radial, and the Jacobian can be computed in terms of the deflection and of the lens density.
A new spherical lens, therefore, just need to define a few functions to have all lensing methods correctly defined.
Gravity.SphericalSingleLens — Type
Abstract type SphericalSingleLens, the parent of all (single) axisymmetric lens models.
An axisymmetric lens has a mass density that is invariant upon rotations; similarly, the potential, the deflection, and the jacobian will transform as a scalar, a vector, or a 2-tensor with respect to rotations.
Members
Any concrete type must define at least the following members
z: the lens redshiftx₀: the lens center
Methods
When defining a concrete type, it is necessary to provide the following methods:
potential2(lens, r²): the lens potential, having as parameter the square of the distance from the lens centerdeflection2(lens, r²): the scalar deflection, having as parameter the square of the distance from the lens centerconvergence2(lens, r²): the lens mass density, having as parameter the square of the distance from the lens center
The standard methods potential, deflection, jacobian, and convergence, for all SphericalSingleLens'es, are automatically defined in terms of the methods in the items above.
Additionally, if one wants to use [EllipticalSingleLens], it is necessary to define also
convergence2′(lens, r²): the derivative ofconvergence2innerslope(lens): the lens logaritmic inner slope
For non-fused broadcasting, one can define the fast_deflection2 method: this is identical to deflection2, but it is used in broadcasting operations and, depending on the equations involved, it might be wise to provide a more efficient implementation.
Optionally, one can also define einsteinradius to compute the Einstein radius at the fiducial distance; if not provided, the Einstein radius is obtained using Newton's method to solve the equation x = deflection(lens, x).
Elliptical lenses
More frequently, lenses have an elliptical symmetry: the mass density depends on an elliptical radius, defined as
\[\xi \equiv \sqrt{\tilde x/q^2 + \tilde y} \; ,\]
where $(\tilde x, \tilde y)$ are Cartesian coordinates in a coordinate system aligned with the lens axes, and $q = b / a$ is the ratio of semi-axes of the ellipse describing the lens distribution. The lensing properties of lenses with elliptical mass distribution, in some cases, can be computed analytically; more often, one needs to perform one-dimensional integrals (see Gravity.EllipticalLens).
Alternatively, one can use models characterized by elliptical symmetry in their potential functions (see Gravity.EllipticalPotLens): in this case, assuming that the corresponding spherical model has known deflection, no integration is needed. Note, however, that these models generally produce mass distributions that are only approximately elliptical, and only for small ellipticities.
In either case, lenses with elliptical symmetry can be better defined in a coordinate system aligned with the lens axes. This technique is implemented using the following abstract class.
Gravity.EllipticalSingleLens — Type
Abstract type EllipticalSingleLens, the parent of all (single) elliptical lens models.
An elliptical lens has a mass density or a lens potential that depends on an elliptical coordinate.
Members
Any concrete type must define at least the following members
z: the lens redshiftx₀: the lens centerq: the axis ratio, computed as $q = b / a$, whereaandbare the two semi-axes (with a ≥ b)θₑ: the lens position angle, computed clockwise from the top (North to West in astronomical sense)cosθₑ: cosine ofθₑsinθₑ: sine ofθₑ
Methods
When defining a concrete type, it is recommended to provide a set of methods that compute the relevant lensing quantities in the rotated coordinate system. In other terms, these methods will accept a parameter δ̃ defined as the vector from the lens center to a point in a coordinate system that is aligned with the lens position angle.
potential∟(lens, δ̃): the lens potential in the lens natural coordinate systemdeflection∟(lens, δ̃): the lens deflection in the lens natural coordinate systemjacobian∟(lens, δ̃): the lens jacobian in the lens natural coordinate systemdeflectionjacobian∟(lens, δ̃): the lens deflection and jacobian in the lens natural coordinate system.convergence∟(lens, , δ̃): the lens convergence in the lens natural coordinate system
The standard methods potential, deflection, jacobian, and convergence, for all EllipticalSingleLens'es, are then automatically defined in terms of the methods in the items above.
Lensing methods
A new lens model should provide various lensing-related methods needed for the various calculations.
Gravity.potential — Function
potential(lens, point)Compute the lens potential at a given point.
The same method can be applied to a LensPlane, in which case the potential is the sum of the potentials of all lenses in the plane.
Gravity.deflection — Function
deflection(lens, point)Compute the lens deflection (that is, the quantity α̂) at a given point.
The same method can be applied to a LensPlane, in which case the deflection is the sum of the deflections of all lenses in the plane.
The result should have the same type as the point type.
Gravity.deflection! — Function
deflection!(α̂s, lens, points [, d=1]; use_threads=false)Compute in place the lens deflection (that is, the quantity α̂) for an array of points.
The same method can be applied to a LensPlane, in which case the deflection is the sum of the deflections of all lenses in the plane.
The result is added to the α̂s array, which should have the type of the points array. Optionally, the result is rescaled by d before being added.
If use_threads is true and length(points) > 1000, the computation will be parallelized using threads: no threads are used for less than 1000 points, as the overhead of thread creation is generally not worth it in this case.
Gravity.jacobian — Function
jacobian(lens, point)Compute the lens jacobian (that is, the partial derivative of α̂) at a given point.
The same method can be applied to a LensPlane, in which case the jacobian is the sum of the jacobians of all lenses in the plane.
The result will generally be an SPointMap, i.e. a static 2×2 matrix.
Gravity.convergence — Function
convergence(lens, point)Compute the lens mass density at a given point.
The convergence is in dimensionless units and is referred to a critical density in the limit where the observer-source distance is equal to the lens-source distance. One can convert this quantity into a dimensional one by multiplying it by the critical density, which can be computed with Gravity.critical_density, or by using Gravity.massdensity.
The same method can be applied to a LensPlane, in which case the convergence is the sum of the convergences of all lenses in the plane.
In some cases, it might be faster to compute at the same time the deflection and the Jacobian of the deflection, instead of calling two separate functions. A lens model can therefore implement this as a single method:
Gravity.deflectionjacobian — Function
deflectionjacobian(lens, point)Compute the lens deflection and jacobian at a given point.
The result will be an SPointDiff, i.e. a static 6-vector: the first two elements represent the deflection α, and the last four its jacobian ∂α/∂θ.
The same method can be applied to a LensPlane, in which case the deflection and jacobian are the sums of the deflections and jacobians of all lenses in the plane.
Lensing methods for spherical lenses
Spherical lenses can be made elliptical using the technique described in Gravity.EllipticalSingleLens. To this purpose, it is necessary to define some additional lensing methods.
Gravity.potential2 — Function
potential2(lens, r²)Return the lensing potential of a SphericalSingleLens.
This function uses r² as an argument. This is often more efficient than using just r and is what is needed by EllipticalLens.
See also: deflection2, convergence2, convergence2′.
Gravity.deflection2 — Function
deflection2(lens, r²)Return the scalar deflection of a SphericalSingleLens.
The computation is carried out over r using r² as an argument: that is,
\[ \mathtt{deflection2}(r^2) = \alpha(\sqrt{r^2}) / \sqrt{r^2} \; .\]
This is often more efficient than using just r and is what is needed by EllipticalLens.
See also: potential2, convergence2, deflectionconvergence2, convergence2′.
Gravity.fast_deflection2 — Function
fast_deflection2(lens, r²)Return the scalar deflection of a SphericalSingleLens.
This method is used only during broadcasting operations by deflection! and (indirectly) lensmap!. Should be defined so to use, whenever possible, SIMD instructions.
If undefined for a specific concrete SphericalSingleLens, then deflection2 is used instead.
Gravity.convergence2 — Function
convergence2(lens, r²)Return the mass density of a SphericalSingleLens.
This function uses r² as an argument. This is often more efficient than using just r and is what is needed by EllipticalLens.
See also: potential2, deflection2, deflectionconvergence2, convergence2′.
Gravity.deflectionconvergence2 — Function
deflectionconvergence2(lens, r²)Return the both the deflection and the mass density of a SphericalSingleLens.
This function uses r² as an argument and return a 2-tuple. It is equivalent to (deflection2(lens, r²), convergence2(r²)), but is usually faster than computing both methods.
See also: potential2, deflection2, convergence2, convergence2′.
Gravity.convergence2′ — Function
convergence2′(lens, r²)Return the derivative of the mass density of a SphericalSingleLens.
This function uses r² as an argument. This is often more efficient than using just r and is what is needed by EllipticalLens.
This method must return $\partial \kappa(r^2) / \partial r^2$, that is the derivative of the mass density with respect to the square of the radius.
See also: potential2, deflection2, convergence2.
Gravity.innerslope — Function
innerslope(lens)Return the logarimic slope of the [massdistribution] of a SphericalSingleLens at its center: that is
\[\alpha = \left. \frac{\partial \log \kappa(r) }{\partial \log r} \right|_{r = 0} \; ,\]
or, equivalently, $\kappa (r) \approx r^\alpha$ for $r \rightarrow 0$.
The result must be returned encapsulated inside a Val object, as in Val(-1).
The definition of these methods for a given spherical lens type, ensures that we are able to compute the lensing potential, vector deflection, Jacobian, and mass density without the need to define the corresponding methods Gravity.potential, Gravity.deflection, and Gravity.jacobian.
For spherical lenses we can also compute the Einstein radius at infinity or at a given redshift: this can be useful to have an idea of the lensing area. The associated method can be implemented for each lens type, but a default implementation that solves the lens equation is provided.
Gravity.einsteinradius — Function
einsteinradius(lens [; guess=(ϵ, 60), tolerance=√ϵ, maxiters=100])
einsteinradius(lens, [zₒ=0,] zₛ [; guess=(ϵ, 60),..., cosmo=defaulltcosmology])
einsteinradius(lens, -f [; guess=(ϵ, 30),...])Return the Einstein radius of a spherical or elliptical lens.
The Einstein radius, defined as the solution of the equation $x = deflection(lens, x)$, is a relevant scale in the image plane. Each SphericalSingleLens can either define this method (returning, for example, the result of an analytical calculation), or can rely on a fallback implementation that solves the equation above iteratively.
For elliptical lenses, the Einstein radius is computed as the geometric mean of the result of einsteinradii (i.e., of the two radii of the associated critical curve).
The first version computes the Einstein radius assuming that f = Dₗₛ / Dₛ = 1. The second version computes the Einstein radius for a source at redshift zₛ as observed by a source at zₒ, in the provided cosmology. Finally, if the second argument is negative, it is taken directly as -f, the ratio between the lens-source distance and the observer-source distance.
Keywords
guess: for spherical lenses, the starting point for the algorithm iteration (default: 1 arcsec); for elliptical lenses, the interval(min, max)used for the braketing search (default: from machine precision to 60 arcsec), or just the upper range if a single number is providedtolerance: the tolerance used to search for zeros of the magnification matrixmaxiters: the maximum numbers of iteration for the search of zeros
Lensing methods for elliptical lenses
The lensing properties of elliptical lenses are often more easily described in their natural coordinate system, aligned with the lens axes. For this reason, it is recommended to define the following methods for all elliptical lens models:
Gravity.potential∟ — Function
potential∟(lens, point)Return the potential of an elliptical lens in its natural coordinate system.
potential uses by default this function for an EllipticalSingleLens.
Gravity.deflection∟ — Function
deflection∟(lens, point)Return the deflection for an elliptical lens in its natural coordinate system.
deflection uses by default this function for an EllipticalSingleLens.
Gravity.jacobian∟ — Function
jacobian∟(lens, point)Return the jacobian for an elliptical lens in its natural coordinate system.
jacobian uses by default this function for an EllipticalSingleLens.
Gravity.deflectionjacobian∟ — Function
deflectionjacobian∟(lens, point)Return the deflection and its jacobian for an elliptical lens in its natural coordinate system.
deflectionjacobian uses by default this function for an EllipticalSingleLens.
Gravity.convergence∟ — Function
convergence∟(lens, point)Return the mass density for an elliptical lens in its natural coordinate system.
convergence uses by default this function for an EllipticalSingleLens.
Gravity.einsteinradii — Function
einsteinradii(lens [; guess=(ϵ, 30), tolerance=√ϵ, maxiters=100])
einsteinradii(lens, [zₒ=0,] zₛ [; guess=(ϵ, 30),..., cosmo=defaulltcosmology])
einsteinradii(lens, -f [; guess=(ϵ, 30),...])Return the size of the critical lines around a spherical or elliptical lens.
This function computes the distance of the critical lines along the semimajor and semiminor axes of the lens, and returns them as a tuple. For a circular lens, the function returns two identical values, equal to the einsteinradius of the lens.
The first version computes the Einstein radii assuming that f = Dₗₛ / Dₛ = 1. The second version computes the Einstein radii for a source at redshift zₛ as observed by a source at zₒ, in the provided cosmology. Finally, if the second argument is negative, it is taken directly as -f, the ratio between the lens-source distance and the observer-source distance.
Keywords
guess: an interval used for the braketing search of the Einsting ring (default: from machine prcecision to 60 arcsec); if a single number is provided, it is taken to be the upper rangetolerance: the tolerance used to search for zeros of the magnification matrixmaxiters: the maximum numbers of iteration for the search of zeros
The function uses Roots.find_zeros as root-finder routine, and all other keywords are directly forwarded to this method.
The corresponding standard methods (ending without the ∟ symbol) are then automatically defined.
Other lensing methods
For debugging purposes it is useful to define a few other methods:
Gravity.redshift — Function
Return the redshift of a lens.
Gravity.moveto — Method
moveto(lens, x₀)Move a lens center to the location x₀.
Gravity.tiltto — Method
tiltto(lens, θₑ)Change the lens position angle to θₑ and return a new object.
Base.rand — Method
Random.rand([rng,] lenstype)Generate a new random lens of type lenstype.
Each LensType should define a way to generate a random lens of that type by defining
Random.rand(rng::AbstractRNG, ::Random.SamplerType{LensType})Example
Random.rand(rng::AbstractRNG, ::Random.SamplerType{PointLens}) = PointLens(rand(rng), randn(rng, SPoint), rand(rng))Lens trait
Gravity.jl also implements a special trait, LensTrait, used to dispatch deflection and Jacobian computations in cases where these functions are not defined. For example, one could have a lens model that does not have a Jacobian method that is simple-enough to be implemented: in this case, the Jacobian is calculated from the deflection using automatic differentiation.
Gravity.LensTrait — Type
Trait type for lens methods.
The three subtypes LensHasPotential, LensHasDeflection, and LensHasJacobian define singletons that can be used to specify the methods implemented for a specific lens model. This allows Julia to infer the correct procedure statically to use from the type of the lens at compile time, with no performance penalty.
Gravity.LensHasPotential — Type
Singleton for lenses that only provide the lensing potential.
All other quantities (deflection angle and jacobian) are computed with forward automatic differentiation using ForwardDiff.
Gravity.LensHasDeflection — Type
Singleton for lenses that provide the lensing potential and the deflection angle.
The jacobian of the deflection angle is computed with forward automatic differentiation using ForwardDiff.
Gravity.LensHasJacobian — Type
Singleton for lenses that provide the lensing potential, the deflection angle, and its jacobian.
These lenses will not need to use any automatic differentiation: this typically makes computations faster.